14.5 Estimation with Custom Objective Function

이제 사용자가 만든 objective function을 이용하여 추정을 하면 다음과 같다.

14.5.1 Estimation with Custom ML Objective Function

ML 방법이 조금 더 간단하므로 먼저 제시한다.

14.5.1.1 Initialization

먼저 자료를 적합한 형태로 준비한다.

nRec = nrow(Orthodont)
IDs = unique(Orthodont$Subject) ; IDs
 [1] M01 M02 M03 M04 M05 M06 M07 M08 M09 M10 M11 M12 M13 M14 M15 M16 F01 F02 F03 F04 F05
[22] F06 F07 F08 F09 F10 F11
27 Levels: M16 < M05 < M02 < M11 < M07 < M08 < M03 < M12 < M13 < M14 < M09 < ... < F11
nID = length(IDs) ; nID
[1] 27

14.5.1.2 Design matrices

설계행렬은 다음과 같이 구할 수 있다.

Y = Orthodont[,"distance"]
X = model.matrix(distance ~ age, Orthodont)
Z = X
head(X)
  (Intercept) age
1           1   8
2           1  10
3           1  12
4           1  14
5           1   8
6           1  10
tail(X)
    (Intercept) age
103           1  12
104           1  14
105           1   8
106           1  10
107           1  12
108           1  14
nParF = ncol(X) # number of beta parameters for fixed effect
nEta = ncol(Z)  # 3 Omega (G matrix)
nEps = 1        # 1 Sigma (R matrix)
nParR = nEta*(nEta + 1)/2 + nEps # number of parameters for random effects
nParA = nParF + nParR # number of all parameters

AIC, BIC 등을 구하기 위해 parameter 개수를 count하였다. 그런데, 이 방식이 SAS와 nlme가 일부 다르다.

14.5.1.3 Record indices

먼저 계산의 편의를 위해 각 subject별 dataset을 list를 이용하여 구분하여 준다. 또한 개인별 V matrix와 inverse V matrix를 저장할 공간을 미리 마련해 두었다.

ind = list()
Vi = list()
iVi = list()
iRec = vector(length=nID)
for (i in 1:nID) {
  ind[[i]] = Orthodont$Subject == IDs[i]
  iRec[i] = sum(ind[[i]])
  Vi[[i]] = matrix(NA, nrow=iRec[i], ncol=iRec[i])
  iVi[[i]] = matrix(NA, nrow=iRec[i], ncol=iRec[i])
}

14.5.1.4 Matrix for EBE (bi)

개인별 realized u를 저장할 공간을 bi 라는 이름의 matrix로 준비해 둔다. 한 사람당 한 행을 차지한다.

bi = matrix(nrow=nID, ncol=ncol(Z))
rownames(bi) = IDs

14.5.1.5 Vector for individual OFV

개인별 목적함수 값(objective function value, OFV)를 저장할 vector도 마련해 둔다.

Oi = rep(NA, nID)

전체 objective function value (Total OFV)는 이 개인별 OFV의 단순합이다.

14.5.1.6 Custom objective function for ML

다음은 R로 작성한 ML method용 objective function이다. 수식을 참조한 reference는 source code 내에 포함되어 있다. (수식 포함 예정)

