I want to create a model matrix for a quadratic regression in R. I have a matrix $X$, where each column is a covariate. Currently, I'm using the code:
Model = model.matrix(~ 1 + X[,1] + X[,2] + X[,3] + I(X[,1]^2/2) + I(X[,2]^2/2) + I(X[,3]^2/2) + I(X[,1]*X[,2]) + I(X[,1]*X[,3])+I(X[,2]*X[,3]))
Then if I'm dealing with a problem with 4 or more covariates, I have to add the new terms by hand. I tried to create the model matrix myself using:
Model = cbind(1,X,X^2/2,mixed_terms)
I don't remember exactly how I was calculating the mixed terms (I think I used a permutation of the columns of $X$), but in any case my code was much slower in the
model.matrix, while I need the model matrix to be created fast because it goes inside a loop in a simulation routine. Thanks!