为什么 poly(raw=T) 与 poly() 的结果大不相同?

机器算法验证 r lme4-nlme 多项式
2022-03-04 06:27:05

我想对两个不同的时间变量进行建模,其中一些在我的数据中高度共线(年龄 + 队列 = 时期)。这样做我遇到了一些麻烦与lmer和交互poly(),但它可能不限于lmer,我得到了与nlmeIIRC 相同的结果。

显然,我对 poly() 函数的作用缺乏了解。我理解什么poly(x,d,raw=T),我认为没有raw=T它会产生正交多项式(我不能说我真的理解这意味着什么),这使得拟合更容易,但不会让你直接解释系数。
读到是因为我使用的是预测函数,所以预测应该是一样的。

但它们不是,即使模型正常收敛。我正在使用中心变量,我首先认为正交多项式可能会导致与共线交互项具有更高的固定效应相关性,但它似乎具有可比性。在这里粘贴了两个模型摘要

这些图有望说明差异的程度。我使用了仅在开发中可用的预测功能。lme4 的版本(在这里听说过),但固定效果在 CRAN 版本中是相同的(而且它们本身也似乎不一样,例如,当我的 DV 的范围为 0-4 时,交互约为 5)。

lmer 电话是

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

预测仅对假数据(所有其他预测变量 = 0)具有固定效应,其中我将原始数据中存在的范围标记为外推 = F。

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

如果需要,我可以提供更多上下文(我没有设法轻松地制作出可重现的示例,但当然可以更加努力),但我认为这是一个更基本的请求:poly()请向我解释一下这个功能。

原始多项式

原始多项式

正交多项式(在Imgur被裁剪,非裁剪)

正交多项式

1个回答

我认为这是 predict 函数中的一个错误(因此是我的错),实际上 nlme 并不共享编辑:应该在最新的 R-forge 版本中修复lme4。)请参阅下面的示例...

我认为您对正交多项式的理解可能还不错。如果您尝试为一类模型编写预测方法,您需要了解的棘手的事情是正交多项式的基础是基于给定的数据集定义的,所以如果您天真(就像我一样! ) 用于model.matrix尝试为一组新数据生成设计矩阵,你会得到一个新的基础——这对于旧参数不再有意义。在我解决这个问题之前,我可能需要设置一个陷阱,告诉人们predict不能使用正交多项式基(或样条基,它们具有相同的属性)。

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)