R/plot.logistf.profile.r

In logistf: Firth's Bias-Reduced Logistic Regression

#' \code{plot} Method for \code{logistf} Likelihood Profiles
#' 
#' Provides the plot method for objects created by \code{profile.logistf} or \code{CLIP.profile}
#' 
#' The plot method provides three types of plots (profile, CDF, and density representation of a profile 
#' likelihood). For objects generated by CLIP.profile, it also allows to show the completed-data 
#' profiles along with the pooled profile.
#'
#' @param x A \code{profile.logistf} object
#' @param type Type of plot: one of c("profile", "cdf", "density")
#' @param max1 If \code{type="density"}, normalizes density to maximum 1
#' @param colmain Color for main profile line
#' @param colimp color for completed-data profile lines (for \code{logistf.profile} objects that also 
#' carry the \code{CLIP.profile} class attribute)
#' @param plotmain if \code{FALSE}, suppresses the main profile line (for \code{logistf.profile} objects that 
#' also carry the \code{CLIP.profile} class attribute)
#' @param ylim Limits for the y-axis
#' @param ... Further arguments to be passed to \code{plot()}.
#'
#' @return The function is called for its side effects
#' 
#' @author Georg Heinze und Meinhard Ploner
#' @references Heinze G, Ploner M, Beyea J (2013). Confidence intervals after multiple imputation: combining 
#' profile likelihood information from logistic regressions. Statistics in Medicine, to appear.
#' @export
#'
#' @examples
#' 
#' data(sex2) 
#' fit<-logistf(case ~ age+oc+vic+vicl+vis+dia, data=sex2)
#' plot(profile(fit,variable="dia"))
#' plot(profile(fit,variable="dia"), "cdf")
#' plot(profile(fit,variable="dia"), "density")
#' 
#' #generate data set with NAs
#' freq=c(5,2,2,7,5,4)
#' y<-c(rep(1,freq[1]+freq[2]), rep(0,freq[3]+freq[4]), rep(1,freq[5]), rep(0,freq[6]))
#' x<-c(rep(1,freq[1]), rep(0,freq[2]), rep(1,freq[3]), rep(0,freq[4]), rep(NA,freq[5]),
#'    rep(NA,freq[6]))
#' toy<-data.frame(x=x,y=y)
#' 
#' # impute data set 5 times
#' set.seed(169)
#' toymi<-list(0)
#' for(i in 1:5){
#'    toymi[[i]]<-toy
#'    y1<-toymi[[i]]$y==1 & is.na(toymi[[i]]$x)
#'    y0<-toymi[[i]]$y==0 & is.na(toymi[[i]]$x)
#'    xnew1<-rbinom(sum(y1),1,freq[1]/(freq[1]+freq[2]))
#'    xnew0<-rbinom(sum(y0),1,freq[3]/(freq[3]+freq[4]))
#'    toymi[[i]]$x[y1==TRUE]<-xnew1
#'    toymi[[i]]$x[y0==TRUE]<-xnew0
#'  }
#'  
#' # logistf analyses of each imputed data set
#' fit.list<-lapply(1:5, function(X) logistf(data=toymi[[X]], y~x, pl=TRUE))
#' 
#' # CLIP profile 
#' xprof<-CLIP.profile(obj=fit.list, variable="x", data=toymi, keep=TRUE)
#' plot(xprof)
#' 
#' #plot as CDF
#' plot(xprof, "cdf")
#' 
#' #plot as density
#' plot(xprof, "density")
#' 
#' @exportS3Method plot logistf.profile
plot.logistf.profile <-
function(x, type="profile", max1=TRUE, colmain="black", colimp="gray", plotmain=T, ylim=NULL, ...){
  obj<-x
 if(any(attr(obj,"class")[1]=="logistf.CLIP.profile")){
   x<-obj$beta
   if(type=="profile") {
     y<-obj$profile
     y.imp<-obj$profile.matrix
     y.label<-"Relative log profile penalized likelihood"
     } 
   else if(type=="cdf") {
      y<-obj$cdf
      y.imp<-obj$cdf.matrix
      y.label<-"Cumulative distribution function"
     }   else
  if(type=="density") {
      x<-obj$beta[-length(obj$beta)]+diff(obj$beta)/2
      xy<-density(x, weights=diff(obj$cdf)/sum(diff(obj$cdf)))
      x<-xy$x
      y<-xy$y
      #y<-y/(sum(y)*sum(diff(obj$beta)))
      # norm to AUC=1
      auc<-function(y,x){
        y.mean<-(y[-length(y)]+y[-1])/2
        x.diff<-diff(x)
        res <- y.mean %*% x.diff
        return(res[1,1])
      }
      if(max1) y<-y/max(y)
      else y<-y/auc(y,x)
      ab<-t(apply(obj$cdf.matrix,1,diff))
      ab[ab<0]<-0
      ab<-ab/rowSums(ab)
      y.imp<-unlist(lapply(1:nrow(obj$cdf.matrix),function(Z) density(obj$beta[-length(obj$beta)]+diff(obj$beta)/2,weights=ab[Z,])$y))
      imputations<-nrow(obj$cdf.matrix)
      y.imp<-t(matrix(y.imp,length(y.imp)/imputations,imputations))
      if(max1) y.imp<-t(apply(y.imp,1,function(Z) Z/max(Z)))
      else y.imp<-t(apply(y.imp,1,function(Z) Z/auc(Z,x)))
      y.label<-"Posterior density"
  }
  plot(x=x, y=y, xlab=expression(gamma), ylab=y.label, type="l", lwd=2, col=colmain, ylim=ylim)
  if(type=="profile") abline(-3.841,0, lty=3)
  else if(type=="cdf") {
     abline(0.025,0, lty=3)
     abline(0.975,0, lty=3)
  }
  if(!is.null(y.imp)) {
     for(i in 1:nrow(y.imp)) lines(x=x, y=y.imp[i,], lwd=0.5, lty=2, col=colimp)
     if(plotmain) lines(x=x, y=y, lty=1, lwd=2, col=colmain)
     }
  } 
  else if(any(attr(obj,"class")=="logistf.profile")) {
  x<-obj$beta
  if(type=="profile") {
     y<-obj$profile
     y.label<-"Relative log profile penalized likelihood"
  } 
  else if(type=="cdf") {
    y<-obj$cdf
    y.label<-"Cumulative distribution function"
  } else if(type=="density") {
    x<-obj$beta[-length(obj$beta)]+diff(obj$beta)/2
    xy<-density(x, weights=diff(obj$cdf)/sum(diff(obj$cdf)))
    x<-xy$x
    y<-xy$y
#       y<-y/(sum(y)*sum(diff(obj$beta)))
# norm to AUC=1
    auc<-function(y,x){
      y.mean<-(y[-length(y)]+y[-1])/2
      x.diff<-diff(x)
      res <- y.mean %*% x.diff
      return(res[1,1])
    }
      y<-y/auc(y,x)
      y.label<-"Posterior density"
  }
  plot(x=x, y=y, xlab=expression(gamma), ylab=y.label, type="l", lwd=2, col=colmain)
  if(type=="profile") abline(max(y)-3.841,0, lty=3)
  else if (type=="cdf") {
    abline(0.025,0, lty=3)
    abline(0.975,0, lty=3)
  }
  } 
  else stop("Object must be of class profile or pool.profile.\n")
}