#Chapter 9
#9.2.1
library(mma)
a0=0
a1=1
b0=0
b1=1
b2=1
b3=1
b4=1
b5=1
rho=0
mean1=c(0,0)
n=200

#direct moderation: simulation 1
var1=matrix(c(1,rho,rho,1),2,2)
set.seed(1)
d1=mult.norm(mean1,var1,n)

#d1=matrix(runif(2*n),n,2)
c1=rnorm(n)
pred=d1[,1]
mo=d1[,2]
m1=a0+a1*pred+rnorm(n)
y=b0+b1*pred+b2*m1+b3*c1+b4*mo*pred+rnorm(n)
x=cbind(m1,mo,c1)

#nonlinear method
data1<-data.org(x,y,pred=pred,mediator=c(1,3),alpha=0.05,alpha2=0.05)
med1=med(data=data1,n=2,nonlinear=T,D=5)
bootmed1<-boot.med(data=data1,n=10,n2=100,nonlinear=TRUE,
                   all.model=T,nu=0.001,D=5)
bootmode1=boot.mod(bootmed1,vari="mo",continuous.resolution=5,n=10)
summary(bootmode1)
summary(bootmode1,bymed=T)

#linear method
pred=data.frame(pred)
x2=form.interaction(x,pred,inter.cov="mo")
x2=cbind(x,x2)
data2<-data.org(x2,y,pred=pred,mediator=c(1,3),alpha=0.05,alpha2=0.05)
med2=med(data=data2,n=2,xmod="mo")
bootmed2<-boot.med(data=data2,n=10,n2=100,all.model=T,xmod="mo")
bootmode2=boot.mod(bootmed2,vari="mo",continuous.resolution=5,n=10)
summary(bootmode2)
summary(bootmode2,bymed=T)

#nonlinear method with the pred-moderator interaction term as a covariate
bootmed1.2<-boot.med(data=data2,n=10,n2=100,nonlinear=TRUE,
                     all.model=T,nu=0.001,D=5)
bootmode1.2=boot.mod(bootmed1.2,vari="mo",continuous.resolution=5,
                     n=10,margin=0.15)
summary(bootmode1.2,bymed=T)

plot2.mma(bootmode2,vari="m1",moderator="mo")
plot2.mma(bootmode1,vari="m1",moderator="mo")

#9.2.2
#exposue-moderated TVE: simulation 2
set.seed(2)
d1=mult.norm(mean1,var1,n)

#d1=matrix(runif(2*n),n,2)
c1=rnorm(n)
pred=d1[,1]
mo=d1[,2]
m1=a0+a1*pred+b4*mo*pred+rnorm(n)
y=b0+b1*pred+b2*m1+b3*c1+rnorm(n)
x=cbind(m1,mo,c1)
pred=data.frame(pred)

cova=cbind(mo,form.interaction(x,pred,inter.cov="mo"))
data3<-data.org(x,y,pred=pred,mediator=c(1,3),alpha=0.05,
                alpha2=0.05,cova=cova)
med3=med(data=data3,n=2,nonlinear=T,cova=cova)

#nonlinear method
bootmed3<-boot.med(data=data3,n=10,n2=100,nonlinear=TRUE,
                   all.model=T,nu=0.05,cova=cova)
bootmode3=boot.mod(bootmed3,vari="mo",continuous.resolution=5,
                   n=10,cova=cova,margin=0.4)
summary(bootmode3,bymed=T)

#linear method
med4=med(data=data3,n=2,xmod="mo",cova=cova)
bootmed4<-boot.med(data=data3,n=10,n2=100,all.model=T,
                   xmod="mo",cova=cova)
bootmode4=boot.mod(bootmed4,vari="mo",continuous.resolution=5,
                   n=10,cova=cova)
summary(bootmode4)
summary(bootmode4,bymed=T)

plot2.mma(bootmode4,vari="m1",moderator="mo")
plot2.mma(bootmode3,vari="m1",moderator="mo")

#9.2.3
#3rd variable-moderated TVE: simulation 3
set.seed(3)
d1=mult.norm(mean1,var1,n)

c1=rnorm(n)
pred=d1[,1]
mo=d1[,2]
m1=a0+a1*pred+rnorm(n)
y=b0+b1*pred+b2*m1+b3*c1+b4*mo*m1+rnorm(n)
x=cbind(m1,mo,c1)
#nonlinear method
data5<-data.org(x,y,pred=pred,mediator=c(1,3),alpha=0.05,alpha2=0.05)
med5=med(data=data5,n=2,nonlinear=T,D=5)
bootmed5<-boot.med(data=data5,n=10,n2=100,nonlinear=TRUE,
                   all.model=T,nu=0.001,D=5)
bootmode5=boot.mod(bootmed5,vari="mo",continuous.resolution=5,
                   n=10,margin=0.5)
summary(bootmode5,bymed=T)

#linear method
m1=as.matrix(m1)
colnames(m1)="m1"
x2=form.interaction(x,m1,inter.cov="mo")
x2=cbind(x,x2)
data6<-data.org(x2,y,pred=pred,mediator=c(1,3),alpha=0.05,alpha2=0.05)
med6=med(data=data6,n=2,xmod="mo")
bootmed6<-boot.med(data=data6,n=10,n2=100,all.model=T,xmod="mo")
bootmode6=boot.mod(bootmed6,vari="mo",continuous.resolution=5,n=10)
summary(bootmode6)
summary(bootmode6,bymed=T)

#nonlinear method with the interaction term covariate
bootmed5.2<-boot.med(data=data6,n=10,n2=100,all.model=T,
                     xmod="mo",nonlinear=T)
