## Simulation Code for Example 2 in McGowan, Nix, Murphy, Bierman, & CPPRG ## Written by Herle McGowan ## We set the standardized values of the parameters and then iterate to find the corresponding unstandardized values ## The unstandardized values are then used to genertate variables that, when standardized, would return the standardized parameters we set ## Parameters are standardized according to the following relations (see endnote 6 in the paper cited above): # For parameters of a continuous predictor of a continuous response (e.g. delta, phi, beta): standardized parameter = unstandardized parameter * (sd(predictor) / sd(response)) # For parameters of a continuous predictor of a binary response (e.g. gamma): standardized parameter = unstandardized parameter * sd(predictor) # Parameters of a binary predictor of a binary response (e.g. theta) are not standardized ################################################# ## Initial set-up and assignments of variables ## ################################################# # Inverse logit (expit) function written by Julian Faraway # Available in library(faraway) ilogit <- function (x) { if (any(omit <- is.na(x))) { lv <- x lv[omit] <- NA if (any(!omit)) lv[!omit] <- Recall(x[!omit]) return(lv) } exp(x)/(1 + exp(x)) } # Set seed so that results can be reproduced set.seed(1023) # Set number of time points to be examined n.time <- 22 # Set sample size to be used for finding unstandardize parameters n <- 10000 # Intitialize matricies to store results of variable generation c.temp<- matrix(0, nrow=n, ncol=(n.time+1)) p.temp<- matrix(0, nrow=n, ncol=(n.time+1)) d.temp<- matrix(0, nrow=n, ncol=(n.time+1)) # Define vectors to store results of unstandardized parameters gamma <- rep(0, n.time) # parameters from confounder to dose delta <- rep(0, n.time) # parameters from confounder to confounder # Define vectors to store set values of standardized parameters for each simulation ## Note that only the code for Simulation A when the size of the counfounding is -0.02 will run with the current parameter values ## To obtain results for other simulations, change the parameter values according to Table 1 in the paper sigma <- 1 # Standard deviation of random error for each variable s.gamma <- rep(c(0,0.14), c(1,n.time)) # Standardized parameters from confounder to dose s.delta <- rep(c(0,1), c(1,n.time)) # Standardized parameters from confounder to confounder s.theta <- 1 # Standardized parameters from dose to dose s.phi <- -0.14 # Standardized parameters from confounder to s.beta <- 0 # Standardized parameters from dose to y ########################################################## ## Iterate to calcuate unstandardized parameter values ## ## According to formula given in header ## ########################################################## # Generate initial pretreatment confounder c.temp[,1] <- rnorm(n, 0, sigma) # Calculate unstandardized parameter from confounder to dose gamma[1] <- s.gamma[1] / sd(c.temp[,1]) # Generate initial first dose p.temp[,2] <- ilogit(gamma[1]*c.temp[,1]) d.temp[,2] <- rbinom(n, 1, p.temp[,2]) # Calculate initial unstandardized parameter from confounder to confounder delta[1] <- s.delta[1] * sigma / sd(c.temp[,1]) # Iterate to refine unstandardized parameter delta repeat{ # Generate initial confounder at time 1 c.temp[,2] <- delta[1]*c.temp[,1] + rnorm(n,0,sigma) # See what standardized parameters would be based on this data delta01.std <- delta[1] * sd(c.temp[,1]) / sd(c.temp[,2]) # Check if these are what we want standardized parameters to be if((delta01.std - s.delta[1])^2 < 1e-4) { break } # If not, update unstandardized parameters and continue delta[1] <- s.delta[1] * sd(c.temp[,2]) / sd(c.temp[,1]) } # End repeat # Loop to generate initial variables for time 2 to end for(i in 3:(n.time+1)){ # Calculate unstandardized parameter from confounder to dose gamma[(i-1)] <- s.gamma[(i-1)] / sd(c.temp[,(i-1)]) # Set unstandardized parameter from dose to dose theta <- s.theta # Generate initial dose at time i p.temp[,i] <- ilogit(gamma[(i-1)]*c.temp[,(i-1)] + theta*d.temp[,(i-1)]) d.temp[,i] <- rbinom(n, 1, p.temp[,i]) # Calculate initial unstandardized parameters from confounder to confounder delta[(i-1)] <- s.delta[(i-1)] * sigma / sd(c.temp[,(i-1)]) # Iterate to refine initial unstandardized delta repeat{ # Generate confounder at time i c.temp[,i] <- delta[(i-1)]*c.temp[,(i-1)] + rnorm(n,0,sigma) # See what standardized parameter would be based on this data delta.std <- delta[(i-1)] * sd(c.temp[,(i-1)]) / sd(c.temp[,i]) # Check if this is what we want standardized parameter to be if((delta.std - s.delta[(i-1)])^2 < 1e-4) {break} # If not, update unstandardized parameter and continue delta[(i-1)] <- s.delta[(i-1)] * sd(c.temp[,i]) / sd(c.temp[,(i-1)]) } # End repeat } # End for loop # Calculate initial cummulative dose across all time points cum.d <- apply(d.temp[,-1], 1, sum) # Calculate initial unstandardized parameter confounder to y phi <- s.phi * sigma / sd(c.temp[,n.time]) # Calculate initial unstandardized parameter dose to y beta <- s.beta * sigma / sd(cum.d) # Iterate to refine initial unstandardized parameters phi and beta repeat{ # Generate y y <- phi*c.temp[,n.time] + beta*cum.d + rnorm(n,0,sigma) # See what standardized parameters would be based on this data phi.std <- phi * sd(c.temp[,n.time]) / sd(y) beta.std <- beta * sd(cum.d) / sd(y) # Check if these are what we want standardized parameters to be if(((phi.std - s.phi)^2 + (beta.std - s.beta)^2) < 1e-4) { break } # If not, update unstandardized parameters and continue phi <- s.phi * sd(y) / sd(c.temp[,n.time]) beta <- s.beta * sd(y) / sd(cum.d) } # End repeat ## End calcuate unstandardized parameter values ################################################## ## Generate data and Estimate regression models ## ################################################## # Set number of datasets to be created nn <- 1000 # Set sample size for each dataset n <- 400 # Intitialize matricies to store data conf <- matrix(0, nrow=n, ncol=(n.time+1)) prob <- matrix(0, nrow=n, ncol=(n.time+1)) dose <- matrix(0, nrow=n, ncol=(n.time+1)) # Define list to save regression results reg <- list() ## Loop to create multiple datasets and estimate regression models for(j in 1:nn){ ####################### ## Generate Datasets ## ####################### # Generate pretreatment confounder conf[,1] <- rnorm(n, 0, sigma) # Generate first dose prob[,2] <- ilogit(gamma[1]*conf[,1]) dose[,2] <- rbinom(n, 1, prob[,2]) # Generate confounder at time 1 conf[,2] <- delta[1]*conf[,1] + rnorm(n,0,sigma) # Loop to generate dose and confounder at remaining time points for(i in 3:(n.time+1)){ # Generate dose at time i prob[,i] <- ilogit(gamma[(i-1)]*conf[,(i-1)] + theta*dose[,(i-1)]) dose[,i] <- rbinom(n, 1, prob[,i]) # Generate confounder at time i conf[,i] <- delta[(i-1)]*conf[,(i-1)] + rnorm(n,0,sigma) } # End for loop # Calculate cummulative dose across all time points cum.d <- apply(dose[,-1], 1, sum) # Generate y y <- phi*conf[,n.time] + beta*cum.d + rnorm(n,0,sigma) ## End generate data ######################################## ## Estimate regression models ## ## Ignoring time-varying confounders ## ######################################## # Estimate linear regression of y on cum.d or regression of y on pretreatment confounder and cum.d, where appropriate ifelse(s.gamma[1]==0,fit <- lm(y~cum.d, x=T),fit <- lm(y~conf[,1]+cum.d, x=T)) # Save estimated betas and calculate se(beta), t-statistics beta.hat <- fit\$coeff se <- sqrt(diag(solve(t(fit\$x)%*%fit\$x)*summary(fit)\$sigma^2)) s.beta.hat <- beta.hat * sd(cum.d) / sd(y) t <- beta.hat / se # Save each estimated model reg[[j]] <- list(beta.hat=beta.hat,se=se,t=t,s.beta.hat=s.beta.hat) ## End estimate models ignoring time-varying confounders } ## End loop to create multiple datasets and estimate regression models ############################################# ## Access estimates from regression models ## ############################################# # Format saved results as dataframe for easier access all <- data.frame(t(mapply(as.matrix,lapply(reg,as.data.frame)))) # Access estimates when there is no effect of the pretreatment confounder (Simulations A and B) if(s.gamma[1]==0){ print("Average estimate of standardized beta") print("Calculated as the mean of standardized betas across 1000 data sets") print(round(mean(all[,8]),3)) print("-----") if(s.beta==0){ print("Type 1 Error Rate for Simulation A (true beta = 0)") print("Calculated as the proportion of |t| > 1.96 across 1000 data sets") print(mean(ifelse(abs(all[,6])>1.96,1,0))) } # End if if(s.beta!=0){ print("Power for Simulation B (true beta > 0)") print("Calculated as the proportion of t > 1.645 across 1000 data sets") print(mean(ifelse(all[,6]>1.645,1,0))) } # End if } # End if statement # Access estimates there is an effect of the pretreatment confounder (Simulation C) if(s.gamma[1]!=0){ print("Average estimate of standardized beta") print("Calculated as the mean of standardized betas across 1000 data sets") print(round(mean(all[,12]),3)) print("-----") if(s.beta==0){ print("Type 1 Error Rate for Simulation C (true beta = 0)") print("Calculated as the proportion of |t| > 1.96 across 1000 data sets") print(mean(ifelse(abs(all[,9])>1.96,1,0))) } # End if } # End if statement ## End access estimates