Say I have the classic mediation scenario, with predictor X, potential mediator M, and outcome Y. My thought to test whether mediation is present would be to do the following:

Set up one structural equation model with only the direct effect of X on M (i.e., setting the mediating coefficients from X to M and M to Y equal to 0). Then set up another model freeing those coefficients (i.e., a standard mediation structural equation model). Then perform a likelihood ratio test between the two models (which are clearly nested) to arrive at a p-value for the presence of mediation.

Does this method make sense? If so, is it mathematically equivalent to any of the other commonly used methods to test for mediation? What would be the advantages or disadvantages of using this method if it is valid?


I don't think it makes sense, because this is a 2 df test, and you're only interested in testing one parameter - that is the mediated effect.

The mediation effect is equal to (X to M) * (M to Y). If that is equal to zero, in the population, then you don't have mediation. This parameter will be zero if either (X to M) OR (M to Y) is zero. If (X to M) is non-zero, and (M to Y) is zero, your test will be statistically significant, which indicates that there is mediation, when there is no mediation.

  • 2
    $\begingroup$ Great answer! Only one of the paths needs to be 0 to reject mediation. My method simultaneously tested both paths. $\endgroup$
    – Noah
    Mar 2 '17 at 5:23

I would offer a different perspective: I believe ML estimation CAN be used to test for mediation. This is a trick I thought of which gets Sobel's test statistic from a linear model. The trick is the following: suppose you have a dataset with an outcome Y, an exposure of interest X, and a mediating variable M. These variables are related via a conditional linear model that satisfies: $E[Y|X,M] = \beta_0 + \beta_1 X + \beta_2 M$, as well as the marginal linear model $E[Y|X] = \gamma_0 + \gamma_1 X$. The assumption of "no mediation" is specifically that $\beta_1 = \gamma_1$.

In R, you can get an approximate result with LM by duplicating the dataset and using an indicator for fitting the conditional and marginal model. For instance, using the mediation example from lavaan:

set.seed(1234) X <- rnorm(100) M <- 0.5*X + rnorm(100) Y <- 0.7*M + rnorm(100)

Obtains the following result for the derived parameter ab (the product-path coefficient ala Baron and Kenny);

Defined parameters: Estimate Std.err Z-value P(>|z|) ab 0.374 0.092 4.059 0.000

(really the p should be <0.001)

The linear model can be obtained:

Data2 <- rbind(Data, Data)
Data2$ind <- rep(1:0, each=nrow(Data))
fit <- lm(Y ~ X*ind + M:ind, data=Data2)

by not including the main effect for the M:ind product term, the product term x:ind is the difference of coefficients between the marginal and conditional model (reversed, but that is correctable).

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.19287    0.11350   1.699 0.090857 .  
X            0.40995    0.11220   3.654 0.000332 ***
ind         -0.02929    0.16057  -0.182 0.855448    
X:ind       -0.37360    0.16690  -2.238 0.026326 *  
ind:M        0.78832    0.10922   7.218 1.15e-11 ***

You'll notice however the inference is quite different, still significant but much less precise. This is because the standard errors are the average of the two models. We know theoretically that the conditional model must have smaller standard errors if the design matrix is non-singular. To permit this, we must fit two regression models, obtain their log likelihoods and add them across, then fit the "reduced" model where the $\beta_1$ and $\gamma_1$ are a common coefficient. Estimating that common coefficient is tricky, but if you aim to do it with ML, there may be estimating equations but I'm too lazy to write them down and did it numerically.

negLogLik1 <- function(parms) {
  resid1 <- Y - parms[1] - parms[2]*X - parms[3]*M
  resid2 <- Y - parms[4] - parms[5]*X
  s1<- {sum(resid1^2)/(n-3)}^.5
  s2<- {sum(resid2^2)/(n-2)}^.5

  -sum(dnorm(resid1, 0, s1, log=T))-sum(dnorm(resid2, 0, s2, log=T))

ml1 <- nlm(negLogLik1, c(0,0,0,0,0))

#parms: b0, b1, b2, g0, g1, s1, s2: constrain g1=b1
negLogLik2 <- function(parms) {
  resid1 <- Y - parms[1] - parms[2]*X - parms[3]*M
  resid2 <- Y - parms[4] - parms[2]*X
  s1<- {sum(resid1^2)/(n-3)}^.5
  s2<- {sum(resid2^2)/(n-2)}^.5

  -sum(dnorm(resid1, 0, s1, log=T))-sum(dnorm(resid2, 0, s2, log=T))

ml2 <- nlm(negLogLik2, c(0,0,0,0))

pchisq(- 2* (ml1$minimum - ml2$minimum), df=1, lower.tail = FALSE)

The test statistic then is 5.2 which follows a chi-square distribution with 1 degree of freedom. Comparing this to the Z value for the lavaan test, which was a normal test statistic of 4, we see the SEM is much more powerful here.

This got me thinking about a whole other way of testing linear models where a referent dataset is duplicated, data are simulated under the specific set of causal assumptions, and a likelihood ratio test can be produced under a wide set of conditions... but I digress.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.