study27_osu5115_analysis.R

# Analysing effect of adding/removing data points
# -----------------------------------------------------------------------------
# Enhancing the algorithm paper by determining the loss of accuracy/precision
#   when reducing the number of sample times. Aiming to prove that less time
#   points are needed.
# Data is from OSU6003; papers below:

# Blum W, Klisovic RB, Becker H, et al. Dose Escalation of Lenalidomide in 
#   Relapsed or Refractory Acute Leukemias. J Clin Oncol. 2010;28(33):4919-4925.
# Maddocks K, Ruppert AS, Browning R, et al. A dose escalation feasibility study
#   of lenalidomide for treatment of symptomatic, relapsed chronic lymphocytic 
#  leukemia. Leuk Res. 2014;38(9):1025-1029.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Load packages
  library(GA)  # genetic algorithm (ga)
  library(plyr)  # iterative functions (ddply)

# Load data visualisation package (set theme, define palette)
  library(ggplot2)  # plots (ggsave etc.)
  theme_bw2 <- theme_set(theme_bw(base_size = 14))
  theme_update(plot.title = element_text(hjust = 0.5))
  cbPalette <- data.frame(
		grey = "#999999",
		orange = "#E69F00",
		skyblue = "#56B4E9",
		green = "#009E73",
		yellow = "#F0E442",
		blue = "#0072B2",
		red = "#D55E00",
		pink = "#CC79A7",
		stringsAsFactors = F
	)
  
# Source algorithm functions
  source("fn_diag/fix_functions.R")

# Define additional functions
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Last observation carried forward
  # Finds an NA and carries forward the previous value
  locf <- function(x) {
    good <- !is.na(x)
    positions <- seq(length(x))
    good.positions <- good * positions
    last.good.position <- cummax(good.positions)
    last.good.position[last.good.position == 0] <- NA
    x[last.good.position]
  }
  
# Determine AUC (regular)
  auc.regular <- function(times, C, log = F) {
    auc <- c(NULL)
    for (i in 2:length(C)) {
      h <- times[i] - times[i-1]
      dC <- C[i-1] - C[i]
      if (log & dC > 0) auc[i-1] <- dC*h/log(C[i-1]/C[i])
      else auc[i-1] <- (C[i-1] + C[i])*h/2
    }
    return(sum(auc))
  }
  auc.regular.lam <- function(t, dv) {
    lambz <- try(-pred.lambdaz(dv, t)[2])
    if (class(lambz) == "try-error") return(0)
    else return(tail(dv, 1)/lambz)
  }
  
# 2 comp. w/ abs sumexp
  pred.d2a <- function(x, p) {
    exp(p[1]*x + p[4]) + exp(p[2]*x + p[5]) - exp(p[3]*x + log(sum(exp(p[4]), exp(p[5]))))
  }
  
 #90% confidence interval functions
  CI90lo <- function(x) quantile(x, probs = 0.05)
  CI90hi <- function(x) quantile(x, probs = 0.95)
  
# Load/Tidy/Prepare the Data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Source lenalidomide data
  nmprep <- read.csv(stringsAsFactors = F, 
    "C:/Users/Jim Hughes/Documents/GitRepos/LEN_PK/Data/data.csv"
  )
# Subset the data to look at OSU6003, without MDVs (excluding time = 0)
# Additionally the weekly troughs taken after the first 48 hours
  clindata <- nmprep[
    nmprep$STUDY == 5115 & 
    (nmprep$MDV == 0 | nmprep$TIME == 0) &
    nmprep$TIME <= 48,
  ]
  # Check the AMT column
  table(clindata$AMT)
  # needs to be filled using last observation carried forward
  clindata$AMT[clindata$AMT == 0] <- NA
  clindata$DOSE <- locf(clindata$AMT)
  clindata <- clindata[
    clindata$DOSE == 25,
  ]
  
# Check the times in the dataset
  table(clindata$TIME)
  # times are not uniform, must be binned to make mean pooled data
# Define nominal sample times
  samp.times <- c(0, 0.5, 1, 2, 3, 4, 6, 8, 24)
  
# Create breaks to use with cut function
  # First find the mean of each pair of values to choose a break half way 
  # between each sample time
  inner.breaks <- rowMeans(cbind(head(samp.times, -1), tail(samp.times, -1)))
  # Then create first and last break points by extending them before and after
  # the final times by the first and final inner break points
  cut.breaks <- c(head(samp.times, 1) - head(inner.breaks, 1), 
    inner.breaks, 
    tail(samp.times, 1) + tail(inner.breaks, 1))
  # and now use cut
  bin.groups <- cut(clindata$TIME, cut.breaks)
  
# Use these bins to define the nominal time after dose 
# Remembering to convert from factor to numeric
  levels(bin.groups) <- samp.times
  clindata$TADNOM <- as.numeric(paste(bin.groups))
  
