Правильный способ масштабирования для многоуровневой регрессии с использованием lmer [R]

Question

Скажем, у меня есть данные для 300 фирм (уровень 1) в 10 странах (уровень 2). Переменные уровня 1 - это PQ и размер. Переменная уровня 2 - это ВВП на душу населения.

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,: Model не сходится с max | grad | = 1.77081 (tol = 0.002, компонент 1) 3: В checkConv (attr (opt, «производные»), opt $ par, ctrl = control $ checkConv,:
Модель почти не идентифицируется: очень большое собственное значение - Перемасштабировать переменные ?; Модель почти не идентифицируется: большое отношение собственных значений - переменные масштабирования?

В самом деле, когда я запустить модель с исходными данными я получаю дополнительное предупреждение:

Модель не сходится: вырожденный Гессе с 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 почти идентичны.)

От 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_. Является коэффициентом масштабирования, применяемым к данный предиктор). Когда они не являются ПД, это становится сложно.