2017-01-20 34 views
0

让我们说我有10个国家(2级)的300家公司(1级)的数据。 1级变量是PQ和大小。二级变量是人均国内生产总值。使用lmer进行多级回归缩放的正确方法[R]

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(选择 “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,:
模型几乎不可识别:非常大的特征值 - 重新定义变量?模型几乎不可识别:较大的特征值比率 - 重新定义变量?

事实上,当我运行的原始数据,我得到一个额外的警告信息模型:

模型没有收敛:堕落黑森州与1个负本征值

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

+0

缩放所有变量'as.data.frame(scale(MLdata))'后,模型正确拟合。 –

+0

谢谢。这在理论上是有效的在多层次上缩放这样的数据吗?由于缩放,每个级别解释的结果和方差都有显着变化。 –

回答

2

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- sizeGDP的统计数字几乎相同。)

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 

非常相似。

关于方差 - 协方差矩阵:

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时,它变得棘手。