Simulate Moderated Regression with Cholesky Decomposition

Question

I want to simulate a moderated regression where the slopes are standardized (i.e., can be interpreted like correlations), and I am wondering how to do this with Cholesky decomposition.

My initial approach was to generate three random variables that represented the outcome, "y," and two predictors, "x1" and "x2."

# Create initial dataset of random variables for predictors and dependent variable
raw.variables <- matrix( cbind(rnorm(1000), rnorm(1000), rnorm(1000)), ncol=3, dimnames=list(NULL, c("y","x1","x2")))

Next, I multiplied x1 and x2 to create the interaction term, and I added this to the raw variable matrix.

# Create interaction term and attach it to raw variable matrix
interaction <- apply(raw.variables[,c('x1','x2')], 1, prod) 
raw.variables <- cbind( raw.variables, interaction )

Next, I defined a correlation matrix that represented the standardized slopes and the relationships between all the predictors that I desired.

# Define desired regression parameters
b1 <- .3
b2 <- 0
b3 <- .5

# Define correlations between predictors
cor_x1_x2 <- 0 # No correlation between x1 and x2
cor_x1_interaction <- 0 # No correlation between x1 and interaction
cor_x2_interaction <- 0 # No correlation between x2 and interaction

# Create correlation matrix
correlation.matrix <- matrix( c( 1, b1, b2, b3,
                                 b1, 1, cor_x1_x2, cor_x1_interaction,
                                 b2, cor_x1_x2, 1, cor_x2_interaction,
                                 b3, cor_x1_interaction, cor_x2_interaction, 1 ),
                              nrow=4 )
print( correlation.matrix )

#     [,1] [,2] [,3] [,4]
#[1,]  1.0  0.3    0  0.5
#[2,]  0.3  1.0    0  0.0
#[3,]  0.0  0.0    1  0.0
#[4,]  0.5  0.0    0  1.0

As a last step, I applied the Cholesky decomposition to the correlation matrix and multiplied it together with the raw variable matrix.

# Multiply matrices to generate correlated data
model.data <- data.frame( t( t( chol( correlation.matrix ) ) %*% t( raw.variables ) ) )
names( model.data ) <- c("y","x1","x2","interaction")

Finally, run the regression and ask for standardized regression coefficients.

# Test code with standardized regression
library( QuantPsyc ) # Need function lm.beta()
lm.beta( lm( y ~ x1 + x2 + eval( x1 * x2 ), data=model.data ) ) # You have to manually multiply x1 and x2 together (i.e., instead of using x1*x2) or lm.beta doesn't estimate the interaction slope correctly

#           x1            x2 eval(x1 * x2) 
# 0.2967814978 -0.0121392943  0.0005289024 

In the standardized slope output, you can see that the main effects for the predictors have the correct slopes (plus some noise), but the interaction term is way off. This is because the interaction term in model.data is not the product of x1 and x2 from model.data:

# Compare generated interaction variable to the product of x1 and x2
sum(apply( model.data[,c("x1","x2")], 1, prod ) - model.data$interaction)==0

# [1] FALSE

The reason it isn't working is that the Cholesky decomposition is creating a multivariate dataset that meets the criteria of my correlation matrix, without "knowing" that the fourth variable must be kept as the product of the middle two variables. I assume I am just totally missing something about the relationship between products and their component variables. Is there any way to simulate data for a moderated regression that would yield an interaction term with a specific effect size using Cholesky decomposition?

For the record, this question is specifically about how to simulate moderated regression with Cholesky decomposition, because you could of course also do it like so:

# (Less scalable) Alternative
y <- with(data.frame(raw.variables), b1*x1 + b2*x2 + b3*x1*x2 + rnorm(1000, 0, sd=sqrt(1-sum( b1**2, b2**2, b3**2 ))))
lm.beta( lm( y ~ x1 + x2 + eval( x1 * x2 ), data=data.frame(raw.variables[,c("x1","x2")]) ) )

