R/psnCdd.r

In PKgraph: Model diagnostics for population pharmacokinetic models

## PsN R script for plotting results of cdd
## Justin Wilkins
## November 2005
## Lars Lindbom
## August 2006


## set basic options
## autogenerated code begins
## we can set other things, like colours, line types, layout, etc here if we wish

## Code modified for PKgraph
psn.cdd <- function(file1, file2)
{
    min.failed    <- FALSE     # do we want to keep minimization failed runs?
    cov.failed    <- FALSE      # do we want to keep covariance failed runs?
    cov.warnings  <- TRUE      # do we want to keep covariance warning runs?
    boundary      <- TRUE      # do we want to keep boundary runs?

    # We have three options for the markers:
    # 1. as circles
    # 2. as text (taken from the column identified by the -case_column
    #    option to the PsN cdd command)
    # 3. as circles and text. The subset of runs which fall outside of 2
    #    standard deviations around the center of the cook-scores and
    #    covariance ratios as calculated using principal components
    #    analysis are marked as in 2.

    markeropt <- 2

    ## autogenerated code ends

    ## read files
    cdd.data <- read.csv(file1) #read.csv("raw_results1.csv")
    cdd.inds <- read.csv(file2, header=F) #read.csv("skipped_individuals1.csv", header=F)

    if (nrow(cdd.data)!= (nrow(cdd.inds) + 1)) return(NULL)
    check.col <- c("minimization_successful", "covariance_step_successful",
                   "covariance_step_warnings", "estimate_near_boundary")
    if (!all(check.col %in% colnames(cdd.data))) return(NULL)

    #cdd.data    <- subset(cdd.data, !is.na(ofv))
    #cdd.skipped <- subset(cdd.data, is.na(ofv))

    names(cdd.inds)[1] <- "ID"

    ## skip first record (no cases deleted)
    cdd.data <- cdd.data[-1,]

    cdd.data <- cbind(cdd.data,cdd.inds[1])

    ## create data frame for current parameter
    p1 <- cdd.data

    #names(p1)[1] <- "minimization.successful"
    #names(p1)[2] <- "covariance.step.successful"
    #names(p1)[3] <- "covariance.step.warnings"
    #names(p1)[4] <- "estimate.near.boundary"
    colnames(p1) <- gsub("_", ".", colnames(p1))

    if (min.failed) {
      mf  <- subset(p1, minimization.successful == 0)
      p1  <- subset(p1, minimization.successful == 1)
    }
    if (cov.failed) {
      cf <- subset(p1, covariance.step.successful == 0)
      p1 <- subset(p1, covariance.step.successful == 1)
    }
    if (!cov.warnings) {
      cw <- subset(p1, covariance.step.warnings == 1)
      p1 <- subset(p1, covariance.step.warnings == 0)
    }
    if (!boundary) {
      nb <- subset(p1, estimate.near.boundary == 1)
      p1 <- subset(p1, estimate.near.boundary == 0)
    }

    if (exists("mf")) {
      cdd.warn <- mf
    }
    if (exists("cf")) {
      if (exists("cdd.warn")) {
        cdd.warn <- data.frame(cdd.warn, cf)
      } else {
        cdd.warn <- cf
      }
    }
    if (exists("cw")) {
      if (exists("cdd.warn")) {
        cdd.warn <- data.frame(cdd.warn, cw)
      } else {
        cdd.warn <- cw
      }
    }
    if (exists("nb")) {
      if (exists("cdd.warn")) {
        cdd.warn <- data.frame(cdd.warn, nb)
      } else {
        cdd.warn <- nb
      }
    }

    if( markeropt == 1 ) {
      cdd.txt  <- subset(p1, FALSE)
      cdd.pt   <- subset(p1, TRUE)
    }

    if( markeropt == 2 ) {
      cdd.txt  <- subset(p1, TRUE)
      cdd.pt   <- subset(p1, FALSE)
    }

    if( markeropt == 3 ) {
      cdd.txt  <- subset(p1, outside.n.sd == 1)
      cdd.pt   <- subset(p1, outside.n.sd == 0)
    }
    
    #plot (cdd.pt$cov.ratios, cdd.pt$cook.scores,
     # type="p",
     # xlab="Cook score",
     # ylab="Covariance ratio",
     # main="Case-deletion diagnostics",
     # xlim=c(0,max(cdd.data$cook.scores, na.rm=T)),
     # ylim=c(0,max(cdd.data$cov.ratio, na.rm=T))
    #)

    #text(cdd.txt$cook.scores, cdd.txt$cov.ratios, labels=as.character(cdd.txt$ID),cex=.8, col=2)
    #if (exists("cdd.warn")) {
    #  text(cdd.warn$cook.scores, cdd.warn$cov.ratios, labels=as.character(cdd.warn$ID),cex=.8, col=4)
    #}

    return(list(ID=cdd.txt$ID, cook.scores=cdd.txt$cook.scores, cov.ratios=cdd.txt$cov.ratios))


}