ObjML = function(TH)
{
  eTH = exp(TH)
  rho = (eTH[3] - 1)/(eTH[3] + 1)
  Cov = rho*eTH[1]*eTH[2]
  OM = matrix(c(eTH[1]*eTH[1], Cov, Cov, eTH[2]*eTH[2]), nrow=2)
  SG = eTH[4]*eTH[4]        # R matrix

  for (i in 1:nID) {
    Zi = Z[ind[[i]],]
    Vi[[i]] <<- Zi %*% OM %*% t(Zi) + diag(SG, nrow=iRec[i])
    iVi[[i]] <<- solve(Vi[[i]])
  }

  S1 = matrix(rep(0, 4), nrow=2)
  S2 = matrix(rep(0, 2), nrow=2)
  for (i in 1:nID) {
    XiV = t(X[ind[[i]],]) %*% iVi[[i]]
    S1 = S1 + XiV %*% X[ind[[i]],]
    S2 = S2 + XiV %*% Y[ind[[i]]]
  }
  FIM <<- S1                # Davidian p78 eq 3.15
  Beta <<- solve(S1) %*% S2 # Davidian p78 eq 3.13, Bonate 2e p195 eq 35

  for (i in 1:nID) {
    qi = Y[ind[[i]]] - X[ind[[i]],] %*% Beta
    iVq = iVi[[i]] %*% qi
    bi[i,] <<- OM %*% t(Z[ind[[i]],]) %*%  iVq
    # Linear Mixed-Effects Models Using R p253 or Bonate p195 or Davidian p78
    Oi[i] <<- determinant(Vi[[i]], logarithm=TRUE)$modulus[[1]] + t(qi) %*% iVq
  }
  OFV = (nRec*log(2*pi) + sum(Oi))/2 # minus loglikelihood
  return(OFV)
}

위 함수의 첫 5줄은 transformation (reparametrization)을 위한 것인데, nlme와 동일한 결과를 쉽게 얻기 위해 동일한 변환 방법을 사용했다. 변환 방법에 대한 좀 더 자세한 설명이 필요하면 Bates의 책을 참고한다.

특이한 점은 ObjML 함수의 argument인 TH에는 G matrix와 sigma만 들어있다. 또한, G matrix는 개개의 요소가 그대로 있는 것이 아니고, 공분산 대신 상관계수로 parametrization되어 있다. \(\beta\) 는 도중에 행렬계산으로 바로 계산한다.

Minimization에는 불필요하지만, Total OFV가 온전한 -log likelihood이 되도록 nReclog(2pi) 가 더해져 있고, 2로도 나누어져 있다.

14.5.1.7 Minimization using optim

앞의 목적함수를 최적화(최소화)하는 데는 optim함수가 편리하다.

t1 = Sys.time()
r1 = optim(c(0.1, 0.1, 0.1, 0.1), ObjML, method="L-BFGS-B", hessian=T)
Sys.time() - t1
Time difference of 0.3034739 secs

수행에 걸린 시간을 확인하기 위해 Sys.time()을 두 번 사용하였다.

14.5.1.8 Minimization result

최소화 결과는 다음과 같다.

r1 # notice $value as -loglihood
$par
[1]  0.7857364 -1.5374861 -1.3293299  0.2700672

$value
[1] 219.6058

$counts
function gradient 
      30       30 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$hessian
          [,1]      [,2]      [,3]        [,4]
[1,] 18.149266  1.601628 7.2382947   8.3685844
[2,]  1.601628 27.740057 9.8832865   6.8324504
[3,]  7.238295  9.883287 6.9589462   0.6108312
[4,]  8.368584  6.832450 0.6108312 136.5027320

먼저 $par reparametrization된 상태에서 추정 값 4개(G matrix 요소 3개, sigma 1개)가 나온다.

두 번째인 $value 가total OFV이다.

$count는 iteration 횟수이고, $covergence는 minimization 성공여부를 알려준다.

$message는 mimization을 종료할 때 나온 종료 사유를 나타내는 message이다.

더 자세한 설명은 optim함수를 참고한다.

마지막의 $hessian이 추정된 parameter (여기서는 G의 요소들과 sigma)의 hessian으로서, standard error를 구하는데 유용하게 사용된다.

앞의 solve(Orth.ML$apVar)와 거의 같은 값이어야 한다. 동일하지는 않은 이유는 minization algorithm이 nlme와 다르기 때문이다. 이는 비선형 함수의 최적화 과정에서 흔히 보는 현상이다.

14.5.1.9 -2 log likelihood (SAS output)

OFV에 2를 곱해주면 된다.

m2LL = 2*r1$value ; m2LL
[1] 439.2116

ML에서 AIC와 AICc는 SAS나 nlme나 같은 계산 방식을 따른다.

14.5.1.10 AIC

2*nParA + m2LL # nlme style = SAS style for ML
[1] 451.2116

14.5.1.11 AICc (SAS output)

