R/plot.R

In multiRec: Analysis of Multi-Type Recurrent Events

#' Plot profile likelihoods
#'
#' @param fit list, the fit object from multiRec. The fit object must be
#'   generated by calling multiRec with fitDetails=TRUE (the default is FALSE)
#' @param showInitial logical, if TRUE the initial parameter estimates are shown
#'   on the figures as a hollow dot.
#' @param maxiter integer, the number of iterations used to identify the
#'   x axis range on each figure.
#' @param tolerance numeric, a multiplicative factor used the set the y axis
#'   range on each figure: the axis will range from approximately the maximum
#'   of the log likelihood to (1+tolerance) times the maximum
#' @param n integer, the number of points on the profile likelihood plot
#'
#' @return No return value; called for side effects (creates profile likelihood plots).
#' @details
#' The function generates a series of plots, one for each parameter in the
#' model. Each plot shows the log likelihood on the vertical axis and a range of
#' values for the parameter on the horizontal axis. The plot for a specific
#' parameter is obtained by fixing the remaining parameters at their maximum
#' likelihood value, and varying the plotted parameter over a range of values.
#' The range of values is selected in such a way that the resulting likelihood
#' ranges from 90% of the maximum to 100% of the maximum.
#'
#' Plotting profile likelihoods is particularly important when using the
#' identity link as it tends to have ill behaved likelihoods. When looking at
#' profile plots, looks for multiple maxima, especially ones for which the
#' corresponding likelihoods are similar in magnitude: if they exist, the
#' interpretation of the fit returned by the model is questionable as there
#' are other solutions that fit the data (almost) as well.  Also look for
#' discontinuities and sharp changes in slope (i.e discontinuities in the first
#' derivatives) as these will often cause convergence problems.
#'
#' @export
#'
#' @importFrom graphics abline title box grconvertX lines par points text xinch
#' @importFrom grDevices gray
#' @importFrom stats approx uniroot
#' @examples
#' fit = multiRec(hf   ~ nPrior.afib(~male),
#'                afib ~ age + nPrior.afib() + nPrior.hf(),
#'                data=multiRecCVD2,
#'                idVar='id',
#'                link='log',
#'                fitDetails=TRUE)
#' plotProfiles(fit)
#'
plotProfiles = function(fit, n=100, showInitial=FALSE, maxiter=50, tolerance=0.1) {
  if (is.null(fit$fitDetails)) halt('Likelihood profile plots require fitDetails=TRUE in the call to multiRec ')

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar), add = TRUE)

  init_val   = fit$fitDetails$init
  mod_obj    = fit$fitDetails$modelObj
  parameters = fit$fitDetails$param

  for (i in seq_along(parameters$hazard)) {
    # Produce a profile for the ith parameter
    estimate = parameters$estimate[i]
    se       = parameters$se[i]
    param    = parameters$estimate
    y        = rep(NA,n)

    yLow = (1+tolerance)*fit$loglik
    f = function(x) {
      param[i] = x
      likelihood(param,mod_obj)-yLow
    }

    if (se==0) se = abs(estimate/3)
    interval = try(list(lo=uniroot(f,interval=c(estimate-se,estimate),maxiter=maxiter,extendInt='upX')$root,
                        hi = uniroot(f,interval=c(estimate,estimate+se),maxiter=maxiter,extendInt='downX')$root),
                   silent=TRUE)
    if (inherits(interval,'try-error')) {
      value = try(uniroot(f,interval=c(estimate-se,estimate+se),maxiter=maxiter,extendInt='yes')$root,
              silent=TRUE)
      if (inherits(value,'try-error')) {
        if (estimate==0) interval=list(lo=-1,hi=1)
        else interval=list(lo=estimate-se,hi=estimate+se)
      }
      else {
        delta = abs(value-estimate)
        interval=list(lo=estimate-delta,hi=estimate+delta)
      }
    }

    x = seq(interval$lo,interval$hi,length.out=n)

    for (j in seq_along(x)) {
      param[i] = x[j]
      y[j]     = likelihood(param,mod_obj)
    }

    # Plotting the MLE
    pts = approx(x,y,c(init_val[i],parameters$estimate[i]))
    par(las=1)
    plot(x,y,type = "n", yaxt='n', xlab='', ylab='')
    abline(h=pts$y[2],col=gray(0.8),lty='dashed')
    lines(x,y)
    box()

    top = fit$loglik
    text(grconvertX(0,from='npc')-xinch(0.1),
         top,
         signif(top,4),
         adj=c(1,0.5),
         xpd=NA)

    bot = min(y)
    text(grconvertX(0,from='npc')-xinch(0.1),
         bot,
         signif(bot,4),
         adj=c(1,0.5),
         xpd=NA)

    if (showInitial) points(pts$x[1],pts$y[1],pch=1)
    points(pts$x[2],pts$y[2],pch=19)

    title(main = paste0(parameters$hazard[i],parameters$component[i]))
  }
  invisible(NULL)

}