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