14.4 Estimation with nlme::lme
nlme에서 linear mixed effects model을 분석할 수 있게 해주는 함수는 lme 함수이다.
lme함수의 default 추정방법은 restricted maximum likelihood (REML) 방법이지만, 그냥 maximum likelihood (ML) 방법이 더 간단한 것이므로 먼저 설명한다.
14.4.1 ML Estimation with nlme::lme
lme함수로 ML 추정을 하려면 method=“ML”이라는 option을 추가해 주어야 한다. (REML을 하려면 없어도 된다.)
lme의 결과를 Orth.ML이라고 이름붙인 object로 저장하고 summary함수로 결과를 보면 다음과 같다.
Orth.ML = lme(distance~age, random=~age|Subject, data=Orthodont, method="ML")
summary(Orth.ML)Linear mixed-effects model fit by maximum likelihood
Data: Orthodont
AIC BIC logLik
451.2116 467.3044 -219.6058
Random effects:
Formula: ~age | Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 2.1940996 (Intr)
age 0.2149244 -0.581
Residual 1.3100399
Fixed effects: distance ~ age
Value Std.Error DF t-value p-value
(Intercept) 16.761111 0.7678975 80 21.827278 0
age 0.660185 0.0705779 80 9.353997 0
Correlation:
(Intr)
age -0.848
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.305969305 -0.487429621 0.007597971 0.482237071 3.922789749
Number of Observations: 108
Number of Groups: 27 위에서 먼저 AIC, BIC, logLik와 같이 minimization 결과로 objective function관련 값들이 어떻게 나왔는지 보여준다.
그 다음은 어떤 random effect 모형을 사용했으며(random intercept와 random slope가 모두 있는 모형이다.), G matrix의 구조가 어떠한지 알려준다. (지정하지 않았으므로 가장 일반적인 general positive-definite matrix가 된다.)
Fixed effect의 결과는 lm의 결과와 유사한데, \(\beta\) 추정값(estimate)의 correlation 도 보여준다.
그 다음 standardized residual의 five number를 보여준다.
마지막으로 관측값 개수(여기서는 108개)와 group 개수(여기서는 사람수 27명)을 보여준다.
다음의 함수로 신뢰구간을 볼 수 있다.
intervals(Orth.ML)Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 15.2471623 16.7611111 18.275060
age 0.5210373 0.6601852 0.799333
Random Effects:
Level: Subject
lower est. upper
sd((Intercept)) 0.83405072 2.1940996 5.7719191
sd(age) 0.09262286 0.2149244 0.4987162
cor((Intercept),age) -0.94284432 -0.5814877 0.4083670
Within-group standard error:
lower est. upper
1.084786 1.310040 1.582068 Fiexd effect(\(\beta\) 추정값)뿐 아니라 random effect(G matrix 추정값)의 신뢰구간을 보여준다. Fixed effect의 신뢰구간을 거의 항상 나오지만, random effect의 신뢰구간은 경우에 따라 못 구할 때도 있다.
마지막으로 G matrix 추정값들의 Hessian matrix는 다음과 같이 볼 수 있다.
solve(Orth.ML$apVar) # hessian matrix for random effects estimates reStruct.Subject1 reStruct.Subject2 reStruct.Subject3 lSigma
reStruct.Subject1 18.151065 1.597546 7.2418501 8.3681397
reStruct.Subject2 1.597546 27.740842 9.8850963 6.8305611
reStruct.Subject3 7.241850 9.885096 6.9565436 0.6196145
lSigma 8.368140 6.830561 0.6196145 136.5142266위 결과에서 re가 붙은 것은 random effect의 약어이며, Sigma 앞의 소문자 l은 log의 약어이다. 목적함수 (주로 -log likelihood 함수 또는 그것의 변형)를 minimization하기 위해서는 변수변환(transformation)을 하여야 하는데(이를 reparametrization이라고도 부른다), 이러한 변환을 하는 방법에는 여러가지가 있으며, nlme는 sigma의 경우 log변환을 한다.
SAS의 경우에는 G, R, V matrix를 간단한 option으로 보여주지만, nlme에서는 getVarCov 함수를 이용해서 봐야 한다. 이때 type option을 적절히 지정해주면 된다.
getVarCov(Orth.ML, type="random.effects") # G matrix, var-cov of random effectRandom effects variance covariance matrix
(Intercept) age
(Intercept) 4.81410 -0.274210
age -0.27421 0.046193
Standard Deviations: 2.1941 0.21492 getVarCov(Orth.ML, type="conditional") # R matrix, conditional var-cov of YSubject M01
Conditional variance covariance matrix
1 2 3 4
1 1.7162 0.0000 0.0000 0.0000
2 0.0000 1.7162 0.0000 0.0000
3 0.0000 0.0000 1.7162 0.0000
4 0.0000 0.0000 0.0000 1.7162
Standard Deviations: 1.31 1.31 1.31 1.31 getVarCov(Orth.ML, type="marginal") # V matrix, marginal var-cov of YSubject M01
Marginal variance covariance matrix
1 2 3 4
1 5.0992 3.5737 3.7644 3.9550
2 3.5737 5.6653 4.3246 4.7000
3 3.7644 4.3246 6.6010 5.4450
4 3.9550 4.7000 5.4450 7.9061
Standard Deviations: 2.2582 2.3802 2.5692 2.8118 VarCorr함수를 이용하여 G matrix의 correlation 값을 볼 수도 있다.
VarCorr(Orth.ML) # Correlation of G matrixSubject = pdLogChol(age)
Variance StdDev Corr
(Intercept) 4.81407327 2.1940996 (Intr)
age 0.04619252 0.2149244 -0.581
Residual 1.71620466 1.3100399 추정하는 과정에서 \(\beta\)값이 산출되며, 확률변수 u의 현실화된 값(realized value, EBE, MAP)도 나온다.
추정된 \(\beta\)값은 fixef 함수로 볼 수 있고, u의 현실화된 값(realized value)은 ranef함수로 볼 수 있다.
fixef(Orth.ML)(Intercept) age
16.7611111 0.6601852 ranef(Orth.ML)[IDs, ] (Intercept) age
M01 1.0712996 0.212833572
M02 -0.6819769 0.010649777
M03 -0.1513173 0.033776865
M04 2.8105092 -0.049757544
M05 -1.1022642 0.019284244
M06 1.1802856 0.085821300
M07 -0.5616968 0.030864762
M08 0.7784275 -0.087411200
M09 -0.3706756 0.129027582
M10 2.5826295 0.215813570
M11 0.7982431 -0.110504341
M12 -0.9013352 0.105900494
M13 -3.7514030 0.379970285
M14 0.8491686 -0.009463365
M15 -0.4301225 0.198307003
M16 -0.2022428 -0.067264111
F01 -0.5234520 -0.174095346
F02 -0.9522607 0.004859518
F03 -0.7117004 0.045289488
F04 0.9991722 -0.023888091
F05 0.4274949 -0.159602725
F06 -0.6536399 -0.182763760
F07 -0.2022428 -0.067264111
F08 1.1180659 -0.162446933
F09 -0.3536328 -0.211613212
F10 -2.2456187 -0.252145024
F11 1.1802856 0.085821300개인별 intercept와 slope는 위의 fixef와 ranef의 조합(여기서는 덧셈)으로 알 수 있는데, coef 함수로 볼 수 있다. 그리고, 개인별 intercept와 slope를 이용하여 fitted graph를 그리자면 augPred 함수를 사용한 뒤 plot하면 된다.
plot(augPred(Orth.ML))
Figure 14.2: Augmented(Individual) prediction plot by ML estimation
coef(Orth.ML)[IDs, ] (Intercept) age
M01 17.83241 0.8730188
M02 16.07913 0.6708350
M03 16.60979 0.6939621
M04 19.57162 0.6104276
M05 15.65885 0.6794694
M06 17.94140 0.7460065
M07 16.19941 0.6910499
M08 17.53954 0.5727740
M09 16.39044 0.7892128
M10 19.34374 0.8759988
M11 17.55935 0.5496808
M12 15.85978 0.7660857
M13 13.00971 1.0401555
M14 17.61028 0.6507218
M15 16.33099 0.8584922
M16 16.55887 0.5929211
F01 16.23766 0.4860898
F02 15.80885 0.6650447
F03 16.04941 0.7054747
F04 17.76028 0.6362971
F05 17.18861 0.5005825
F06 16.10747 0.4774214
F07 16.55887 0.5929211
F08 17.87918 0.4977383
F09 16.40748 0.4485720
F10 14.51549 0.4080402
F11 17.94140 0.7460065이 그림이나 개인별 intercept와 slope는 standard two-stage때와 유사하지만, 다르다. 어떤 것이 더 진실에 가까울 지는 경우에 따라 다르므로, 자료를 분석하는 연구자는 심사숙고해야 한다. 다만 standard two-stage method에 비해 mixed effects model이 추정한 parameter갯수는 훨씬 적다고 본다. 계산하자면 standard two-stage에서는 한 사람당 2개씩 27명에 대해 추정하였으므로 54개의 fixed effect parameter를 추정하였으며(그 만큼의 자유도 감소가 일어났다. 단, 이 경우 분산 추정은 \(\beta\)추정 뒤에 일의적으로 결정되는 것으로 자유도를 감소시키지 않는다.), mixed effects model에서는 \(\beta\) 2개와 G matrix의 원소 3개와 R matrix의 원소 1개를 추정하였으므로 fixed effect parameter 2개와 random effect paramter 4개를 추정한 것이고, 총 6개를 추정한 것이다. (추정된 EBE는 앞의 6개 파라미터가 정해지면 역시 일의적으로 결정되는 것이어서 자유도를 감소시키지 않는다.) 그런데, 대부분의 경우 더 정확히 추정해야 하는 것은 fixed effect (\(\beta\))이며, random effect (G matrix)는 fixed effect를 더 잘 추정하는데 도움이 된다면 자신의 유의성 여부(신뢰구간이 0을 포함하는지 여부)와 상관없이 추가하는 것을 바람직하다고 본다.
14.4.2 REML Estimation with nlme::lme
REML의 경우에는 분산의 크기가 과소추정되는 것을 방지하기 위해서 objective function만 더 복잡하게 되었는데, fixed effect (\(\beta\))의 추정값은 동일하고, 분산 term들의 추정값이 조금씩 더 크게 되었다.
아래에서 method=“REML”은 생략가능하며, 나머지 함수들의 사용볍과 해석방법은 동일하다.
Orth.REML = lme(distance~age, random=~age|Subject, data=Orthodont, method="REML")
summary(Orth.REML)Linear mixed-effects model fit by REML
Data: Orthodont
AIC BIC logLik
454.6367 470.6173 -221.3183
Random effects:
Formula: ~age | Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 2.3270341 (Intr)
age 0.2264278 -0.609
Residual 1.3100397
Fixed effects: distance ~ age
Value Std.Error DF t-value p-value
(Intercept) 16.761111 0.7752460 80 21.620377 0
age 0.660185 0.0712533 80 9.265333 0
Correlation:
(Intr)
age -0.848
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.223106060 -0.493761143 0.007316631 0.472151118 3.916033216
Number of Observations: 108
Number of Groups: 27 intervals(Orth.REML)Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 15.2183224 16.7611111 18.3038999
age 0.5183867 0.6601852 0.8019837
Random Effects:
Level: Subject
lower est. upper
sd((Intercept)) 0.9485605 2.3270341 5.7087424
sd(age) 0.1025090 0.2264278 0.5001467
cor((Intercept),age) -0.9382505 -0.6093329 0.2981686
Within-group standard error:
lower est. upper
1.084977 1.310040 1.581788 solve(Orth.REML$apVar) # hessian matrix for random effects estimates reStruct.Subject1 reStruct.Subject2 reStruct.Subject3 lSigma
reStruct.Subject1 19.3769133 0.2681486 7.1856664 8.3494422
reStruct.Subject2 0.2681486 28.4894419 9.6908977 6.7147935
reStruct.Subject3 7.1856664 9.6908977 7.0001137 0.5879244
lSigma 8.3494422 6.7147935 0.5879244 133.4748628getVarCov(Orth.REML, type="random.effects") # G matrixRandom effects variance covariance matrix
(Intercept) age
(Intercept) 5.41510 -0.32106
age -0.32106 0.05127
Standard Deviations: 2.327 0.22643 getVarCov(Orth.REML, type="conditional") # R matrixSubject M01
Conditional variance covariance matrix
1 2 3 4
1 1.7162 0.0000 0.0000 0.0000
2 0.0000 1.7162 0.0000 0.0000
3 0.0000 0.0000 1.7162 0.0000
4 0.0000 0.0000 0.0000 1.7162
Standard Deviations: 1.31 1.31 1.31 1.31 getVarCov(Orth.REML, type="marginal") # V matrixSubject M01
Marginal variance covariance matrix
1 2 3 4
1 5.2756 3.7376 3.9158 4.0939
2 3.7376 5.8370 4.5041 4.8874
3 3.9158 4.5041 6.8087 5.6808
4 4.0939 4.8874 5.6808 8.1904
Standard Deviations: 2.2969 2.416 2.6093 2.8619 VarCorr(Orth.REML)Subject = pdLogChol(age)
Variance StdDev Corr
(Intercept) 5.41508758 2.3270341 (Intr)
age 0.05126955 0.2264278 -0.609
Residual 1.71620400 1.3100397 R, V matrix는 개인(subject)별로 다를 수 있다. (예를 들어, 측정 횟수가 현재는 개인별로 모두 4개로 동일하지만, 실제 상황에서는 1 - 4개로 다양할 수 있다. 이 경우 X와 Z matrix는 개인별로 달라진다.) 따라서, 따로 지정하지 않으면 첫 번째 subject의 R과 V가 출력된다.
fixef(Orth.REML)(Intercept) age
16.7611111 0.6601852 ranef(Orth.REML)[IDs, ] (Intercept) age
M01 1.0515830 0.215684556
M02 -0.7275013 0.014507808
M03 -0.1739024 0.035852350
M04 3.0004527 -0.065884768
M05 -1.1766673 0.025600299
M06 1.2176432 0.083191283
M07 -0.6079720 0.034899998
M08 0.8602971 -0.094736217
M09 -0.4443954 0.135908591
M10 2.6532361 0.211146620
M11 0.8904899 -0.118825903
M12 -0.9979943 0.114564049
M13 -4.1295437 0.413668503
M14 0.9043445 -0.014119813
M15 -0.5349736 0.208177648
M16 -0.1877570 -0.068853740
F01 -0.4859593 -0.178209679
F02 -1.0118489 0.009857959
F03 -0.7727904 0.050642337
F04 1.0691629 -0.029862153
F05 0.5168057 -0.167957627
F06 -0.6205849 -0.186557025
F07 -0.1877570 -0.068853740
F08 1.2503194 -0.174400268
F09 -0.2909481 -0.218041704
F10 -2.2813817 -0.250590650
F11 1.2176432 0.083191283plot(augPred(Orth.REML))
Figure 14.3: Augmented(Individual) prediction plot by REML estimation
coef(Orth.REML)[IDs, ] (Intercept) age
M01 17.81269 0.8758697
M02 16.03361 0.6746930
M03 16.58721 0.6960375
M04 19.76156 0.5943004
M05 15.58444 0.6857855
M06 17.97875 0.7433765
M07 16.15314 0.6950852
M08 17.62141 0.5654490
M09 16.31672 0.7960938
M10 19.41435 0.8713318
M11 17.65160 0.5413593
M12 15.76312 0.7747492
M13 12.63157 1.0738537
M14 17.66546 0.6460654
M15 16.22614 0.8683628
M16 16.57335 0.5913314
F01 16.27515 0.4819755
F02 15.74926 0.6700431
F03 15.98832 0.7108275
F04 17.83027 0.6303230
F05 17.27792 0.4922276
F06 16.14053 0.4736282
F07 16.57335 0.5913314
F08 18.01143 0.4857849
F09 16.47016 0.4421435
F10 14.47973 0.4095945
F11 17.97875 0.7433765여기까지의 사용법이나 해석방법은 ML방법 때와 동일하다.
14.4.3 SAS Type 3 Tests of Fixed Effects
SAS 결과에는 Type 3 tests of fixed effects 부분이 나오지만, nlme에서는 나오지 않는다.
만약 SAS와 같은 Type 3 test결과를 보고 싶으면 다음과 같이 sasLM함수를 사용한다.
일반적으로 많이 사용하는 anova, aov, car::Anova 함수로는 동일한 결과가 나오지 않는 경우가 많다.
require(sasLM)
T3test(distance ~ age + Subject + age:Subject, Orthodont, "age:Subject")$age
Error: 1*age:Subject + 0*MSE
Df Sum Sq Mean Sq F value Pr(>F)
age 1 235.356 235.356 85.846 1.013e-09 ***
Error 26 71.281 2.742
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`The Rest Terms`
Df Sum Sq Mean Sq F value Pr(>F)
Subject 26 66.969 2.5757 1.5008 0.10395
age:Subject 26 71.281 2.7416 1.5975 0.07345 .
RESIDUALS 54 92.675 1.7162
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1위에서 마지막 옵션은 그 term을 분모로하여 test할 수 있는 것은 그렇게 하라는 것이다.
aov3(distance ~ age + Subject + age:Subject, Orthodont)Response : distance
Df Sum Sq Mean Sq F value Pr(>F)
MODEL 53 825.02 15.566 9.0702 6.573e-14 ***
age 1 235.36 235.356 137.1376 2.220e-16 ***
Subject 26 66.97 2.576 1.5008 0.10395
age:Subject 26 71.28 2.742 1.5975 0.07345 .
RESIDUALS 54 92.67 1.716
CORRECTED TOTAL 107 917.69
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1다음은 R 내장 함수인 anova와 aov를 사용한 결과이다.
anova(Orth.ML) numDF denDF F-value p-value
(Intercept) 1 80 3156.0380 <.0001
age 1 80 87.4973 <.0001anova(Orth.REML) numDF denDF F-value p-value
(Intercept) 1 80 3096.4871 <.0001
age 1 80 85.8464 <.0001summary(aov(distance ~ age + Error(age:Subject), Orthodont))
Error: age:Subject
Df Sum Sq Mean Sq F value Pr(>F)
age 1 235.4 235.4 11.71 0.00207 **
Residuals 26 522.7 20.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 80 159.6 1.996 summary(aov(distance ~ age + Subject + Error(age:Subject), Orthodont))
Error: age:Subject
Df Sum Sq Mean Sq
age 1 235.4 235.4
Subject 26 522.7 20.1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Subject 26 66.97 2.576 1.501 0.104
Residuals 54 92.67 1.716 car::Anova 함수도 별로 만족스럽지 못하다.
require(car)
options(contrasts = c("contr.sum", "contr.poly"))
Anova(Orth.ML, type=3, test.statistic="F", singular.ok=TRUE)Analysis of Deviance Table (Type III tests)
Response: distance
Chisq Df Pr(>Chisq)
(Intercept) 485.419 1 < 2.2e-16 ***
age 89.148 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(Orth.REML, type=3, test.statistic="F", singular.ok=TRUE)Analysis of Deviance Table (Type III tests)
Response: distance
Chisq Df Pr(>Chisq)
(Intercept) 467.441 1 < 2.2e-16 ***
age 85.846 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1options(contrasts = c("contr.sum", "contr.poly"))
r1 = lm(distance ~ age + Subject + age:Subject, Orthodont)
Anova(r1, type=3, test.statistic="F", singular.ok=TRUE)Anova Table (Type III tests)
Response: distance
Sum Sq Df F value Pr(>F)
(Intercept) 1204.01 1 701.5522 < 2e-16 ***
age 235.36 1 137.1376 < 2e-16 ***
Subject 66.97 26 1.5008 0.10395
age:Subject 71.28 26 1.5975 0.07345 .
Residuals 92.67 54
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1