### Sample program to compute the SE, Confidence Interval and p-value for a paramter ### ### estimate based on calibrated log total energy intake, using a booststrap resampling ### ### procedure. This method correctly allows for the sampling variation in the ### ### calibration coefficient estimate. ### ### In this example, 'datin' is the analytic dataset. An initial step is not included ### ### in this code which is to compute the variable called 'lgffqenergyC', which is the ### ### calibrated estimate of the log FFQ total energy (i.e. the calibration equation was ### ### applied to each observation's log FFQ total energy. ### Note: this program is intended as an example only; modifications specific to each ### ### project are necessary. ### library(survival) library(gee) ### Define all macros ### #Macro to compute CI and p-value calCI<-function(Est,SE,alpha=0.05,lowbound=NULL,highbound=NULL){ level<-1-alpha/2 low<-Est-qnorm(level)*SE if (!missing(lowbound)) low<-ifelse(lowhighbound,highbound,high) pvalue<-2*(1-pnorm(abs(Est/SE))) return(list(Est=Est,low=low,high=high,pvalue=pvalue)) } #Macro to run Cox model coxreg<-function(datain,t2ev,ev,expvar){ modfit <-coxph(Surv(t2ev,ev) ~ expvar + ethnr + educrx + smokstat + brcarel1x + ealone + ep + whidt + strata(ageffq5y,hrtarm,osflag,whixstrat),data=datain,ties="breslow") est<-c(modfit$coef[1]) #return parameter estimate for exposure variable only return(est) } #Macro to carry out random sampling with replacement, separately in cases and non-cases getsamp<-function(datin,varCase){ dat<-datin datinca <- dat[which(varCase==1),] datinco <- dat[which(varCase==0),] caseid <- datinco$id contid <- datinca$id sampcaid <- sample(1:length(caseid),replace=T) sampcoid <- sample(1:length(contid),replace=T) newdat<- rbind(datinca[sampcaid,],datinco[sampcoid,]) return(newdat) } #Macro runs the calibration model on the nutrient biomarker dataset, saves the coefficients # and calculates calibrated total energy in analytic dataset. #Assumes the calibration equation has previously been defined calRCEnergy <- function(nutdat,anadat){ ## For this example, compute values to center covariates from nutrient biomarker dataset## Calories.M<-mean(nutdat$logffqenergy,na.rm=T) Bmi.M<-mean(nutdat$bmi,na.rm=T) Age.M<-mean(nutdat$agev1,na.rm=T) nutdat$logffqenergy.C<-nutdat$logffqenergy-Calories.M nutdat$bmi.C<-nutdat$bmi-Bmi.M nutdat$age.C<-nutdat$agev1-Age.M ## computed centered covariates in analytic dataset ## anadat$logcalories.FFQ.C<-anadat$lgf60enrgy-Calories.M anadat$bmi.C<-anadat$bmi-Bmi.M anadat$age.C<-anadat$age-Age.M ## run calibration equation and use coefficients to compute calibrated log FFQ energy in analytic dataset ## logfit.FFQ<-gee(logtee~logffqenergy.C +bmi.C+age.C+black+hispanic+othrace+inclt20+inc20to35+inc50to75+incge75+texpwk,data=nutdat) calibcoeff<-coef(logfit.FFQ) ## Calculate calibrated log total energy, called zhat.FFQ V<-cbind(rep(1,length(anadat$age.C)),anadat$logcalories.FFQ.C,anadat$bmi.C,anadata$age.C,anadat$black,anadat$hispanic,anadat$othrace,anadat$inclt20,anadat$inc20to35,anadat$inc50to75,anadat$incge75,anadat$texpwk) anadat$zhat.FFQ<-t(matrix(as.numeric(calibcoeff), nrow=1) %*% t(V)) return(anadat) } #read in analytic dataset - must include all variables used in calibration equation datin <- read.csv(file="S://...//R//dmcosffqS.csv") names(datin) #below specific to each project datin$hrtarm <-as.factor(datin$hrtarm) datin$ageffq5y <-as.factor(datin$ageffq5y) datin$smokstat <-as.factor(datin$smokstat) datin$brcarel1x <-as.factor(datin$brcarel1x) datin$ethnr <-as.factor(datin$ethnr) datin$educrx <-as.factor(datin$educrx) datin$whix<-as.factor(datin$whix) # Analytic sample combines DMC and OS cohorts - create separate datasets for use in resampling# datinos <- datin[datin$osflag==1,] datindm <- datin[datin$osflag==0,] #read in nutrient biomarker dataset (NBS or NPAAS) # in this example, the nutrient biomarker dataset has one record per participant nutdat<-read.csv(file="S://...//R//nbs_inv.csv",na.strings=c("NA","Missing","MISSING",""," ",".")) names(nutdat) #create log-transformed variables nutdat$logffqenergy<-log(nutdat$ffqenergy) nutdat$logtee<-log(nutdat$bioenergy_TEE) ### Bootstrap procedure - uses 10000 iterations ### ### Resampling takes place in both the nutrient biomarker and analytic datasets ### In this example, there are two Cox models of interest- one each for ER+ and ER- breast cancer set.seed(1) # set up output datasets outallerpos <- outallerneg <- array(NA,dim=c(10000,2)) for (k in 1:10000) { #resample with replacement nutrient biomarker dataset nut <- nutdat nut.id <- nut$id nutresamp.id <- sample(1:length(nut.id),replace=T) nut.boot <- nut[nutresamp.id,] # ER pos #resample analytic dataset, stratified by DM/OS and case/non-case newdatos1 <- getsamp(datinos,datinos$invbrcaerp) newdatdm1 <- getsamp(datindm,datindm$invbrcaerp) newdatall1 <- rbind(newdatdm1,newdatos1) newdatall1rc <- calRCEnergy(nut.boot,newdatall1) outallerpos[k,] <-coxreg(newdatall1rc,newdatall1rc$t2invbrca,newdatall1rc$invbrcaerp,newdatall1rc$zhat.FFQ) # ER neg newdatos2 <- getsamp(datinos,datinos$invbrcaern) newdatdm2 <- getsamp(datindm,datindm$invbrcaern) newdatall2 <- rbind(newdatdm2,newdatos2) newdatall2rc <- calRCEnergy(nut.boot,newdatall2) outallerneg[k,] <-coxreg(newdatall2rc,newdatall2rc$t2invbrca,newdatall2rc$invbrcaern,newdatall2rc$zhat.FFQ) ### include all other models run as part of this analysis ### } # Compute 95% CI and p-value for the Hazard Ratio based on calibrated total energy # This example uses a normal approximation method to get a measure of spread, # in this case the SE of the 10000 coefficients from the bootstrap models, # an alternative is the percentile method ## ER pos sd.boot.allerpos<-apply(outallerpos,2,sd) # get measure of spread # #get HR estimate from actual model # outcoxerpos <-coxreg(datin,datin$t2invbrca,datin$invbrcaerp,datin$lgffqenergyC) #get CI and pvalues based on bootstrap pval1 <- calCI(outcoxerpos,sd.boot.allerpos)$pvalue lowCI1 <-calCI(outcoxerpos,sd.boot.allerpos)$low hiCI1 <-calCI(outcoxerpos,sd.boot.allerpos)$high ResultAllerpos <- c(exp(outcoxerpos[1]),exp(lowCI1),exp(hiCI1),pval1) colnames(ResultAllerpos)<-c("HR","LowCI","UpCI","pval") write.csv(ResultAllerpos ,file="S://...//Results//dmcosBootCIERPOS.csv") ## ER neg sd.boot.allerpos<-apply(outallerneg,2,sd) # gets measure of spread # #get estimates from actual model - All DMC/OS # outcoxerneg <-coxreg(datin,datin$t2invbrca,datin$invbrcaern,datin$lgffqenergyC) #get CI and pvalues based on bootstrap pval2 <- calCI(outcoxerneg,sd.boot.allerneg)$pvalue lowCI2 <-calCI(outcoxerneg,sd.boot.allerneg)$low hiCI2 <-calCI(outcoxerneg,sd.boot.allerneg)$high ResultAllerneg <- c(exp(outcoxerneg[1]),exp(lowCI2),exp(hiCI2),pval2) colnames(ResultAllerneg)<-c("HR","LowCI","UpCI","pval") write.csv(ResultAllerneg ,file="S://...//Results//dmcosBootCIERNEG.csv")