6
votes
$\begingroup$

My data mydata consists of columns of x1,x2,..,x100,y in R. But I am thinking a linear model with second order terms such as y ~ x1^2 + x2^2 + x1*x2 + ... how do I achieve that within formula or in any other way in R?

When I tried above, my model pls ignored all second order terms. Do I have to manually create those columns?

$\endgroup$
0

3 Answers 3

Reset to default
9
votes
$\begingroup$

The formula documentation for R shows how to do this. In short, you use poly(). For example, make some quadratic data:

x <- rnorm(100)
y <- x + x**2 * 0.5 + rnorm(100)

Now fit this using a second order polynomial (i.e x and x**2) like this

mod <- lm(y ~ poly(x, 2))

Note that this will fit an orthogonal polynomial, so it won't recover 1 and 0.5 as the coefficients in the generating distribution. If for some reason you want that, use poly(x, 2, raw=TRUE). In general you don't for stability reasons, so stick with the cooked version.

There is also polym as in: lm(y ~ polym(x, z, degree=2) for a model with a full set of crossed variables, which is a bit more trouble to interpret, but that's presumably not important with hundreds of variables.

$\endgroup$
8
  • $\begingroup$ Thanks. But I have 100 predictors and I already tried poly. It's not working well with multiple predictors - I get a series of errors when I use multiple predictors with poly. $\endgroup$
    – Tae-Sung Shin
    Apr 6, 2012 at 13:27
  • $\begingroup$ what does 'not working well' mean? $\endgroup$
    – conjugateprior
    Apr 6, 2012 at 13:28
  • $\begingroup$ As I said, I get a series of errors. $\endgroup$
    – Tae-Sung Shin
    Apr 6, 2012 at 13:30
  • 1
    $\begingroup$ To automate separate polys for all the variables, you could construct the formula. This is ugly but seems to work: forma <- eval(paste("y ~", paste(paste('poly(', names(data), ', 2)', sep=''), collapse=" + "))); mod <- lm(forma) $\endgroup$
    – conjugateprior
    Apr 6, 2012 at 13:39
  • 1
    $\begingroup$ It's easier if you do it by pasteing then use as.formula on the resulting character string.... $\endgroup$
    – Ari B. Friedman
    Mar 11, 2014 at 23:13
8
votes
$\begingroup$

Type : 

lm(y ~ x1 + x2 + I(x1*x2) + I(x1^2) + ...)
$\endgroup$
2
  • $\begingroup$ Thanks a lot. But one more thing: I have 100 variables to add like what you describe. Is there any way to automatically generate those terms? Thanks again. $\endgroup$
    – Tae-Sung Shin
    Apr 6, 2012 at 13:24
  • $\begingroup$ The I() operator tells R to treat it as math instead of treating it as an interaction and linear terms. $\endgroup$
    – thomasrive
    Feb 29, 2016 at 11:27
0
votes
$\begingroup$

Here's how to do it in principle, illustrated on a smaller dataset with only 10 predictors:

# Make fake data
mydata = as.data.frame(matrix(rnorm(1100), 100))
names(mydata) = c(paste0("x", 1:10), "y")

# Form a matrix containing all predictor columns but not y
x = as.matrix(mydata[, 1:10])
# Use poly() to form all 2-way interactions and 2nd order terms
x2 = poly(x, degree = 2, raw = TRUE)
# Resave as a data frame including y
mydata2 = as.data.frame(cbind(x2, y = mydata$y))
# Fit the complete linear model
lm2 = lm(y ~ ., data = mydata2)

However, you have 100 predictors. In my experience, with more than 10-15 predictors, R usually cannot allocate enough memory for the matrix containing every 2-way interaction. You will get unhelpful errors or R will simply crash.

If so, consider whether you really need all 2-way interactions. Maybe just a subset would make sense. For instance, you could use poly() as above to form all 2-way interactions within one subset of x's, then again to form interactions between another subset of x's, but not have any interactions across those subsets.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.