#Chapter 6 codes are at the appendix
library(mma)
library(MASS)
a0=0.2
a1=0.1
a=0.773
b=0.701
n=50
nSims=100
result1=NULL
result2=NULL
set.seed(13)
for (rho in c(-0.9, -0.5, -0.1, 0, 0.1, 0.5, 0.9))
  for (c1 in c(0, 0.518, 1.402, 2.150))
    for (cy in c(0, 0.259, 0.701, 1.705))
      for (i in 1:nSims){
        e1=rnorm(n,0,1)
        e2=rnorm(n,0,1)
        sigma=matrix(c(1,rho,rho,1),2,2)
        G=mvrnorm(n=50,rep(0,2),sigma)
        X1=G[,1]
        C=G[,2]
        pred=X1
        X=a0+a*pred+e1
        Y=a1+b*X+c1*pred+cy*C+e2
        Xo=data.frame(X)
        X2=data.frame(X,C)
        Yo=data.frame(Y)
        pred1=data.frame(X1)
        #linear results
        temp2<-mma(Xo,Yo,pred=pred1,mediator=1,jointm=list(n=1,
                                                           j1=1),alpha=0.1,alpha2=0.1)
        #M is forced in as a third-variable
        D=summary(temp2,plot=F)
        result1=rbind(result1,c(a,b,c1,rho,cy,i,D$results$total.effect,
                                D$results$direct.effect,D$results$indirect.effect[[1]][,2]))
        #nonlinear results
        temp2<-mma(Xo,Yo,pred=pred1,mediator=1,jointm=list(n=1,j1=1),
                   alpha=0.1,alpha2=0.1,nonlinear=T)
        D=summary(temp2,plot=F)
        result2=rbind(result2,c(a,b,c1,rho,cy,i,D$results$total.effect,
                                D$results$direct.effect,D$results$indirect.effect[[1]][,2]))
      }
# The right model if C is treated as a covariate
temp2<-mma(X2,Yo,pred=pred1,mediator=1,jointm=list(n=1,
                                                   j1=1),alpha=0.1,alpha2=0.1)
# The right model if C is treated as another third variable
temp2<-mma(X2,Yo,pred=pred1,mediator=1:2,jointm=list(n=1,
                                                     j1=1),alpha=0.1,alpha2=0.1)
# Codes to generate data for Assumption 2
b=0.701
c1=1.402
n=50
nSims=100
result1=NULL
result2=NULL
set.seed(14)
for (rho in c(-0.9, -0.5, -0.1, 0, 0.1, 0.5, 0.9))
  for (a in c(0,0.286,0.773,1.185))
    for (cm in c(0, 0.259, 0.701, 1.705))
      for (i in 1:nSims){
        e1=rnorm(n,0,1)
        e2=rnorm(n,0,1)
        sigma=matrix(c(1,rho,rho,1),2,2)
        G=mvrnorm(n=50,rep(0,2),sigma)
        X1=G[,1]
        C=G[,2]
        pred=X1
        X=a0+a*pred+cm*C+e1
        Y=a1+b*X+c1*pred+e2
        Xo=data.frame(X)
        Yo=data.frame(Y)
        pred1=data.frame(pred)
        Co=data.frame(C)
        #linear results
        temp2<-mma(Xo,Yo,pred=pred1,mediator=1,jointm=list(n=1,j1=1),
                   alpha=0.1,alpha2=0.1)
        D=summary(temp2,plot=F)
        result1=rbind(result1,c(a,b,c1,rho,cm,i,D$results$total.effect,
                                D$results$direct.effect,
                                D$results$indirect.effect[[1]][,2]))
        #nonlinear results
        temp2<-mma(Xo,Yo,pred=pred1,mediator=1,jointm=list(n=1,j1=1),
                   alpha=0.1,alpha2=0.1,nonlinear=T)
        D=summary(temp2,plot=F)
        result2=rbind(result2,c(a,b,c1,rho,cm,i,D$results$total.effect,
                                D$results$direct.effect,
                                D$results$indirect.effect[[1]][,2]))
      }
# The right model
temp2<-mma(Xo,Yo,pred=pred1,mediator=1, cova=data.frame(C),
           jointm=list(n=1,j1=1),alpha=0.1,alpha2=0.1)        
        
        
#Codes to generate data for Assumption 3
a0=0.2        
a1=0.1
a=0.773
c1=0.701
n=50
nSims=100
result1=NULL
result2=NULL
set.seed(2)
for (rho in c(-0.9, -0.5, -0.1, 0, 0.1, 0.5, 0.9))
  for (b in c(0,0.518,1.402,2.150))
    for (cy in c(0,0.259,0.701,1.705))
      for (i in 1:nSims){
        e1=rnorm(n,0,1)
        e2=rnorm(n,0,1)
        sigma=matrix(c(1,rho,rho,1),2,2)
        G=mvrnorm(n=50,rep(0,2),sigma)
        e3=G[,1]
        C=G[,2]
        pred=e1
        X=a0+a*pred+e3
        Y=a1+b*X+c1*pred+cy*C+e2
        Xo=data.frame(X)
        X2=data.frame(X,C)
        Yo=data.frame(Y)
        pred1=data.frame(pred)
        #linear results
        temp2<-mma(Xo,Yo,pred=pred1,mediator=1,jointm=list(n=1,j1=1),
                   alpha=0.1,alpha2=0.1)
        D=summary(temp2,plot=F)
        result1=rbind(result1,c(a,b,c1,rho,cy,i,D$results$total.effect,
                                D$results$direct.effect,
                                D$results$indirect.effect[[1]][,2]))
        #nonlinear results        
        temp2<-mma(Xo,Yo,pred=pred1,mediator=1,jointm=list(n=1,j1=1),
                   alpha=0.1,alpha2=0.1,nonlinear=T)
        D=summary(temp2,plot=F)
        result2=rbind(result2,c(a,b,c1,rho,cy,i,D$results$total.effect,
                                D$results$direct.effect,
                                D$results$indirect.effect[[1]][,2]))
      }
#codes to fit the right model
temp2<-mma(X2,Yo,pred=pred1,mediator=1,jointm=list(n=1,j1=1),
           alpha=0.1,alpha2=0.1)

#Codes to generate data for Assumption 4
a1=0.773
b1=0.701
c1=1.250
n=50
nSims=100
result1=NULL
set.seed(5)
for (cm in c(-0.9, -0.5, -0.1, 0, 0.1, 0.5,0.9))
  for (a2 in c(0,0.286,0.773,1.185))
    for (b2 in c(0,0.259,0.701,1.705))
      for (i in 1:nSims){
        pred=rnorm(n,0,1)
        e2=rnorm(n,0,1)
        e3=rnorm(n,0,1)
        e4=rnorm(n,0,1)
        M2=a2*pred+e2
        M1=a1*pred+cm*M2+e3
        Y=b1*M1+c1*pred+b2*M2+e4
        Xo=data.frame(M1,M2)
        Yo=data.frame(Y)
        pred1=data.frame(pred)
        #linear results
        temp1<-mma(Xo,Yo,pred=pred1,mediator=1:2,jointm=list(n=1,j1=1:2),
                   alpha=0.1,alpha2=0.1)
        D1=summary(temp1,plot=F)
        result1=rbind(result1,c(a,b,c1,d,cm,cy,i,D1$results$total.effect,
                                D1$results$direct.effect,
                                D1$results$indirect.effect[[1]][,2],
                                D1$results$indirect.effect[[1]][,3],
                                D1$results$indirect.effect[[1]][,4]))
      }