13.2 이항 로지스틱 회귀 (Binomial Logistic Regression)

이분형 로지스틱 회귀를 위한 자료는 다음과 같은 beetle 사망 자료가 있다. 이황화탄소의 로그 농도(x)에 따라 beetle의 사망 개체수(y)를 나타낸 것이다.

beetle

       x  n  y
1 1.6907 59  6
2 1.7242 60 13
3 1.7552 62 18
4 1.7842 56 28
5 1.8113 63 52
6 1.8369 59 53
7 1.8610 62 61
8 1.8839 60 60

R에서 기본자료 분석 script는 다음과 같다.

r5 = glm(y/n ~ x, weight=n, "binomial", beetle)
summary(r5)

Call:
glm(formula = y/n ~ x, family = "binomial", data = beetle, weights = n)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5941  -0.3944   0.8329   1.2592   1.5940  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -60.717      5.181  -11.72   <2e-16 ***
x             34.270      2.912   11.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 284.202  on 7  degrees of freedom
Residual deviance:  11.232  on 6  degrees of freedom
AIC: 41.43

Number of Fisher Scoring iterations: 4

logLik(r5)

'log Lik.' -18.71513 (df=2)

vcov(r5)

            (Intercept)          x
(Intercept)    26.83966 -15.082090
x             -15.08209   8.480525

위의 계산 과정은 다음과 같다.

먼저 fitting을 위한 design matrix와 y 값들의 vector를 준비한다.

X = model.matrix(~x, data=beetle) ; X

  (Intercept)      x
1           1 1.6907
2           1 1.7242
3           1 1.7552
4           1 1.7842
5           1 1.8113
6           1 1.8369
7           1 1.8610
8           1 1.8839
attr(,"assign")
[1] 0 1

n = beetle$n
y = beetle$y
sumlogny = sum(lchoose(beetle$n, beetle$y)) # to reduce calculation amount

Binomial distribution을 사용하는 경우 목적함수는 다음과 같이 할 수 있다.

LLb = function(b) # -log likelihood function for binomial distribution
{
  Xb = X %*% b
  LogLik = sum(y*Xb  - n*log(1 + exp(Xb))) + sumlogny
  return(-LogLik)
}

마지막 항 sumlogny는 optimize 과정에는 불필요하다. 따라서, 이 함수에서는 없애고 나중에 log likelihood를 계산할 때 다시 넣어 주는 것이 더 효율적이다. 하지만, 목적함수값이 -log likelihood가 아니면 생길 수 있는 혼란을 줄이기 위해 포함시켰다.

Fitting 결과는 다음과 같다.

r6 = optim(c(-1, 1), LLb, method="L-BFGS-B", hessian=T)
r6

$par
[1] -60.71753  34.27037

$value
[1] 18.71513

$counts
function gradient 
      20       20 

$convergence
[1] 0

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

$hessian
          [,1]     [,2]
[1,]  58.48413 104.0104
[2,] 104.01040 185.0940

-r6$value # log likelihood

[1] -18.71513

Vcov0 = solve(r6$hessian) ; Vcov0

          [,1]       [,2]
[1,]  26.83972 -15.082124
[2,] -15.08212   8.480544

SE0 = sqrt(diag(Vcov0)) ; SE0

[1] 5.180707 2.912137

r6$par / SE0 # Wald statistics

[1] -11.71993  11.76811

\(\hat{y}\)는 다음과 같다.

expXbhat = exp(X %*% r6$par)
pihat = expXbhat/(1 + expXbhat)
yhat = n*pihat ; t(yhat)

            1        2        3        4        5        6        7        8
[1,] 3.457448 9.841653 22.45136 33.89764 50.09584 53.29093 59.22217 58.74297

Pearson residual을 구해보고 R 내장 함수의 결과와 비교해 보면 다음과 같다.

ri = (y - yhat)/sqrt(n*pihat*(1 - pihat)) ; t(ri) # Pearson residual

            1        2         3         4         5          6        7        8
[1,] 1.409305 1.101108 -1.176256 -1.612382 0.5944412 -0.1281148 1.091419 1.133108

resid(r5, "pearson")

         1          2          3          4          5          6          7          8 
 1.4092960  1.1011003 -1.1762596 -1.6123815  0.5944454 -0.1281090  1.0914228  1.1331102 

이 모형의 deviance를 구하고 summary(r6)의 결과와 비교해보자.

Di = 2*(y*log(y/yhat) + ifelse(n - y > 0, (n - y)*log((n - y)/(n - yhat)), 0))
sum(Di) # deviance

[1] 11.23223

Deviance residual을 구하고 R 내장 함수 결과와 비교하면 다음과 같다.

di = sign(y - yhat)*sqrt(Di) ; t(di) # deviance residual

            1        2         3         4         5          6        7        8
[1,] 1.283685 1.059697 -1.196109 -1.594125 0.6061361 -0.1271641 1.251066 1.593981

resid(r5, "deviance")

         1          2          3          4          5          6          7          8 
 1.2836777  1.0596900 -1.1961123 -1.5941244  0.6061405 -0.1271584  1.2510711  1.5939850 

Null model의 log likelihood와 deviance는 다음과 같이 구할 수 있다.

X = matrix(1, nrow=nrow(beetle), ncol=1) # intercept only
r6n = optim(c(1), LLb, method="L-BFGS-B")
r6n

$par
[1] 0.4262992

$value
[1] 155.2002

$counts
function gradient 
       6        6 

$convergence
[1] 0

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

r6n$par                               # intercept = 0.4263

[1] 0.4262992

-r6n$value                            # log likelihood of null model

[1] -155.2002

D1b = 2*(-r6$value + r6n$value) ; D1b # deviance from null model

[1] 272.9702