2*nParA + m2LL + 2*nParA*(nParA + 1)/(nRec - nParA - 1) # nlme = SAS for ML
[1] 452.0433

14.5.1.12 BIC

하지만, BIC의 계산 방식이 다음과 같이 조금 다르다.

nParA*log(nRec) + m2LL # nlme style for ML
[1] 467.3044
nParA*log(nID) + m2LL  # SAS style for ML
[1] 458.9866

14.5.1.13 logLik, log likelihood, log lihood

log likelihood는 OFV의 부호만 바꾸어 주면 된다.

-r1$value
[1] -219.6058

14.5.1.14 Estimates of random effects

Reparametrization된 것을 original scale로 바꾸기 위해서는 objective function 첫 부분의 계산을 역으로 해주면 된다.

SDs = exp(r1$par)[-3] ; SDs
[1] 2.1940220 0.2149207 1.3100524
rho = (exp(r1$par)[3] - 1)/(exp(r1$par)[3] + 1) ; rho
[1] -0.5814595

SAS나 NONMEM의 경우에는 variance를 그대로 사용하여 minization 하지만, nlme는 standard deviation(sd, \(\sigma\))과 correlation(\(\rho\)) scale에서 계산한다.

14.5.1.15 Estimates of fixed effects

\(\beta\) 추정값의 신뢰구간을 구하기 위해서는 추정값의 분산(에 대한 추정값)을 알아야 한다.

Var1 = solve(FIM)                    # SAS, NONMEM style
Var2 = Var1*nRec/(nRec - qr(X)$rank) # nlme style

그런데, 위에서 보는 것과 같이 SAS와 nlme가 다르다. SAS는 ML이론에 충실하게 추정된 Fisher’s Information Matrix(FIM)을 그대로 사용하지만, nlme는 REML과 유사하게 variance의 bias를 correction해서 약간 더 크게 산출한다.

나머지 계산은 단순선형회귀 때와 동일하다.

SE = sqrt(diag(Var2))  # Note: This used Var2 not Var1
Df = nRec - qr(X)$rank - (nID - 1) # Df 80 # nlme
tVal = Beta/SE
fixedPE = cbind(Beta, SE, Df, tVal, 2*(1 - pt(abs(tVal), Df)))
colnames(fixedPE) = c("Value", "Std.Error", "DF", "t-value", "p-value")
fixedPE
                 Value  Std.Error DF   t-value      p-value
(Intercept) 16.7611111 0.76683287 80 21.857581 0.000000e+00
age          0.6601852 0.07048624 80  9.366157 1.687539e-14
cov2cor(Var2)[lower.tri(Var2)] # correlation between intercept and age estimates
[1] -0.8477951

14.5.1.16 SAS Output

SAS의 경우에는 \(\beta\) 추정값의 분산이 조금 다르므로 다음과 같이 계산한다.

SE1 = sqrt(diag(Var1)) # SAS OUTPUT
Df1 = nID - 1 # Df1 26 # SAS
tVal1 = Beta/SE1
fixedSol = cbind(Beta, SE1, Df1, tVal1, 2*(1 - pt(abs(tVal1), Df1)))
colnames(fixedSol) = c("Value", "Std.Error", "DF", "t-value", "p-value")
fixedSol
                 Value  Std.Error DF   t-value      p-value
(Intercept) 16.7611111 0.75969938 26 22.062821 0.000000e+00
age          0.6601852 0.06983054 26  9.454104 6.733913e-10

14.5.1.17 Standardized within-group residuals

잔차를 구하기 위해서는 individual prediction값이 필요하고, 이를 위해서는 또 개인별 intercept와 slope가 필요하다. 개인간의 interindividual random variability를 고려하지 않은 (즉, u 또는 bi를 모두 0으로 고정한) 경우의 예측값을 typical prediction(전형적인 예측)이라 부르기로 하자. 다음은 이를 계산하는 과정이다.

