Quadratic models with R. The use of poly(..) and I(..) functions (R-language)

Question

What causes the different results below?

var1 = c(0.04875,0.13725,0.28350,0.50975,0.77425,0.94700,0.05325,0.14050,0.29725,0.51525,0.79000,0.95400,0.04625,0.15250,0.29000,0.53300,0.79825,0.95225,0.05025,0.14625,0.28800,0.52625,0.78200,0.95925,0.04700,0.14225,0.30325,0.53500,0.79325,0.95875,0.04775,0.13850,0.28675,0.54250,0.78300,0.95175,0.05150,0.12725,0.30175,0.54725,0.79475,0.96275,0.05375,0.14100,0.30050,0.53275,0.78100,0.96175,0.05450,0.15300,0.29650,0.52850,0.80100,0.95675,0.05425,0.13975,0.30875,0.56025,0.80575,0.96100,0.05100,0.15350,0.31175,0.53300,0.78900,0.96000,0.04650,0.13525,0.29600,0.53625,0.78475,0.96375,0.05375,0.13900,0.29600,0.53725,0.78700,0.95800,0.05075,0.14350,0.29225,0.54525,0.80275,0.95800,0.05050,0.13200,0.29850,0.52700,0.80525,0.96150,0.05150,0.14050,0.29450,0.54375,0.79450,0.96375,0.05375,0.13525,0.30475,0.55250,0.79425,0.96025,0.04950,0.14500,0.29425,0.52250,0.78475,0.95650,0.05225,0.14425,0.29225,0.53150,0.80425,0.95375)

var2 = c(1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6)

var3 = c(4,4,4,4,4,4,6,6,6,6,6,6,8,8,8,8,8,8,10,10,10,10,10,10,12,12,12,12,12,12,14,14,14,14,14,14,16,16,16,16,16,16,18,18,18,18,18,18,20,20,20,20,20,20,22,22,22,22,22,22,24,24,24,24,24,24,26,26,26,26,26,26,28,28,28,28,28,28,30,30,30,30,30,30,32,32,32,32,32,32,34,34,34,34,34,34,36,36,36,36,36,36,38,38,38,38,38,38,40,40,40,40,40,40)


summary(aov(var1~as.factor(var2)*var3))
summary(aov(var1~as.factor(poly(var2,1))*poly(var3,1)))

I might not understand well how poly(..) and I(..) work. Could you give me some hints about these two functions?

(I aim to create models of different degree (quadratic, cubic, etc..) and compare their BIC (or AIC)).

Answer 1

You have the first model all sorts of wrong; that model

var1 ~ var2 * var3

says the variance in var1 is explained by the main effects of factor var2 and continuous covariate var3 and their interaction. In other words, the model is one where each level of var2 has a separate intercept and slope for the regressions lines fitted. There are no polynomials here.

The second model is wrong also, but what you actually want is unclear from the description given. That model

var1 ~ poly(var2,1) * poly(var3,1)

where you cast poly(var2,1) as a factor is effectively the same as the first, just with extra effort.

poly() generates orthogonal (by default) polynomials of its first argument of degree specified by the second argument. Hence the first order polynomial of 1:10 is

> poly(1:10, 1)
                1
 [1,] -0.49543369
 [2,] -0.38533732
 [3,] -0.27524094
 [4,] -0.16514456
 [5,] -0.05504819
 [6,]  0.05504819
 [7,]  0.16514456
 [8,]  0.27524094
 [9,]  0.38533732
[10,]  0.49543369
attr(,"degree")
[1] 1
attr(,"coefs")
attr(,"coefs")$alpha
[1] 5.5

attr(,"coefs")$norm2
[1]  1.0 10.0 82.5

attr(,"class")
[1] "poly"   "matrix"

The second orthogonal polynomial of the vector 1:10 is essentially 1:10 and (1:10) * (1:10) but done in a way as to make the two new vectors orthogonal (or uncorrelated)

> poly(1:10, 2)
                1           2
 [1,] -0.49543369  0.52223297
 [2,] -0.38533732  0.17407766
 [3,] -0.27524094 -0.08703883
 [4,] -0.16514456 -0.26111648
 [5,] -0.05504819 -0.34815531
 [6,]  0.05504819 -0.34815531
 [7,]  0.16514456 -0.26111648
 [8,]  0.27524094 -0.08703883
 [9,]  0.38533732  0.17407766
[10,]  0.49543369  0.52223297
attr(,"degree")
[1] 1 2
attr(,"coefs")
attr(,"coefs")$alpha
[1] 5.5 5.5

attr(,"coefs")$norm2
[1]   1.0  10.0  82.5 528.0

attr(,"class")
[1] "poly"   "matrix"

I'm not clear what you want, but if you want to explore models for different polynomials of var2 and var3 then just use

var1 ~ poly(var2, 2) + poly(var3, 2)

for main effects of quadratic polynomials of var2 and var3. Or more complex

var1 ~ poly(var2, 2) * poly(var3, 3)

which is the main effects of a quadratic in var2 and a cubic in var3, plus their interaction.

I( ) isolates or insulates the contents in the parentheses from R's formula parsing code. For example, you might commonly see

var1 ~ var2 + var2^2 + var3 + var3^2

which is the same as var1 ~ poly(var2, 2) + poly(var3, 2) (except the polynomials are not orthogonal), or it should be. Unfortunately, ^ in a formula means ordered terms, i.e. itself plus its interaction, because ^ has special meaning. To stop R interpreting ^ incorrectly, you wrap those terms in I( ). i.e.

var1 ~ var2 + I(var2^2) + var3 + I(var3^2)

However, do note that var2 and I(var2^2) will be correlated (likewise for var3 and I(var3^2)) and correlated variables in a model can cause issues. Hence the use of poly() which produces orthogonal polynomials, as discussed above.

Note also that poly() can give you the usual raw polynomials via use of raw = TRUE. Hence this might be more what you were expecting for the quadratic of the vector 1:10

> poly(1:10, 2, raw = TRUE)
       1   2
 [1,]  1   1
 [2,]  2   4
 [3,]  3   9
 [4,]  4  16
 [5,]  5  25
 [6,]  6  36
 [7,]  7  49
 [8,]  8  64
 [9,]  9  81
[10,] 10 100
attr(,"degree")
[1] 1 2
attr(,"class")
[1] "poly"   "matrix"

But poly(1:10, 2) would be better in a model.