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 #  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