d1 = Orthodont
d1$PRED = NA
d1$IPRE = NA
for (i in 1:nID) {
  d1[ind[[i]], "PRED"] =  X[ind[[i]],] %*% Beta  # Typical prediction
  d1[ind[[i]], "IPRE"] =  X[ind[[i]],] %*% Beta + Z[ind[[i]],] %*% bi[i,]
}
d1$RES = d1$distance - d1$PRED
d1$IRES = d1$distance - d1$IPRE
d1$IWRE = d1$IRES/SDs[3]
quantile(d1$IWRE, c(0, 0.25, 0.5, 0.75, 1)) # fivenum(d1$IWRE)
          0%          25%          50%          75%         100% 
-3.303582775 -0.487855750  0.007694123  0.481942430  3.922339927

위에서 PRED에 typical prediction이 담겨있고, IPRE에 individual prediction 값이 담겨있다. 따라서, residual도 두 가지 종류가 발생하며, 후자에 의한 residual이 individual residual이다. 이를 다시 자신의 sigma로 나누어 준 것을 nlme는 standardized within-group residual이라 부른다.

14.5.1.18 Intervals of fixed effects

\(\beta\)의 신뢰구간은 다음과 같다.

alpha = 0.05
Len = qt(1 - alpha/2, Df)*sqrt(diag(Var1)) # Note: This used Var1
CIbeta = cbind(Beta - Len, Beta, Beta + Len)
colnames(CIbeta) = c("lower", "est", "upper")
CIbeta
                lower        est      upper
(Intercept) 15.249261 16.7611111 18.2729611
age          0.521218  0.6601852  0.7991524

14.5.1.19 Intervals of random effects

분산 term들의 신뢰구간을 구하는 것은 reparametrization이 되었기 때문에 좀 더 복잡한데, 다행히 변환에 따른 계산 수식들이 알려져 있다. 만약 이 변환에 따라 추정값들의 variance-covariance matrix가 어떻게 변하는지 알려져 있지 않다면 (NONMEM의 경우), final estimate에서 original scale로 hessian matrix만 다시 구하여 standard error를 계산하여야 한다.

vPE = r1$par
Var = solve(r1$hessian)
vB = qnorm(1 - alpha/2)*sqrt(diag(Var))

vLL = exp(vPE[1:2] - vB[1:2])
vUL = exp(vPE[1:2] + vB[1:2])

rhoLL = (exp(vPE[3] - vB[3]) - 1)/(exp(vPE[3] - vB[3]) + 1) # rho LL
rhoUL = (exp(vPE[3] + vB[3]) - 1)/(exp(vPE[3] + vB[3]) + 1) # rho UL

CIomega = cbind(c(vLL, rhoLL), c(SDs[1:2], rho), c(vUL, rhoUL))
rownames(CIomega) = c("sd((Intercept))", "sd(age)", "cor((Intercept),age)")
colnames(CIomega) = c("lower", "est.", "upper")

CIsigma = SDs[3]*exp(c(-1, 0, 1)*vB[4])
names(CIsigma) = c("lower", "est.", "upper")
  • Approximate confidence interval of random effects: Subject level
CIomega
                           lower       est.     upper
sd((Intercept))       0.83680006  2.1940220 5.7525479
sd(age)               0.09288575  0.2149207 0.4972873
cor((Intercept),age) -0.94235079 -0.5814595 0.4047436
  • Approximate confidence interval of random effects: Within group
CIsigma
   lower     est.    upper 
1.084877 1.310052 1.581965

14.5.1.20 Varaiance-Covariance Matrix

SAS의 경우는 앞과 달리 G, R, V matrix를 그대로 보여준다. SAS style로 변환하는 과정은 다음과 같다.

Cov = rho*SDs[1]*SDs[2]
OM = matrix(c(SDs[1]*SDs[1], Cov, Cov, SDs[2]*SDs[2]), nrow=2) ; OM # G matrix
           [,1]        [,2]
[1,]  4.8137325 -0.27418187
[2,] -0.2741819  0.04619091
SG = SDs[3]*SDs[3] ; SG     # R matrix
[1] 1.716237
Z[ind[[1]],] %*% OM %*% t(Z[ind[[1]],]) + diag(nrow(Z[ind[[1]],]))*SG # V matrix
         1        2        3        4
