Umediation<-function(n=100,Atype="D",Mtype="C",Ytype="C",Ctype="C",Utype="C",
interact=FALSE,muC=0,varC=1,muU=0,varU=1,
gamma0=0,gammaC=0,gammaU=0,varA=1,
alpha0=0,alphaA=0,alphaC=0,alphaU=0,varM=1,
beta0=0,betaA=0,betaM=0,betaI=0,betaC=0,betaU=0,varY=1,
alpha=0.05,nSim=250,nBoot=400,seed=1,atreat=1,acontrol=0){
library(mediation) #load the mediation library
library(car) #load the car library
#######################################
# Check Input for Errors
#######################################
ErrorCheck(n,Atype,Mtype,Ytype,Ctype,Utype,interact,muC,varC,muU,varU,gamma0,gammaC,gammaU,varA,alpha0,alphaA,alphaC,alphaU,varM,beta0,betaA,betaM,betaI,betaC,betaU,varY,alpha,nSim,nBoot,seed,atreat,acontrol)
#######################################
# Create Results Matrix
#######################################
Results<-matrix(0,nrow=12,ncol=1)
rownames(Results)<-c("Prop. of simulations w/ significant ACME excluding U","Prop. of simulations w/ significant ACME including U","Prop. of simulations where conclusions based on ACME match","Average ACME excluding U","Average ACME including U","Average absolute difference of ACME including U minus ACME excluding U","Prop. of simulations w/ significant ADE excluding U","Prop. of simulations w/ significant ADE including U","Prop. of simulations where conclusions based on ADE match","Average ADE excluding U","Average ADE including U","Average absolute difference of ADE including U minus ADE excluding U")
#correlation
corA<-matrix(0,nrow=(3+length(Ctype)+length(Utype)),ncol=(3+length(Ctype)+length(Utype)))
set.seed(seed)
#######################################
# Loop through all simulations
#######################################
for(si in 1:nSim){ #loop through all simulations
#set.seed(seed+(si-1)) #set the seed to ensure reproducibility
simStep<-5
if(floor(si/simStep)==ceiling(si/simStep)){print(paste("simulation",si, "of",nSim))}
#######################################
# Generate the unmeasured confounder,
# measured confounder C, exposure A,
# mediator M, and outcome Y
#######################################
# generate unmeasured confounder U
U<-matrix(0,nrow=n,ncol=length(Utype))
for(iu in 1:length(Utype)){
if(Utype[iu]=="C"){
if(varU[iu]==0|varU[iu]<0){stop(paste("Error: Variance of U",iu," (varU) is 0 or negative.",sep=""))}
U[,iu]<-rnorm(n,muU[iu],sqrt(varU[iu]))
}
if(Utype[iu]=="D"){
if(muU[iu]<0|muU[iu]==0|muU[iu]==1|muU[iu]>1){stop(paste("Error: Prob(U",iu,"=1) (i.e. muU",iu,") is less than or equal to zero or greater than or equal to 1.",sep=""))}
U[,iu]<-rbinom(n,1,muU[iu])
}
colnames(U)<-paste("U",length(Utype),Utype,sep="")
}
# generate measured confounder C
CC<-matrix(0,nrow=n,ncol=length(Ctype))
for(ic in 1:length(Ctype)){
if(Ctype[ic]=="C"){
if(varC[ic]==0|varC[ic]<0){stop(paste("Error: Variance of C",ic," (varC) is 0 or negative.",sep=""))}
CC[,ic]<-rnorm(n,muC[ic],sqrt(varC[ic]))
}
if(Ctype[ic]=="D"){
if(muC[ic]<0|muC[ic]==0|muC[ic]==1|muC[ic]>1){stop(paste("Error: Prob(C",ic,"=1) is less than or equal to zero or greater than or equal to 1.",sep=""))}
CC[,ic]<-rbinom(n,1,muC[ic])
}
colnames(CC)<-paste("C",length(Ctype),Ctype,sep="")
}
# generate exposure A as a function of C and U
# logit(P(A=1)) or E[A]=gamma0+gammaU*U+gammaC*C
muA<-gamma0+CC%*%gammaC+U%*%gammaU
if(Atype=="C"){
A<-rnorm(n,muA,sqrt(varA)) # normally distributed mediator
if(varA==0|varA<0){stop("Error: Variance of M (varM) is 0 or negative.")}
}
if(Atype=="D"){
Ap<-exp(muA)/(1+exp(muA))
if(max(Ap)=="NaN"|is.na(max(Ap))){stop("Error: Prob(A=1) is too close to 0 or 1. This can occur if the absolute value of gamma0 is too large.")}
A<-rbinom(n,1,Ap)# binary mediator
if(length(A[A==1])/length(A)<(2/n)|length(A[A==1])/length(A)>(1-2/n)){
stop("Error: there was no enough varaibality in A (i.e. Prob(A=1)=1 or Prob(A=1)=0). Considering increasing the sample size n or using mean centered confounders (i.e. muC=0 and muU=0).")
}
}
# generate mediator M
# Logit(P(M=1)) or E[M]=alpha0+alphaA*A+alphaC*C+alphaU*U)
muM<-alpha0+alphaA*A+CC%*%alphaC+U%*%alphaU
if(Mtype=="C"){
M<-rnorm(n,muM,sqrt(varM)) # normally distributed mediator
if(varM==0|varM<0){stop("Error: Variance of M (varM) is 0 or negative.")}
}
if(Mtype=="D"){
Mp<-exp(muM)/(1+exp(muM))
if(max(Mp)=="NaN"|is.na(max(Mp))){stop("Error: Prob(M=1) is too close to 0 or 1. This can occur if the absolute value of alpha0 is too large.")}
M<-rbinom(n,1,Mp)# binary mediator
if(length(M[M==1])/length(M)<(5/n)|length(M[M==1])/length(M)>(1-5/n)){
stop("Error: there was no enough varaibality in M (i.e. Prob(M=1)=1 or Prob(M=1)=0). Considering increasing the sample size n or using mean centered confounders (i.e. muC=0 and muU=0).")
}
}
# generate outcome Y
# logit(P(Y=1)) or E[Y]=beta0+betaA*A+betaM*M+betaC*C+betaU*U
if(interact==TRUE){ muY<-beta0+betaA*A+betaM*M+betaI*A*M+CC%*%betaC+U%*%betaU}
if(interact==FALSE){ muY<-beta0+betaA*A+betaM*M+CC%*%betaC+U%*%betaU}
if(Ytype=="C"){
if(varY==0|varY<0){stop("Error: Variance of Y (varY) is 0 or negative.")}
Y<-rnorm(n,muY,sqrt(varY)) # normally distributed outcome
}
if(Ytype=="D"){
Yp<-exp(muY)/(1+exp(muY))
if(max(Yp)=="NaN"|is.na(max(Yp))){stop("Error: Prob(Y=1) is too close to 0 or 1. This can occur if the absolute value of beta0 is too large.")}
Y<-rbinom(n,1,Yp) # binary outcome
if(length(Y[Y==1])/length(Y)<(5/n)|length(Y[Y==1])/length(Y)>(1-5/n)){
stop("Error: there was no enough varaibality in Y (i.e. Prob(Y=1)=1 or Prob(Y=1)=0). Considering increasing the sample size n or using mean centered confounders (i.e. muC=0 and muU=0).")
}
}
#######################################
# Check for collinearity
#######################################
if(length(Utype)>1|length(Ctype)>1){
if(max(vif(lm(Y~A+M+CC+U))[,3])>10 & interact==FALSE ){stop("Error: There is strong evidence of collinearity (i.e. VIF>10) when A,M,C,U are regressed on Y. Consider reducing the absolute value of large gammas, alphas, or betas.")}
if(max(vif(lm(Y~A+M+A*M+CC+U))[,3])>10 & interact==TRUE ){stop("Error: There is strong evidence of collinearity (i.e. VIF>10) when A,M,A*M,C,U are regressed on Y. Consider reducing the absolute value of large gammas, alphas, or betas.")}
if(max(vif(lm(M~A+CC+U))[,3])>10){stop("Error: There is strong evidence of collinearity (i.e. VIF>10) when A,C,U are regressed on M. Consider reducing the absolute value of large gammas or alphas.")}
}
if(length(Utype)==1&length(Ctype)==1){
if(max(vif(lm(Y~A+M+CC+U)))>10 & interact==FALSE ){stop("Error: There is strong evidence of collinearity (i.e. VIF>10) when A,M,C,U are regressed on Y. Consider reducing the absolute value of large gammas, alphas, or betas.")}
if(max(vif(lm(Y~A+M+A*M+CC+U)))>10 & interact==TRUE ){stop("Error: There is strong evidence of collinearity (i.e. VIF>10) when A,M,A*M,C,U are regressed on Y. Consider reducing the absolute value of large gammas, alphas, or betas.")}
if(max(vif(lm(M~A+CC+U)))>10){stop("Error: There is strong evidence of collinearity (i.e. VIF>10) when A,C,U are regressed on M. Consider reducing the absolute value of large gammas or alphas.")}
}
#######################################
# Correlation for U,C,A,M,Y
#######################################
matA<-cbind(A,M,Y,CC,U)
corA<-corA+cor(matA)
#print(cor(matA))
#######################################
# Mediation analysis w/ and w/out U
#######################################
# models without U
if(Mtype=="C"){med.fit <-lm(M~A+CC)} #fit the model for the mediator
if(Mtype=="D"){med.fit <-glm(M~A+CC,family=binomial(link = "logit"))}
if(interact==FALSE){
if(Ytype=="C"){out.fit<-lm(Y~M+A+CC)} #fit the model for the outcome
if(Ytype=="D"){out.fit<-glm(Y~M+A+CC,family=binomial(link = "logit"))}
}
if(interact==TRUE){
if(Ytype=="C"){out.fit<-lm(Y~M+A+CC+A*M)} #fit the model for the outcome
if(Ytype=="D"){out.fit<-glm(Y~M+A+CC+A*M,family=binomial(link = "logit"))}
}
med.out <- mediate(med.fit, out.fit, treat = "A", mediator = "M",sim=nBoot)
# models with U
if(Mtype=="C"){med.fitU <-lm(M~A+CC+U)} #fit the model for the mediator
if(Mtype=="D"){med.fitU <-glm(M~A+CC+U,family=binomial(link = "logit"))}
if(interact==FALSE){
if(Ytype=="C"){out.fitU<-lm(Y~M+A+CC+U)} #fit the model for the outcome
if(Ytype=="D"){out.fitU<-glm(Y~M+A+CC+U,family=binomial(link = "logit"))}
}
if(interact==TRUE){
if(Ytype=="C"){out.fitU<-lm(Y~M+A+CC+U+A*M)} #fit the model for the outcome
if(Ytype=="D"){out.fitU<-glm(Y~M+A+CC+U+A*M,family=binomial(link = "logit"))}
}
med.outU <- mediate(med.fitU, out.fitU, treat = "A", mediator = "M",sim=nBoot)
#Results
Results["Average ACME excluding U",1]<-Results["Average ACME excluding U",1]+summary(med.out)$d.avg
Results["Average ACME including U",1]<-Results["Average ACME including U",1]+summary(med.outU)$d.avg
Results["Average absolute difference of ACME including U minus ACME excluding U",1]<-Results["Average absolute difference of ACME including U minus ACME excluding U",1]+abs(summary(med.out)$d.avg-summary(med.outU)$d.avg)
Results["Average ADE excluding U",1]<-Results["Average ADE excluding U",1]+summary(med.out)$z.avg
Results["Average ADE including U",1]<-Results["Average ADE including U",1]+summary(med.outU)$z.avg
Results["Average absolute difference of ADE including U minus ADE excluding U",1]<-Results["Average absolute difference of ADE including U minus ADE excluding U",1]+abs(summary(med.out)$z.avg-summary(med.outU)$z.avg)
if(summary(med.out)$d.avg.p<alpha){Results["Prop. of simulations w/ significant ACME excluding U",1]<- Results["Prop. of simulations w/ significant ACME excluding U",1]+1}
if(summary(med.outU)$d.avg.p<alpha){Results["Prop. of simulations w/ significant ACME including U",1]<- Results["Prop. of simulations w/ significant ACME including U",1]+1}
if(((summary(med.out)$d.avg.p<alpha)&(summary(med.outU)$d.avg.p<alpha))|((summary(med.out)$d.avg.p>alpha)&(summary(med.outU)$d.avg.p>alpha))){Results["Prop. of simulations where conclusions based on ACME match",1]<- Results["Prop. of simulations where conclusions based on ACME match",1]+1}
if(summary(med.out)$z.avg.p<alpha){Results["Prop. of simulations w/ significant ADE excluding U",1]<- Results["Prop. of simulations w/ significant ADE excluding U",1]+1}
if(summary(med.outU)$z.avg.p<alpha){Results["Prop. of simulations w/ significant ADE including U",1]<- Results["Prop. of simulations w/ significant ADE including U",1]+1}
if(((summary(med.out)$z.avg.p<alpha)&(summary(med.outU)$z.avg.p<alpha))|((summary(med.out)$z.avg.p>alpha)&(summary(med.outU)$z.avg.p>alpha))){Results["Prop. of simulations where conclusions based on ADE match",1]<- Results["Prop. of simulations where conclusions based on ADE match",1]+1}
} #end of simulation loop
#######################################
# Results
#######################################
matC<-round(corA/nSim,digits=2)
colnames(matC)<-c("A","M","Y",paste("C",c(1:length(Ctype)),sep=""),paste("U",c(1:length(Utype)),sep=""))
rownames(matC)<-c("A","M","Y",paste("C",c(1:length(Ctype)),sep=""),paste("U",c(1:length(Utype)),sep=""))
matR<-Results/nSim
listA<-list(matR,matC,paste("Warning: correlations are only valid if at least one of the variables is normally distributed."))
names(listA)<-c("Results","Correlations_Between_Variables","Warning")
listA
} #end of function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.