반복 측정에 문제가있는 일부 데이터로 작업하고 있습니다. 사이에 그렇게 나는 매우 다른 행동을 발견 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의 무작위 효과로 키에 대한 체중의 회귀를 실행했습니다.nlmelmer()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()

왜이 모델들이 그렇게 많이 다른가?

답변

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

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

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()

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