Simulating data for mediation model

Question

There is a prior question here, but it doesn't quite cover what I need.

I am trying to simulate a scenario with 3 variables in a mediation setup and with measurement error. All variables must be standardized (mean ≈ 0, sd ≈ 1).

The arrows from nowhere are residuals/errors, round corners = latent variables, rectangles = manifest (observed) variables.

Normally when simulating data with measurement error, I do the following:

n = 10e4
x_measure_beta = .8
x_true = rnorm(n)
x_orb = x_true * x_measure_beta + rnorm(n) * get_error(x_measure_beta)
d = data.frame(x_true, x_orb)

get_error() is:

get_error = function(x) {
  x_sq = x^2
  remain_var = 1 - x_sq
  error_beta = remain_var %>% sqrt
  
  return(error_beta)
}

In other words, get error takes a (standardized) beta, squares it to calculate the proportion of variance explained (R2), then subjects that from 1 to get the remaining variance, then takes the square root to get back to a beta coefficient.

Note that the variance of the predictor betas sum to 1:

> c(x_measure_beta, get_error(x_measure_beta))^2 %>% sum %>% sqrt
[1] 1

In the above case, the output is as expected:

> cor(d)
          x_true     x_orb
x_true 1.0000000 0.7993681
x_orb  0.7993681 1.0000000
> sapply(d, sd)
   x_true     x_orb 
0.9994189 0.9987206 

However, in the mediation scenario, this doesn't work. First some settings:

# mediation scenario ------------------------------------------------------
library(stringr)

n = 10e4

#mediation
med_xym = .5
beta_xy_total = .8

#measurement error -- score x true score beta
beta_xT = .9
beta_mT = .9
beta_yT = .9

#betas
beta_xy = (beta_xy_total^2 * (1-med_xym)) %>% sqrt
beta_xm = (beta_xy_total^2 * med_xym) %>% sqrt %>% sqrt
beta_my = beta_xm
beta_ey = get_error(beta_xy_total)

So, the total effect of x on y is .80, and it is 50% mediated by m. For measurement error, I use a score x true score of .90.

Betas are then: x→y is half of the .80, which I calculated by first converting to proportion variance, multiply by .50 (1-med_xym) because mediation is 50%, and then square root to get the beta, .57.

I stipulate that paths x→m = m→y, thus betas x→m and m→y are the square root of the indirect beta x→m→y (.57) giving them both .75. This is because the link goes thru both of them to get to y.

The residual for y is then the square of the total x→y path, subtracted from 1, and then square rooted, .60.

Based on my understanding, the number should thus be:

Then generate:

generate

#generate
x_true = rnorm(n)
x_orb = x_true * beta_xT + rnorm(n) * get_error(beta_xT)

m_true = x_true * beta_xm + rnorm(n) * get_error(beta_xm)
m_orb = m_true * beta_mT + rnorm(n) * get_error(beta_mT)

y_true = x_true * beta_xy + m_true * beta_my + rnorm(n) * beta_ey
y_orb = y_true * beta_yT + rnorm(n) * get_error(beta_yT)

d_t = data.frame(x_true, m_true, y_true)

However, result is incorrect:

> cor(d_t) %>% round(2)
       x_true m_true y_true
x_true   1.00   0.75   0.82
m_true   0.75   1.00   0.86
y_true   0.82   0.86   1.00
> sapply(d_t, sd)
  x_true   m_true   y_true 
0.998325 1.000613 1.371431 

The correlation is a bit too large (should be .80). It is not a coincidence because I have tried with larger samples. Upon examining the variables, I see that the sd of the y is too large. This baffled me because betas seem to be right in terms of variance:

> c(beta_xy, beta_xm*beta_my)^2 %>% sum %>% sqrt
[1] 0.8
> c(beta_xy, beta_xm*beta_my, beta_ey)^2 %>% sum
[1] 1

Thinking it may be due to the non-zero correlation between x and m, I tried with some random variables:

# how does sd increase when adding uncorrelated variables ------------------------------
> t = rnorm(n) + rnorm(n)
> sd(t)
[1] 1.411864
> t = rnorm(n) + rnorm(n) + rnorm(n)
> sd(t)
[1] 1.732708
> t = rnorm(n) + rnorm(n) + rnorm(n) + rnorm(n)
> sd(t)
[1] 2.000496

We see that for random variables, the effect of summing them is that the new sd is the square root of the sum of the former, i.e. sqrt(1+1)=1.41, sqrt(1+1+1)=1.73 etc.

However, what happens with correlated variables:

# how does sd increase when adding correlated variables -------------------
#same as above but with correlated variables
library(MASS)

#simulate
mat_length = 5
r_mat = rep(.5, mat_length*mat_length) %>% matrix(nrow=mat_length, ncol=mat_length)
diag(r_mat) = 1
t = mvrnorm(n = n, mu = rep(0, mat_length), Sigma = r_mat)
t2 = t[, 1] + t[, 2]
t3 = t[, 1] + t[, 2] + t[, 3]
t4 = t[, 1] + t[, 2] + t[, 3] + t[, 4]

> data.frame(t, t2, t3, t4) %>% sapply(sd)
       X1        X2        X3        X4        X5        t2        t3        t4 
0.9960690 0.9998767 1.0001114 0.9996555 0.9979987 1.7251341 2.4412936 3.1531223 

They have higher variances than before (1.41 vs. 1.73, 1.71 vs. 2.44, 2 vs. 2.44). The intuitive explanation is that because they correlate, cases that are extreme in the same direction get more extreme while cases that are extreme in opposite directions are reduced. Thus, more inequality among cases is created by the summing across variables.

---

How do I simulate a mediation scenario as specified above? I'd like to do it manually so that I can understand what is going on, not use software that deals with the complexity for me.

Answer 1

The question has been answered by Jake Westfall, but, as I find the question (and the answer) interesting, I decided to develop further.

Let say you have the mediation model : $x \rightarrow m \rightarrow y $ and you know the correlations $a=\rho(x,m)$, $b=\rho(m,y|x) $ and $ c=\rho(x,y)$. You then know the indirect effect $ab=a*b$ and the direct effect $c' = c-ab$.

To generate data corresponding to this pattern, you can enter the following code:

> n = 100000
> a = .75
> b = .75
> c = .8
> ab = a*b
> cp= c-ab #cp = c'
> x = rnorm(n)
> m = a*x + sqrt(1-a^2)*rnorm(n)
> ey = 1-(cp^2 +b^2 + 2*a*cp*b)
> y = cp*x + b*m + ey*rnorm(n)

where effectively the error of $y$ (ey) is the residual variance of $var(x)+var(m)+2cov(x,m)$.

> Results #------------------------
> d = data.frame(x,m,y)
> round(cor(d),2)
     x    m    y
x 1.00 0.75 0.80
m 0.75 1.00 0.93
y 0.80 0.93 1.00
> sapply(d, sd)
        x         m         y 
0.9961686 0.9973353 0.9961831 
> lm(y ~ x + m)

Call:
lm(formula = y ~ x + m)

Coefficients:
(Intercept)            x            m  
 -0.0008164    0.2393110    0.7491112  
> lm(m~x)

Call:
lm(formula = m ~ x)

Coefficients:
(Intercept)            x  
  0.0007106    0.7478490  

which is exactly what you asked for : variance are all approximately equal to $1$, correlation $a$ and $c$ are as expected (respectively $.75$ and $.80$) as well as parameter estimates.