bootmode5.2=boot.mod(bootmed5.2,vari="mo",continuous.resolution=5,
                     n=10,margin=0.5)
summary(bootmode5.2,bymed=T)

plot2.mma(bootmode5,vari="m1",moderator="mo")
plot2.mma(bootmode6,vari="m1",moderator="mo")

#relationship plot: Figure 9.5
x=bootmode6
full.model=x$a.contx$model$model[[1]]
b1<-predict(full.model,se.fit=T,type=x$a.contx$model$type)$fit
coef<-full.model$coefficients[grep(vari,
                                   names(full.model$coefficients))] 
#plot the straight line instead of the loess line
location<-"Figure 5.eps"
postscript(location,width=16,height=8)
par(mfrow=c(5,2),mar=c(6,6,1,1),oma=c(3,2,2,4))
for(q1 in 1:5){ 
  temp.all=x$a.contx$moder.level$levels[,q1]
  b<-marg.den(x$a.contx$data$x[temp.all,"m1"],b1[temp.all],
              x$a.contx$data$w[temp.all]) 
  plot(x$a.contx$data$x[,"m1"],b1,type="n",xlab="m",ylab=paste("f(","m",")",sep=
                                                                 ""),main=paste("mean(mo)=",round(x$a.contx$moder.level$moder.level[q1],2)))
  
  points(b,col=q1)
  if(length(coef)>1){#browser()
    b3=ifelse(is.factor(x$a.contx$moder.level$moder.level),1,
              x$a.contx$moder.level$moder.level[q1])
    b2=coef[2]*b3
    b2.1=coef[2]*b3+coef[1]*mean(b[,1])
    abline(a=mean(b[,2],na.rm=T)-b2.1,b=b2,col=q1)} #
  axis(1,at=x$a.contx$data$x[temp.all,"m1"],labels=F)
  
  a<-marg.den(x$a.contx$pred.new[temp.all,],x$a.contx$data$x[temp.all,"m1"],
              x$a.contx$data$w[temp.all])   
  scatter.smooth(a[,1],a[,2],family="gaussian", xlab=colnames(x$data$dirx)
                 [1],ylim=xlim,ylab=paste("Mean","m",sep="."),main=
                   paste("mean(mo)=",round(x$a.contx$moder.level$moder.level[q1],2)))}
dev.off(2)

#Figure 9.3: Figure of indirect effect
location<-"Figure 3.eps"
postscript(location,width=8,height=8)
par(mfrow=c(2,2),mar=c(6,6,1,1),oma=c(3,2,2,4))
x=summary(bootmode3,bymed=T,plot=F)
re<-sapply(x$a.contx$results$indirect.effect, "[[", 2)
upper<-sapply(x$a.contx$results$indirect.effect, "[[", 6) #use the quantile interval
lower<-sapply(x$a.contx$results$indirect.effect, "[[", 7)
name1<-c(-1.46,-0.29, 0.02, 0.64, 1.52)
bp <- barplot2(re, horiz = TRUE, main=paste("IE from MART for Sim 2"), 
               names.arg=name1,plot.ci = TRUE, ci.u = upper, ci.l = lower,
               cex.names=0.9,beside=FALSE,cex.axis=0.9,las=1,xlim=range(c(upper,
                                                                          lower)),col = rainbow(length(re), start = 3/6, end = 4/6))
x=summary(bootmode4,bymed=T,plot=F)
re<-sapply(x$a.contx$results$indirect.effect, "[[", 2)
upper<-sapply(x$a.contx$results$indirect.effect, "[[", 6) #use the quantile interval
lower<-sapply(x$a.contx$results$indirect.effect, "[[", 7)
#name1<-colnames(x$a.contx$results$direct.effect)
bp <- barplot2(re, horiz = TRUE, main=paste("IE from Linear Models for Sim 2"), 
               names.arg=name1,plot.ci = TRUE, ci.u = upper, ci.l = lower,
               cex.names=0.9,beside=FALSE,cex.axis=0.9,las=1,xlim=range(c(upper,
                                                                          lower)),col = rainbow(length(re), start = 3/6, end = 4/6))
x=summary(bootmode5.2,bymed=T,plot=F)
re<-sapply(x$a.contx$results$indirect.effect, "[[", 2)
upper<-sapply(x$a.contx$results$indirect.effect, "[[", 6) #use the quantile interval
lower<-sapply(x$a.contx$results$indirect.effect, "[[", 7)
name1<-c(-1.32,-0.60, 0.03, 0.58, 1.38)
bp <- barplot2(re, horiz = TRUE, main=paste("IE from MART for Sim 3"), 
               names.arg=name1,plot.ci = TRUE, ci.u = upper, ci.l = lower,
               cex.names=0.9,beside=FALSE,cex.axis=0.9,las=1,xlim=range(c(upper,
                                                                          lower)),col = rainbow(length(re), start = 3/6, end = 4/6))
x=summary(bootmode6,bymed=T,plot=F)
re<-sapply(x$a.contx$results$indirect.effect, "[[", 2)
upper<-sapply(x$a.contx$results$indirect.effect, "[[", 6) #use the quantile interval
lower<-sapply(x$a.contx$results$indirect.effect, "[[", 7)
#name1<-colnames(x$a.contx$results$direct.effect)
bp <- barplot2(re, horiz = TRUE, main=paste("IE from Linear Models for Sim 3"), 
               names.arg=name1,plot.ci = TRUE, ci.u = upper, ci.l = lower,
               cex.names=0.9,beside=FALSE,cex.axis=0.9,las=1,xlim=range(c(upper,
                                                                          lower)),col = rainbow(length(re), start = 3/6, end = 4/6))
dev.off(2)