1
$\begingroup$

i am an absolute noob when it comes to R programming but, alas, I have a college course I have to attend (and pass) and i am really, coming from social sciences, struggling. we're doing this assignment where we have to do a parallel mediation model/analysis with more than 2 mediators, dataset of choice (i have a dataset but in class we've only done it with one mediator)... any nice soul that could help out with an "for dummies" r script? i have this script from the lectures but struggling to replicate it with more mediators... I get the basic condition influences result, condition influences mediator than by doing so the result but not with two more mediators. much obliged to anyone willing to help. i suck at this, i know

read in "pmi.sav" data and

## assign a name to that dataset as "sample.data2" sample.data2 <- read_spss("pmi.sav") sample.data2 <- data.table(sample.data2)

## regression predicting reported reactions model1 <- lm(reaction ~ cond, data = sample.data2) summary(model1)

model.M <- lm(pmi ~ cond, data = sample.data2) model.Y <- lm(reaction ~ cond + pmi, data = sample.data2) summary(model.M) summary(model.Y)

## sobel test ## model.M = regression model predicting mediator, with "apath" being the name of the IV ## model.Y = regression model predicting DV, with "bpath" being the name of the mediator sobel.test <- function(model.M, apath, model.Y, bpath) { a <- coef(model.M)[apath] sa <- summary(model.M)$coef[apath, 2] b <- coef(model.Y)[bpath] sb <- summary(model.Y)$coef[bpath, 2] indirect <- a*b

z <- indirect/sqrt(b^2*sa^2 + a^2*sb^2 + sa^2*sb^2)
pval <- 2*pnorm(-abs(z));

return(c(indirect = indirect, z = z, pval = pval))

}

sobel.test(model.M, "cond", model.Y, "pmi")

## 10,000 resample example ## function to get empirical indirect effect indirect.effect <- function(model.M, apath, model.Y, bpath) { a <- coef(model.M)[apath] b <- coef(model.Y)[bpath]

indirect <- a*b
indirect

}

indirect.effect(model.M, "cond", model.Y, "pmi")

rep.indirect <- sapply(1:10000, function(i) { resample <- sample.data2[sample(1:123, replace = T), ] model.M.resample <- lm(formula(model.M), data = resample) model.Y.resample <- lm(formula(model.Y), data = resample) indirect.effect(model.M.resample, "cond", model.Y.resample, "pmi") }, simplify = "array")

quantile(rep.indirect, c(0.025, 0.975))

## plot for 10,000 replication of indirect effects hist(rep.indirect, breaks = 100, col = "gray", main = "", xlab = "Indirect effects") abline(v = quantile(rep.indirect, 0.025), col = "red", lty = 2) ## red dotted (lty = 2) line, lower 95th percentile abline(v = quantile(rep.indirect, 0.975), col = "red", lty = 2) ## red dotted (lty = 2) line, upper 95th percentile

## ------------------------ ## ## indirect effect est in R ## ## ------------------------ ##

## function to get empirical indirect effect indirect.effect <- function(model.M, apath, model.Y, bpath) { a <- coef(model.M)[apath] b <- coef(model.Y)[bpath]

indirect <- a*b
indirect

}

indirect.effect(model.M, "cond", model.Y, "pmi")

## function to make inferecences for indirect effect boot.indirect <- function(data, i, model.M, apath, model.Y, bpath) {

## resample data n times
resample <- data[i,]

## refit regression model on resampled data
model.M.resample <- lm(formula(model.M), data = resample)
model.Y.resample <- lm(formula(model.Y), data = resample)

a <- coef(model.M.resample)[apath] ## a path
b <- coef(model.Y.resample)[bpath] ## b path
c <- coef(model.Y.resample)[apath] ## c path

## effect decomposition
## indirect effect
indirect.boot <- a*b
direct.boot <- c
return(c(direct.boot = direct.boot, indirect.boot = indirect.boot))

}

require(boot) boot_indirect_effect <- boot(data = sample.data2, statistic = boot.indirect, model.M = "model.M", apath = "cond", model.Y = "model.Y", bpath = "pmi", R = 5000, parallel = "multicore", ncpus = 8)

## inference for direct and indirect effect boot.ci(boot_indirect_effect, index = 1, type = "perc") boot.ci(boot_indirect_effect, index = 2, type = "perc")

## more complex model ## Model of M1 model.m1 <- lm(pmi ~ cond + gender + age, data = sample.data2) summary(model.m1)

## Model of M2 model.m2 <- lm(import ~ cond + gender + age, data = sample.data2) summary(model.m2)

## Model of Y model.Y <- lm(reaction ~ cond + pmi + import + gender + age, data = sample.data2) summary(model.Y)

## custom function for inference boot.indirect2 <- function(data, i) {

## resample data based on index i
resample.data <- data[i, ]

## re-estimate the regression models 
model.m1.resample <- lm(formula(model.m1), data = resample.data)
model.m2.resample <- lm(formula(model.m2), data = resample.data)
model.Y.resample <- lm(formula(model.Y), data = resample.data)

## extract coefficients and compute indirect effects
a1 <- coef(model.m1.resample)["cond"]
a2 <- coef(model.m2.resample)["cond"]
b1 <- coef(model.Y.resample)["pmi"]
b2 <- coef(model.Y.resample)["import"]
c <- coef(model.Y.resample)["cond"]