1 5.099278 3.573732 3.764422 3.955113
2 3.573732 5.665423 4.324641 4.700095
3 3.764422 4.324641 6.601096 5.445077
4 3.955113 4.700095 5.445077 7.906296

R matrix는 가장 단순한 모형을 사용하였으므로 개인별 차이가 없고, V matrix는 (역시 개인별 측정횟수가 동일해서 차이가 없지만) 첫 번째 subject것만 출력하였다.

14.5.1.21 EBE, MAP, realized Eta, post-hoc Eta

EBE는 minimization과정에서 매 iteration마다 update하여 bi에 저장하였으므로 그 값을 보면 된다.

bi
          [,1]         [,2]
M01  1.0710140  0.212967314
M02 -0.6833646  0.010755330
M03 -0.1519626  0.033841277
M04  2.8162911 -0.050200198
M05 -1.1045294  0.019457550
M06  1.1815675  0.085774658
M07 -0.5630770  0.030980528
M08  0.7808197 -0.087629880
M09 -0.3727535  0.129245766
M10  2.5851336  0.215748081
M11  0.8009204 -0.110755818
M12 -0.9041555  0.106159819
M13 -3.7624882  0.380970281
M14  0.8508557 -0.009589840
M15 -0.4330554  0.198623577
M16 -0.2018980 -0.067324701
F01 -0.5225595 -0.174252168
F02 -0.9540908  0.004993841
F03 -0.7135156  0.045444236
F04  1.0012943 -0.024053548
F05  0.4300073 -0.159868441
F06 -0.6528974 -0.182914397
F07 -0.2018980 -0.067324701
F08  1.1218981 -0.162809172
F09 -0.3520203 -0.211841814
F10 -2.2471034 -0.252172236
F11  1.1815675  0.085774658

14.5.1.22 Coefficients of Individuals

개인별 intercept와 slope는 fixed effect(\(\beta\))와 bi를 조합(여기서는 더하기)하여 구할 수 있다.

iCoef = bi
for (i in 1:nID) iCoef[i, ] = iCoef[i,] + t(Beta)
iCoef
        [,1]      [,2]
M01 17.83213 0.8731525
M02 16.07775 0.6709405
M03 16.60915 0.6940265
M04 19.57740 0.6099850
M05 15.65658 0.6796427
M06 17.94268 0.7459598
M07 16.19803 0.6911657
M08 17.54193 0.5725553
M09 16.38836 0.7894310
M10 19.34624 0.8759333
M11 17.56203 0.5494294
M12 15.85696 0.7663450
M13 12.99862 1.0411555
M14 17.61197 0.6505953
M15 16.32806 0.8588088
M16 16.55921 0.5928605
F01 16.23855 0.4859330
F02 15.80702 0.6651790
F03 16.04760 0.7056294
F04 17.76241 0.6361316
F05 17.19112 0.5003167
F06 16.10821 0.4772708
F07 16.55921 0.5928605
F08 17.88301 0.4973760
F09 16.40909 0.4483434
F10 14.51401 0.4080129
F11 17.94268 0.7459598

14.5.2 Estimation with Custom REML Objective Function

REML방법의 경우에는objective function이 다르다. 그런데, AIC, AICc, BIC 등의 계산법이 SAS와 nlme간에 차이가 난다.

14.5.2.1 Custom Objective Function for REML

아래와 같이 약간 더 복잡하게 바뀌었음을 알 수 있다.

