lme () и lmer () дают противоречивые результаты

Я работал с некоторыми данными, которые имеют некоторые проблемы с повторными измерениями. При этом я заметил очень различное поведение lme()и lmer()использование моих тестовых данных и хочу знать почему.

Поддельный набор данных, который я создал, содержит измерения роста и веса для 10 предметов, взятых дважды каждый. Я настроил данные так, чтобы между субъектами была положительная связь между ростом и весом, но отрицательная связь между повторными измерениями в каждом человеке.

set.seed(21)
Height=1:10; Height=Height+runif(10,min=0,max=3) #First height measurement
Weight=1:10; Weight=Weight+runif(10,min=0,max=3) #First weight measurement

Height2=Height+runif(10,min=0,max=1) #second height measurement
Weight2=Weight-runif(10,min=0,max=1) #second weight measurement

Height=c(Height,Height2) #combine height and wight measurements
Weight=c(Weight,Weight2)

DF=data.frame(Height,Weight) #generate data frame
DF$ID=as.factor(rep(1:10,2)) #add subject ID
DF$Number=as.factor(c(rep(1,10),rep(2,10))) #differentiate between first and second measurement

Вот график данных, с линиями, соединяющими два измерения от каждого человека.

Так что я провел две модели, одна с lme()из nlmeпакета и один с lmer()от lme4. В обоих случаях я проводил регрессию веса против роста со случайным эффектом ID, чтобы контролировать повторные измерения каждого человека.

library(nlme)
Mlme=lme(Height~Weight,random=~1|ID,data=DF)
library(lme4)
Mlmer=lmer(Height~Weight+(1|ID),data=DF)

Эти две модели часто (хотя и не всегда в зависимости от семян) давали совершенно разные результаты. Я видел, где они генерируют немного разные оценки дисперсии, вычисляют различные степени свободы и т. Д., Но здесь коэффициенты находятся в противоположных направлениях.

coef(Mlme)
#   (Intercept)    Weight
#1   1.57102183 0.7477639
#2  -0.08765784 0.7477639
#3   3.33128509 0.7477639
#4   1.09639883 0.7477639
#5   4.08969282 0.7477639
#6   4.48649982 0.7477639
#7   1.37824171 0.7477639
#8   2.54690995 0.7477639
#9   4.43051687 0.7477639
#10  4.04812243 0.7477639

coef(Mlmer)
#   (Intercept)    Weight
#1     4.689264 -0.516824
#2     5.427231 -0.516824
#3     6.943274 -0.516824
#4     7.832617 -0.516824
#5    10.656164 -0.516824
#6    12.256954 -0.516824
#7    11.963619 -0.516824
#8    13.304242 -0.516824
#9    17.637284 -0.516824
#10   18.883624 -0.516824

Для наглядной иллюстрации модель с lme()

И модель с lmer()

Почему эти модели так сильно расходятся?

tl; dr, если вы измените оптимизатор на "nloptwrap", я думаю, что он избежит этих проблем (вероятно).

Поздравляем, вы нашли один из простейших примеров множественных оптимумов в задаче статистической оценки! Параметр, который lme4использует внутренне (таким образом, удобный для иллюстрации), представляет собой масштабированное стандартное отклонение случайных эффектов, то есть стандартное отклонение между группами, деленное на остаточное стандартное отклонение .

Извлеките эти значения для оригинала lmeи lmerподходит:

(sd1 <- sqrt(getVarCov(Mlme)[[1]])/sigma(Mlme))
## 2.332469
(sd2 <- getME(Mlmer,"theta")) ## 14.48926

Переустановите с помощью другого оптимизатора (это, вероятно, будет по умолчанию в следующей версии lme4):

Mlmer2 <- update(Mlmer,
  control=lmerControl(optimizer="nloptwrap"))
sd3 <- getME(Mlmer2,"theta")   ## 2.33247

Матчи lme... посмотрим, что происходит. Функция отклонения (-2 * log правдоподобия), или в этом случае аналогичная функция критерия REML, для LMM с одним случайным эффектом принимает только один аргумент, поскольку параметры фиксированного эффекта профилируются ; они могут быть вычислены автоматически для заданного значения стандартного отклонения RE.

ff <- as.function(Mlmer)
tvec <- seq(0,20,length=101)
Lvec <- sapply(tvec,ff)
png("CV38425.png")
par(bty="l",las=1)
plot(tvec,Lvec,type="l",
     ylab="REML criterion",
     xlab="scaled random effects standard deviation")
abline(v=1,lty=2)
points(sd1,ff(sd1),pch=16,col=1)
points(sd2,ff(sd2),pch=16,col=2)
points(sd3,ff(sd3),pch=1,col=4)
dev.off()

Я продолжал дальше зацикливаться на этом и побежал припадки для случайных семян от 1 до 1000, установки lme, lmerи lmer+ nloptwrap для каждого случая. Вот числа из 1000, где данный метод получает ответы, которые по крайней мере на 0,001 единицы отклонения хуже, чем другой ...

          lme.dev lmer.dev lmer2.dev
lme.dev         0       64        61
lmer.dev      369        0       326
lmer2.dev      43        3         0

Другими словами, (1) нет метода, который всегда работает лучше всего; (2) lmerс оптимизатором по умолчанию наихудший (дает сбой примерно в 1/3 времени); (3) lmerс «nloptwrap» лучше (хуже, чем в lme4% случаев, редко хуже lmer).

Чтобы быть немного обнадеживающим, я думаю, что эта ситуация, вероятно, будет наихудшей для небольших, неправильно определенных случаев (т. Е. Остаточная ошибка здесь является равномерной, а не нормальной). Было бы интересно изучить это более систематически, хотя ...