使用 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)

当我运行收到警告消息时

Warning messages: 1: Some predictor variables are on very different scales: consider rescaling 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 1.77081 (tol = 0.002, component 1) 3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue - Rescale variables?; Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?

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

Model failed to converge: degenerate Hessian with 1 negative eigenvalues

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

Answer 1

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

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

使用您的示例和 re-running 缩放参数：

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)

比较log-likelihoods：

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

这表明模型非常接近，但缩放模型的效果稍好（尽管 0.03 log-likelihood 个单位的改进非常小）。

固定效应看起来不同，但我要证明它们是等价的：

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(.))，您会发现 size 和 GDP 的 t-statistics 几乎相同。）

来自this question

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

非常相似。

关于 variance-covariance 矩阵：

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

variance-covariance 矩阵都不是正定的；第一个具有各国截距的变化与各国斜率的变化完全相关，而第二个将零方差分配给各国的截距。此外，在任一情况下，among-country 相对于剩余方差的变化都非常小。如果两个矩阵都是正定的，我们可以努力找到从一种情况扩展到另一种情况的变换（我认为这只是将每个元素乘以（s_i s_j)，其中 s_. 是应用于给定预测变量的比例因子）。当他们不是 PD 时，它就变得棘手了。

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

Correct way to scale for multilevel regression using lmer [R]

scaling

r

lme4

multi-level