ObjREML = function(TH)
{
  eTH = exp(TH)
  rho = (eTH[3] - 1)/(eTH[3] + 1)
  Cov = rho*eTH[1]*eTH[2]
  OM = matrix(c(eTH[1]*eTH[1], Cov, Cov, eTH[2]*eTH[2]), nrow=2)
  SG = eTH[4]*eTH[4]        # R matrix

  for (i in 1:nID) {
    Zi = Z[ind[[i]],]
    Vi[[i]] <<- Zi %*% OM %*% t(Zi) + diag(SG, nrow=iRec[i])
    iVi[[i]] <<- solve(Vi[[i]])
  }

  S1 = matrix(rep(0, 4), nrow=2)
  S2 = matrix(rep(0, 2), nrow=2)
  S3 = 0
  for (i in 1:nID) {
    XiV = t(X[ind[[i]],]) %*% iVi[[i]]
    S1 = S1 + XiV %*% X[ind[[i]],]
    S2 = S2 + XiV %*% Y[ind[[i]]]
    S3 = S3 + determinant(Vi[[i]], logarithm=TRUE)$modulus[[1]]
  }
  FIM <<- S1
  Beta <<- solve(S1) %*% S2

  S4 = 0
  for (i in 1:nID) {
    qi = Y[ind[[i]]] - X[ind[[i]],] %*% Beta
    iVq = iVi[[i]] %*% qi
    bi[i,] <<- OM %*% t(Z[ind[[i]],]) %*% iVq
    # Linear Mixed-Effects Models Using R p253 or Bonate p195 or Davidian p78
    S4 = S4 + t(qi) %*% iVq
  }
  S0 = (nRec - nParF)*log(2*pi) # Constant part
  OFV = (S0 + S3 + S4 + determinant(S1, logarithm=TRUE)$modulus[[1]])/2
  return(OFV) # minus loglikelihood
}

수식은 위 R script내에 표기한 문헌을 참고한다.

14.5.2.2 Minimization using optim or nlminb

최적화 과정은 ML의 경우와 동일하다.

t2 = Sys.time()
r2 = optim(c(0.1, 0.1, 0.1, 0.1), ObjREML, method="L-BFGS-B", hessian=T)
Sys.time() - t2
Time difference of 0.2967291 secs

14.5.2.3 Minimization result

r2 # notice $value as -loglihood
$par
[1]  0.8445663 -1.4854348 -1.4155778  0.2700548

$value
[1] 221.3183

$counts
function gradient 
      31       31 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$hessian
           [,1]       [,2]      [,3]        [,4]
[1,] 19.3736963  0.2706394 7.1843146   8.3596208
[2,]  0.2706394 28.4826002 9.6901327   6.7164627
[3,]  7.1843146  9.6901327 7.0006189   0.5788023
[4,]  8.3596208  6.7164627 0.5788023 133.4552270

결과 해석도 ML의 경우와 동일하다.

14.5.2.4 -2 log likelihood (SAS output)

m2LL = 2*r2$value ; m2LL
[1] 442.6367

14.5.2.5 AIC

AIC의 계산법이 nlme는 ML이나 REML이나 동일한 반면, SAS는 다르게 되어 있다. 즉, random effect의 parameter 개수로 보정해주고 있다. 이는 minimization 과정의 인자가 random effect들 뿐이라고 간주하기 때문인 것으로 생각된다.

2*nParA + m2LL # nlme style
[1] 454.6367
2*nParR + m2LL # SAS style for REML
[1] 450.6367

14.5.2.6 AICc (SAS output)

AICc의 계산법도 SAS에서는 random effect의 parameter 개수만 count하는데, 이 역시 minimization이 꼭 필요한 것만 count하기 때문으로 보인다.

2*nParA + m2LL + 2*nParA*(nParA + 1)/(nRec - nParA - 1) # nlme style
[1] 455.4684
2*nParR + m2LL + 2*nParR*(nParR + 1)/(nRec - nParR - 1) # SAS style
[1] 451.025

14.5.2.7 BIC

BIC계산법도 다른데, SAS에서는 대상자 수(subject count)로 하는 반면, nlme에서는 총 관측값 개수, total parameter 개수, random effect parameter 개수로 계산한다.

nParA*log(nRec - nParF) + m2LL # nlme style for REML
[1] 470.6173
nParR*log(nID) + m2LL          # SAS style for REML
[1] 455.82

14.5.2.8 Log likelihood

-r2$value
[1] -221.3183

14.5.2.9 Estimates of random effects

Reparametrization 방법은 ML때와 같다.

SDs = exp(r2$par)[-3] ; SDs
[1] 2.3269685 0.2264039 1.3100362
rho = (exp(r2$par)[3] - 1)/(exp(r2$par)[3] + 1) ; rho
[1] -0.6092885

