使用 lmer 缩放多级回归的正确方法 [R]

Question

假设我有 10 个国家（第 2 级）的 300 家公司（第 1 级）的数据。1 级变量是 PQ 和大小。水平 2 变量是人均 GDP。

library(lme4)

set.seed(1)
PQ <- runif(300,7,21)
id <- (1:300)
country <- sample(1:10,300,replace=T)
size <- sample(1:25000,300,replace=T)
GDP <- sample(800:40000,10,replace=T)
Country1 <- 1:10
L1data <- as.data.frame(cbind(id,country,PQ,size))
L2data <- as.data.frame(cbind(Country1,GDP))
MLdata <- merge(L1data,L2data, by.x = "country", by.y = "Country1")
dummymodel <- lmer(PQ ~ size + GDP + (size|country), data = MLdata, REML = F)

当我运行它时，我收到警告消息

警告消息： 1：一些预测变量的尺度非常不同：考虑重新调整 2：在 checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, 中：模型未能与 max|grad| 收敛= 1.77081 (tol = 0.002, component 1) 3: 在 checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
模型几乎无法识别：非常大的特征值 - 重新调整变量？；模型是几乎无法识别：大特征值比 - 重新调整变量？

事实上，当我使用原始数据运行模型时，我会收到一条额外的警告消息：

模型未能收敛：退化 Hessian 具有 1 个负特征值

我想我需要缩放自变量来解决这个问题。在这样的多级回归中进行缩放的正确方法是什么？这个问题很重要，因为多级模型的结果取决于缩放。或者还有其他一些我无法找到的问题吗？

Answer 1

tl;dr模型具有几乎相同的拟合优度；缩放模型稍微好一点，更可靠；固定效应几乎相等；两种模型都估计奇异随机效应方差协方差矩阵，这使得比较更加困难，这意味着在任何情况下您都不应该依赖这些模型来得出关于方差的结论......

模型在居中和重新缩放后应该具有相同的含义。正如我将在下面展示的，固定效应本质上是等效的。我发现很难让自己相信方差-协方差估计是等效的，但无论如何该模型是奇异的（即，没有足够的信息来拟合正定方差-协方差矩阵）。

使用您的示例并使用缩放参数重新运行：

MLdata <- transform(MLdata,
    size_cs=scale(size),
    GDP_cs=scale(GDP))
m2 <- lmer(PQ ~ size_cs + GDP_cs + (size_cs|country), data = MLdata,
                   REML = FALSE)

比较对数似然：

logLik(dummymodel)  ## -828.4349
logLik(m2)          ## -828.4067

这表明模型非常接近，但缩放模型的表现稍好一些（尽管 0.03 对数似然单位的改进非常小）。

固定效果看起来不同，但我将证明它们是等效的：

fixef(dummymodel)
##   (Intercept)          size           GDP 
##  1.345754e+01  3.706393e-05 -5.464550e-06 
##     fixef(m2)
## (Intercept)     size_cs      GDP_cs 
## 13.86155940  0.27300688 -0.05914308 

（如果你看一下coef(summary(.))两个模型，你会看到，T-统计size和GDP几乎是一样的。）

从这个问题

rescale.coefs <- function(beta,mu,sigma) {
   beta2 <- beta ## inherit names etc.
   ## parameters other than intercept:
   beta2[-1] <- sigma[1]*beta[-1]/sigma[-1]
   ## intercept:
   beta2[1]  <- sigma[1]*beta[1]+mu[1]-sum(beta2[-1]*mu[-1])
   return(beta2)
}     
cm <- sapply(MLdata[c("size","GDP")],mean)
csd <- sapply(MLdata[c("size","GDP")],sd)

rescaled <- rescale.coefs(fixef(m2),mu=c(0,cm),sigma=c(1,csd))
all.equal(rescaled,fixef(m2))
cbind(fixef(dummymodel),rescaled)
##                                rescaled
## (Intercept)  1.345754e+01  1.345833e+01
## size         3.706393e-05  3.713062e-05
## GDP         -5.464550e-06 -5.435151e-06

非常相似。

关于方差-协方差矩阵：

VarCorr(dummymodel)
##  Groups   Name        Std.Dev.   Corr  
##  country  (Intercept) 2.3420e-01       
##           size        1.5874e-05 -1.000
##  Residual             3.8267e+00

VarCorr(m2)
##  Groups   Name        Std.Dev.   Corr 
##  country  (Intercept) 0.0000e+00      
##           size_cs     4.7463e-08   NaN
##  Residual             3.8283e+00      

方差-协方差矩阵都不是正定的；第一个将不同国家间截距的变化与不同国家斜率的变化完全相关，而第二个将不同国家间的截距分配零方差。此外，在任何一种情况下，国家间的变化相对于残差都非常小。如果两个矩阵都是正定的，我们就可以找到从一种情况扩展到另一种情况的变换（我认为这只是将每个元素乘以 (s_i s_j)，其中 s_. 是应用于给定的缩放因子预测器）。当他们不是 PD 时，事情就变得棘手了。