## Created by EmpowerStats @ Mon, 15 Sep 25 17:05:04 +0800## 
#******************** Regarding ALL Following R Code   ******************** 
#*** COPYRIGHT (c) 2010, 2021 X&Y Solutions, ALL RIGHT RESERVED *********** 
#******************* www.EmpowerStats.com ********************************* 
#************************************************************************** 
Sys.setlocale(category = 'LC_ALL', locale = 'English_United States.1252'); 
.libPaths(file.path(R.home(),'library')); 
library(doBy); 
options(timeout=600); 
library(plotrix); 
library(stringi); 
library(stringr); 
library(survival); 
library(rms); 
library(nnet); 
library(car); 
library(mgcv); 
pdfwd<-6; pdfht<-6; 
load('D:/EmpowerRCH/Analysis/zfg2/longitudinal 2010.Rdata'); 
if (length(which(ls()=='EmpowerStatsR'))==0) EmpowerStatsR<-get(ls()[1]); 
names(EmpowerStatsR)<-toupper(names(EmpowerStatsR)); 
originalVNAME<-names(EmpowerStatsR); 
ofname<-'zfg2_179_tbl'; 
attach(EmpowerStatsR); 
sink(paste(ofname,'_datastep.lst',sep='')); 
print('Creating new variable: HDL1.LOG'); 
HDL1.LOG<- ifelse(HDL1>0, log(HDL1), NA); 
summary(HDL1.LOG); 
EmpowerStatsR<-cbind(EmpowerStatsR,HDL1.LOG); 
print('Creating new variable: HDL1.LOG.1'); 
HDL1.LOG.1<-HDL1.LOG; 
HDL1.LOG.1[(HDL1.LOG<0)]<-NA; 
HDL1.LOG.1[(HDL1.LOG==0)]<-NA; 
summary(HDL1.LOG.1); 
EmpowerStatsR<-cbind(EmpowerStatsR,HDL1.LOG.1); 
print('Creating new variable: A'); 
A<-(2*GLU1+TG1)*BMI; 
summary(A); 
EmpowerStatsR<-cbind(EmpowerStatsR,A); 
print('Creating new variable: A.LOG'); 
A.LOG<- ifelse(A>0, log(A), NA); 
summary(A.LOG); 
EmpowerStatsR<-cbind(EmpowerStatsR,A.LOG); 
print('Creating new variable: METS'); 
METS<-A.LOG/HDL1.LOG.1; 
summary(METS); 
EmpowerStatsR<-cbind(EmpowerStatsR,METS); 
print('Creating new variable: BMI2'); 
BMI2<-BMI*BMI; 
summary(BMI2); 
EmpowerStatsR<-cbind(EmpowerStatsR,BMI2); 
print('Creating new variable: CUNBAE'); 
CUNBAE<-rep(NA, times=nrow(EmpowerStatsR)); 
tmp <- is.na(CUNBAE) & (SEX==0); 
tmp <- (!is.na(tmp) & tmp); 
CUNBAE[tmp]<-((0.503*AGE)+10.689+(3.172*BMI)-(0.026*BMI2)+(0.181*BMI)-(0.02*BMI*AGE)-(0.005*BMI2)+(0.00021*BMI2*AGE)-44.988)[tmp]; 
tmp <- is.na(CUNBAE) & (SEX==1); 
tmp <- (!is.na(tmp) & tmp); 
CUNBAE[tmp]<-((0.503*AGE)+(3.172*BMI)-(0.026*BMI2)-(0.02*BMI*AGE)+(0.00021*BMI2*AGE)-44.988)[tmp]; 
summary(CUNBAE); 
EmpowerStatsR<-cbind(EmpowerStatsR,CUNBAE); 
rm(HDL1.LOG,HDL1.LOG.1,A,A.LOG,METS,BMI2,CUNBAE); 
sink(); 
vname<-c(NA,'SEX','SEX.0','SEX.1','AGE','NAFLD','NAFLD.0','NAFLD.1','GLU1','TG1','HDL1','BMI','HDL1.LOG','HDL1.LOG.1','A','A.LOG','METS','BMI2','CUNBAE')[-1]; 
vlabel<-c(NA,'SEX','  0','  1','AGE','NAFLD','  0','  1','GLU1','TG1','HDL1','BMI','HDL1 Log','HDL1 Log','A','A Log','METS','BMI2','CUNBAE')[-1]; 
varused4this <- c('SEX','AGE','NAFLD','GLU1','TG1','HDL1','BMI','HDL1.LOG','HDL1.LOG.1','A','A.LOG','METS','BMI2','CUNBAE'); 
pkgs<-c('quantreg','mgcv','mediation'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(quantreg); 
library(mgcv); 
library(mediation); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'因果中介效应分析'; 
wd.subset=''; 
weights.var<- NA; 
colvname<-'NAFLD'; colvlv<- 2; 
par1<-2; 
timevar<- NA; 
vname.start<- NA; 
par2<- NA; 
subjvname<-'METS'; subjlv<- 0; 
par3<-1; 
chk<- 0; 
cox<- 0; 
bvar<-'CUNBAE'; blv<- 0; 
xvname<-c(); 
dec<-4; 
##R package## quantreg mgcv mediation ##R package##;
pvformat<-function(p,dec) {
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p)); colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  lw <- paste("<",substr("0.00000000000",1,dec+1),"1",sep="");
  pp[as.numeric(p)<(1/10^dec)]<-lw
  return(pp)
}
numfmt<-function(p,dec) {
  if (is.list(p)) p<-as.matrix(p)
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p));colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  pp[as.numeric(p)>10000000]<- "inf."
  pp[is.na(p) | gsub(" ","",p)==""]<- ""
  pp[p=="-Inf"]<-"-Inf"
  pp[p=="Inf"]<-"Inf"
  return(pp)
}
mat2htmltable<-function(mat) {
  t1<- apply(mat,1,function(z) paste(z,collapse="</td><td>"))
  t2<- paste("<tr><td>",t1,"</td></tr>")
  return(paste(t2,collapse=" "))
}
factor2dummy<-function(v, v.name,keepNA=FALSE) {
  vl<-levels(factor(v)); nc<-length(vl);nr<-length(v)
  new.name<-paste(v.name,vl,sep=".")
  d.out<-matrix(0,nrow=nr,ncol=nc)
  for (i in 1:nc) d.out[v==vl[i],i]<-1
  if (sum(is.na(v))>0) {
     if (keepNA) {
       d.out[is.na(v),]<-NA;
     } else {
       d.out<-cbind(d.out,rep(0,nr))
       d.out[is.na(v),nc+1]<-1
       new.name<-c(new.name,paste(v.name,"NA",sep="."))
     }
  }
  colnames(d.out)<-new.name; return(d.out);
}
createXvars<-function(xvname,xlv,WD) {
  nx<-length(xvname); xx.mat<-NA; cname<-c() 
  for (i in 1:nx) {
    v.name<-xvname[i]; v<-WD[,v.name]
    if (xlv[i]>0 & xlv[i]<10) {
      tmp<-factor2dummy(v,v.name,FALSE)
      xx.mat<-cbind(xx.mat,tmp[,-1])
      cname<-c(cname,colnames(tmp)[-1])
    } else {
      xx.mat<-cbind(xx.mat,v)
      cname<-c(cname,v.name)
    }
  }
  out<-xx.mat[,-1];
  if (!is.matrix(out)) out<-matrix(out,ncol=1)
  colnames(out)<-cname
  return(out)
}
dec<-6
vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN];
allname<-c(colvname,timevar,bvar,subjvname,xvname)
allname<-allname[!is.na(allname)]
WD<-WD[apply(is.na(WD[,allname]),1,sum)==0,]
if (!is.na(weights.var)) {weights<-WD[,weights.var];} else {weights<-rep(1,nrow(WD));}
colvb<-vlabelV[vnameV==colvname];
subjb<-vlabelV[vnameV==subjvname];
expob<-vlabelV[vnameV==bvar];
nx<-length(xvname); xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)]
if (nx>0) {
  xx<-createXvars(xvname,xlv,WD)
  WD<-cbind(WD[ , !(names(WD) %in% xvname)],xx)
  xvars<-colnames(xx)
  xvstr<-paste("+", paste(xvars,collapse="+"))
} else {xvstr<-""};
if (blv==2) {
  blvl <- levels(factor(WD[,bvar]))
  if (blvl[1]!="0" | blvl[2]!="1") {
    tmpx<-rep(NA,nrow(WD))
    tmpx<-(WD[,bvar]==blvl[2])*1
    tmpx[is.na(WD[,bvar])]<-NA	
    WD[,bvar]<-tmpx
	rm(tmpx)
  }	
  control<-(WD[,bvar]==0)*1
  control[is.na(WD[,bvar])]<-NA
  WD<-cbind(WD,control)
  control.value <- 0; treat.value <- 1
} else {
  ct.value<- quantile(WD[,bvar],probs=c(0.25,0.75))
  control.value <- ct.value[1]; treat.value <- ct.value[2]
}
if (is.na(par2)) {tau<-0;
} else if (par2==1) {tau<-0.5;
} else if (par2<1 & par2>0) {tau<-par2;
}
if (is.na(timevar) & par1==3) par1<-2
if (par1!=3) timevar<-NA
subj.lv<-levels(factor(WD[,subjvname]))
if (length(subj.lv)==2) {
   if (subj.lv[1]!="0" | subj.lv[2]!="1") {
     tmps<-rep(NA,nrow(WD))
	 tmps<-(WD[,subjvname]==subj.lv[2])*1
	 tmps[is.na(WD[,subjvname])]<-NA
	 WD[,subjvname]<-tmps
     rm(tmps)
   }
   par3<-2
} else {if (length(subj.lv)>20) {par3<-1;} else {par3<-3;}
}
if (par3!=1) chk<-0
fml.m<-paste(subjvname,"~",bvar,xvstr)
if (is.na(timevar)) {
  fml.y<-paste(colvname,"~",bvar,"+",subjvname,xvstr,sep="")
  fml.yxm<-paste(colvname,"~",bvar,"+",subjvname,"+",paste(bvar,subjvname,sep=":"),xvstr,sep="")
  fml.ysm<-paste(colvname,"~",bvar,"+s(",subjvname,",bs='cr')",xvstr,collapse="")
  fml.yxsm<-paste(colvname,"~",bvar,"+s(",subjvname,",by=",bvar,")+s(",subjvname,",by=control)",xvstr,sep="")
} else if (is.na(vname.start)) {
  fml.y<-paste("Surv(",timevar,",",colvname,")~",bvar,"+",subjvname,xvstr,sep="")
  fml.yxm<-paste("Surv(",timevar,",",colvname,")~",bvar,"+",subjvname,"+",paste(bvar,subjvname,sep=":"),xvstr,sep="")
  fml.ysm<-paste("Surv(",timevar,",",colvname,")~",bvar,"+rcs(",subjvname,")",xvstr,sep="")
} else {
  fml.y<-paste("Surv(",vname.start,",",timevar,",",colvname,")~",bvar,"+",subjvname,xvstr,sep="")
  fml.yxm<-paste("Surv(",vname.start,",",timevar,",",colvname,")~",bvar,"+",subjvname,"+",paste(bvar,subjvname,sep=":"),xvstr,sep="")
  fml.ysm<-paste("Surv(",vname.start,",",timevar,",",colvname,")~",bvar,"+ rcs(",subjvname,")",xvstr,sep="")
}
if (par3==1) {
  mdl.m<-lm(formula(fml.m),data=WD, weights=weights) 
} else if (par3==2) {
  mdl.m<-glm(formula(fml.m),family=binomial(link="probit"),data=WD, weights=weights)
} else if (par3==3) {
  WD[,subjvname]<-as.factor(WD[,subjvname])
  mdl.m<-polr(formula(fml.m),method="probit",Hess=TRUE,data=WD, weights=weights)
}
fml4y<-fml.y
if (par1==1) { 
  if (cox==0) {
    if (tau>0) {         mdl.y<-rq(formula(fml.y),tau=tau,data=WD, weights=weights);
    } else if (chk==1) { mdl.y<-gam(formula(fml.ysm), data=WD, weights=weights);          fml4y<-fml.ysm;
    } else {             mdl.y<-lm(formula(fml.y),data=WD, weights=weights)
	}
  } else {
    if (tau>0) {         mdl.y<-rq(formula(fml.yxm),tau=tau,data=WD, weights=weights);    fml4y<-fml.yxm;
    } else if (chk==1) { mdl.y<-gam(formula(fml.yxsm), data=WD, weights=weights);         fml4y<-fml.yxsm;
    } else {             mdl.y<-lm(formula(fml.yxm),data=WD, weights=weights);            fml4y<-fml.yxm;
	}    
  }	
} else if (par1==2) {
  if (cox==0) {
    if (chk==1) { mdl.y<-gam(formula(fml.ysm),family=binomial(link="probit"), data=WD, weights=weights);    fml4y<-fml.ysm;
    } else {      mdl.y<-glm(formula(fml.y),  family=binomial(link="probit"), data=WD, weights=weights)
	}
  } else {
    if (chk==1) { mdl.y<-gam(formula(fml.yxsm),family=binomial(link="probit"), data=WD, weights=weights);   fml4y<-fml.yxsm;
    } else {      mdl.y<-glm(formula(fml.yxm), family=binomial(link="probit"), data=WD, weights=weights);   fml4y<-fml.yxm;
	}   
  }	
} else if (par1==3) {
  if (cox==0) {
    if (chk==1) { mdl.y<-survreg(formula(fml.ysm),dist="weibull",data=WD, weights=weights);    fml4y<-fml.ysm;
    } else {      mdl.y<-survreg(formula(fml.y),  dist="weibull",data=WD, weights=weights)
    } 
  } else {        mdl.y<-survreg(formula(fml.yxm),dist="weibull",data=WD, weights=weights);    fml4y<-fml.yxm;
  }	
}
set.seed(123456)
if (par1==3) {
  med<-try(mediate(mdl.m,mdl.y,treat=bvar,mediator=subjvname,outcome=colvname,control.value=control.value,treat.value=treat.value,weights=weights))
} else {
  if (chk==1 & cox==1) {
    med<-try(mediate(mdl.m,mdl.y,boot=TRUE,treat=bvar,control="control",mediator=subjvname,control.value=control.value,treat.value=treat.value,weights=weights))
  } else {
    med<-try(mediate(mdl.m,mdl.y,boot=TRUE,treat=bvar,mediator=subjvname,control.value=control.value,treat.value=treat.value,weights=weights))
  }
} 
if (substr(med[[1]][1],0,5) != "Error") {
 if (cox==1 | par1==2 | par1==3) {
  es<-rbind(c(med$tau.coef,med$tau.ci),c(med$d0,med$d0.ci),c(med$d1,med$d1.ci),c(med$z0,med$z0.ci),c(med$z1,med$z1.ci),
  c(med$n0,med$n0.ci),c(med$n1,med$n1.ci),c(med$d.avg,med$d.avg.ci),c(med$z.avg,med$z.avg.ci),c(med$n.avg,med$n.avg.ci))
  pv<-c(med$tau.p,med$d0.p,med$d1.p,med$z0.p,med$z1.p,med$n0.p,med$n1.p,med$d.avg.p,med$z.avg.p,med$n.avg.p)
  oo<-cbind(c("Total effect","Mediation effect (control)","Mediation effect (treated)","Direct effect (control)",
      "Direct effect (treated)","Propotion mediated (control)","Propotion mediated (treated)", 
      "Mediation effect (average)","Direct effect (average)","Propotion mediated (average)"),numfmt(es,dec),pvformat(pv,4))
 } else {
  es<-rbind(c(med$tau.coef,med$tau.ci),c(med$d.avg,med$d.avg.ci),c(med$z.avg,med$z.avg.ci),c(med$n.avg,med$n.avg.ci))
  pv<-c(med$tau.p,med$d.avg.p,med$z.avg.p,med$n.avg.p)
  oo<-cbind(c("Total effect","Mediation effect","Direct effect","Propotion mediated"),numfmt(es,dec),pvformat(pv,4))
 }
 if (par1<3) {resid.cor <- numfmt(cor(mdl.y$residuals, mdl.m$residuals), 6);} else {resid.cor <- NA;}
 if (cox==1) {
  tmx <- test.TMint(med)
  tmi <- paste(numfmt(tmx$statistic,dec), " (",paste(numfmt(tmx$conf.int,dec),collapse=", "),"), p = ",pvformat(tmx$p.value,4),sep="")
 }
 oo<-rbind(c("","Estimate","95% CI lower","95% CI upper","P-value"),oo)
 png(paste(ofname,"png",sep="."),width=720,height=560);plot(med);dev.off()
}
sink(paste(ofname,".lst",sep="")); 
print(fml.m); summary(mdl.m); print(fml4y); summary(mdl.y);summary(med);
sink()
w<-"<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>"
w<-c(w, "<h2>Causal Mediation Analysis</h2>")
w<-c(w, "Outcome variable: ", colvb)
w<-c(w, "</br>Exposure variable: ", expob)
if (blv != 2) w<-c(w, "</br>Control value = ", control.value, "; Treat value = ", treat.value)
w<-c(w, "</br>Mediator variable: ", subjb)
if (chk==1) w<-c(w, "</br>Generalized additive (GAM) model was used to smooth mediator effect on outcome") 
if (tau>0) w<-c(w, "</br>Using uantile regression (tau=", tau, ") for outcome")
if (par1==3) w<-c(w, "</br>Survival regression for outcome, time variable:", timevar)
if (par1==3 & !is.na(vname.start)) w<-c(w, ", start time:", vname.start)
if (par1==2) w<-c(w, "</br>Generalized linear model for outcome with link function = probit")
w<-c(w, "</br>Adjusting variables:", paste(xb, collapse="; "))
if (substr(med[[1]][1],0,5) != "Error") {
if (par1!=3) {
   w<-c(w, "</br>Nonparametric Bootstrap (resamples=1000) Confidence Intervals with the Percentile Method");
 } else {w<-c(w, "</br>Quasi-Bayesian Confidence Intervals");
 }
 w<-c(w, "</br></br><table border=3>", mat2htmltable(oo), "</table>")
 if (cox==1) {
   w<-c(w, "</br>Test of X*M interaction: difference of mediation effect between X=1 and X=0 (95%CI): ")
   w<-c(w, "</br> ", tmi, "</br>")
 }  
 if (!is.na(resid.cor)) w<-c(w, "</br>Correlation of residuals from mediator model with residuals from outcome model: ", resid.cor, "</br>")  
 w<-c(w, "</br>Sample size used: ", med$nobs, "</br>")
}
w<-c(w, wd.subset)
w<-c(w,"</br>Reference:")
w<-c(w,"</br>Imai, K., Keele, L., and Tingley, D. (2010b). A general approach to causal mediation analysis.
Psychological Methods, 15(4), 309–334.")
w<-c(w,"</br>https://cran.r-project.org/web/packages/mediation/vignettes/mediation-old.pdf")
w<-c(w,paste("</br></br>Created by EmpowerStats (www.empowerstats.com) and R on",Sys.Date()))
w<-c(w, "</body></html>")
fileConn<-file(paste(ofname,".htm",sep="")); writeLines(w, fileConn)