14.5.2.10 Estimates of fixed effects

nlme에서는 다음과 같이 fixed effect에 대한 검정을 한다. 그런데, 자유도 계산법이 SAS와 다르다.

Var1 = solve(FIM)  # no need to further correction
SE = sqrt(diag(Var1))
Df = nRec - qr(X)$rank - (nID - 1) # Df # 80
tVal = Beta/SE
fixedPE = cbind(Beta, SE, Df, tVal, 2*(1 - pt(abs(tVal), Df)))
colnames(fixedPE) = c("Value", "Std.Error", "DF", "t-value", "p-value")
fixedPE
                 Value  Std.Error DF   t-value      p-value
(Intercept) 16.7611111 0.77420564 80 21.649430 0.000000e+00
age          0.6601852 0.07116124 80  9.277314 2.509104e-14
cov2cor(Var1)[lower.tri(Var1)] # correlation between intercept and age estimates
[1] -0.8478078

14.5.2.11 SAS Output

SAS에서는 nlme와 달리 subject수에서 1을 뺀 것을 자유도로 간주하였다.

Df1 = nID - 1 # Df # 26
fixedSol = cbind(Beta, SE, Df1, tVal, 2*(1 - pt(abs(tVal), Df1)))
colnames(fixedSol) = c("Value", "Std.Error", "DF", "t-value", "p-value")
fixedSol
                 Value  Std.Error DF   t-value      p-value
(Intercept) 16.7611111 0.77420564 26 21.649430 0.000000e+00
age          0.6601852 0.07116124 26  9.277314 9.871619e-10

14.5.2.12 Standardized within-group residuals

Standardized individual residual은 다음과 같이 구할 수 있다.

d1 = Orthodont
d1$PRED = NA
d1$IPRE = NA
for (i in 1:nID) {
  d1[ind[[i]], "PRED"] =  X[ind[[i]],] %*% Beta
  d1[ind[[i]], "IPRE"] =  X[ind[[i]],] %*% Beta + Z[ind[[i]],] %*% bi[i,]
}
d1$RES = d1$distance - d1$PRED
d1$IRES = d1$distance - d1$IPRE
d1$IWRE = d1$IRES/SDs[3]
quantile(d1$IWRE, c(0, 0.25, 0.5, 0.75, 1)) # fivenum(d1$IWRE)
          0%          25%          50%          75%         100% 
-3.220629484 -0.494171775  0.007408944  0.471845267  3.915650761

14.5.2.13 Intervals of fixed effects

\(\beta\)의 신뢰구간은 다음과 같이 구할 수 있다.

alpha = 0.05
Len = qt(1 - alpha/2, Df)*sqrt(diag(Var1))
CIbeta = cbind(Beta - Len, Beta, Beta + Len)
colnames(CIbeta) = c("lower", "est", "upper")
CIbeta
                 lower        est      upper
(Intercept) 15.2203928 16.7611111 18.3018294
age          0.5185698  0.6601852  0.8018006

14.5.2.14 Intervals of random effects

G matrix와 R matrix의 요소들에 대한 신뢰구간은 다음과 같이 구할 수 있다. ML때와 마찬가지의 역변환이 필요하다.

vPE = r2$par
Var = solve(r2$hessian)
vB = qnorm(1 - alpha/2)*sqrt(diag(Var))

vLL = exp(vPE[1:2] - vB[1:2])
vUL = exp(vPE[1:2] + vB[1:2])

rhoLL = (exp(vPE[3] - vB[3]) - 1)/(exp(vPE[3] - vB[3]) + 1) # rho LL
rhoUL = (exp(vPE[3] + vB[3]) - 1)/(exp(vPE[3] + vB[3]) + 1) # rho UL

CIomega = cbind(c(vLL, rhoLL), c(SDs[1:2], rho), c(vUL, rhoUL))
rownames(CIomega) = c("sd((Intercept))", "sd(age)", "cor((Intercept),age)")
colnames(CIomega) = c("lower", "est.", "upper")

