In rchendrickson/pmsimstats: Personalized medicine simulated statistics

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(pmsimstats)
library(data.table)
library(lme4)
library(lmerTest)
library(corpcor)
library(ggplot2)
library(MASS)
library(svMisc)
library(tictoc)
library(ggpubr)
library(gridExtra)

pmsimstats package

This is the second part of a pair of vignettes designed to walk through the steps used to generate the results published in [___]. The steps fall into four basic steps:

Define your trial designs
Define the parameter spaces you want to explore
Run the simulations and save the results
Look at the results

The first 3 steps are in the previous vignette, Produce_Publication_Results_1_generate_data. Here, we look at the results

Look at the trajectories of each factor on each path, all pre-censoring

data(results_trajectories)
simresults<-results_trajectories
data<-simresults$rawdata$precensor # this gives a list of length (total reps done)
tinfo<-simresults$results # dim[1] of tinfo gives the parameter space that went in
tds<-simresults$parameterselections$trialdesigns
mergeacrossreps<-TRUE
plots1<-vector(mode="list",length=length(simresults$parameterselections$trialdesigns))
for(ip in 1:length(plots1)){
  param2hold<-data.table(carryover_t1half=0,trialdesign=ip)
  plots1[[ip]]<-plotfactortrajectories(data,tinfo,tds,param2hold,mergeacrossreps)
}
legend<-as_ggplot(get_legend(plots1[[1]]))
# Add titles
for(ip in 1:(length(plots1))){
  plots1[[ip]]<-plots1[[ip]]+theme(legend.position="none",text=element_text(size=10))+
#    ggtitle(simresults$parameterselections$trialdesigns[[ip]]$metadata$name_longform)
    ggtitle(paste(c("","A","B","C","D","E","F","G")[ip],
                  simresults$parameterselections$trialdesigns[[ip]]$metadata$name_longform,
                  sep=". "))+
    theme(plot.title=element_text(size=8))
}
# Repeat with carryover term
plots2<-vector(mode="list",length=length(simresults$parameterselections$trialdesigns))
for(ip in 1:length(plots2)){
  param2hold<-data.table(carryover_t1half=.5,trialdesign=ip)
  plots2[[ip]]<-plotfactortrajectories(data,tinfo,tds,param2hold,mergeacrossreps)
}
# Add titles
for(ip in 1:(length(plots2))){
  plots2[[ip]]<-plots2[[ip]]+theme(legend.position="none",text=element_text(size=10))+
#    ggtitle(simresults$parameterselections$trialdesigns[[ip]]$metadata$name_longform)
    ggtitle(paste(c("","E","F","G","H")[ip],
                  simresults$parameterselections$trialdesigns[[ip]]$metadata$name_longform,
                  sep=". "))+
    theme(plot.title=element_text(size=8))
}

# adjust ylab and ylab and legends
for(ip in c(2,3,4)){
  plots1[[ip]]<-plots1[[ip]]+ylab("")
  plots2[[ip]]<-plots2[[ip]]+ylab("")
}
for(ip in c(1,2,3,4)){
  plots1[[ip]]<-plots1[[ip]]+xlab("")
}
plots1[[1]]<-plots1[[1]]+ylab("NO CARRYOVER \n Improvement in symptoms")
plots2[[1]]<-plots2[[1]]+ylab("WITH CARRYOVER \n Improvement in symptoms")

lay <- rbind(c(1,1,2,2,2,3,3,3,4,4,4,NA),
             c(5,5,6,6,6,7,7,7,8,8,8,9))
plots<-c(plots1,plots2,list(legend))
grid.arrange(grobs = plots, layout_matrix = lay)
#ml <- arrangeGrob(grobs = plots, layout_matrix = lay)
#setwd(sourcedir)
#ggsave("Fig3_factortrajectories.pdf",ml,width=9,height=5,units="in",dp=600)

Look at how our core response parameters affect power

# Our core parameters
data(results_core)
simresults<-results_core
param2vary<-c("trialdesign","N","c.bm","censorparamset") 
param2hold<-data.table(blparamset=1,respparamset=1,carryover_t1half=0)
param2nothold<-c("c.cfct","modelparamset")
param2plot<-"power" # Options: power, mean_frac_NA, mean_beta, mean_betaSE, sd_beta...
p1<-PlotModelingResults(simresults,param2plot,param2vary,param2hold,param2nothold)

param2vary<-c("trialdesign","N","carryover_t1half","censorparamset") 
param2hold<-data.table(blparamset=1,respparamset=1,c.bm=.6)
param2nothold<-c("c.cfct","modelparamset")
param2plot<-"power" # Options: power, mean_frac_NA, mean_beta, mean_betaSE, sd_beta...
p2<-PlotModelingResults(simresults,param2plot,param2vary,param2hold,param2nothold)

p1<-p1+ggtitle("A. Effect of trial design, censoring and biomarker effect on power")+
  theme(plot.title=element_text(hjust=0))
p2<-p2+ggtitle("B. Effect of nonzero carryover term on power")+
  theme(plot.title=element_text(hjust=0))

lay <- rbind(c(1),
             c(2))
plots<-list(p1,p2)
grid.arrange(grobs = plots, layout_matrix = lay)
#ml <- arrangeGrob(grobs = plots, layout_matrix = lay)
#setwd(sourcedir)
#ggsave("Fig4_CorePowerResults.pdf",ml,width=9,height=9,units="in",dp=600)

# Resp parameters:

