13.1 포아송 회귀 (Poisson Regression)
다음과 같은 영국 의사들을 상대로 한 자료가 있다.
참고로 R에서 “tbl_df” class는reshape이나 subsetting시 문제가 발생할 수 있으므로 아래와 같이 class를 “data.frame”으로 해두는 것이 좋다.
class(doctors) = "data.frame" # tbl_df class has problem with reshape or subsetting
class(beetle) = "data.frame" # tbl_df class has problem with reshape or subsetting
class(Cars) = "data.frame" # tbl_df class has problem with reshape or subsettingdoctors age smoking deaths person-years
1 35 to 44 smoker 32 52407
2 45 to 54 smoker 104 43248
3 55 to 64 smoker 206 28612
4 65 to 74 smoker 186 12663
5 75 to 84 smoker 102 5317
6 35 to 44 non-smoker 2 18790
7 45 to 54 non-smoker 12 10673
8 55 to 64 non-smoker 28 5710
9 65 to 74 non-smoker 28 2585
10 75 to 84 non-smoker 31 1462흡연 여부가 사망과 어떤 연관이 있는지 연령과 함께 모형화해 보고자 한다. 그러기 위해 우선 연령군을 1 ~ 5의 숫자로 코딩하고, smoker를 1로, non-smoker를 0으로 코딩한다.
그래프를 그려보면 연령에 따라 직선적으로 증가하는 것이 아니라, 약간 아래로 꺽이는 것을 볼 수 있다. 따라서, age의 제곱항도 넣어서 모형화하면 다음과 같다.
doctors[,"Age"] = as.numeric(as.ordered(doctors$age))
doctors[,"Smoke"] = as.numeric(doctors$smoking == "smoker")
doctors[,"ni"] = doctors$'person-years'
doctors age smoking deaths person-years Age Smoke ni
1 35 to 44 smoker 32 52407 1 1 52407
2 45 to 54 smoker 104 43248 2 1 43248
3 55 to 64 smoker 206 28612 3 1 28612
4 65 to 74 smoker 186 12663 4 1 12663
5 75 to 84 smoker 102 5317 5 1 5317
6 35 to 44 non-smoker 2 18790 1 0 18790
7 45 to 54 non-smoker 12 10673 2 0 10673
8 55 to 64 non-smoker 28 5710 3 0 5710
9 65 to 74 non-smoker 28 2585 4 0 2585
10 75 to 84 non-smoker 31 1462 5 0 1462r1 = glm(deaths ~ Smoke + Age + I(Age*Age) + offset(log(ni)), "poisson", doctors)
summary(r1)
Call:
glm(formula = deaths ~ Smoke + Age + I(Age * Age) + offset(log(ni)),
family = "poisson", data = doctors)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.14815 -0.71275 0.02995 0.33549 1.92253
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.83396 0.29516 -33.317 < 2e-16 ***
Smoke 0.35452 0.10737 3.302 0.000961 ***
Age 2.09471 0.18119 11.561 < 2e-16 ***
I(Age * Age) -0.19438 0.02715 -7.159 8.14e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 935.067 on 9 degrees of freedom
Residual deviance: 12.176 on 6 degrees of freedom
AIC: 75.243
Number of Fisher Scoring iterations: 4일반화 선형 모형의 경우 잔차의 종류가 다양하다.
r1$fitted.values # y hat 1 2 3 4 5 6 7 8 9 10
26.781 100.199 203.739 187.868 111.413 6.736 17.347 28.523 26.904 21.491 residuals(r1, type="deviance") 1 2 3 4 5 6 7 8 9 10
0.9782 0.3773 0.1581 -0.1365 -0.9048 -2.1482 -1.3600 -0.0982 0.2100 1.9225 residuals(r1, type="pearson") 1 2 3 4 5 6 7 8 9 10
1.0085 0.3797 0.1584 -0.1363 -0.8918 -1.8248 -1.2838 -0.0979 0.2114 2.0513 rstandard(r1, type="deviance") 1 2 3 4 5 6 7 8 9 10
1.4361 0.4872 0.2348 -0.1815 -1.8414 -2.3738 -1.5304 -0.1196 0.2499 2.3218 rstandard(r1, type="pearson") # preferred, use this 1 2 3 4 5 6 7 8 9 10
1.4807 0.4902 0.2352 -0.1812 -1.8149 -2.0164 -1.4446 -0.1192 0.2515 2.4773 위의 결과는 다음과 같이 계산된다.
Objective function을 계산하기 편하게 하기 위해, y 값, n 값, sum of log factorial (yi)를 다음과 같이 미리 준비해 둔다.
y = doctors[,"deaths"]
ni = doctors[,"ni"]
sumlfacty = sum(lfactorial(y))Design matrix는 다음과 같다.
X = model.matrix(~ Smoke + Age + I(Age*Age), data=doctors) ; X (Intercept) Smoke Age I(Age * Age)
1 1 1 1 1
2 1 1 2 4
3 1 1 3 9
4 1 1 4 16
5 1 1 5 25
6 1 0 1 1
7 1 0 2 4
8 1 0 3 9
9 1 0 4 16
10 1 0 5 25
attr(,"assign")
[1] 0 1 2 3Objective function은 -log likelihood로 다음과 같다.
LLp = function(Beta) # -log likelihood function for poisson distribution
{
LogLam = log(ni) + X %*% Beta
LogLik = sum(y*LogLam - exp(LogLam)) - sumlfacty
return(-LogLik)
}마지막 항 sumlfacty는 optimize 과정에는 불필요하다. 따라서, 이 함수에서는 없애고 나중에 log likelihood를 계산할 때 다시 넣어 주는 것이 더 효율적이다. 하지만, 목적함수값이 -log likelihood가 아니면 생길 수 있는 혼란을 줄이기 위해 포함시켰다.
optim 함수를 이용하여 적합하면 다음과 같다.
r2 = optim(c(1, 1, 1, 1), LLp, method="L-BFGS-B", hessian=TRUE) ; r2$par
[1] -9.8372498 0.3544376 2.0970299 -0.1947360
$value
[1] 33.62183
$counts
function gradient
84 84
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
$hessian
[,1] [,2] [,3] [,4]
[1,] 731.0001 630.0001 2489.000 9353.001
[2,] 630.0001 630.0001 2146.930 8052.424
[3,] 2489.0003 2146.9295 9352.531 37601.080
[4,] 9353.0012 8052.4235 37601.080 158764.770COV2 = solve(r2$hessian) ; COV2 [,1] [,2] [,3] [,4]
[1,] 0.087170563 -0.008217121 -0.048971494 0.006879627
[2,] -0.008217121 0.011528811 -0.001150242 0.000171765
[3,] -0.048971494 -0.001150242 0.032839515 -0.004834250
[4,] 0.006879627 0.000171765 -0.004834250 0.000737221SE2 = sqrt(diag(COV2)) ; SE2[1] 0.29524661 0.10737230 0.18121676 0.02715181Saturated model의 log likelihood와 그 deviance는 다음과 같다.
LLs = sum(y*log(y) - y - lfactorial(y)) ; LLs # LL of saturated model[1] -27.533972*(LLs + r2$value) # Deviance[1] 12.17572먼저 적합시킨 모형의 deviance를 이용하여 p-value를 구하면 다음과 같다.
1 - pchisq(2*(LLs + r2$value), length(y) - length(r2$par))[1] 0.05816164optim에 의한 적합값과 glm에 의한 적합값이 약간 다름을 알 수 있다. 이는 miminization algorithm에 따라 다른 것이다.
이제 glm에 의한 값을 이용해서 나머지를 계산하면 다음과 같다.
b = coef(r1) # beta hat
pihat = exp(X %*% b) # pi hat
yhat = ni*pihat ; t(yhat) # y hat 1 2 3 4 5 6 7 8 9
[1,] 26.78081 100.1992 203.7388 187.868 111.4132 6.735889 17.34677 28.52302 26.9036
10
[1,] 21.49073Pearson residual은 다음과 같다.
pr = (y - yhat)/sqrt(yhat) ; t(pr) # Pearson residual 1 2 3 4 5 6 7 8
[1,] 1.008536 0.3796996 0.1584157 -0.1362827 -0.8918021 -1.824753 -1.283755 -0.09793023
9 10
[1,] 0.2113798 2.051265모형의 deviance와 deviance residual은 다음과 같이 구할 수 있다.
Di = 2*(y*log(y/yhat) - (y - yhat)) ; sum(Di) # deviance[1] 12.17555dr = sign(y - yhat)*sqrt(Di) ; t(dr) # deviance residual 1 2 3 4 5 6 7 8
[1,] 0.9781824 0.3773363 0.158124 -0.1365095 -0.9048239 -2.148154 -1.359981 -0.09823183
9 10
[1,] 0.2099679 1.922531NULL model (intercept만 있는 모형)의 deviance는 다음과 같이도 구할 수 있다.
X = model.matrix(~ 1, data=doctors)
r3 = optim(1, LLp, method="L-BFGS-B") ; r3$par
[1] -5.514416
$value
[1] 495.0676
$counts
function gradient
17 17
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"2*(LLs + r3$value) # null deviance of r0[1] 935.0673그런데, 자세히 보면 age에 따른 curve의 기울기가 흡연 여부에 따라 다른 것을 볼 수 있다. 따라서, 이에 대한 상호작용항(interaction term)을 추가하면 다음과 같다.
r4 = glm(deaths ~ Smoke + Age + I(Age*Age) + Age:Smoke + offset(log(ni)), "poisson", doctors)
summary(r4)
Call:
glm(formula = deaths ~ Smoke + Age + I(Age * Age) + Age:Smoke +
offset(log(ni)), family = "poisson", data = doctors)
Deviance Residuals:
1 2 3 4 5 6 7 8
0.43820 -0.27329 -0.15265 0.23393 -0.05700 -0.83049 0.13404 0.64107
9 10
-0.41058 -0.01275
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -10.79176 0.45008 -23.978 < 2e-16 ***
Smoke 1.44097 0.37220 3.872 0.000108 ***
Age 2.37648 0.20795 11.428 < 2e-16 ***
I(Age * Age) -0.19768 0.02737 -7.223 5.08e-13 ***
Smoke:Age -0.30755 0.09704 -3.169 0.001528 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 935.0673 on 9 degrees of freedom
Residual deviance: 1.6354 on 5 degrees of freedom
AIC: 66.703
Number of Fisher Scoring iterations: 4추정된 \(\beta\) 계수들로 odds ratio를 구하면 다음과 같다.
exp(r4$coef[-1]) Smoke Age I(Age * Age) Smoke:Age
4.2247998 10.7669184 0.8206353 0.7352475 Interaction term이 있는 모형과 없는 모형은 다음과 같이 비교할 수 있다.
c1 = anova(r1, r4) ; c1Analysis of Deviance Table
Model 1: deaths ~ Smoke + Age + I(Age * Age) + offset(log(ni))
Model 2: deaths ~ Smoke + Age + I(Age * Age) + Age:Smoke + offset(log(ni))
Resid. Df Resid. Dev Df Deviance
1 6 12.1755
2 5 1.6354 1 10.541 - pchisq(c1[2, "Deviance"], c1[2, "Df"])[1] 0.001168073통계적으로 유의한 차이가 있으므로 더 복잡한 모형인 상호작용이 있는 모형을 선택한다.