Why do I get wildly different results for poly(raw=T) vs. poly()?

I want to model two different time variables, some of which are heavily collinear in my data (age + cohort = period). Doing this I ran into some trouble with lmer and and interactions of poly(), but it's probably not limited to lmer, I got the same results with nlme IIRC.

Obviously, my understanding of what the poly() function does is lacking. I understand what poly(x,d,raw=T) does and I thought without raw=T it makes orthogonal polynomials (I can't say I really understand what that means), which makes fitting easier, but doesn't let you interpret the coefficients directly.
I read that because I'm using the predict function, the predictions should be the same.

But they aren't, even when the models converge normally. I'm using centered variables and I first thought that maybe the orthogonal polynomial leads to higher fixed effect correlation with the collinear interaction term, but it seems comparable. I've pasted two model summaries over here.

These plots hopefully illustrate the extent of the difference. I used the predict-function which is only available in the dev. version of lme4 (heard about it here), but the fixed effects are the same in the CRAN version (and they also seem off by themselves, e.g. ~ 5 for the interaction when my DV has a range of 0-4).

The lmer call was

cohort2_age =lmer(churchattendance ~
poly(cohort_c,2,raw=T) * age_c +
ctd_c + dropoutalive + obs_c + (1+ age_c |PERSNR), data=long.kg)


The prediction was fixed effects only, on fake data (all other predictors=0) where I marked the range present in the original data as extrapolation=F.

predict(cohort2_age,REform=NA,newdata=cohort.moderates.age)


I can provide more context if need be (I didn't manage to produce a reproducible example easily, but can of course try harder), but I think this is a more basic plea: explain the poly() function to me, pretty please.

Orthogonal polynomials (clipped, nonclipped at Imgur)

-

I think this is a bug in the predict function (and hence my fault), which in fact nlme does not share. (Edit: should be fixed in most recent R-forge version of lme4.) See below for an example ...

I think your understanding of orthogonal polynomials is probably just fine. The tricky thing you need to know about them if you are trying to write a predict method for a class of models is that the basis for the orthogonal polynomials is defined based on a given set of data, so if you naively (like I did!) use model.matrix to try to generate the design matrix for a new set of data, you get a new basis -- which no longer makes sense with the old parameters. Until I get this fixed, I may need to put a trap in that tells people that predict doesn't work with orthogonal polynomial bases (or spline bases, which have the same property).

d <- expand.grid(x=seq(0,1,length=50),f=LETTERS[1:10])
set.seed(1001)
u.int <- rnorm(10,sd=0.5)
u.slope <- rnorm(10,sd=0.2)
d <- transform(d,
ypred = (1+u.int[f])+
(2+u.slope[f])*x-
d$y <- rnorm(nrow(d),mean=d$ypred,sd=0.2)
ggplot(d,aes(x=x,y=y,colour=f))+geom_line()+
geom_line(aes(y=ypred),linetype=2)

library(lme4)
fm1 <- lmer(y~poly(x,2,raw=TRUE)+(1|f)+(0+x|f)+(0+I(x^2)|f),
data=d)

fm2 <- lmer(y~poly(x,2)+(1|f)+(0+x|f)+(0+I(x^2)|f),
data=d)
newdat <- data.frame(x=unique(d$x)) plot(predict(fm1,newdata=newdat,REform=NA)) lines(predict(fm2,newdata=newdat,REform=NA),col=2) detach("package:lme4") library(nlme) fm3 <- lme(y~poly(x,2,raw=TRUE), random=list(~1|f,~0+x|f,~0+I(x^2)|f), data=d) VarCorr(fm3) fm4 <- lme(y~poly(x,2), random=list(~1|f,~0+x|f,~0+I(x^2)|f), data=d) newdat <- data.frame(x=unique(d$x))
lines(predict(fm3,newdata=newdat,level=0),col=4)
lines(predict(fm4,newdata=newdat,level=0),col=5)

-
Thank you, that's reassuring. To reiterate: I read that you can't take the orthogonal polynomial fixed effects at face value, but sometimes they seem insanely large. For example if I run an interaction of two cubic polynomials, I get fixed effects for the polynomials and their interactions in the range from -22 to -127400. That just seems way off to me, especially considering that all fixed effects are negative. Would a revised predict function make sense of these fixed effects or did the models falsely converge or is their something wrong in lmer after all? –  Ruben Jun 27 '12 at 19:41
Again, I suspect (but obviously don't know for sure) that everything's fine. Orth. polynomials are good for numerical stability and hypothesis testing, but (as you're figuring out) the actual parameter values can be harder to interpret. The current version of lme4-devel (I just posted a version that should pass tests, it might take ~24 hours to rebuild on r-forge, unless you can build from SVN yourself) should give you matching predictions between raw/ortho polynomials. An alternative is to center and scale continuous predictors à la Schielzeth 2010 Methods in Ecology & Evolution ... –  Ben Bolker Jun 28 '12 at 10:14
Yes, the two polynomials agree perfectly well now. Thanks a lot! I had scaled and centered my predictors, but some models didn't fit with raw polynomials. –  Ruben Jul 2 '12 at 16:20