# Check nominal times in the dataset
  table(clindata$TADNOM)
  
# Create mean pooled data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Use ddply to mean the normalised DV for each of the bins
  meanpool <- ddply(clindata, .(TADNOM), function(df) {
    out <- data.frame(MEAN = mean(df$DV))
  })
  
# Visualisation of mean pooled data
  p1 <- NULL
  p1 <- ggplot()
  p1 <- p1 + geom_line(aes(x = TADNOM, y = MEAN), data = meanpool, 
    col = cbPalette$red, size = 0.8)
  p1 <- p1 + geom_point(aes(x = TIME, y = DV), data = clindata,
    col = cbPalette$blue, size = 1, shape = 1, alpha = 0.7)
  p1  # looks pretty decent
  
# Fit sum of exponentials to mean pooled data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Optimise sum of exponential equation, then choose one with best aic
  sumexp.res <- best.sumexp.aic(optim.sumexp.new(meanpool, oral = T, nexp = 2))
  opt.sumexp <- sumexp.res$sumexp
  
# Visualisation of sum of exponentials
  # Create dataset
  sumexpdata <- data.frame(
    TIME = seq(0, 24, 0.1), 
    PRED = pred.sumexp(opt.sumexp, seq(0, 24, 0.1))
  )
  # Plot dataset against mean pooled
  p2 <- p1 + geom_line(aes(x = TIME, y = PRED), sumexpdata,
    col = cbPalette$blue, size = 0.8, alpha = 0.8)
  p2  # looks good!
  
# Determine optimal sample times and compare AUC
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Optimise time intervals for various numbers of times (6-13)
# Output is list due to variable lengths of time vectors
  interv.res <- llply(6:9, function(i) {
    metric <- llply(1:1000, function(j) {
      res <- optim.interv.dtmax(opt.sumexp, seq(0, 24, length.out = i))
      out.times <- res$times
      auc024 <- auc.interv.sumexp(out.times, opt.sumexp)
      pred <- pred.sumexp(opt.sumexp, out.times)
      list(
        times = rbind(out.times),
        auc024 = auc024,
        auc0inf = auc024 + auc.interv.lam(opt.sumexp, out.times),
        cmax = max(pred),
        tmax = out.times[which(pred == max(pred))]
      )
    })
    out.metrics <- cbind(
      colwise(mean)(ldply(metric, function(x) {x$auc024})),
      colwise(CI90lo)(ldply(metric, function(x) {x$auc024})),
      colwise(CI90hi)(ldply(metric, function(x) {x$auc024})),
      colwise(mean)(ldply(metric, function(x) {x$auc0inf})),
      colwise(CI90lo)(ldply(metric, function(x) {x$auc0inf})),
      colwise(CI90hi)(ldply(metric, function(x) {x$auc0inf})),
      colwise(mean)(ldply(metric, function(x) {x$cmax})),
      colwise(CI90lo)(ldply(metric, function(x) {x$cmax})),
      colwise(CI90hi)(ldply(metric, function(x) {x$cmax})),
      colwise(mean)(ldply(metric, function(x) {x$tmax})),
      colwise(CI90lo)(ldply(metric, function(x) {x$tmax})),
      colwise(CI90hi)(ldply(metric, function(x) {x$tmax}))
    )
    names(out.metrics) <- c(
      "mean.auc024", "CI90lo.auc024", "CI90hi.auc024",
      "mean.aucinf", "CI90lo.aucinf", "CI90hi.aucinf",
      "mean.cmax", "CI90lo.cmax", "CI90hi.cmax",
      "mean.tmax", "CI90lo.tmax", "CI90hi.tmax"
    )
    list(
      mean.times = colwise(mean)(ldply(metric, function(x) {x$times})),
      metrics = out.metrics
    )
  })

# Determine AUC0-24 of original mean pooled data
  out <- list(
    auc024 = data.frame(
      data = with(meanpool, auc.regular(TADNOM, MEAN)),
      sumexp = integrate(pred.d2a, 0, 24, p = opt.sumexp)$value
    ),
    aucinf = data.frame(
      data = with(meanpool, 
        auc.regular(TADNOM, MEAN) + auc.regular.lam(TADNOM, MEAN)
      ),
      sumexp = integrate(pred.d2a, 0, 168, p = opt.sumexp)$value),
    cmax = data.frame(
      data = with(meanpool, max(MEAN)),
      sumexp = pred.sumexp(opt.sumexp, tmax.sumexp(opt.sumexp))),
    tmax = data.frame(
      data = with(meanpool, TADNOM[which(MEAN == max(MEAN))]),
      sumexp = tmax.sumexp(opt.sumexp))
  )
  
  saveRDS(interv.res, "osu5115_optim.rds")
  saveRDS(out, "osu5115_simandobs.rds")