Generate simulated data with predefined partial correlation structure

Question

I teach statistics (undergraduate level). It's an introduction course to statistics and very basic. Something that I find very time-consuming is to create training data for the students.

I use the formula $X_3= (r \times X_1)+(\sqrt{(1 - r^2)} \times X_2$ to create data with different degrees of correlation. Using this equation I will get two variables $X_3$ and $X_2$ that has a correlations coefficient is $r$ (Choice by me). $X_1$ and $X_2$ are two randomly generated variables (same mean, standard deviation). So far no problem.

Does anyone have any idea how to easily create a data where the correlation between $X$ and $Y$ can be controlled away by adding variable $Z$? Usually I want to give students a task like; first make a regression analysis between $X$ and $Y$, interpret the results. Then add $Z$ (multivariate analysis). Does the impact of $X$ on $Y$ change when controlling for $Z$. ($Z$ is typically an mediating or confounding factor)

What I want to accomplish is a data where the relationship between $X$ and $Y$ disappears (completely or partially) when the students add $Z$. Someone who has a smart solution in an easy way to create such data?

Answer 1

Why not just create data following your causal model? If the model is $X \leftarrow Z \rightarrow Y$, then just create random $Z$, and build $X$ as some function of $Z$ plus noise (same for $Y$). (I left error SD=1, but just set it higher to get more realistic significance.)

> set.seed(100)
> z = rnorm(500)
> x = z + rnorm(500, 0, 1)
> y = z + rnorm(500, 0, 1)

In the uncontrolled model, you will get fake association between $X$ and $Y$, as requested:

> summary(lm(y ~ x))    
Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1534 -0.8674  0.0423  0.8893  4.0521 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.06574    0.05513  -1.192    0.234    
x            0.49819    0.03809  13.081   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.233 on 498 degrees of freedom
Multiple R-squared:  0.2557,    Adjusted R-squared:  0.2542 
F-statistic: 171.1 on 1 and 498 DF,  p-value: < 2.2e-16

Association with $X$ disappears upon proper controlling:

> summary(lm(y ~ x + z))

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

Residuals:
     Min       1Q   Median       3Q      Max 
-3.03113 -0.70603  0.04694  0.62110  2.89843 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01238    0.04464  -0.277    0.782    
x            0.02431    0.04228   0.575    0.566    
z            0.99553    0.06095  16.334   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9952 on 497 degrees of freedom
Multiple R-squared:  0.5157,    Adjusted R-squared:  0.5138 
F-statistic: 264.6 on 2 and 497 DF,  p-value: < 2.2e-16