机器算法验证 - 随机效应逻辑回归中的高收缩率 - 吾爱随笔录

考虑以下简单示例：

library( rms )
library( lme4 )
params <- structure(list(Ns = c(181L, 191L, 147L, 190L, 243L, 164L, 83L, 
                            383L, 134L, 238L, 528L, 288L, 214L, 502L, 307L, 302L, 199L, 156L, 
                            183L), means = c(0.09, 0.05, 0.03, 0.06, 0.07, 0.07, 0.1, 0.1, 
                                             0.06, 0.11, 0.1, 0.11, 0.07, 0.11, 0.1, 0.09, 0.1, 0.09, 0.08
                            )), .Names = c("Ns", "means"), row.names = c(NA, -19L), class = "data.frame")
SimData <- data.frame( ID = as.factor( rep( 1:nrow( params ), params$Ns ) ),
                   Res = do.call( c, apply( params, 1, function( x ) c( rep( 0, x[ 1 ]-round( x[ 1 ]*x[ 2 ] ) ),
                                                                        rep( 1, round( x[ 1 ]*x[ 2 ] ) ) ) ) ) )
tapply( SimData$Res, SimData$ID, mean )
dd <- datadist( SimData )
options( datadist = "dd" )
fitFE <- lrm( Res ~  ID, data = SimData )
fitRE <- glmer( Res ~ ( 1|ID ), data = SimData, family = binomial( link = logit ), nAGQ = 50 )

即，我们为同一个非常简单的问题（逻辑回归，仅截距）给出了固定效应和随机效应模型。

这就是固定效应模型的样子：

plot( summary( fitFE ) )

这就是随机效应：

dotplot( ranef( fitRE, condVar = TRUE ) )

收缩本身并不令人惊讶，但其程度却令人惊讶。这是一个更直接的比较：

xyplot( plogis(fe)~plogis(re),
    data = data.frame( re = coef( fitRE )$ID[ , 1 ],
                       fe = c( 0, coef( fitFE )[ -1  ] )+coef( fitFE )[ 1 ] ),
    abline = c( 0, 1 ) )

固定效应估计值范围从小于 3% 到大于 11，但随机效应在 7.5% 到 9.5% 之间。（包含协变量使这更加极端。）

我不是逻辑回归随机效应方面的专家，但从线性回归来看，我的印象是，只有非常小的群体规模才会出现如此大幅度的收缩。然而，在这里，即使是最小的组也有近百个观察值，样本量超过 500 个。

这是什么原因？还是我忽略了什么……？

编辑（2017 年 7 月 28 日）。根据@Ben Bolker 的建议，我尝试了如果响应是连续的会发生什么（这样我们就可以消除有关有效样本量的问题，这是特定于二项式数据的）。

因此新SimData的

SimData <- data.frame( ID = as.factor( rep( 1:nrow( params ), params$Ns ) ),
                   Res = do.call( c, apply( params, 1, function( x ) c( rep( 0, x[ 1 ]-round( x[ 1 ]*x[ 2 ] ) ),
                                                                        rep( 1, round( x[ 1 ]*x[ 2 ] ) ) ) ) ),
                   Res2 = do.call( c, apply( params, 1, function( x ) rnorm( x[1], x[2], 0.1 ) ) ) )
data.frame( params, Res = tapply( SimData$Res, SimData$ID, mean ), Res2 = tapply( SimData$Res2, SimData$ID, mean ) )

新模型是

fitFE2 <- ols( Res2 ~ ID, data = SimData )
fitRE2 <- lmer( Res2 ~ ( 1|ID ), data = SimData )

结果与

xyplot( fe~re, data = data.frame( re = coef( fitRE2 )$ID[ , 1 ],
                       fe = c( 0, coef( fitFE2 )[ -1  ] )+coef( fitFE2 )[ 1 ] ),
    abline = c( 0, 1 ) )

是

到目前为止，一切都很好！

然而，我决定再次检查以验证 Ben 的想法，但结果却非常奇怪。我决定以另一种方式检查理论：我返回二元结果，但增加均值以使有效样本量变大。我只是运行params$means <- params$means + 0.5然后重试原始示例，结果如下：

尽管最小（有效）样本量确实急剧增加......

> summary(with(SimData,tapply(Res,list(ID),
+                             function(x) min(sum(x==0),sum(x==1)))))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
33.0    72.5    86.0   100.3   117.5   211.0

...收缩实际上增加了！（成为总数，估计方差为零。）