CIsigma = SDs[3]*exp(c(-1, 0, 1)*vB[4])
names(CIsigma) = c("lower", "est.", "upper")
  • Approximate confidence interval of random effects: Subject level
CIomega
                          lower       est.     upper
sd((Intercept))       0.9486210  2.3269685 5.7080564
sd(age)               0.1025006  0.2264039 0.5000823
cor((Intercept),age) -0.9382154 -0.6092885 0.2980301
  • Approximate confidence interval of random effects: Within group
CIsigma
   lower     est.    upper 
1.084865 1.310036 1.581944

14.5.2.15 Variance-Covariance Matrix

SAS style의 결과물을 보고 싶으면 다음과 같이 추가로 계산해주어야 한다.

Cov = rho*SDs[1]*SDs[2]
OM = matrix(c(SDs[1]*SDs[1], Cov, Cov, SDs[2]*SDs[2]), nrow=2) ; OM # G matrix
           [,1]        [,2]
[1,]  5.4147822 -0.32099428
[2,] -0.3209943  0.05125871
SG = SDs[3]*SDs[3] ; SG     # R matrix
[1] 1.716195
Z[ind[[1]],] %*% OM %*% t(Z[ind[[1]],]) + diag(nrow(Z[ind[[1]],]))*SG # V matrix
         1        2        3        4
1 5.275626 3.737582 3.915733 4.093884
2 3.737582 5.836963 4.503954 4.887139
3 3.915733 4.503954 6.808369 5.680395
4 4.093884 4.887139 5.680395 8.189845

14.5.2.16 EBE, MAP, realized Eta, post-hoc Eta

EBE (random variable u의 realized value)는 miminization과정에서 이미 매 iteration마다 update하여 저장해두었다.

bi
          [,1]        [,2]
M01  1.0515429  0.21579176
M02 -0.7289931  0.01462285
M03 -0.1745595  0.03591728
M04  3.0066615 -0.06636624
M05 -1.1791000  0.02578884
M06  1.2191425  0.08312293
M07 -0.6094301  0.03502246
M08  0.8627813 -0.09496192
M09 -0.4464858  0.13612569
M10  2.6562259  0.21103221
M11  0.8932539 -0.11908352
M12 -1.0009194  0.11483126
M13 -4.1410867  0.41470442
M14  0.9061626 -0.01425830
M15 -0.5379038  0.20849050
M16 -0.1874682 -0.06890795
F01 -0.4852118 -0.17834998
F02 -1.0138280  0.01000604
F03 -0.7747021  0.05080526
F04  1.0714346 -0.03004110
F05  0.5193286 -0.16822154
F06 -0.6200111 -0.18668879
F07 -0.1874682 -0.06890795
F08  1.2542705 -0.17477072
F09 -0.2894672 -0.21825438
F10 -2.2833118 -0.25057206
F11  1.2191425  0.08312293

14.5.2.17 Coefficients of Individuals

개인별 intercept와 slope는 \(\beta\)와 EBE를 조합(이 model에서는 더하기)으로 구할 수 있다.

iCoef = bi
for (i in 1:nID) iCoef[i, ] = iCoef[i,] + t(Beta)
iCoef
        [,1]      [,2]
M01 17.81265 0.8759769
M02 16.03212 0.6748080
M03 16.58655 0.6961025
M04 19.76777 0.5938189
M05 15.58201 0.6859740
M06 17.98025 0.7433081
M07 16.15168 0.6952076
M08 17.62389 0.5652233
M09 16.31463 0.7963109
M10 19.41734 0.8712174
M11 17.65437 0.5411017
M12 15.76019 0.7750164
M13 12.62002 1.0748896
M14 17.66727 0.6459269
M15 16.22321 0.8686757
M16 16.57364 0.5912772
F01 16.27590 0.4818352
F02 15.74728 0.6701912
F03 15.98641 0.7109904
F04 17.83255 0.6301441
F05 17.28044 0.4919636
F06 16.14110 0.4734964
F07 16.57364 0.5912772
F08 18.01538 0.4854145
F09 16.47164 0.4419308
F10 14.47780 0.4096131
F11 17.98025 0.7433081