In eniclas/FAIRsimulator: Simulation Platform for the FAIR Trial
Model 1
library(FAIRsimulator)
set.seed(3243)
theme_set(theme_gray(base_size=14))
Treatments
- Cell 0 (SoC)
- Counseling
- 6 month HAZ effect: 0
- Cell 1
- Counseling
- Food: 7.5% PER, 10% ASP
- Iron + Vitamin A + Zinc
- 6 month HAZ effect: 0.0633
- Cell 2
- Counseling
- Food: 15% PER, 30% ASP
- Multiple micronutrient supplementation
- 6 month HAZ effect: 0.1037
- Cell 3
- Counseling
- Food: 15% PER, 30% ASP
- Multiple micronutrient supplementation
- Antibiotics + Probiotics
- 6 month HAZ effect: 0.1574
- Cell 4
- Counseling
- Food: 15% PER, 30% ASP
- Multiple micronutrient supplementation
- Antibiotics + Probiotics
- WASH Intensified
- Psychosocial care and support
- 6 month HAZ effect: 0.1687
Study design settings
These study design options were used for the current set of simulations:
- 3 recruitment cycles: 0, 6, and 12 months (study time).
- 2 recruitment age groups: 6 and 12 months (age).
- Instantaneous recruitment (all children recruited on day 1).
- The initial randomization set to 1:1:1:1:1.
- The minimum randomization probability for SoC is 20%.
- Treatment arms with a probability less than 10% before adjusting for the minimum SoC randomization probability is dropped.
- All available data up and including the timepoint for the interim analyses is used to update the probabilities.
- The dropout rate is approximately 20%/6 months.
- Bi-monthly HAZ observations.
- Interim analyses every 6 months.
- Square root probability updates (small probabilities are inflated)
## One adapted cohort
AddNewSixMonthCohortEvent<-function(StudyObj) {
#Note - assuming birth cohort is first element in StudyDesignSettings
if (!is.null(StudyObj$CohortList)) {
for (i in 1:length(StudyObj$CohortList)) {
Cohort<-StudyObj$CohortList[[i]]
if ((Cohort$Active==FALSE &&
Cohort$RandomizationAgeRange[1]==6*30 &&
Cohort$CycleNum==1 &&
Cohort$CohortStartTime+Cohort$CurrentTime==StudyObj$CurrentTime &&
StudyObj$CurrentTime<=StudyObj$StudyDesignSettings$LatestTimeForNewBirthCohorts) || #If this cohort just ended and its time to add a new cohort
(Cohort$Active==TRUE &&
Cohort$RandomizationAgeRange[1]==6*30 &&
Cohort$CycleNum==1 &&
Cohort$CohortStartTime+6*30==StudyObj$CurrentTime &&
StudyObj$CurrentTime<=StudyObj$StudyDesignSettings$LatestTimeForNewBirthCohorts))
{
BirthCohortIndex<-which(unlist(lapply(StudyObj$StudyDesignSettings$CohortAgeRange,function(x) {x[1]==6*30}))==TRUE) #Get the birth cohort index (i.e. lower AgeRangeAtRandomization=0)
NewChildCohort<-NewCohort(StudyObj,CohortNum=BirthCohortIndex)
NewChildCohort$StartNum <- max(StudyObj$CohortList %listmap% "StartNum")+1
NewChildCohort$Name<-paste0("C-",NewChildCohort$StartNum,"-",1," [",Cohort$RandomizationAgeRange[1]/30,"-",Cohort$RandomizationAgeRange[2]/30,"m @ rand]")
DebugPrint(paste0("Create new 6 month cohort ",NewChildCohort$Name," based on probabilities in ",Cohort$Name," at time: ",StudyObj$CurrentTime),1,StudyObj)
NewChildCohort$RandomizationProbabilities<-Cohort$UpdateProbabilities #Use the latest probabilities from previous birth cohort
NewChildCohort$UnWeightedRandomizationProbabilities<-Cohort$UnWeightedUpdateProbabilities
StudyObj$CohortList[[(length(StudyObj$CohortList)+1)]]<-NewChildCohort
break
}
}
}
return(StudyObj)
}
## One adapted cohort
AddNewTwelveMonthCohortEvent<-function(StudyObj) {
#Note - assuming birth cohort is first element in StudyDesignSettings
if (!is.null(StudyObj$CohortList)) {
for (i in 1:length(StudyObj$CohortList)) {
Cohort<-StudyObj$CohortList[[i]]
if (Cohort$Active==FALSE && Cohort$RandomizationAgeRange[1]==12*30 && Cohort$CycleNum==1 && Cohort$CohortStartTime+Cohort$CurrentTime==StudyObj$CurrentTime && StudyObj$CurrentTime<=StudyObj$StudyDesignSettings$LatestTimeForNewBirthCohorts) #if this cohort just ended and its time to add a new cohort
{
BirthCohortIndex<-which(unlist(lapply(StudyObj$StudyDesignSettings$CohortAgeRange,function(x) {x[1]==12*30}))==TRUE) #Get the birth cohort index (i.e. lower AgeRangeAtRandomization=0)
NewChildCohort<-NewCohort(StudyObj,CohortNum=BirthCohortIndex)
NewChildCohort$StartNum <- max(StudyObj$CohortList %listmap% "StartNum")+1
NewChildCohort$Name<-paste0("C-",NewChildCohort$StartNum,"-",1," [",Cohort$RandomizationAgeRange[1]/30,"-",Cohort$RandomizationAgeRange[2]/30,"m @ rand]")
DebugPrint(paste0("Create new 12 month cohort ",NewChildCohort$Name," based on probabilities in ",Cohort$Name," at time: ",StudyObj$CurrentTime),1,StudyObj)
NewChildCohort$RandomizationProbabilities<-Cohort$UpdateProbabilities #Use the latest probabilities from previous birth cohort
NewChildCohort$UnWeightedRandomizationProbabilities<-Cohort$UnWeightedUpdateProbabilities
StudyObj$CohortList[[(length(StudyObj$CohortList)+1)]]<-NewChildCohort
break
}
}
}
return(StudyObj)
}
probTemperation <- function(probs) {
probs <- sqrt(probs)/sum(sqrt(probs))
return(probs)
}
InterimAnalyzesTime<-function(Cohort,StudyObj) {
DebugPrint(paste0("Check if cohort ",Cohort$Name," is about to have an interim analysis",StudyObj$CurrentTime),3,StudyObj)
TimeToPerformInterim<-FALSE
tmp<-sum(unlist(lapply(Cohort$SubjectList,function(Subject,LastSample){
return(as.numeric(Subject$Status==0 || (LastSample %in% Subject$SubjectSampleTime)))},max(Cohort$SamplingDesign))))/Cohort$MaxNumberOfSubjects
if (length(tmp)==0) tmp<-0
if (Cohort$RandomizationAgeRange[1]==6*30 && (StudyObj$CurrentTime==6*30-1 || StudyObj$CurrentTime==18*30-1 )) {
DebugPrint(paste0("Time to do an interim analyses for cohort ",Cohort$Name," at time: ",StudyObj$CurrentTime," (",round(tmp*100,1)," % subjects completed)"),1,StudyObj)
TimeToPerformInterim<-TRUE
}
return(TimeToPerformInterim)
}
FinalAnalysesEvent<-function(StudyObj) {
DebugPrint(paste0("Check if time to do a final analyses at study time: ",StudyObj$CurrentTime),3,StudyObj)
if (StudyObj$CurrentTime==StudyObj$StudyDesignSettings$StudyStopTime-1) {
DebugPrint(paste0("Time to do a final analyses with all data at time: ",StudyObj$CurrentTime),1,StudyObj)
df<-getAllSubjectData(StudyObj,prevTreatment = TRUE) #Get all data that is neede for final analyses
## Rename, select and order columns
if(length(names(df)[grep(names(df),pattern = "PTRT")]) >0) {
df <- df %>%
rename(ID=StudyID,DATA=HAZ,AGE=Age,TRT=TreatmentIndex,TRTS=Treatment) %>%
select(ID,DATA,AGE,TRT,TRTS,one_of(StudyObj$StudyDesignSettings$Covariates),Level,CohortName,matches("PTRT"))
for (cstr in names(df)[grep(names(df),pattern = "PTRT")]) { #Set no previous treatment to "0"
df[,cstr]<-ifelse(is.na(df[,cstr]),0,df[,cstr])
}
} else {
df <- df %>%
rename(ID=StudyID,DATA=HAZ,AGE=Age,TRT=TreatmentIndex,TRTS=Treatment) %>%
select(ID,DATA,AGE,TRT,TRTS,one_of(StudyObj$StudyDesignSettings$Covariates),Level,CohortName)
}
df[df==-99]<-NA #Set -99 to missing
df<-ImputeCovariates(df,StudyObj,method=StudyObj$StudyDesignSettings$ImpMethod) #Impute missing covariates
#Sort to get correct TRT order
df <- df[order(df["TRT"],df["ID"], df["AGE"]),]
##### Make some covariate factors
df$TRT<-as.factor(df$TRT)
df$SEXN<-as.factor(df$SEXN)
df$SANITATN<-as.factor(df$SANITATN)
df$AGE<-df$AGE/(12*30) #Rescale time to years
### If we have previous treatments
ptrti <- grep("PTRT",names(df))
if(length(ptrti)>0) {
for(i in 1:length(ptrti)) {
df[,ptrti[i]] <- as.factor(df[,ptrti[i]])
}
myPTRTs <- names(df)[ptrti]
}
cohortlevels<-unique(StudyObj$CohortList %listmap% "Level")
for (l in 1:length(cohortlevels)) {
dflevel<-subset(df,Level==cohortlevels[l])
#### Perform LME estimation based on some covariates and treatment effects for each level
if(length(ptrti)>0) {
lmeFormula <- paste0("DATA~1 + AGE + AGE:TRT + (AGE|ID) +",paste0(c(StudyObj$StudyDesignSettings$Covariates,myPTRTs),collapse = " + "))
} else {
lmeFormula <- paste0("DATA~1 + AGE + AGE:TRT + (AGE|ID) +",paste0(StudyObj$StudyDesignSettings$Covariates,collapse = " + "))
}
lmefit <- lmer(lmeFormula,data=dflevel,REML=FALSE) # IIV on baseline only
##### Calculate new probabilites based on another cohort LME results
lmecoef<-summary(lmefit)$coefficients[,1] #Get coefficicents from LME
lmese<-summary(lmefit)$coefficients[,2] #Get SE from LME
lmecoef<-lmecoef[regexpr('AGE:TRT.*',names(lmecoef))==1]
lmese<-lmese[regexpr('AGE:TRT.*',names(lmese))==1]
print(lmefit)
### Get first cohort with level = l
if (!is.null(StudyObj$CohortList)) {
for (i in 1:length(StudyObj$CohortList)) {
if (StudyObj$CohortList[[i]]$Level==l) {
Cohort<-StudyObj$CohortList[[i]]
break
}
}
if ((length(lmecoef)+1)!=length(Cohort$RandomizationProbabilities)) {
lmecoefnew<-rep(0,length(Cohort$RandomizationProbabilities)-1)
lmesenew<-rep(0,length(Cohort$RandomizationProbabilities)-1)
for (i in 1:(length(Cohort$RandomizationProbabilities)-1)) {
iIndex<-which(names(lmecoef)==paste0("AGE:TRT",i+1))
if (length(iIndex)!=0) {
lmecoefnew[i]<-lmecoef[iIndex]
lmesenew[i]<-lmese[iIndex]
}
}
names(lmecoefnew)<-paste0("AGE:TRT",2:length(Cohort$RandomizationProbabilities))
names(lmesenew)<-paste0("AGE:TRT",2:length(Cohort$RandomizationProbabilities))
lmecoef<-lmecoefnew
lmese<-lmesenew
}
}
#Get probability of beeing best
probs <- GetNewRandomizationProbabilities(trtcoeff=lmecoef,trtse=lmese,
StudyObj$StudyDesignSettings$iNumPosteriorSamples) #Calculate randomization probs based on posterior distribution
## Apply any probability temperation function, e.g. sqrt
probstemp <- StudyObj$StudyDesignSettings$probTemperation(probs)
## Save finala analyses output
FinalAnalyses<-list()
FinalAnalyses$Level = l
FinalAnalyses$Time = StudyObj$CurrentTime
FinalAnalyses$LMEFit<-lmefit
FinalAnalyses$UnWeightedUpdateProbabilities<-probs
FinalAnalyses$UnWeightedTemperatedUpdateProbabilities<-probstemp
FinalAnalyses$LMECoeff<-lmecoef
FinalAnalyses$LMESE<-lmese
StudyObj$FinalAnalysesList[[length(StudyObj$FinalAnalysesList)+1]]<-FinalAnalyses
## Also add probabilities and coeffs on the StudyObj level
StudyObj$FinalUnWeightedUpdateProbabilities<-probs
StudyObj$LMECoeff<-lmecoef
}
}
return(StudyObj)
}
StudyObjIni <- createStudy(
nCohorts = 2,
recruitmentAges = list(c(6,7)*30,c(12,13)*30),
nSubjects = c(500,500),
cohortStartTimes = c(0*30,0*30),
newCohortLink = list(NULL,NULL),
Recruitmentfunction = function(...) {return(5000)},
samplingDesign = list(seq(0,12,by=2)*30,seq(0,6,by=2)*30),
studyStopTime = 18*30+5,
latestTimeForNewBirthCohorts = 14*30+1,
treatments = list(c("SoC-1","Cell 1","Cell 2"," Cell 3"," Cell 4"),c("SoC-1","Cell 1","Cell 2"," Cell 3"," Cell 4")),
effSizes = list(c(0,0.0633,0.1037,0.1574,0.1687),c(0,0.0633,0.1037,0.1574,0.1687)),
randomizationProbabilities = list(rep(0.20,5),rep(0.20,5)),
minAllocationProbabilities = list(c(0.2,rep(0,4)),c(0.2,rep(0,4))),
probTemperationFunction = probTemperation,
AddNewBirthCohortEventFunction = AddNewSixMonthCohortEvent,
interimAnalyzesTimeFunction = InterimAnalyzesTime,
accumulatedData = TRUE
)
StudyObjIni$EventList[[length(StudyObjIni$EventList)+1]]<-AddNewTwelveMonthCohortEvent
StudyObjIni$EventList[[length(StudyObjIni$EventList)+1]]<-FinalAnalysesEvent
StudyObj <- AdaptiveStudy(StudyObjIni)
p <- plotStudyCohorts(StudyObj,plotAnaTimes = T,shiftWithinLevel = 1)
p + scale_y_continuous(breaks=c(7.5,13.5),labels=c("6-18 months","12-18 months")) +
ylab(NULL) +
theme(legend.position="none")
Results
The simulated HAZ data.
allData <- getAllSubjectData(StudyObj,scalarItems="AgeAtRand",covariates=NULL)
allData$Cohort2 <- as.character(allData$Cohort)
allData$Cohort2 <- factor(allData$Cohort2,levels=c("Cohort 2","Cohort 4","Cohort 6","Cohort 1","Cohort 3","Cohort 5"))
plotHAZ(StudyObj,data=allData) + aes(color=Cohort) + facet_wrap(~Cohort2) + theme(legend.position="none") + scale_x_continuous(breaks=c(6,12,18))
The randomization probabilties before* adjusting for the minimum SoC randomization probability.
plotProbs(StudyObj,cohortAgeNames=c("6-18 months","12-18 months"),strProb="UnWeightedRandomizationProbabilities")
The randomization probabilties after adjusting for the minimum SoC randomization probability.
plotProbs(StudyObj,cohortAgeNames=c("6-18 months","12-18 months"),strProb="RandomizationProbabilities")
The final analysis of all the data gives the following probabilities:
finalProbDf <- data.frame(Treatment=StudyObj$StudyDesignSettings$Treatments[[1]],
Col1= StudyObj$FinalAnalysesList[[1]]$UnWeightedUpdateProbabilities,
Col2= StudyObj$FinalAnalysesList[[2]]$UnWeightedUpdateProbabilities)
kable(finalProbDf,dig=2,col.names=c("Treatment","Prob 6-18 months","Prob 12-18 months"))
The final analysis of all the data gives the following estimated treatment effects (change in HAZ/6 months):
finalCoefDf <- data.frame(Treatment=(StudyObj$StudyDesignSettings$Treatments[[1]])[-1],
Col1= StudyObj$FinalAnalysesList[[1]]$LMECoeff/2,
Col2= StudyObj$FinalAnalysesList[[2]]$LMECoeff/2)
kable(finalCoefDf,dig=2,col.names=c("Treatment","Eff 6-18 months","Eff 12-18 months"),row.names=FALSE)
Repeating the study multiple times
iter <- 2
ncores <- 7
myMultStud <- runMultiSim(StudyObjIni,iter=iter,ncores=ncores,cohortAgeNames=c("6-18 months","12-18 months"),strProb="UnWeightedRandomizationProbabilities")
probDfUnWeightedRandProb <- myMultStud[[2]]
probDfRandProb <- getMultiProbList(myMultStud[[1]],cohortAgeNames=c("6-18 months","12-18 months"),strProb="RandomizationProbabilities",ncore=ncores)
Running the study multiple time (in this case r iter times) allows for the computation of performance statistics.
The randomization probabilties after adjusting for the minimum SoC randomization probability. The shaded areas are the 95% prediction intervals.
plotMultiProb(probDfUnWeightedRandProb)
The randomization probabilties after adjusting for the minimum SoC randomization probability. The shaded areas are the 95% prediction intervals.
plotMultiProb(probDfRandProb)
Performance metrics
To compare design choices quantitatively it is useful with performance metrics. Two suggestions are the randomization probabilities based on an analysis of all data (i.e. at the end of the trial) and the estimated treatment effects based on all data.
trtNams1 <- c("SoC-1","Cell1", "Cell2","Cell3","Cell4")
nStud <- length(myMultStud[[1]])
dfProbs1 <- data.frame(matrix(0,nrow=nStud,ncol=length(trtNams1)))
dfCoeff1 <- dfProbs1
dfProbs2 <- data.frame(matrix(0,nrow=nStud,ncol=length(trtNams1)))
dfCoeff2 <- dfProbs2
for(i in 1:length(myMultStud[[1]])) {
dfProbs1[i,] <- myMultStud[[1]][[i]]$FinalAnalysesList[[1]]$UnWeightedUpdateProbabilities
dfProbs2[i,] <- myMultStud[[1]][[i]]$FinalAnalysesList[[2]]$UnWeightedUpdateProbabilities
dfCoeff1[i,] <- c(0,myMultStud[[1]][[i]]$FinalAnalysesList[[1]]$LMECoeff)
dfCoeff2[i,] <- c(0,myMultStud[[1]][[i]]$FinalAnalysesList[[2]]$LMECoeff)
}
names(dfProbs1) <- trtNams1
names(dfProbs2) <- trtNams1
names(dfCoeff1) <- trtNams1
names(dfCoeff2) <- trtNams1
dfProbs1$AgeGroup <- "6-18 months"
dfProbs2$AgeGroup <- "12-18 months"
dfCoeff1$AgeGroup <- "6-18 months"
dfCoeff2$AgeGroup <- "12-18 months"
dfProbs <- rbind(dfProbs1,dfProbs2)
dfCoeff <- rbind(dfCoeff1,dfCoeff2)
sumFinProb <-
dfProbs %>% gather("Treatment","Probability",-AgeGroup) %>%
mutate(Treatment = factor(Treatment,levels=trtNams1),AgeGroup=factor(AgeGroup,levels=c("6-18 months","12-18 months"))) %>%
group_by(Treatment,AgeGroup) %>%
summarise(
Mean = mean(Probability),
Low=quantile(Probability,p=0.025),
High=quantile(Probability,p=0.975))
sumFinProbPow <- dfProbs
for(i in 1:nrow(sumFinProbPow)) {
sumFinProbPow$Pow1[i] <- ifelse(sumFinProbPow$Cell1[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0)
sumFinProbPow$Pow2[i] <- ifelse(sumFinProbPow$Cell2[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0)
sumFinProbPow$Pow3[i] <- ifelse(sumFinProbPow$Cell3[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0)
sumFinProbPow$Pow4[i] <- ifelse(sumFinProbPow$Cell4[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0)
}
sumFinProbPow2 <- sumFinProbPow %>%
select(AgeGroup,Pow1:Pow4) %>%
dplyr::rename(Cell1=Pow1,Cell2=Pow2,Cell3=Pow3,Cell4=Pow4) %>%
gather("Treatment","value",-AgeGroup) %>%
group_by(Treatment,AgeGroup) %>%
summarise(Power=sum(value/n())) %>%
filter(Treatment=="Cell4")
sumFinCoeff <-
dfCoeff %>% gather("Treatment","Coeff",-AgeGroup) %>%
mutate(Treatment = factor(Treatment,levels=trtNams1),Effect = Coeff/2,AgeGroup=factor(AgeGroup,levels=c("6-18 months","12-18 months"))) %>%
group_by(Treatment,AgeGroup) %>%
summarise(
Mean = mean(Effect),
Low = Mean - 1.96*sd(Effect)/sqrt(n()),
High = Mean + 1.96*sd(Effect)/sqrt(n()))
The "power" for Cell4 to be the best treatment in the two age groups are:
kable(sumFinProbPow2,dig=2)
The randomization probabilities (the bars are 95% prrdiction intervals):
ggplot(sumFinProb,aes(Treatment,Mean,color=Treatment)) +
geom_point() +
theme(legend.position="none") +
geom_errorbar(aes(ymin=Low,ymax=High),width=0.2) +
facet_wrap(~AgeGroup)
The treatment effects/6 months (the black dots are the true effects an the bars are 95% CI).
trueEff <- data.frame(Treatment=trtNams1,Mean=c(0,0.0633,0.1037,0.1574,0.1687))
ggplot(sumFinCoeff,aes(Treatment,Mean,color=Treatment)) +
geom_point() +
geom_point(data=trueEff,color="black") +
theme(legend.position="none") +
geom_errorbar(aes(ymin=Low,ymax=High),width=0.2) +
facet_wrap(~AgeGroup)
eniclas/FAIRsimulator documentation built on May 16, 2019, 5:12 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Model 1
library(FAIRsimulator) set.seed(3243) theme_set(theme_gray(base_size=14))
Treatments
- Cell 0 (SoC)
- Counseling
- 6 month HAZ effect: 0
- Cell 1
- Counseling
- Food: 7.5% PER, 10% ASP
- Iron + Vitamin A + Zinc
- 6 month HAZ effect: 0.0633
- Cell 2
- Counseling
- Food: 15% PER, 30% ASP
- Multiple micronutrient supplementation
- 6 month HAZ effect: 0.1037
- Cell 3
- Counseling
- Food: 15% PER, 30% ASP
- Multiple micronutrient supplementation
- Antibiotics + Probiotics
- 6 month HAZ effect: 0.1574
- Cell 4
- Counseling
- Food: 15% PER, 30% ASP
- Multiple micronutrient supplementation
- Antibiotics + Probiotics
- WASH Intensified
- Psychosocial care and support
- 6 month HAZ effect: 0.1687
Study design settings
These study design options were used for the current set of simulations:
- 3 recruitment cycles: 0, 6, and 12 months (study time).
- 2 recruitment age groups: 6 and 12 months (age).
- Instantaneous recruitment (all children recruited on day 1).
- The initial randomization set to 1:1:1:1:1.
- The minimum randomization probability for SoC is 20%.
- Treatment arms with a probability less than 10% before adjusting for the minimum SoC randomization probability is dropped.
- All available data up and including the timepoint for the interim analyses is used to update the probabilities.
- The dropout rate is approximately 20%/6 months.
- Bi-monthly HAZ observations.
- Interim analyses every 6 months.
- Square root probability updates (small probabilities are inflated)
## One adapted cohort AddNewSixMonthCohortEvent<-function(StudyObj) { #Note - assuming birth cohort is first element in StudyDesignSettings if (!is.null(StudyObj$CohortList)) { for (i in 1:length(StudyObj$CohortList)) { Cohort<-StudyObj$CohortList[[i]] if ((Cohort$Active==FALSE && Cohort$RandomizationAgeRange[1]==6*30 && Cohort$CycleNum==1 && Cohort$CohortStartTime+Cohort$CurrentTime==StudyObj$CurrentTime && StudyObj$CurrentTime<=StudyObj$StudyDesignSettings$LatestTimeForNewBirthCohorts) || #If this cohort just ended and its time to add a new cohort (Cohort$Active==TRUE && Cohort$RandomizationAgeRange[1]==6*30 && Cohort$CycleNum==1 && Cohort$CohortStartTime+6*30==StudyObj$CurrentTime && StudyObj$CurrentTime<=StudyObj$StudyDesignSettings$LatestTimeForNewBirthCohorts)) { BirthCohortIndex<-which(unlist(lapply(StudyObj$StudyDesignSettings$CohortAgeRange,function(x) {x[1]==6*30}))==TRUE) #Get the birth cohort index (i.e. lower AgeRangeAtRandomization=0) NewChildCohort<-NewCohort(StudyObj,CohortNum=BirthCohortIndex) NewChildCohort$StartNum <- max(StudyObj$CohortList %listmap% "StartNum")+1 NewChildCohort$Name<-paste0("C-",NewChildCohort$StartNum,"-",1," [",Cohort$RandomizationAgeRange[1]/30,"-",Cohort$RandomizationAgeRange[2]/30,"m @ rand]") DebugPrint(paste0("Create new 6 month cohort ",NewChildCohort$Name," based on probabilities in ",Cohort$Name," at time: ",StudyObj$CurrentTime),1,StudyObj) NewChildCohort$RandomizationProbabilities<-Cohort$UpdateProbabilities #Use the latest probabilities from previous birth cohort NewChildCohort$UnWeightedRandomizationProbabilities<-Cohort$UnWeightedUpdateProbabilities StudyObj$CohortList[[(length(StudyObj$CohortList)+1)]]<-NewChildCohort break } } } return(StudyObj) } ## One adapted cohort AddNewTwelveMonthCohortEvent<-function(StudyObj) { #Note - assuming birth cohort is first element in StudyDesignSettings if (!is.null(StudyObj$CohortList)) { for (i in 1:length(StudyObj$CohortList)) { Cohort<-StudyObj$CohortList[[i]] if (Cohort$Active==FALSE && Cohort$RandomizationAgeRange[1]==12*30 && Cohort$CycleNum==1 && Cohort$CohortStartTime+Cohort$CurrentTime==StudyObj$CurrentTime && StudyObj$CurrentTime<=StudyObj$StudyDesignSettings$LatestTimeForNewBirthCohorts) #if this cohort just ended and its time to add a new cohort { BirthCohortIndex<-which(unlist(lapply(StudyObj$StudyDesignSettings$CohortAgeRange,function(x) {x[1]==12*30}))==TRUE) #Get the birth cohort index (i.e. lower AgeRangeAtRandomization=0) NewChildCohort<-NewCohort(StudyObj,CohortNum=BirthCohortIndex) NewChildCohort$StartNum <- max(StudyObj$CohortList %listmap% "StartNum")+1 NewChildCohort$Name<-paste0("C-",NewChildCohort$StartNum,"-",1," [",Cohort$RandomizationAgeRange[1]/30,"-",Cohort$RandomizationAgeRange[2]/30,"m @ rand]") DebugPrint(paste0("Create new 12 month cohort ",NewChildCohort$Name," based on probabilities in ",Cohort$Name," at time: ",StudyObj$CurrentTime),1,StudyObj) NewChildCohort$RandomizationProbabilities<-Cohort$UpdateProbabilities #Use the latest probabilities from previous birth cohort NewChildCohort$UnWeightedRandomizationProbabilities<-Cohort$UnWeightedUpdateProbabilities StudyObj$CohortList[[(length(StudyObj$CohortList)+1)]]<-NewChildCohort break } } } return(StudyObj) } probTemperation <- function(probs) { probs <- sqrt(probs)/sum(sqrt(probs)) return(probs) } InterimAnalyzesTime<-function(Cohort,StudyObj) { DebugPrint(paste0("Check if cohort ",Cohort$Name," is about to have an interim analysis",StudyObj$CurrentTime),3,StudyObj) TimeToPerformInterim<-FALSE tmp<-sum(unlist(lapply(Cohort$SubjectList,function(Subject,LastSample){ return(as.numeric(Subject$Status==0 || (LastSample %in% Subject$SubjectSampleTime)))},max(Cohort$SamplingDesign))))/Cohort$MaxNumberOfSubjects if (length(tmp)==0) tmp<-0 if (Cohort$RandomizationAgeRange[1]==6*30 && (StudyObj$CurrentTime==6*30-1 || StudyObj$CurrentTime==18*30-1 )) { DebugPrint(paste0("Time to do an interim analyses for cohort ",Cohort$Name," at time: ",StudyObj$CurrentTime," (",round(tmp*100,1)," % subjects completed)"),1,StudyObj) TimeToPerformInterim<-TRUE } return(TimeToPerformInterim) } FinalAnalysesEvent<-function(StudyObj) { DebugPrint(paste0("Check if time to do a final analyses at study time: ",StudyObj$CurrentTime),3,StudyObj) if (StudyObj$CurrentTime==StudyObj$StudyDesignSettings$StudyStopTime-1) { DebugPrint(paste0("Time to do a final analyses with all data at time: ",StudyObj$CurrentTime),1,StudyObj) df<-getAllSubjectData(StudyObj,prevTreatment = TRUE) #Get all data that is neede for final analyses ## Rename, select and order columns if(length(names(df)[grep(names(df),pattern = "PTRT")]) >0) { df <- df %>% rename(ID=StudyID,DATA=HAZ,AGE=Age,TRT=TreatmentIndex,TRTS=Treatment) %>% select(ID,DATA,AGE,TRT,TRTS,one_of(StudyObj$StudyDesignSettings$Covariates),Level,CohortName,matches("PTRT")) for (cstr in names(df)[grep(names(df),pattern = "PTRT")]) { #Set no previous treatment to "0" df[,cstr]<-ifelse(is.na(df[,cstr]),0,df[,cstr]) } } else { df <- df %>% rename(ID=StudyID,DATA=HAZ,AGE=Age,TRT=TreatmentIndex,TRTS=Treatment) %>% select(ID,DATA,AGE,TRT,TRTS,one_of(StudyObj$StudyDesignSettings$Covariates),Level,CohortName) } df[df==-99]<-NA #Set -99 to missing df<-ImputeCovariates(df,StudyObj,method=StudyObj$StudyDesignSettings$ImpMethod) #Impute missing covariates #Sort to get correct TRT order df <- df[order(df["TRT"],df["ID"], df["AGE"]),] ##### Make some covariate factors df$TRT<-as.factor(df$TRT) df$SEXN<-as.factor(df$SEXN) df$SANITATN<-as.factor(df$SANITATN) df$AGE<-df$AGE/(12*30) #Rescale time to years ### If we have previous treatments ptrti <- grep("PTRT",names(df)) if(length(ptrti)>0) { for(i in 1:length(ptrti)) { df[,ptrti[i]] <- as.factor(df[,ptrti[i]]) } myPTRTs <- names(df)[ptrti] } cohortlevels<-unique(StudyObj$CohortList %listmap% "Level") for (l in 1:length(cohortlevels)) { dflevel<-subset(df,Level==cohortlevels[l]) #### Perform LME estimation based on some covariates and treatment effects for each level if(length(ptrti)>0) { lmeFormula <- paste0("DATA~1 + AGE + AGE:TRT + (AGE|ID) +",paste0(c(StudyObj$StudyDesignSettings$Covariates,myPTRTs),collapse = " + ")) } else { lmeFormula <- paste0("DATA~1 + AGE + AGE:TRT + (AGE|ID) +",paste0(StudyObj$StudyDesignSettings$Covariates,collapse = " + ")) } lmefit <- lmer(lmeFormula,data=dflevel,REML=FALSE) # IIV on baseline only ##### Calculate new probabilites based on another cohort LME results lmecoef<-summary(lmefit)$coefficients[,1] #Get coefficicents from LME lmese<-summary(lmefit)$coefficients[,2] #Get SE from LME lmecoef<-lmecoef[regexpr('AGE:TRT.*',names(lmecoef))==1] lmese<-lmese[regexpr('AGE:TRT.*',names(lmese))==1] print(lmefit) ### Get first cohort with level = l if (!is.null(StudyObj$CohortList)) { for (i in 1:length(StudyObj$CohortList)) { if (StudyObj$CohortList[[i]]$Level==l) { Cohort<-StudyObj$CohortList[[i]] break } } if ((length(lmecoef)+1)!=length(Cohort$RandomizationProbabilities)) { lmecoefnew<-rep(0,length(Cohort$RandomizationProbabilities)-1) lmesenew<-rep(0,length(Cohort$RandomizationProbabilities)-1) for (i in 1:(length(Cohort$RandomizationProbabilities)-1)) { iIndex<-which(names(lmecoef)==paste0("AGE:TRT",i+1)) if (length(iIndex)!=0) { lmecoefnew[i]<-lmecoef[iIndex] lmesenew[i]<-lmese[iIndex] } } names(lmecoefnew)<-paste0("AGE:TRT",2:length(Cohort$RandomizationProbabilities)) names(lmesenew)<-paste0("AGE:TRT",2:length(Cohort$RandomizationProbabilities)) lmecoef<-lmecoefnew lmese<-lmesenew } } #Get probability of beeing best probs <- GetNewRandomizationProbabilities(trtcoeff=lmecoef,trtse=lmese, StudyObj$StudyDesignSettings$iNumPosteriorSamples) #Calculate randomization probs based on posterior distribution ## Apply any probability temperation function, e.g. sqrt probstemp <- StudyObj$StudyDesignSettings$probTemperation(probs) ## Save finala analyses output FinalAnalyses<-list() FinalAnalyses$Level = l FinalAnalyses$Time = StudyObj$CurrentTime FinalAnalyses$LMEFit<-lmefit FinalAnalyses$UnWeightedUpdateProbabilities<-probs FinalAnalyses$UnWeightedTemperatedUpdateProbabilities<-probstemp FinalAnalyses$LMECoeff<-lmecoef FinalAnalyses$LMESE<-lmese StudyObj$FinalAnalysesList[[length(StudyObj$FinalAnalysesList)+1]]<-FinalAnalyses ## Also add probabilities and coeffs on the StudyObj level StudyObj$FinalUnWeightedUpdateProbabilities<-probs StudyObj$LMECoeff<-lmecoef } } return(StudyObj) } StudyObjIni <- createStudy( nCohorts = 2, recruitmentAges = list(c(6,7)*30,c(12,13)*30), nSubjects = c(500,500), cohortStartTimes = c(0*30,0*30), newCohortLink = list(NULL,NULL), Recruitmentfunction = function(...) {return(5000)}, samplingDesign = list(seq(0,12,by=2)*30,seq(0,6,by=2)*30), studyStopTime = 18*30+5, latestTimeForNewBirthCohorts = 14*30+1, treatments = list(c("SoC-1","Cell 1","Cell 2"," Cell 3"," Cell 4"),c("SoC-1","Cell 1","Cell 2"," Cell 3"," Cell 4")), effSizes = list(c(0,0.0633,0.1037,0.1574,0.1687),c(0,0.0633,0.1037,0.1574,0.1687)), randomizationProbabilities = list(rep(0.20,5),rep(0.20,5)), minAllocationProbabilities = list(c(0.2,rep(0,4)),c(0.2,rep(0,4))), probTemperationFunction = probTemperation, AddNewBirthCohortEventFunction = AddNewSixMonthCohortEvent, interimAnalyzesTimeFunction = InterimAnalyzesTime, accumulatedData = TRUE ) StudyObjIni$EventList[[length(StudyObjIni$EventList)+1]]<-AddNewTwelveMonthCohortEvent StudyObjIni$EventList[[length(StudyObjIni$EventList)+1]]<-FinalAnalysesEvent StudyObj <- AdaptiveStudy(StudyObjIni)
p <- plotStudyCohorts(StudyObj,plotAnaTimes = T,shiftWithinLevel = 1) p + scale_y_continuous(breaks=c(7.5,13.5),labels=c("6-18 months","12-18 months")) + ylab(NULL) + theme(legend.position="none")
Results
The simulated HAZ data.
allData <- getAllSubjectData(StudyObj,scalarItems="AgeAtRand",covariates=NULL) allData$Cohort2 <- as.character(allData$Cohort) allData$Cohort2 <- factor(allData$Cohort2,levels=c("Cohort 2","Cohort 4","Cohort 6","Cohort 1","Cohort 3","Cohort 5")) plotHAZ(StudyObj,data=allData) + aes(color=Cohort) + facet_wrap(~Cohort2) + theme(legend.position="none") + scale_x_continuous(breaks=c(6,12,18))
The randomization probabilties before* adjusting for the minimum SoC randomization probability.
plotProbs(StudyObj,cohortAgeNames=c("6-18 months","12-18 months"),strProb="UnWeightedRandomizationProbabilities")
The randomization probabilties after adjusting for the minimum SoC randomization probability.
plotProbs(StudyObj,cohortAgeNames=c("6-18 months","12-18 months"),strProb="RandomizationProbabilities")
The final analysis of all the data gives the following probabilities:
finalProbDf <- data.frame(Treatment=StudyObj$StudyDesignSettings$Treatments[[1]], Col1= StudyObj$FinalAnalysesList[[1]]$UnWeightedUpdateProbabilities, Col2= StudyObj$FinalAnalysesList[[2]]$UnWeightedUpdateProbabilities) kable(finalProbDf,dig=2,col.names=c("Treatment","Prob 6-18 months","Prob 12-18 months"))
The final analysis of all the data gives the following estimated treatment effects (change in HAZ/6 months):
finalCoefDf <- data.frame(Treatment=(StudyObj$StudyDesignSettings$Treatments[[1]])[-1], Col1= StudyObj$FinalAnalysesList[[1]]$LMECoeff/2, Col2= StudyObj$FinalAnalysesList[[2]]$LMECoeff/2) kable(finalCoefDf,dig=2,col.names=c("Treatment","Eff 6-18 months","Eff 12-18 months"),row.names=FALSE)
Repeating the study multiple times
iter <- 2 ncores <- 7 myMultStud <- runMultiSim(StudyObjIni,iter=iter,ncores=ncores,cohortAgeNames=c("6-18 months","12-18 months"),strProb="UnWeightedRandomizationProbabilities")
probDfUnWeightedRandProb <- myMultStud[[2]] probDfRandProb <- getMultiProbList(myMultStud[[1]],cohortAgeNames=c("6-18 months","12-18 months"),strProb="RandomizationProbabilities",ncore=ncores)
Running the study multiple time (in this case r iter times) allows for the computation of performance statistics.
The randomization probabilties after adjusting for the minimum SoC randomization probability. The shaded areas are the 95% prediction intervals.
plotMultiProb(probDfUnWeightedRandProb)
The randomization probabilties after adjusting for the minimum SoC randomization probability. The shaded areas are the 95% prediction intervals.
plotMultiProb(probDfRandProb)
Performance metrics
To compare design choices quantitatively it is useful with performance metrics. Two suggestions are the randomization probabilities based on an analysis of all data (i.e. at the end of the trial) and the estimated treatment effects based on all data.
trtNams1 <- c("SoC-1","Cell1", "Cell2","Cell3","Cell4") nStud <- length(myMultStud[[1]]) dfProbs1 <- data.frame(matrix(0,nrow=nStud,ncol=length(trtNams1))) dfCoeff1 <- dfProbs1 dfProbs2 <- data.frame(matrix(0,nrow=nStud,ncol=length(trtNams1))) dfCoeff2 <- dfProbs2 for(i in 1:length(myMultStud[[1]])) { dfProbs1[i,] <- myMultStud[[1]][[i]]$FinalAnalysesList[[1]]$UnWeightedUpdateProbabilities dfProbs2[i,] <- myMultStud[[1]][[i]]$FinalAnalysesList[[2]]$UnWeightedUpdateProbabilities dfCoeff1[i,] <- c(0,myMultStud[[1]][[i]]$FinalAnalysesList[[1]]$LMECoeff) dfCoeff2[i,] <- c(0,myMultStud[[1]][[i]]$FinalAnalysesList[[2]]$LMECoeff) } names(dfProbs1) <- trtNams1 names(dfProbs2) <- trtNams1 names(dfCoeff1) <- trtNams1 names(dfCoeff2) <- trtNams1 dfProbs1$AgeGroup <- "6-18 months" dfProbs2$AgeGroup <- "12-18 months" dfCoeff1$AgeGroup <- "6-18 months" dfCoeff2$AgeGroup <- "12-18 months" dfProbs <- rbind(dfProbs1,dfProbs2) dfCoeff <- rbind(dfCoeff1,dfCoeff2) sumFinProb <- dfProbs %>% gather("Treatment","Probability",-AgeGroup) %>% mutate(Treatment = factor(Treatment,levels=trtNams1),AgeGroup=factor(AgeGroup,levels=c("6-18 months","12-18 months"))) %>% group_by(Treatment,AgeGroup) %>% summarise( Mean = mean(Probability), Low=quantile(Probability,p=0.025), High=quantile(Probability,p=0.975)) sumFinProbPow <- dfProbs for(i in 1:nrow(sumFinProbPow)) { sumFinProbPow$Pow1[i] <- ifelse(sumFinProbPow$Cell1[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0) sumFinProbPow$Pow2[i] <- ifelse(sumFinProbPow$Cell2[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0) sumFinProbPow$Pow3[i] <- ifelse(sumFinProbPow$Cell3[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0) sumFinProbPow$Pow4[i] <- ifelse(sumFinProbPow$Cell4[i] == max(sumFinProbPow$Cell1[i],sumFinProbPow$Cell2[i],sumFinProbPow$Cell3[i],sumFinProbPow$Cell4[i]),1,0) } sumFinProbPow2 <- sumFinProbPow %>% select(AgeGroup,Pow1:Pow4) %>% dplyr::rename(Cell1=Pow1,Cell2=Pow2,Cell3=Pow3,Cell4=Pow4) %>% gather("Treatment","value",-AgeGroup) %>% group_by(Treatment,AgeGroup) %>% summarise(Power=sum(value/n())) %>% filter(Treatment=="Cell4") sumFinCoeff <- dfCoeff %>% gather("Treatment","Coeff",-AgeGroup) %>% mutate(Treatment = factor(Treatment,levels=trtNams1),Effect = Coeff/2,AgeGroup=factor(AgeGroup,levels=c("6-18 months","12-18 months"))) %>% group_by(Treatment,AgeGroup) %>% summarise( Mean = mean(Effect), Low = Mean - 1.96*sd(Effect)/sqrt(n()), High = Mean + 1.96*sd(Effect)/sqrt(n()))
The "power" for Cell4 to be the best treatment in the two age groups are:
kable(sumFinProbPow2,dig=2)
The randomization probabilities (the bars are 95% prrdiction intervals):
ggplot(sumFinProb,aes(Treatment,Mean,color=Treatment)) + geom_point() + theme(legend.position="none") + geom_errorbar(aes(ymin=Low,ymax=High),width=0.2) + facet_wrap(~AgeGroup)
The treatment effects/6 months (the black dots are the true effects an the bars are 95% CI).
trueEff <- data.frame(Treatment=trtNams1,Mean=c(0,0.0633,0.1037,0.1574,0.1687)) ggplot(sumFinCoeff,aes(Treatment,Mean,color=Treatment)) + geom_point() + geom_point(data=trueEff,color="black") + theme(legend.position="none") + geom_errorbar(aes(ymin=Low,ymax=High),width=0.2) + facet_wrap(~AgeGroup)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.