반복 측정에 문제가있는 일부 데이터로 작업하고 있습니다. 사이에 그렇게 나는 매우 다른 행동을 발견 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()
에서 하나와 from에서 두 가지 모델을 실행 했습니다 . 두 경우 모두 나는 각 개인의 반복 측정을 제어하기 위해 ID의 무작위 효과로 키에 대한 체중의 회귀를 실행했습니다.nlme
lmer()
lme4
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
Matches lme
… 무슨 일인지 봅시다. 고정 효과 매개 변수가 프로파일 링 되기 때문에 단일 랜덤 효과를 갖는 LMM에 대한 이탈 함수 (-2 * log 우도) 또는이 경우 유사한 REML 기준 함수는 하나의 인수 만 취합니다 . 주어진 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을. 주어진 방법이 다른 방법 보다 적어도 0.001 이탈 단위 보다 나쁜 답변을 얻는 1000 중 숫자는 다음과 같습니다 .
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”을 사용하는 것이 가장 좋습니다 ( lme
시간의 4 % 미만 , lmer
.
약간 안심할 수 있도록,이 상황은 작고 잘못 지정 된 사례 (즉, 여기서 잔류 오차가 정상이 아닌 균일 한 경우)에 최악이라고 생각합니다. 좀 더 체계적으로 탐색하는 것이 흥미로울 것입니다 …