# A - maxes, holding rates steady...
data(results_maxes)
simrsults<-results_maxes
param2vary<-c("trialdesign","tv_max","pb_max","br_max") 
param2hold<-data.table(blparamset=1,censorparamset=0,carryover_t1half=0,c.bm=.6)
param2nothold<-c("c.cfct","modelparamset","respparamset")
param2plot<-"power" # Options: power, mean_frac_NA, mean_beta, mean_betaSE, sd_beta...
p1<-PlotModelingResults(simresults,param2plot,param2vary,param2hold,param2nothold)

# B - rates, holding maxes steady...
data("results_rates")
d1<-results_rates
d1$results<-d1$results[respparamset%in%
                         (setdiff(1:length(d1$parameterselections$respparamsets),c(1:5)))]
param2vary<-c("trialdesign","tv_rate","pb_rate","br_rate") 
param2hold<-data.table(blparamset=1,censorparamset=1,carryover_t1half=0,c.bm=.3)
param2nothold<-c("c.cfct","modelparamset","respparamset")
p2<-PlotModelingResults(d1,param2plot,param2vary,param2hold,param2nothold)

p1<-p1+ggtitle("A. Effect of response parameter maximum values on power")+
  theme(plot.title=element_text(hjust=0))
p2<-p2+ggtitle("B. Effect of response parameter rate values on power")+
  theme(plot.title=element_text(hjust=0))

lay <- rbind(c(1),
             c(2))
plots<-list(p1,p2)
grid.arrange(grobs = plots, layout_matrix = lay)
#ml <- arrangeGrob(grobs = plots, layout_matrix = lay)
#setwd(sourcedir)
#ggsave("Fig5_CoreResponseparams.pdf",ml,width=9,height=9,units="in",dp=600)

Make the two plots having to do with variance and bias

# SE of the betas
data(results_core)
simresults<-results_core
param2vary<-c("trialdesign","respparamset","carryover_t1half","censorparamset") 
param2hold<-data.table(N=70,blparamset=1)
param2nothold<-c("c.cfct","modelparamset")
param2plot<-"mean_betaSE" # Options: power, mean_frac_NA, mean_beta, mean_betaSE, sd_beta...
p<-PlotModelingResults(simresults,param2plot,param2vary,param2hold,param2nothold)
p
#setwd(sourcedir)
#ggsave("Fig6_mean_betaSE.pdf",p,width=9,height=7,units="in",dp=600)

# Bias in the betas
data(results_core)
data<-results_core
d1<-data$results[carryover_t1half==0][respparamset==1]#[c.bm==.6]
d1$censorparamset<-as.factor(d1$censorparamset)
d1$N<-as.factor(d1$N)
# Change the labels on the stuff that will be faceted
bklabels<-vector(mode="character",length=length(data$parameterselections$trialdesigns))
for(iT in 1:length(bklabels)){
  bklabels[iT]<-data$parameterselections$trialdesigns[[iT]]$metadata$name_shortform[[1]]
}
d1$trialdesign<-as.factor(d1$trialdesign)
levels(d1$trialdesign)<-bklabels
bklabels<-vector(mode="character",length=length(data$parameterselections$censorparams))
for(iT in 1:length(bklabels)){
  bklabels[iT]<-data$parameterselections$censorparams$names[[iT]]
}
levels(d1$censorparamset)<-c("none",bklabels)

d1[,bdelta:=(beta-ETbeta)]
d1[,tpval:=as.numeric(NA)]
d1[,test:=as.numeric(NA)]
for(td in unique(d1$trialdesign)){
  for(cp in unique(d1$censorparamset)){
    for(bm in unique(d1$c.bm)){
      if(bm==0){
        tout<-t.test(x=d1[trialdesign==td][censorparamset==cp][c.bm==bm]$beta,mu=0)
      }else{
        tout<-t.test(x=d1[trialdesign==td][censorparamset==cp][c.bm==bm]$bdelta,mu=0)
      }
      d1[trialdesign==td][censorparamset==cp][c.bm==bm]$tpval<-tout$p.value
      d1[trialdesign==td][censorparamset==cp][c.bm==bm]$test<-tout$estimate
    }
  }
}
d2<-unique(d1[,.(trialdesign,censorparamset,c.bm,tpval,test)])
d2$test<-d2$test*1000
p1<-ggplot(data=d2[c.bm==0])+geom_point(aes(y=tpval,x=test,color=censorparamset,shape=censorparamset),size=2)+
  scale_y_continuous(trans='log10')+  geom_vline(aes(xintercept=0))+
  geom_hline(aes(yintercept=.05),linetype=3)+geom_hline(aes(yintercept=.0001),linetype=4)+
  xlab("Estimated bias in model coefficient")+ylab("Signifiance (p-value)")+
  labs(color="Censoring",shape="Censoring")+
  facet_grid(c.bm~trialdesign)+ggtitle("A. Bias with known 0 effect size")
p2<-ggplot(data=d2[c.bm>0])+geom_point(aes(y=tpval,x=test,color=censorparamset,shape=censorparamset),size=2)+
  scale_y_continuous(trans='log10')+  geom_vline(aes(xintercept=0))+
  geom_hline(aes(yintercept=.05),linetype=3)+geom_hline(aes(yintercept=.0001),linetype=4)+
  xlab("Estimated bias in model coefficient")+ylab("Signifiance (p-value)")+
  labs(color="Censoring",shape="Censoring")+
  facet_grid(c.bm~trialdesign)+ggtitle("B. Bias comapered to full data set without censoring")
plots<-list(p1,p2)

lay <- rbind(1,1,1,2,2,2,2)
grid.arrange(grobs = plots, layout_matrix = lay)
#ml <- arrangeGrob(grobs = plots, layout_matrix = lay)
#setwd(sourcedir)
#ggsave("Fig7_Bias.pdf",ml,width=9,height=7,units="in",dp=600)