150 likes | 309 Views
Example 1.0 Indirect Effects and the Test of Mediation using the winBUGS software. Use this tutorial along with the generic example 1.0. One of the many findings: vegetation recovery was a function of the age of the stand that burned. r = -0.35. cover is in proportions. age in years. age of
E N D
Example 1.0Indirect Effects and the Test of Mediation using the winBUGS software Use this tutorial along with the generic example 1.0.
One of the many findings: vegetation recovery was a function of the age of the stand that burned. r = -0.35 cover is in proportions. age in years.
age of stand that burned e1 post-fire vegetation cover Indirect Effects as Causal Tests: Step 1 Model A
winBUGS code for Model A # Title: Simple Regression of cov (cover) on age MODEL simple { # Priors b ~ dnorm(0,0.00001) a ~ dnorm(0,0.00001) tau ~ dgamma(0.5,1) # Likelihoods for(i in 1:N) { covhat[i] <- a + b*age[i] cov[i] ~ dnorm(covhat[i],tau) } # Calculate standardized coefficients covstd <- sd(cov[]) covvar <- covstd*covstd agestd <- sd(age[]) bstd <- b*agestd/covstd # Calculate Prediction Efficiency (Rsqr) covmu <- mean(cov[]) covhatmu <- mean(covhat[]) for (i in 1:N) { covdev[i] <- cov[i]-covmu covhatdev[i] <- covhat[i]-covhatmu } covar <- inprod(covdev[] , covhatdev[])/(N-1) bpred <- covar/(sd(cov[])*sd(covhat[])) #bpred is correlation between predicted and observed Rsqr <- bpred*bpred } CONTINUED ON NEXT PAGE
winBUGS code for Model A (cont.) Initial Values list(a=0, b=1, tau=0.01) Data list(N=90,age=c( 40,25,15,15,23,24,35,14,45,35,45,35,35,30,15,40,21,39,39,29,12,19,20,21,3, 17,40,40,57,52,35,40,5,5,28,33,31,48,55,22,21,13,13,25,15,12,28,16,25,28, 38,10,26,23,16,19,15,27,30,30,3,6,25,28,27,19,25,11,11,15,36,24,24,26,26, 6,31,20,15,15,15,16,20,33,13,20,48,35,60,36), cov=c(1.03879744,0.47759241,0.94893572,1.19490019,1.29818904,1.17348657,0.86158595,0.41906623,0.12851131,0.3062645,0.68247702,0.5310083,0.2957547,0.98468564,0.69869688,1.16663393,0.84068273,1.11949189,0.41204136,0.84594715,1.17590928,0.66100507,0.6517726,0.41172415,0.53282988,0.49739296,0.20203482,0.05557657,0.12124606,0.08297883,0.54920017,0.54646221,1.27310233,0.84805791,0.37941846,0.35841475,0.63674533,0.290822,0.50717877,0.63749067,1.06473979,0.76755405,1.32067202,1.07900425,0.9377264,0.69915774,0.30582847,0.55437438,0.83649961,0.8507124,0.74700881,0.51896384,1.12002246,0.56923477,0.79578273,0.9167156,0.57708428,0.67706467,0.61410302,0.62828578,1.1074659,0.84054394,0.90855822,0.73840393,0.27628054,0.77771328,0.59338111,0.58821905,0.51846397,0.34806766,1.09761107,0.54498063,0.67605775,0.28858364,0.68021606,1.12022933,0.44701614,0.48445559,1.535408,0.92860808,0.46919542,0.45107384,0.51761156,0.5572632,0.5144298,1.16590727,1.08351251,0.56744869,0.5808883,0.37861547))
Indirect Effects as Causal Tests: Step 1 – Results for Model A Statistics: node mean sd MC error 2.5% median 97.5% start sample Rsqr 0.1228 0.0 8.238E-13 0.1228 0.1228 0.1228 1000 19001 a 0.9175 0.08044 5.282E-4 0.7604 0.9169 1.078 1000 19001 b -0.008861 0.002836 1.967E-5 -0.01452 -0.008858 -0.003299 1000 19001 bstd -0.351 0.1123 7.792E-4 -0.5753 -0.3509 -0.1307 1000 19001 0
e1 e2 age of stand that burned fire severity Model B post-fire vegetation cover Indirect Effects as Causal Tests: Step 2
winBUGS code for Model B # Title: Example 1.1 Model B, with cov, fire (fire severity), and age MODEL simple { # Priors b1 ~ dnorm(0,0.00001); b2 ~ dnorm(0,0.00001) a1 ~ dnorm(0,0.00001); a2 ~ dnorm(0,0.00001) tau1 ~ dgamma(0.5,1); tau2 ~ dgamma(0.5,1) # Likelihoods for(i in 1:N) { firehat[i] <- a1 + b1*age[i] fire[i] ~ dnorm(firehat[i],tau1) covhat[i] <- a2 + b2*fire[i] cov[i] ~ dnorm(covhat[i],tau2) } # Calculate standardized coefficients firestd <- sd(fire[]); agestd <- sd(age[]); b1std <- b1*agestd/firestd covstd <- sd(cov[]); b2std <- b2*firestd/covstd CODE CONTINUED ON NEXT PAGE
winBUGS code for Model B # Calculate Prediction Efficiency (Rsqr) firemu <- mean(fire[]) firehatmu <- mean(firehat[]) covmu <- mean(cov[]) covhatmu <- mean(covhat[]) for (i in 1:N) { firedev[i] <- fire[i]-firemu firehatdev[i] <- firehat[i]-firehatmu covdev[i] <- cov[i]-covmu covhatdev[i] <- covhat[i]-covhatmu } covar1 <- inprod(firedev[] , firehatdev[])/(N-1) b1pred <- covar1/(sd(fire[])*sd(firehat[])) #b1pred is cor between predicted and observed Rsqr1 <- b1pred*b1pred covar2 <- inprod(covdev[] , covhatdev[])/(N-1) b2pred <- covar2/(sd(cov[])*sd(covhat[])) #b2pred is cor between predicted and observed Rsqr2 <- b2pred*b2pred } DATA CONTINUED ON NEXT PAGE
winBUGS code for Model B (cont.) Initial Values list(a1=0, a2=0, b1=1, b2=1, tau1=0.01, tau2=0.01) Data list(N=90,age=c( 40,25,15,15,23,24,35,14,45,35,45,35,35,30,15,40,21,39,39,29,12,19,20,21,3, 17,40,40,57,52,35,40,5,5,28,33,31,48,55,22,21,13,13,25,15,12,28,16,25,28, 38,10,26,23,16,19,15,27,30,30,3,6,25,28,27,19,25,11,11,15,36,24,24,26,26, 6,31,20,15,15,15,16,20,33,13,20,48,35,60,36), cov=c(1.03879744,0.47759241,0.94893572,1.19490019,1.29818904,1.17348657,0.86158595,0.41906623,0.12851131,0.3062645,0.68247702,0.5310083,0.2957547,0.98468564,0.69869688,1.16663393,0.84068273,1.11949189,0.41204136,0.84594715,1.17590928,0.66100507,0.6517726,0.41172415,0.53282988,0.49739296,0.20203482,0.05557657,0.12124606,0.08297883,0.54920017,0.54646221,1.27310233,0.84805791,0.37941846,0.35841475,0.63674533,0.290822,0.50717877,0.63749067,1.06473979,0.76755405,1.32067202,1.07900425,0.9377264,0.69915774,0.30582847,0.55437438,0.83649961,0.8507124,0.74700881,0.51896384,1.12002246,0.56923477,0.79578273,0.9167156,0.57708428,0.67706467,0.61410302,0.62828578,1.1074659,0.84054394,0.90855822,0.73840393,0.27628054,0.77771328,0.59338111,0.58821905,0.51846397,0.34806766,1.09761107,0.54498063,0.67605775,0.28858364,0.68021606,1.12022933,0.44701614,0.48445559,1.535408,0.92860808,0.46919542,0.45107384,0.51761156,0.5572632,0.5144298,1.16590727,1.08351251,0.56744869,0.5808883,0.37861547), fire=c(3.5,4.05,2.6,2.9,4.3,4,4.8,4.8,7.25,6.2,8.05,7.55,7.25, 7.3,3.8,5.7,4.3,4.2,4.25,5.1,1.95,4.4,4.45,3,1.2,6.4,8.2,6.85,6.15,7.7,6.1,4.3,3.75,3.7,5.7,5.4,4.85,4.85,4.3,3.25,2.5,3.85,2.55,3.3,4,3.4,3.9,4,7.5,4.9,4.7,1.7,3,6,6.7,2,6.05,3.9,2.9,3.7,1.2,1.3,3.6,4.4,3.7,4.85,4.85,5.6,4.1,5.7,4.2,4.8,7.2,9.2,6.2,2.1,3.85,3.9,4.85,5.3,4.6,3.3,5.4,4.3,3.6,3,3.8,4.5,3.9,4.6))
Indirect Effects as Causal Tests: Step 2 – Results for Model B Statistics: node mean sd MC error 2.5% median 97.5% start sample Rsqr1 0.206 0.0 4.517E-13 0.206 0.206 0.206 1000 49001 Rsqr2 0.1911 0.0 4.517E-13 0.1911 0.1911 0.1911 1000 49001 a1 3.043 0.3595 0.001567 2.338 3.042 3.752 1000 49001 a2 1.074 0.1017 4.759E-4 0.873 1.074 1.271 1000 49001 b1 0.05955 0.01263 5.623E-5 0.03447 0.05955 0.08426 1000 49001 b1std 0.4529 0.09608 4.276E-4 0.2622 0.4529 0.6408 1000 49001 b2 -0.08385 0.02095 9.758E-5 -0.1246 -0.08386 -0.04258 1000 49001 b2std -0.4367 0.1091 5.082E-4 -0.6491 -0.4368 -0.2218 1000 49001
e1 e2 age of stand that burned fire severity post-fire vegetation cover Indirect Effects as Causal Tests: Step 3 Model C
winBUGS code for Model C # Title: Example 1.1 Model C, with cov, fire (fire severity), and age MODEL simple { # Priors b1 ~ dnorm(0,0.00001); b2 ~ dnorm(0,0.00001); b3 ~ dnorm(0,00001) a1 ~ dnorm(0,0.00001); a2 ~ dnorm(0,0.00001) tau1 ~ dgamma(0.5,1); tau2 ~ dgamma(0.5,1) # Likelihoods for(i in 1:N) { firehat[i] <- a1 + b1*age[i] fire[i] ~ dnorm(firehat[i],tau1) covhat[i] <- a2 + b2*fire[i] + b3*age[i] cov[i] ~ dnorm(covhat[i],tau2) } # Calculate standardized coefficients firestd <- sd(fire[]); agestd <- sd(age[]); b1std <- b1*agestd/firestd covstd <- sd(cov[]); b2std <- b2*firestd/covstd b3std <- b3*agestd/covstd THE REST OF THE CODE AND THE DATA ARE THE SAME AS FOR MODEL B
Indirect Effects as Causal Tests: Step 3 – Results for Model C Statistics: node mean sd MC error 2.5% median 97.5% start sample Rsqr1 0.206 0.0 4.517E-13 0.206 0.206 0.206 1000 49001 Rsqr2 0.2084 0.0171 8.172E-5 0.1597 0.215 0.2202 1000 49001 a1 3.039 0.3575 0.001604 2.337 3.041 3.737 1000 49001 a2 1.121 0.1052 4.759E-4 0.9138 1.122 1.327 1000 49001 b1 0.0597 0.01259 5.691E-5 0.03501 0.05964 0.08446 1000 49001 b1std 0.454 0.09578 4.328E-4 0.2663 0.4536 0.6423 1000 49001 b2 -0.06713 0.02323 1.113E-4 -0.1122 -0.06716 -0.02123 1000 49001 b2std -0.3497 0.121 5.796E-4 -0.5842 -0.3498 -0.1106 1000 49001 b3 -0.00485 0.003076 1.287E-5 -0.01085 -0.004854 0.001184 1000 49001 b3std -0.1921 0.1218 5.099E-4 -0.4298 -0.1923 0.04689 1000 49001 DENSITY FUNCTIONS ON NEXT PAGE