indirect_pmi <- a1*b1
indirect_import <- a2*b2
direct <- c

## output vector, name - values
out <- c(indirect_pmi = indirect_pmi,
         indirect_import = indirect_import,
         total_indirect = indirect_pmi + indirect_import,
         direct = direct)
return(out)

}

boot_indirect_effect2 <- boot(data = sample.data2, statistic = boot.indirect2, R = 10000, parallel = "multicore", ncpus = 8)

## inference for direct and indirect effect boot.ci(boot_indirect_effect2, index = 1, type = "perc") boot.ci(boot_indirect_effect2, index = 2, type = "perc") boot.ci(boot_indirect_effect2, index = 3, type = "perc")

## comparing two specific indirect effects boot_indirect_effect2$t0[1] - boot_indirect_effect2$t0[2] test <- boot_indirect_effect2$t[,1] - boot_indirect_effect2$t[,2] quantile(test, c(0.025, 0.975))

## Multicategorical variables ## read in "web.sav" data and ## assign a name to that dataset as "sample.data3" sample.data3 <- read_spss("web.sav")

multicat.model <- lm(attitude ~ d1 + d2, data = sample.data3) summary(multicat.model)

multicat.model.M <- lm(inter ~ d1 + d2, data = sample.data3) multicat.model.Y <- lm(attitude ~ d1 + d2 + inter, data = sample.data3) summary(multicat.model.M) summary(multicat.model.Y)

boot.indirect.multicat <- function(data, i) {

## resample data
resample.data <- data[i, ]

## re-estimate model of M and Y based on resampled data
multicat.M.rspl <- lm(formula(multicat.model.M), data = resample.data)
multicat.Y.rspl <- lm(formula(multicat.model.Y), data = resample.data)

## recover coefficients
a1 <- coef(multicat.M.rspl)["d1"]
a2 <- coef(multicat.M.rspl)["d2"]
b <- coef(multicat.Y.rspl)["inter"]
cp1 <- coef(multicat.Y.rspl)["d1"]
cp2 <- coef(multicat.Y.rspl)["d2"]

## effect decomposition
relative.ind.d1 <- a1*b
relative.ind.d2 <- a2*b
relative.direct.d1 <- cp1
relative.direct.d2 <- cp2

out <- c(relative.ind.d1 = relative.ind.d1,
         relative.ind.d2 = relative.ind.d2,
         relative.direct.d1 = relative.direct.d1,
         relative.direct.d2 = relative.direct.d2)

return(out)

}

boot_indirect_effect3 <- boot(data = sample.data3, statistic = boot.indirect.multicat, R = 10000, parallel = "multicore", ncpus = 8)

## inference for direct and indirect effect boot_indirect_effect3$t0 boot.ci(boot_indirect_effect3, index = 1, type = "perc") boot.ci(boot_indirect_effect3, index = 2, type = "perc") boot.ci(boot_indirect_effect3, index = 3, type = "perc") boot.ci(boot_indirect_effect3, index = 4, type = "perc")

## omnibus test of total effect ## total effect model without no group dummies multicat.model.null <- lm(attitude ~ 1, data = sample.data3) summary(multicat.model.null) ## total effect model with group dummies (we already have it) summary(multicat.model)

## changes in R-sq (as a whole)
summary(multicat.model)$r.squared - summary(multicat.model.null)$r.squared
## anova (whether R-sq change is significant?)
anova(multicat.model, multicat.model.null)

## omnibus test of direct effect
## direct effect model without no group dummies multicat.model.null2 <- lm(attitude ~ inter, data = sample.data3) summary(multicat.model.null2) ## direct effect model with group dummies (we already have it) summary(multicat.model.Y)

## changes in R-sq (as a whole)
summary(multicat.model.Y)$r.squared - summary(multicat.model.null2)$r.squared
## anova (whether R-sq change is significant?)
anova(multicat.model.null2, multicat.model.Y)

## omnibus test of indirect effect
boot.omnibus.multicat <- function(data, i) {

## resample data
resample.data <- data[i, ]

## re-estimate model of M and Y based on resampled data
multicat.M.null.rspl <- lm(inter ~ 1, data = resample.data)
multicat.M.rspl <- lm(formula(multicat.model.M), data = resample.data)
multicat.Y.rspl <- lm(formula(multicat.model.Y), data = resample.data)

## recover quantity of interest
rsq.change <- summary(multicat.M.rspl)$r.squared - summary(multicat.M.null.rspl)$r.squared
b <- coef(multicat.Y.rspl)["inter"]

omnibus.indirect.multicat <- rsq.change*b

return(omnibus.indirect.multicat)

}

boot_omnibus.test <- boot(data = sample.data3, statistic = boot.omnibus.multicat, R = 10000, parallel = "multicore", ncpus = 8)

## inference for direct and indirect effect boot_omnibus.test$t0 boot.ci(boot_omnibus.test, index = 1, type = "perc")

$\endgroup$

closed as off-topic by Michael Chernick, Peter Flom Apr 24 at 11:11

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Michael Chernick, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

Browse other questions tagged or ask your own question.