R/nca.npde.plot.R

Defines functions nca.npde.plot

Documented in nca.npde.plot

# NPDE plot
# roxygen comments
#' Plots population histogram of the NCA metrics selected for model diagnosis.
#'
#' \pkg{nca.npde.plot} plots individual NPDE values and histogram of the NPDE
#' values within a population group
#'
#' \pkg{nca.npde.plot} individual NPDE values and histogram of the NPDE
#' values of NCA metrics within a population group.
#' 
#' @param plotdata Data frame with the values of the NPDE values of each
#'   individual for the NCA metrics
#' @param xvar Name of the independent variable column against which NPDE values
#'   will be plotted
#' @param npdecol Column names or column numbers of containing the NPDE values
#' @param figlbl Figure label based on dose identifier and/or population
#'   stratifier (\strong{NULL})
#' @param cunit Unit for concentration (default is \strong{\code{NULL}})
#' @param tunit Unit for time (default is \strong{\code{NULL}})
#'
#' @return returns a data frame with the mean and SD of population NPDE values
#'   of each NCA metric and two graphical objects created by arrangeGrob
#'   function for the individual and population histogram of the NPDE values
#' @export
#'

nca.npde.plot <- function(plotdata,
                          xvar=NULL,
                          npdecol=NULL,
                          figlbl=NULL,
                          cunit=NULL,
                          tunit=NULL){
  
  "npde" <- "type" <- "mcil" <- "mciu" <- "sdu" <- "sducil" <- "sduciu" <- "..density.." <- "sdl" <- "XVAR" <- "melt" <- "xlab" <- "ylab" <- "theme" <- "element_text" <- "unit" <- "geom_point" <- "facet_wrap" <- "scale_linetype_manual" <- "scale_color_manual" <- "guides" <- "guide_legend" <- "element_rect" <- "geom_histogram" <- "aes" <- "geom_vline" <- "ggplot" <- "labs" <- "..count.." <- "..PANEL.." <- "na.omit" <- "sd" <- "qt" <- "qchisq" <- "scale_y_continuous" <- "percent" <- "sdcil" <- "sdciu" <- NULL
  rm(list=c("npde","type","mcil","mciu","sdu","sducil","sduciu","..density..","sdl","XVAR","melt","xlab","ylab","theme","element_text","unit","geom_point","facet_wrap","scale_linetype_manual","scale_color_manual","guides","guide_legend","element_rect","geom_histogram","aes","geom_vline","ggplot","labs","..count..","..PANEL..","na.omit","sd","qt","qchisq","scale_y_continuous","percent","sdcil","sdciu"))
  
  if (!is.data.frame(plotdata)) stop("plotdata must be a data frame.")
  if (is.null(xvar)){
    stop("Missing x-variable against which the NPDE data will be plotted.")
  }else{
    names(plotdata)[which(names(plotdata)==xvar)] <- "XVAR"
  }
  if (is.null(npdecol)){
    stop("Missing column names/numbers of the data frame containing NPDE data.")
  }else if (!is.null(npdecol) && !is.numeric(npdecol)){
    if (!all(npdecol%in%names(plotdata))) stop("All column names given as NPDE data column must be present in plotdata")
  }else if (!is.null(npdecol) && is.numeric(npdecol)){
    if (any(npdecol <= 0) | (max(npdecol) > ncol(plotdata))) stop("Column number for NPDE data out of range of plotdata.")
    npdecol <- names(plotdata)[npdecol]
  }
  
  plotdata <- subset(plotdata, select=c("XVAR",npdecol))
  longdata <- melt(plotdata, measure.vars = npdecol)
  names(longdata)[c((ncol(longdata)-1):ncol(longdata))] <- c("type","npde")
  longdata <- na.omit(longdata)
  longdata <- subset(longdata, npde!="NaN")
  
  if (nrow(longdata)==0){
    return(list(forestdata=NULL))
  }else{
    npdecol         <- levels(factor(longdata$type))
    longdata$mean   <- numeric(nrow(longdata))
    longdata$sd     <- numeric(nrow(longdata))
    longdata$length <- numeric(nrow(longdata))
    longdata$ci     <- numeric(nrow(longdata))
    longdata$sd1    <- numeric(nrow(longdata))
    longdata$sd2    <- numeric(nrow(longdata))
    longdata$mcil   <- numeric(nrow(longdata))
    longdata$mciu   <- numeric(nrow(longdata))
    longdata$sdl    <- numeric(nrow(longdata))
    longdata$sdu    <- numeric(nrow(longdata))
    longdata$sdlcil <- numeric(nrow(longdata))
    longdata$sdlciu <- numeric(nrow(longdata))
    longdata$sducil <- numeric(nrow(longdata))
    longdata$sduciu <- numeric(nrow(longdata))
    
    for (i in 1:length(npdecol)){
      npdeval <- na.omit(as.numeric(longdata[longdata$type==npdecol[i], "npde"]))
      longdata[longdata$type==npdecol[i],"mean"]   <- mean(npdeval)
      longdata[longdata$type==npdecol[i],"sd"]     <- sd(npdeval)
      longdata[longdata$type==npdecol[i],"length"] <- length(npdeval)
      longdata[longdata$type==npdecol[i],"ci"]     <- abs(qt(0.025,length(npdeval)-1)*(sd(npdeval)/sqrt(length(npdeval))))
      longdata[longdata$type==npdecol[i],"sd1"]    <- abs(sd(npdeval)-(sd(npdeval)*sqrt((length(npdeval)-1)/qchisq(0.975,length(npdeval)-1))))
      longdata[longdata$type==npdecol[i],"sd2"]    <- abs(sd(npdeval)-(sd(npdeval)*sqrt((length(npdeval)-1)/qchisq(0.025,length(npdeval)-1))))
      longdata[longdata$type==npdecol[i],"sdcil"]  <- sd(npdeval)*sqrt((length(npdeval)-1)/qchisq(0.975,length(npdeval)-1))
      longdata[longdata$type==npdecol[i],"sdciu"]  <- sd(npdeval)*sqrt((length(npdeval)-1)/qchisq(0.025,length(npdeval)-1))
    }
    
    longdata$FCT    <- paste0(longdata$type,"\nmean=", signif(longdata$mean,3), "+/-", signif(longdata$ci,3), "\nSD=", signif(longdata$sd,3), "(-", signif(longdata$sd1,3), ",+", signif(longdata$sd2,3),")")
    
    longdata$mcil   <- longdata$mean - longdata$ci
    longdata$mciu   <- longdata$mean + longdata$ci
    longdata$sdl    <- longdata$mean - longdata$sd
    longdata$sdu    <- longdata$mean + longdata$sd
    longdata$sdlcil <- longdata$mean - longdata$sd - longdata$sd1
    longdata$sdlciu <- longdata$mean - longdata$sd + longdata$sd2
    longdata$sducil <- longdata$mean + longdata$sd - longdata$sd1
    longdata$sduciu <- longdata$mean + longdata$sd + longdata$sd2
    
    fctNm <- data.frame()
    npr   <- length(npdecol)
    nc    <- ifelse(npr<2, 1, ifelse(npr>=2 & npr<=6, 2, 3))
    for (p in 1:npr){
      if (npdecol[p] == "npdeAUClast" | npdecol[p] == "npdeAUClower_upper" | npdecol[p] == "npdeAUCINF_obs" | npdecol[p] == "npdeAUCINF_pred"){
        if(is.null(cunit) | is.null(tunit)){
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=gsub("^npde", "", npdecol)[p]))
        }else{
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=paste0(gsub("^npde", "", npdecol)[p]," (",cunit,"*",tunit,")")))
        }
      }else if (npdecol[p] == "npdeAUMClast"){
        if(is.null(cunit) | is.null(tunit)){
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=gsub("^npde", "", npdecol)[p]))
        }else{
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=paste0(gsub("^npde", "", npdecol)[p]," (",cunit,"*",tunit,"^2)")))
        }
      }else if (npdecol[p] == "npdeCmax"){
        if(is.null(cunit)){
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=gsub("^npde", "", npdecol)[p]))
        }else{
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=paste0(gsub("^npde", "", npdecol)[p]," (",cunit,")")))
        }
      }else if (npdecol[p] == "npdeTmax" | npdecol[p] == "npdeHL_Lambda_z"){
        if(is.null(cunit)){
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=gsub("^npde", "", npdecol)[p]))
        }else{
          fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=paste0(gsub("^npde", "", npdecol)[p]," (",tunit,")")))
        }
      }else{
        fctNm <- rbind(fctNm, data.frame(prmNm=npdecol[p],prmUnit=npdecol[p]))
      }
    }
    longdata$type <- factor(longdata$type, levels=fctNm$prmNm, labels=fctNm$prmUnit)
    forestdata    <- subset(longdata[!duplicated(longdata$type),], select=c(type,mean,mcil,mciu,sd,sdcil,sdciu))
    
    # ggplot options for individual NPDE plots
    nc    <- ifelse(length(npdecol)<2, 1, ifelse(length(npdecol)>=2 & length(npdecol)<=6, 2, 3))
    ggOpt_npde <- list(xlab("\nID"), ylab("NPDE\n"),
                       theme(axis.text.x = element_text(angle=45,vjust=1,hjust=1,size=10),
                             axis.text.y = element_text(hjust=0,size=10),
                             strip.text.x = element_text(size=10),
                             legend.text = element_text(size=12),
                             title = element_text(size=14,face="bold"),
                             legend.position = "bottom", legend.direction = "horizontal",
                             legend.background = element_rect(),
                             legend.key.height = unit(1,"cm")),
                       geom_point(size=2), facet_wrap(~type, ncol=nc, scales="free"))
    
    # ggplot options for histogram of NPDE
    ggOpt_histnpde <- list(scale_linetype_manual(name="",values=c("meanNormal"="solid","meanNPDE"="solid","SD"="dashed")),
                           scale_color_manual(name="",values=c("meanNormal"="red","meanNPDE"="blue","SD"="blue")),
                           xlab("\nNPDE"), ylab(""),
                           guides(fill = guide_legend(override.aes = list(linetype = 0 )), shape = guide_legend(override.aes = list(linetype = 0))),
                           theme(axis.text.x = element_text(angle=45,vjust=1,hjust=1,size=10),
                                 axis.text.y = element_text(hjust=0,size=10),
                                 strip.text.x = element_text(size=10),
                                 legend.text = element_text(size=12),
                                 title = element_text(size=14,face="bold"),
                                 legend.position = "bottom", legend.direction = "horizontal",
                                 legend.background = element_rect(),
                                 legend.key.height = unit(1,"cm")),
                           scale_y_continuous(labels = percent),
                           geom_vline(aes(xintercept=0, color="meanNormal", lty="meanNormal"), size=1, show.legend=T),
                           geom_vline(aes(xintercept=as.numeric(mean), color="meanNPDE", lty="meanNPDE"), size=1),
                           geom_vline(aes(xintercept=as.numeric(sdl), color="SD", lty="SD"), size=1),
                           geom_vline(aes(xintercept=as.numeric(sdu), color="SD", lty="SD"), size=1),
                           facet_wrap(~FCT, ncol=nc, scales="free"))
    
    if (is.null(figlbl)){
      ttl   <- paste0("NPDE vs. ",xvar,"\n")
    }else{
      ttl   <- paste0("NPDE vs. ",xvar," (",figlbl,")\n")
    }
    ggnpde <- ggplot(longdata,aes(as.numeric(as.character(XVAR)),as.numeric(as.character(npde)))) + ggOpt_npde + labs(title = ttl) + theme(plot.title = element_text(size = rel(0.85)))
    
    if (is.null(figlbl)){
      ttl   <- paste0("Histogram of NPDE\n")
    }else{
      ttl   <- paste0("Histogram of NPDE (",figlbl,")\n")
    }
    
    bw <- diff(range(as.numeric(longdata$npde)))/(2*IQR(as.numeric(longdata$npde)))/length(as.numeric(longdata$npde))^(1/3)
    gghnpde <- ggplot(longdata,aes(as.numeric(as.character(npde)))) +
      geom_histogram(aes(y=(..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..]), size=0.6, color="black", fill="white", binwidth = bw) +
      ggOpt_histnpde + labs(title = ttl) + theme(plot.title = element_text(size = rel(0.85)))
    return(list(forestdata=forestdata,ggnpde=ggnpde,gghnpde=gghnpde))
  }
}

Try the ncappc package in your browser

Any scripts or data that you put into this service are public.

ncappc documentation built on May 1, 2019, 7:31 p.m.