GLMM: Logistic Regression with a random effect

[안재형]

모델: g[E(Y|bi)] = mu + bi , g()=logit-link, bi~N(0, sigma_b^2)

간단한 데이터로 해보죠.

4명이고, 3번 반복측정한 (0, 1)입니다. 이항분포죠.

> cbind(id,y)
id y
[1,] 1 1
[2,] 1 0
[3,] 1 0
[4,] 2 1
[5,] 2 1
[6,] 2 1
[7,] 3 0
[8,] 3 0
[9,] 3 0
[10,] 4 0
[11,] 4 1
[12,] 4 0

# R로 fit해보면
> library(lme4)
> id = rep(1:4, rep(3,4))
> y = c(1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0)
>
> glmer(y ~ 1 + (1|id), family=binomial, nAGQ=10)
Random effects:
Groups Name Std.Dev.
id (Intercept) 1.572
Number of obs: 12, groups: id, 4
Fixed Effects:
(Intercept)
-0.4549

# 이게 결과입니다.

-------------------
mu = -0.4549
sd(bi) = 1.572

### 이제 데이터를 가지고 Likelihood를 만들고, mu와 sd(bi)의 MLE를 구해보죠.

각 4명이 3번의 반복측정에서 성공횟수는 다음과 같습니다.
id=1, 1 success 
id=2, 3 successes
id=3, 0 success
id=4, 1 success

### Likelihood

각각 4 사람의 이항분포에서 bi의 분포인 정규분포를 곱해서 적분한 겁니다. 맨 위의 수식처럼요.

parm[1]은 mu, parm[2]는 sd(bi) 입니다.

integrate()이 quadrature라는 방법으로 적분합니다.

log-likelihood는 다음과 같이 4명의 log-likelihood의 합입니다.

---------------------------------------------------
lik=function(parm){
loglik=
log(integrate(function(x){dbinom(1,size=3,prob=exp(parm[1]+x)/(1+exp(parm[1]+x)))*dnorm(x,mean=0,sd=parm[2])}, -20, 20)$value)+
log(integrate(function(x){dbinom(3,size=3,prob=exp(parm[1]+x)/(1+exp(parm[1]+x)))*dnorm(x,mean=0,sd=parm[2])}, -20, 20)$value)+
log(integrate(function(x){dbinom(0,size=3,prob=exp(parm[1]+x)/(1+exp(parm[1]+x)))*dnorm(x,mean=0,sd=parm[2])}, -20, 20)$value)+
log(integrate(function(x){dbinom(1,size=3,prob=exp(parm[1]+x)/(1+exp(parm[1]+x)))*dnorm(x,mean=0,sd=parm[2])}, -20, 20)$value)
return(-loglik)
}

# optim()으로 MLE를 구해보면 동일합니다.
---------------------------------------------------
> optim(c(-0.1,3), lik, control=list(reltol=1e-20))$par
[1] -0.4549106 1.5719257
---------------------------------------------------

# SAS PROC GLIMMIX로 해도 거의 같은 결과가 나옵니다.
data glmm;
input id y @@;
datalines;
1 1 1 0 1 0
2 1 2 1 2 1
3 0 3 0 3 0
4 0 4 1 4 0
;

proc glimmix data=glmm method=quad;
class id;
model y = / s link=logit dist=binomial;
random intercept / subject=id;
run;

[강성찬]

감사합니다. 어려운 내용이네요.

[안재형]

자세히 설명 못했는데, 한번 이해하면 그리 어려울 것도 없습니다. MLE를 구하려면 likelihood부터 만들어야 하는데, bi가 random이어서 적분해버려야 합니다. 저 수식이 그 적분입니다. 적분후 그냥 다 곱한게 likelihood이고 (log-likelihood여서 더했지만) 그걸 극대화 시키는 mu와 sd를 구한 겁니다.

[강성찬]

그렇군요, 잘 보겠습니다.