Pourquoi est-ce que j'obtiens des résultats extrêmement différents pour poly (raw = T) et poly ()?

Je souhaite modéliser deux variables temporelles différentes, dont certaines sont fortement colinéaires dans mes données (âge + cohorte = période). Ce faisant, j'ai rencontré des problèmes avec lmeret et les interactions de poly(), mais ce n'est probablement pas limité à lmer, j'ai obtenu les mêmes résultats avec l' nlmeIIRC.

De toute évidence, ma compréhension de ce que fait la fonction poly () fait défaut. Je comprends ce qui poly(x,d,raw=T)fait et je pensais que sans raw=Tcela, cela fait des polynômes orthogonaux (je ne peux pas dire que je comprends vraiment ce que cela signifie), ce qui facilite l'ajustement, mais ne vous permet pas d'interpréter directement les coefficients.
J'ai lu que parce que j'utilise la fonction de prédiction, les prédictions devraient être les mêmes.

Mais ce n'est pas le cas, même lorsque les modèles convergent normalement. J'utilise des variables centrées et j'ai d'abord pensé que peut-être le polynôme orthogonal conduit à une corrélation à effet fixe plus élevée avec le terme d'interaction colinéaire, mais il semble comparable. J'ai collé deux modèles de résumés ici .

Espérons que ces graphiques illustrent l'ampleur de la différence. J'ai utilisé la fonction prédire qui n'est disponible que dans le dev. version de lme4 (entendu à ce sujet ici ), mais les effets fixes sont les mêmes dans la version CRAN (et ils semblent également éteints par eux-mêmes, par exemple ~ 5 pour l'interaction lorsque mon DV a une plage de 0-4).

L'appel lmer était

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

La prédiction était des effets fixes uniquement, sur de fausses données (tous les autres prédicteurs = 0) où j'ai marqué la plage présente dans les données originales comme extrapolation = F.

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

Je peux fournir plus de contexte si nécessaire (je n'ai pas réussi à produire un exemple reproductible facilement, mais je peux bien sûr essayer plus fort), mais je pense que c'est un moyen plus simple: expliquez- poly()moi la fonction, très bien s'il vous plaît.

Polynômes bruts

Polynômes orthogonaux (coupés, non coupés à Imgur )

Je pense que c'est un bug dans la fonction de prédiction (et donc ma faute), qui en fait nlme ne partage pas . ( Modifier : devrait être corrigé dans la version la plus récente de R-forge lme4.) Voir ci-dessous pour un exemple ...

Je pense que votre compréhension des polynômes orthogonaux est probablement très bien. La chose délicate que vous devez savoir à leur sujet si vous essayez d'écrire une méthode de prédiction pour une classe de modèles est que la base des polynômes orthogonaux est définie en fonction d'un ensemble de données donné, donc si vous naïvement (comme je l'ai fait! ) utiliser model.matrixpour essayer de générer la matrice de conception pour un nouvel ensemble de données, vous obtenez une nouvelle base - qui n'a plus de sens avec les anciens paramètres. Jusqu'à ce que ce problème soit résolu, je devrais peut-être mettre un piège qui indique aux personnes qui predictne fonctionnent pas avec des bases polynomiales orthogonales (ou des bases splines, qui ont la même propriété).

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)