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

Question

假设我有10个国家(2级)300家公司(1级)的数据。1级变量是PQ和大小。二级变量是人均国内生产总值。

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：一些预测变量的值非常不同scales:考虑重新缩放2:在checkConv(attr(opt，"derivs")中，opt$par，ctrl=control$checkConv，：模型未能与max | grad |=1.77081(tol=0.002，组分1)3:IncheckConv(attr(opt，"derivs")，opt$par，ctrl=control$checkConv，：
模型几乎是不可识别的：非常大的特征值-重新缩放变量？；模型几乎不可识别：特征值比率大-重新缩放变量？

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

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

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

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

从这个问题

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_is_j)，其中s_。是应用于给定预测器的缩放因子)。当他们不是PD时，这就变得棘手了。