具有相关回归量的线性 OLS v 混合效应模型

机器算法验证 r 线性模型 混合模式
2022-04-17 06:46:45

阅读@gung 的这篇文章让我尝试重现他的精彩插图，并最终导致对我读过或听过的东西提出质疑，但我想更直观地理解：为什么 OLSlm控制相关变量更好（我可以试探性地说“更好”吗？）比具有不同相交和斜率的混合效应模型？

这是玩具示例，再次尝试与上面引用的帖子平行：

先上数据：

set.seed(0)    
x1 <- c(rnorm(10,3,1),rnorm(10,5,1),rnorm(10,7,1))
x2 <- rep(c(1:3),each=10)
    y1 <- 2  -0.8 * x1[1:10] + 8 * x2[1:10] +rnorm(10)
    y2 <- 6  -0.8 * x1[11:20] + 8 * x2[11:20] +rnorm(10)
    y3 <- 8  -0.8 * x1[21:30] + 8 * x2[21:30] +rnorm(10)
y <- c(y1, y2, y3)

以及不同的型号：

library(lme4)

fit1 <- lm(y ~ x1)
fit2 <- lm(y ~ x1 + x2)
fit3 <- lmer(y ~ x1 + (1|x2))
fit4 <- lmer(y ~ x1|x2, REML=F)

比较模型之间的 Akaike 信息准则 (AIC)：

AIC(fit1, fit2, fit3, fit4)

     df      AIC
     df      AIC
fit1  3 184.5330
fit2  4  97.6568
fit3  4 112.0120
fit4  5 114.8401

因此，最好的模型似乎是，鉴于和lm(y ~ x1 + x2)之间的强相关性，我认为这是有道理的。x1x2 cor(x1,x2) [1] 0.8619565

但问题是，这种行为背后的直觉是什么，当具有不同截距和斜率的混合模型似乎产生了与数据完美拟合的系数时？

coef(lmer(y ~ x1|x2))
$x2

      x1       (Intercept)
1 -1.1595730    11.37746
2 -0.2586303    19.38601
3  0.2829754    24.20038

library(lattice)    
xyplot(y ~ x1, groups = x2, pch=19,
           panel = function(x, y,...) {
             panel.xyplot(x, y,...);
             panel.abline(a=coef(fit4)$x2[1,2], b=coef(fit4)$x2[1,1],lty=2,col='blue');
             panel.abline(a=coef(fit4)$x2[2,2], b=coef(fit4)$x2[2,1],lty=2,col='magenta');
             panel.abline(a=coef(fit4)$x2[3,2], b=coef(fit4)$x2[3,1],lty=2,col='green')
           })

在此处输入图像描述

我确实意识到二元 OLS 的图形拟合看起来也不错：

library(scatterplot3d)
plot1 <- scatterplot3d(x1,x2,y, type='h', pch=16,
              highlight.3d=T)
plot1$plane3d(fit2)

在此处输入图像描述

...而且我不知道这是否会使问题无效。

1个回答

与评论中的对话一致，实际上可能非常明显：在混合效应中，的总体值随着的增加而增加，如下所示：yx1

正相关cor(x1,y) [1] 0.7924759, 在由lm(y ~ x1 + x2)（在 OP 图中描绘）或在下图中形成的平面中很好地捕捉到：

library(car)
library(rgl)
scatter3d (y ∼ x1  +  x2)
rgl.snapshot(filename="scat1.png", fmt="png")

另一方面，这只是在lmer混合效应中丢失了。而不是回归中的正斜率y超过x1，所有斜率系数始终为负，因为在由离散值创建的每个分层级别x2( 1,2和3)，对应的x1向下倾斜：

scatter3d (x=x1,y=x2,z=y, groups=as.factor(x2))

（来源在这里）

我倾向于假设（并且喜欢得到确认）这可能对混合效果的想法很重要，并且不能通过修改lmer调用合成器来纠正它。例如lmer(y ~ (x1|x2), REML=F)，我尝试了lmer(y ~ x1 + (x1|x2), REML=F)概念上相似的输出。

撇开一个事实不谈，在 OP 的示例中，每个离散值x2假设作为一个独特的“单元”角色，这是对分层模型思想的强制，其中单元通常是主体，很明显lmer混合效应模型在这种情况下不是最优的，即使在这种情况下的值x2是更“自然”的单位（例如，x2==1对应于John;x2==2是Tom，并且x2==3是 的值），带有虚拟变量Mary的线性 OLS将比混合效应模型更合适。从不同的角度来看，分布lmy沿每个点的值x1由于单元之间和单元内的变化，轴不遵循高斯分布，正如混合效应模型所假设的那样。

如其中一条评论所示，混合效应模型在 {lme4} 数据库中的情况下会更自然sleepstudy，其中不同的个体对睡眠剥夺的反应时间逐渐增加，并且异方差似乎随着时间的推移而增加上：

require(lme4)
attach(sleepstudy)
fit <- lm(Reaction ~ Days)
plot(fit, pch=19, cex=0.5, col='slategray')

在此处输入图像描述

有助于为每个主题拟合不同的截距和斜率：

plot(Reaction ~ Days, xlim=c(0, 9), ylim=c(200, 480), type='n', 
     xlab='Days', ylab='Reaction')

for (i in sleepstudy$Subject){
  points(Reaction ~ Days, sleepstudy[Subject==i,], pch=19, 
         cex = 0.5, col=i)
  fit<-lm(Reaction ~ Days, sleepstudy[sleepstudy$Subject==i,])
  a <- predict(fit)
  lines(x=c(0:9),predict(fit), col=i, lwd=1)
}

在此处输入图像描述

该模型将对应于会导致每个受试者不同的截距和斜率的lme4调用。另一方面，每个主题会产生不同的截距，但所有主题的斜率完全相同：lmer(Reaction ~ Days +(Days|Subject))coef(lmer(Reaction ~ Days +(Days|Subject)))lmer(Reaction ~ Days +(1|Subject))

plot(Reaction ~ Days, xlim=c(0, 9), ylim=c(200, 480), type='n', 
     xlab='Days', ylab='Reaction')

for (i in sleepstudy$Subject){
  points(Reaction ~ Days, sleepstudy[Subject==i,], pch=19, 
         cex = 0.5, col=i)
  a<-coef(lmer(Reaction ~ Days + (1|Subject)))$Subject[i,1]
  b<-coef(lmer(Reaction ~ Days + (1|Subject)))$Subject[1,2]
  abline(a=a, b=b, col=i, lwd=1)
}

在此处输入图像描述

... 显然，模型不如相对 AIC 值似乎也表明 ( AIC(lmer(Reaction ~ Days +(1|Subject)), lmer(Reaction ~ Days +(Days|Subject)) 1794.465v 1755.628)。

其它你可能感兴趣的问题
上一篇在命名实体识别的上下文中，支持向量机和条件随机场模型有什么区别？ 下一篇带滞后算子的除法