message("\n")
message("############################################")
message("#        Mediation analysis example        #")
message("############################################")
message("\n")

message(" + Import packages\n")
library('mediation')
# For sobel filter:
library('multilevel')

message("\n + Read the data")
data = read.csv('ExamAnxiety.csv')

message("\n + Calculate the results")
outcome = with(data, sobel(pred=Revise, med=Anxiety, out=Exam))

message("\n + Overview of the analysis")
outcome[1:3]

message("\n + Get coefficients")
a <- outcome$`Mod3: M~X`[2,1]; b <- outcome$`Mod2: Y~X+M`[3,1]
sa <- outcome$`Mod3: M~X`[2,2]; sb <- outcome$`Mod2: Y~X+M`[3,2]

message("\n + Z Aroian")
(z.aroian <- (a*b)/sqrt(b^2*sa^2 + a^2*sb^2 + sa^2 * sb^2))

message("\n + Aroian p-value")
round(2*(1-pnorm(z.aroian)),5)

message("\n + Mediation percentage")
(VAF <- outcome$Indirect.Effect/outcome$`Mod1: Y~X`[2,1])

message("\n + outcome$Indirect.Effect")
outcome$Indirect.Effect
message("\n + outcome$SE")
outcome$SE
message("\n + outcome$z.value")
outcome$z.value
message("\n + outcome$N")
outcome$N

message("\n\np-value from the from CDF of the normal distribution")
round(2*(1-pnorm(outcome$z.value)),5)

message("\n")
message("############################################")
message("#    End of Mediation analysis example     #")
message("############################################")
message("\n")