# Interpreting a mediation model with binary outcome

I have a mediation model with two continuous mediators (m1; m2), a continuous input variable (x) and a dichotomous output variable (y).

The two mediators are different mechanisms of the input variable. I understand mediation should be useful to determine whether one, both or none of the mediators increases the chances of the outcome variable to be 1.

I have tried to calculate the model with two different methods: a) the lavaan package for R, and b) Baron and Kenny's steps.

One difference in the calculations is that the lavaan modeling approach calculates the input variable x based on regression coefficients of 4 predictors, whereas for the Baron and Kenny approach the variable x is just the mean of these 4 predictors (I just didn't know how to do it).

The outcomes of the two models differ quite a lot as can be seen in the attached figures. However the main problem I have at this moment of time, is the interpretation of the coefficients. In figure 3 and 4 the estimates as they are returned from R are on the respective arrows.

## Baron and Kenny model with estimates

To make some sense out of these estimates, I tried to transform then into a percentage. A percentage increase of the odds ratio of the dependent variable in case of increasing the respective predictor by 1. These are in figures 5 and 6.

## Baron and Kenny model with standardized estimates

To do so I did the following for the regression coefficients with binary outcomes:

round((exp(COEFF)-1)*100,2)


Now if I look at the lavaan model, could I say:

(1) The total effect of x on y is high (11.52%), the direct effect of x on y is rather low (-4.78%).

(2) An increase of 1 in m1, increases the odds of Y to be 1 by 90%.

Generally, is it a bad idea to try some kind of transformation here? I can't really understand how to make sense out of the different coefficients in the indirect effects paths.

#### Minimum Working Example ####

library(QuantPsyc)
library(psych)
library(lavaan)

data <- read.csv("http://dl.dropboxusercontent.com/u/63191123/data.csv", header = TRUE, sep ="\t", stringsAsFactors=TRUE)
sapply(data, mode);sapply(data, class);attach(data);head(data);psych::describe(data)

model <- '
# measurement model
m1 =~ m11 + m12 + m13 + m14 + m15 + m16
m2 =~ m21 + m22 + m23 + m24 + m25
m24 ~~ m25
m13 ~~ m14
m16 ~~ m21
m25 ~~ m22
m23 ~~ m24
m13 ~~ m15
m12 ~~ m13
m24 ~~ m22
'

# fit model
fit <- cfa(model, mimic="Mplus", data=data)

m1 <- (parameterEstimates(fit)[1,4] * m11 + parameterEstimates(fit)[2,4] * m12 + parameterEstimates(fit)[3,4] * m13 + parameterEstimates(fit)[4,4] * m14 + parameterEstimates(fit)[5,4] * m15 + parameterEstimates(fit)[6,4] * m16)/6
m2 <- (parameterEstimates(fit)[7,4] * m21 + parameterEstimates(fit)[8,4] * m22 + parameterEstimates(fit)[9,4] * m23 + parameterEstimates(fit)[10,4] * m24 + parameterEstimates(fit)[11,4] * m25)/5

#### Baron and Kenny ####

#<<<<<< Baron and Kenny - Step 1 >>>>>>#
Baron_c.mod <- glm(Y ~ X, family = binomial(link = "logit"))

#<<<<<< Baron and Kenny - Step 2 >>>>>>#
Baron_a1.mod <- lm(m1 ~ X)
Baron_a2.mod <- lm(m2 ~ X)

#<<<<<< Baron and Kenny - Step 3 >>>>>>#
Baron_b.mod <- (glm(Y ~ X + m1 + m2, family = binomial(link = "logit")))

# Baron and Kenny Summary:
Path_a1 <- summary(Baron_a1.mod)$coefficients[2,] Path_a2 <- summary(Baron_a2.mod)$coefficients[2,]
Path_b1 <- summary(Baron_b.mod)$coefficients[3,] Path_b2 <- summary(Baron_b.mod)$coefficients[4,]
Path_c <- summary(Baron_c.mod)$coefficients[2,] Path_c_ <- summary(Baron_b.mod)$coefficients[2,]

Baron <- rbind(Path_a1, Path_a2, Path_b1, Path_b2, Path_c, Path_c_)
print(Baron, digits=3)

#Baron.std <- rbind(lm.beta(Baron_a1.mod), lm.beta(Baron_a2.mod), lm.beta(Baron_b.mod)[2], lm.beta(Baron_b.mod)[3], lm.beta(Baron_c.mod)[1], lm.beta(Baron_b.mod)[1])
#Baron.std <- Baron.std[,1]
#Baron <- data.frame(Baron, Baron.std)

################
#### Lavaan ####
################

data$m1 <- m1 data$m2 <- m2

mediation <- '
X =~ x1 + x2 + x3 + x4
Y ~ b1*m1 + b2*m2 + c*X
m1 ~ a1*X
m2 ~ a2*X
indirect_m1 := a1*b1
indirect_m2 := a2*b2
total := c + (a1*b1) + (a2*b2)
'

mediation.fit <- sem(mediation, std.lv=FALSE, ordered = "Y", data=data) #
summary(mediation.fit)

###### Standardized coefficients for logit regressions:
# e.g. for the logistic regressions:
# Baron and Kenny Model
round((exp(summary(Baron_b.mod)$coefficients[3,1])-1)*100,2) # b1 round((exp(summary(Baron_b.mod)$coefficients[4,1])-1)*100,2) # b2
round((exp(summary(Baron_c.mod)$coefficients[2,1])-1)*100,2) # c round((exp(summary(Baron_b.mod)$coefficients[2,1])-1)*100,2) # c'

# Lavaan Model
round((exp(parameterEstimates(mediation.fit)[5,5])-1)*100,2) #b1
round((exp(parameterEstimates(mediation.fit)[6,5])-1)*100,2) #b2
round((exp(parameterEstimates(mediation.fit)[28,5])-1)*100,2) #c
round((exp(parameterEstimates(mediation.fit)[7,5])-1)*100,2) #c'

• I'm confused by a couple of things - your estimates are 0.109 (from glm?) and 0.172 (from lavaan?). Why aren't these the same? Second, are the percentage of variance? I'm not sure this is appropriate with a dichotomous variable. – Jeremy Miles Oct 1 '13 at 18:34
• Also, could you post some code with a toy example? – Jeremy Miles Oct 1 '13 at 19:35
• Thanks for your comments Jeremy. I think the differences from lavaan to glm come from my different approach in how I calculated X. I am not sure about the percentages, or the estimates in general. I can't quite understand how to find whether there is mediation or not when my mediating paths (a1*b1 and a2*b2) consist of different types of estimates (log odds and betas from lm). – Roman Oct 2 '13 at 12:44
• You can multiple the two values, and this gives you a log odds. – Jeremy Miles Oct 2 '13 at 13:39