2017-01-20 66 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時,它變得棘手。