Solved – get wildly different results for poly(raw=T) vs. poly()

lme4-nlmepolynomialr

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.

Raw polynomials

Raw polynomials

Orthogonal polynomials (clipped, nonclipped at Imgur)

Orthogonal polynomials

Best Answer

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)
u.quad <- rnorm(10,sd=0.1)
d <- transform(d,
               ypred = (1+u.int[f])+
               (2+u.slope[f])*x-
               (1+u.quad[f])*x^2)
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)
Related Question