Por que obtenho resultados muito diferentes para poli (bruto = T) vs. poli ()?

10

Quero modelar duas variáveis ​​de tempo diferentes, algumas das quais são fortemente colineares nos meus dados (idade + coorte = período). Ao fazer isso, tive alguns problemas com ee lmerinterações poly(), mas provavelmente não está limitado a lmer, obtive os mesmos resultados com o nlmeIIRC.

Obviamente, falta entender o que a função poly () faz. Entendo o que poly(x,d,raw=T)faz e pensei que sem raw=Tpolinômios ortogonais (não posso dizer que realmente entendo o que isso significa), o que facilita o ajuste, mas não permite que você interprete os coeficientes diretamente.
Eu li isso porque estou usando a função de previsão, as previsões devem ser as mesmas.

Mas eles não são, mesmo quando os modelos convergem normalmente. Estou usando variáveis ​​centralizadas e pensei primeiro que talvez o polinômio ortogonal leve a uma correlação de efeito fixo mais alta com o termo de interação colinear, mas parece comparável. Eu colei dois resumos de modelos por aqui .

Esperamos que esses gráficos ilustrem a extensão da diferença. Eu usei a função de previsão, que está disponível apenas no desenvolvedor. versão do lme4 (ouvi falar sobre isso aqui ), mas os efeitos fixos são os mesmos na versão CRAN (e também parecem isolados, por exemplo, ~ 5 para a interação quando meu DV tem um intervalo de 0 a 4).

A ligação mais recente foi

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

A previsão foi apenas de efeitos fixos, em dados falsos (todos os outros preditores = 0) em que marquei o intervalo presente nos dados originais como extrapolação = F.

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

Posso fornecer mais contexto, se necessário (não consegui produzir um exemplo reproduzível facilmente, mas é claro que posso me esforçar mais), mas acho que esse é um apelo mais básico: explique a poly()função para mim, por favor.

Polinômios brutos

Polinômios brutos

Polinômios ortogonais (cortados, não cortados em Imgur )

Polinômios ortogonais

Ruben
fonte

Respostas:

10

Eu acho que isso é um bug na função de previsão (e, portanto, minha culpa), que na verdade o nlme não compartilha. ( Editar : deve ser corrigido na versão mais recente do R-forge lme4.) Veja abaixo um exemplo ...

Acho que sua compreensão dos polinômios ortogonais provavelmente está boa. A coisa complicada que você precisa saber sobre eles se estiver tentando escrever um método de previsão para uma classe de modelos é que a base para os polinômios ortogonais é definida com base em um determinado conjunto de dados, portanto, se você ingenuamente (como eu fiz! ) use model.matrixpara tentar gerar a matriz de design para um novo conjunto de dados, você obtém uma nova base - que não faz mais sentido com os parâmetros antigos. Até que eu conserte isso, talvez seja necessário colocar uma armadilha para informar às pessoas que predictnão trabalham com bases polinomiais ortogonais (ou bases spline, que têm a mesma propriedade).

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)
Ben Bolker
fonte
Obrigado, isso é reconfortante. Para reiterar: li que você não pode considerar os efeitos fixos polinomiais ortogonais pelo valor de face, mas às vezes eles parecem insanamente grandes. Por exemplo, se eu executar uma interação de dois polinômios cúbicos, obtenho efeitos fixos para os polinômios e suas interações no intervalo de -22 a -127400. Isso me parece muito distante, especialmente considerando que todos os efeitos fixos são negativos. Uma função de previsão revisada entenderia esses efeitos fixos ou os modelos convergiriam falsamente ou, afinal, é algo errado?
Ruben
Mais uma vez, suspeito (mas obviamente não sei ao certo) que está tudo bem. Orth. os polinômios são bons para estabilidade numérica e teste de hipóteses, mas (como você está descobrindo) os valores reais dos parâmetros podem ser mais difíceis de interpretar. A versão atual do lme4-devel (eu acabei de publicar uma versão que deve passar nos testes, pode levar ~ 24 horas para ser reconstruída no r-forge, a menos que você possa construir a partir do SVN), fornecendo previsões correspondentes entre polinômios brutos / orto. Uma alternativa é centro e escala preditores contínuas à la Schielzeth 2010 Métodos em Ecologia e Evolução ...
Ben Bolker
Sim, os dois polinômios concordam perfeitamente bem agora. Muito obrigado! Eu havia escalado e centralizado meus preditores, mas alguns modelos não se encaixavam em polinômios brutos.
Ruben