########################################## ### data simulation ########################################## # set the size of the sample n=1000 # set seed set.seed(123) # simulate x (normal distribution) X <- rnorm(n, 5, 4) # calculate the mean of X mean_x=mean(X) # simulate a residual for M residual_1<-rnorm(n,0,1) M_1<-0.3+0.5*X+residual_1 residual_2<-rnorm(n,0,1.1) M_2<-0.1+0.7*X+0.5*M_1+residual_2 # mu for Poisson regression via a log link mu_1 <- exp(0.2 + 0.3*M_1+0.3*M_2+0.01*X) # use rpois to generate Y Y <- rpois(n, lambda=mu_1) # combine into a dataframe and print out the first 6 rows data <- data.frame(X=X, M_1=M_1,M_2=M_2, Y=Y) head(data) ################################################## ### Serial Mediation Analysis for Count data ################################################## library(boot) Serial_Mediation_poisson<-function(data_used,i,x_predetermined=0) { # Sample a data data_temp=data_used[i,] # Deciding which X value to use if(x_predetermined==0){x_predetermined=mean(data_temp\$X)} else if (x_predetermined==-1){x_predetermined=mean(data_temp\$X)-sd(data_temp\$X)} else(x_predetermined=mean(data_temp\$X)+sd(data_temp\$X)) # a path result_a<-lm(M_1~X, data = data_temp) a_0<-result_a\$coefficients[1] a_1<-result_a\$coefficients[2] # d path result_d<-lm(M_2~X+M_1, data = data_temp) d_0<-result_d\$coefficients[1] d_1<-result_d\$coefficients[2] d_2<-result_d\$coefficients[3] # b path result_b<-glm(Y~M_1+M_2+X, data = data_temp, family = quasipoisson) b_0<-result_b\$coefficients[1] b_1<-result_b\$coefficients[2] b_2<-result_b\$coefficients[3] c_1_apostrophe<-result_b\$coefficients[4] #calculating the indirect effect M_1_estimated=a_0+a_1*x_predetermined M_2_estimated=d_0+d_1*x_predetermined+d_2*M_1_estimated indirect_effect<-a_1*d_2*b_2*exp(b_0+b_1*M_1_estimated+b_2*M_2_estimated+c_1_apostrophe*x_predetermined) return(indirect_effect) } # use boot() to do bootstrapping mediation analysis boot_mediation <- boot(data, Serial_Mediation_poisson, R=1000) boot_mediation # plot the 1000 indirect effects plot(boot_mediation) # print out confidence intervals boot.ci(boot.out = boot_mediation, type = c("norm", "perc")) ################################################## ### calculate the theoretical value ################################################## x_mean=mean(data\$X) M_1=0.3+0.5*x_mean M_2=0.1+0.7*x_mean+0.5*M_1 0.5*0.5*0.3*exp(0.2 + 0.3*M_1+0.3*M_2+0.01*x_mean) ################################################## ### calculate the original indirect value ################################################## # a path result_a<-lm(M_1~X, data = data) a_0<-result_a\$coefficients[1] a_1<-result_a\$coefficients[2] a_0 a_1 # d path result_d<-lm(M_2~X+M_1, data = data) d_0<-result_d\$coefficients[1] d_1<-result_d\$coefficients[2] d_2<-result_d\$coefficients[3] d_0 d_1 d_2 # b path result_b<-glm(Y~M_1+M_2+X, data = data, family = quasipoisson) b_0<-result_b\$coefficients[1] b_1<-result_b\$coefficients[2] b_2<-result_b\$coefficients[3] c_1_apostrophe<-result_b\$coefficients[4] b_0 b_1 b_2 c_1_apostrophe #calculating the indirect effect M_1_estimated=a_0+a_1*x_mean M_2_estimated=d_0+d_1*x_mean+d_2*M_1_estimated indirect_effect<-a_1*d_2*b_2*exp(b_0+b_1*M_1_estimated+b_2*M_2_estimated+c_1_apostrophe*x_mean) indirect_effect