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

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)).

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.

• This is agreat answer to a not-so-easy-to-answer question! thanks GavinSimpson! I got all the stuff about I() and poly(). But I'm not sure about the effect that as.factor has. And did you actually make a mistake on the very first grey box of your answer? didn't you mean to use as.factor on var2? When would someone use a factor in the formula? When the thing it applies on is considered as a covariate? – Remi.b Aug 1 '13 at 18:49
• @Remi.b Thanks. I removed as.factor intentionally and stated that I assumed that var2 was a factor already - to simplify things. as.factor(var2) would create a factor and thus it would be expanded to $l - 1$ columns in the model matrix (assuming the default contrasts and inclusion of an intercept). I wouldn't use as.factor() in a formula - you should separate the model fitting from the data analysis code, but people do and it is valid as model.frame arranges for the function to be called as it is creating the data from which to build the model matrix. – Gavin Simpson Aug 1 '13 at 19:02
• OK! But still, I don't really get what it means that "as.factor(var2) would create a factor and thus it would be expanded to l-1 columns in the model matrix" Why would someone use such a model? Maybe giving a practical example of experimental design where someone would use a factor migh help me! Many Thks – Remi.b Aug 1 '13 at 19:11
• R might have read in the values c(1,2,3,1,2,3,1,2,2,2,2,1,2) as an integer vector. If these are classes, the wrong model, using 1 df for this variable, would be fitted. The $l-1$ columns represent, using default contrasts, the difference in means between the reference level (which the intercept represents) and each other level/class in the factor. Consider a variable Treatment taking values c(0, 5, 10, 50) where 0 is the control. R will most likely not read that in as a factor and if you plug that into aov() it would use a single df and treat Treatment as a linear covariate. – Gavin Simpson Aug 1 '13 at 20:34
• Also it should be noted that the orthogonal polynomials generated by R add extra wrinkles in terms of interpreting your regression coefficients. If you want to know the size of the effect in real terms you'd have to back out the transformation that poly (raw=FALSE, the default) does. I'm give to understand that you can get most of the way to the advantage of poly by mean centering your variables prior to calculating your interaction terms. – russellpierce Aug 2 '13 at 0:07