#           x1            x2 eval(x1 * x2) 
#   0.28163137    0.02982228    0.49765757 

Answer 1

This derivation may help explain what went wrong.

Let's say you have a random vector $\mathbf{z} \sim \mathcal{N}\left( \mathbf{0}, I \right) $ (i.e. mean zero, uncorrelated, variance 1).

From this you'd like to apply a linear transformation to $z_1$ and $z_2$ to get $x_1$ and $x_2$ such that:

$$ \mathrm{Cov}\left( \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} \right) = \begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix} $$

One way to do this is to multiply $\begin{bmatrix}z_1\\z_2 \end{bmatrix}$ by the Cholesky decomposition of the above covariance matrix. For $|\rho| < 1$ the Cholesky decomposition is given by $\begin{bmatrix}1& 0\\\rho&\sqrt{1 - \rho^2} \end{bmatrix}$ Hence multiplying $\mathbf{z}$ by this matrix gives:

$$\begin{align*} x_1 &= z_1\\ \\ x_2 &= \rho z_1 + \left(\sqrt{1 - \rho^2}\right)z_2\\ \\ \end{align*}$$ If you want $x_3$ to be the interaction of $x_1$ and $x_2$ then: $$\begin{align*} x_3 &= x_1x_2 \\ &= z_1\left(\rho z_1 + \left(\sqrt{1 - \rho^2}\right)z_2 \right) \\ &= \rho z_1^2 + \left( \sqrt{1 - \rho^2}\right)z_1z_2 \end{align*} $$

Observe that not only do you have a $z_1z_2$ term, you have a $z_1^2$ term!

Now let's use $z_3$ to make the error term have variance $\sigma^2$, hence error term $\epsilon = \sigma z_3$. Hence $ y = \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \epsilon$ can be written with the variables $\mathbf{z}$ as:

$$ y = \left( \beta_1 + \beta_2 \rho \right) z_1 + \beta_2 \left( \sqrt{1 - \rho^2} \right)z_2 + \beta_3 \rho z_1^2 + \left(\beta_3 \sqrt{1 - \rho^2}\right)z_1z_2 + \sigma z_3 $$

Hence if you constructed a random vector $\mathbf{w} = \begin{bmatrix} z_1 \\ z_2 \\ z_1^2 \\ z_1z_2 \\ z_3 \end{bmatrix} $ then: $$ \begin{bmatrix}y \\x_1\\x_2\\x_3 \end{bmatrix} = \begin{bmatrix} \beta_1 + \beta_2 \rho & \beta_2 \sqrt{1 - \rho^2} & \beta_3 \rho &\beta_3 \sqrt{1 - \rho^2} & \sigma \\ 1 & 0 & 0 & 0 & 0\\ \rho & \sqrt{ 1 - \rho^2} & 0 & 0 & 0\\ 0 & 0 & \rho & \sqrt{ 1 - \rho^2} & 0 & \end{bmatrix} \mathbf{w}$$

Anyway, the non-linearities of an interaction term severely complicates matters.

Also observe that $\mathrm{Cov} \left( x_1 , x_1x_2\right)$ isn't some free parameter you get to specify (i.e. you can't set cov_x1_interaction like that.)

$$ \begin{align*}\mathrm{Cov} \left( x_1 , x_1x_2\right) &= \mathrm{E}\left[x_1^2x_2 \right] \\ &= \mathrm{E}\left[\rho z_1^3 + \sqrt{1 - \rho^2}z_1^2z_2 \right]\\ &= 0 \end{align*} $$ (If all the variables weren't mean zero, this would be different.)

If I were doing this, I would:

1. Generate random variables $z$
2. Do the standard Cholesky stuff to generate $x_1$ and $x_2$ using $z_1$ and $z_2$.
3. Compute $x_3 = x_1x_2$
4. Compute $y = \beta_1 x_1 + \beta_2 x_2 + + \beta_3 x_3 + \epsilon$

This is basically your last, "of course you could always do it like so." What's wrong with that?