일반화 선형모형(Generalized Linear Model)
회귀분석이나 분산분석은 종속변수가 정규분포되어 있는 연속형 변수이다. 하지만 많은 경우에 있어서 종속변수가 정규분포되어 있다는 가정을 할 수 없는 경우도 있으며 범주형 변수가 종속변수인 경우도 있다. 다음과 같은 경우에 일반화 선형모형을 사용한다.
- 종속변수가 범주형변수인 경우 : 이항변수( 0 또는 1, 합격/불합격, 사망.생존 등)인 경우도 있으며 다항변수(예를 들어 poor/good/excellent 또는 공화당/민주당/무소속 등)인 경우 정규분포 하지 않는다.
- 종속변수가 count(예를 들면 한 주간 교통사고 발생 건수, 하루에 마시는 물이 몇잔인지 등)인 경우. 이들 값은 매우 제한적이며 음수가 되지 않고 평균과 분산이 밀접하게 관련되어 있고 정규분포하지 않는다.
일반화 선형 모형은 종속변수가 정규분포하지 않는 경우를 포함하는 선형모형의 확장이며 glm()함수를 사용한다. 이 장에서는 대표적으로 로지스틱회귀(Logistic regression)와 포아송회귀(Poisson regression)를 다룬다.
glm() 함수
일반화선형모형은 glm()함수를 사용한다. glm() 함수의 사용방법은 lm()함수와 유사하나 추가로 family라는 인수를 지정해준다. family에 따라 연결된 함수가 달라지는데 사용법은 다음과 같다.
glm(formula, family=family(link=function), data)
family는 종속변수의 분포에 따라 다음과 같은 것들을 사용할 수 있다. 종속변수의 분포가 정규분포인 경우 gaussian, 이항분포인 경우 binomial, 포아송분포인 경우 poisson, 역정규분포인 경우 inverse.gaussian, 감마분포인 경우 gamma, 그리고 응답분포가 확실하지 않은 때를 위한 유사가능도 모형인 경우 quasi를 사용할 수 있다. glm()함수의 결과를 anova()와 조합하면 분산분석표를 생성할 수 있고 summary()에 넣어서 잔차와 추정값 등을 얻을 수 있다. coef() 함수를 사용하여 모형 인수들의 절편과 기울기 등을 얻을수 있으며 residual()함수는 잔차를 얻을 수 있다. plot()함수를 사용하여 회귀진단 plot을 얻을 수 있고 회귀모형을 사용하여 predict() 함수로 새로운 데이타에 대한 예측치를 추정할 수 있다.
로지스틱 회귀(Logistic regression)
로지스틱 회귀는 종속변수가 이항변수( 0 또는 1, 합격/불합격, 사망.생존 등)인 경우로 의학연구에 많이 사용된다. 로지스틱 회귀분석을 할때에는 family=binomial로 지정해주어야 한다. 다음의 예를 보자.
대장암 데이타
survival 패키지에 있는 colon데이타는 stage B/C의 colon cancer 1858명의 데이타이다. 자세한 내용은 ?survival::colon을 입력해 보면 알수 있다. survival 패키지를 설치하지 않은 경우는 install.packages(“survival”) 을 이용하여 설치한다.
require(survival)Loading required package: survival
Loading required package: splinesstr(colon)'data.frame': 1858 obs. of 16 variables:
$ id : num 1 1 2 2 3 3 4 4 5 5 ...
$ study : num 1 1 1 1 1 1 1 1 1 1 ...
$ rx : Factor w/ 3 levels "Obs","Lev","Lev+5FU": 3 3 3 3 1 1 3 3 1 1 ...
$ sex : num 1 1 1 1 0 0 0 0 1 1 ...
$ age : num 43 43 63 63 71 71 66 66 69 69 ...
$ obstruct: num 0 0 0 0 0 0 1 1 0 0 ...
$ perfor : num 0 0 0 0 0 0 0 0 0 0 ...
$ adhere : num 0 0 0 0 1 1 0 0 0 0 ...
$ nodes : num 5 5 1 1 7 7 6 6 22 22 ...
$ status : num 1 1 0 0 1 1 1 1 1 1 ...
$ differ : num 2 2 2 2 2 2 2 2 2 2 ...
$ extent : num 3 3 3 3 2 2 3 3 3 3 ...
$ surg : num 0 0 0 0 0 0 1 1 1 1 ...
$ node4 : num 1 1 0 0 1 1 1 1 1 1 ...
$ time : num 1521 968 3087 3087 963 ...
$ etype : num 2 1 2 1 2 1 2 1 2 1 ...환자의 치료 방법에 따라 rx는 Obs(ervation), Lev(amisole), Lev(amisole)+5-FU의 세군으로 나뉘고 obstruct는 종양에 의한 장의 폐쇄(obstruction), perfor은 장의 천공(perforation), adhere는 인접장기와의 유착(adherence), nodes는 암세포가 확인된 림프절의 수, differ는 암세포의 조직학적인 분화의 정도(1=well, 2=moderate, 3=poor) 를 나타내며 extent는 암세포가 침습한 깊이(1=submucosa, 2=muscle, 3=serosa, 4=인접장기)를 나타내며 surg는 수술 후 등록 까지의 시간으로 0=short, 1=long을 나타낸다. status는 재발 또는 사망인 경우 1로 되어있다.이 중 조직학적인 분화의 정도가 알려져 있지 않은 46례 전이된 림프절의 수가 누락된 36례 를 제외한 데이타로 로지스틱 회귀분석을 해보면 다음과 같다.
colon1<-na.omit(colon)
result<-glm(status~rx+sex+age+obstruct+perfor+adhere+nodes+differ+
extent+surg,family=binomial,data=colon1)
summary(result)
Call:
glm(formula = status ~ rx + sex + age + obstruct + perfor + adhere +
nodes + differ + extent + surg, family = binomial, data = colon1)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.575 -1.046 -0.584 1.119 2.070
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.430926 0.478301 -5.082 3.73e-07 ***
rxLev -0.069553 0.122490 -0.568 0.570156
rxLev+5FU -0.585606 0.124579 -4.701 2.59e-06 ***
sex -0.086161 0.101614 -0.848 0.396481
age 0.001896 0.004322 0.439 0.660933
obstruct 0.219995 0.128234 1.716 0.086240 .
perfor 0.085831 0.298339 0.288 0.773578
adhere 0.373527 0.147164 2.538 0.011144 *
nodes 0.185245 0.018873 9.815 < 2e-16 ***
differ 0.031839 0.100757 0.316 0.752003
extent 0.563617 0.116837 4.824 1.41e-06 ***
surg 0.388068 0.113840 3.409 0.000652 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2461.7 on 1775 degrees of freedom
Residual deviance: 2240.4 on 1764 degrees of freedom
AIC: 2264.4
Number of Fisher Scoring iterations: 4제 13 장 회귀분석에서 배운 backward elimination방법으로 stepwise logistic regression 을 해보면 다음과 같다.
reduced.model=step(result)Start: AIC=2264.43
status ~ rx + sex + age + obstruct + perfor + adhere + nodes +
differ + extent + surg
Df Deviance AIC
- perfor 1 2240.5 2262.5
- differ 1 2240.5 2262.5
- age 1 2240.6 2262.6
- sex 1 2241.2 2263.2
<none> 2240.4 2264.4
- obstruct 1 2243.4 2265.4
- adhere 1 2246.9 2268.9
- surg 1 2252.1 2274.1
- rx 2 2266.7 2286.7
- extent 1 2265.5 2287.5
- nodes 1 2363.5 2385.5
Step: AIC=2262.52
status ~ rx + sex + age + obstruct + adhere + nodes + differ +
extent + surg
Df Deviance AIC
- differ 1 2240.6 2260.6
- age 1 2240.7 2260.7
- sex 1 2241.2 2261.2
<none> 2240.5 2262.5
- obstruct 1 2243.6 2263.6
- adhere 1 2247.4 2267.4
- surg 1 2252.2 2272.2
- rx 2 2266.8 2284.8
- extent 1 2265.8 2285.8
- nodes 1 2363.7 2383.7
Step: AIC=2260.61
status ~ rx + sex + age + obstruct + adhere + nodes + extent +
surg
Df Deviance AIC
- age 1 2240.8 2258.8
- sex 1 2241.3 2259.3
<none> 2240.6 2260.6
- obstruct 1 2243.7 2261.7
- adhere 1 2247.6 2265.6
- surg 1 2252.4 2270.4
- rx 2 2266.8 2282.8
- extent 1 2266.2 2284.2
- nodes 1 2367.7 2385.7
Step: AIC=2258.8
status ~ rx + sex + obstruct + adhere + nodes + extent + surg
Df Deviance AIC
- sex 1 2241.5 2257.5
<none> 2240.8 2258.8
- obstruct 1 2243.7 2259.7
- adhere 1 2247.9 2263.9
- surg 1 2252.7 2268.7
- rx 2 2266.9 2280.9
- extent 1 2266.4 2282.4
- nodes 1 2368.5 2384.5
Step: AIC=2257.49
status ~ rx + obstruct + adhere + nodes + extent + surg
Df Deviance AIC
<none> 2241.5 2257.5
- obstruct 1 2244.5 2258.5
- adhere 1 2248.8 2262.8
- surg 1 2253.3 2267.3
- rx 2 2267.1 2279.1
- extent 1 2266.9 2280.9
- nodes 1 2369.7 2383.7최종적으로 치료(rx), 장폐쇄(obstruct), 유착(adhere), 전이된 림프절수(nodes), 침습깊이(extent), 수술후 항암치료까지 기간(surg)을 포함하는 모형이 선택되었다. 최종적으로 선택된 reduced.model의 예측 인자들의 회귀계수를 보려면 다음과 같이 한다.
summary(reduced.model)
Call:
glm(formula = status ~ rx + obstruct + adhere + nodes + extent +
surg, family = binomial, data = colon1)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5583 -1.0490 -0.5884 1.1213 2.0393
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.30406 0.35138 -6.557 5.49e-11 ***
rxLev -0.07214 0.12221 -0.590 0.554978
rxLev+5FU -0.57807 0.12428 -4.651 3.30e-06 ***
obstruct 0.22148 0.12700 1.744 0.081179 .
adhere 0.38929 0.14498 2.685 0.007251 **
nodes 0.18556 0.01850 10.030 < 2e-16 ***
extent 0.56510 0.11643 4.854 1.21e-06 ***
surg 0.38989 0.11371 3.429 0.000606 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2461.7 on 1775 degrees of freedom
Residual deviance: 2241.5 on 1768 degrees of freedom
AIC: 2257.5
Number of Fisher Scoring iterations: 4이들 각각의 예측 인자들의 odds ratios를 구하려면 exp(coef(result))와 같이 구해야 하고 95%신뢰구간을 구하려면 exp(confint(result))와 같이 구해야 하며 p value 는 summary(result)$coefficient[,4]와 같이 구해야 한다. 번거로움을 피하기 위해 다음과 같은 함수를 만들어 사용하면 좋을 것이다.
ORtable=function(x,digits=2){
suppressMessages(a<-confint(x))
result=data.frame(exp(coef(x)),exp(a))
result=round(result,digits)
result=cbind(result,round(summary(x)$coefficient[,4],3))
colnames(result)=c("OR","2.5%","97.5%","p")
result
}실제 사용 예는 다음과 같다.
ORtable(reduced.model) OR 2.5% 97.5% p
(Intercept) 0.10 0.05 0.20 0.000
rxLev 0.93 0.73 1.18 0.555
rxLev+5FU 0.56 0.44 0.72 0.000
obstruct 1.25 0.97 1.60 0.081
adhere 1.48 1.11 1.96 0.007
nodes 1.20 1.16 1.25 0.000
extent 1.76 1.41 2.22 0.000
surg 1.48 1.18 1.85 0.001과산포(Overdispersion)
이항분포를 따르는 데이타의 기대 분산은 \(\sigma^2 = n \pi (1- \pi)\) 인데 이때 \(n\)은 관찰 수이고 \(\pi\) 는 Y=1 군에 속할 가능성이다. 과산포는 종속변수의 실제 분산이 이항분포에서 기대되는 예측분산보다 클 때 생길수 있으며 과산ㅊ포는 검정의 표준오차를 뒤틀리게 만들어 검정을 부정확하게 만들수 있다. 과산포가 있는 경우에도 glm()함수로 로지스틱 회귀분석을 할 수 있는데 이 때에는 family=binomial 대신 family=quasibinomial을 사용하여야 한다. 과산포를 발견할 수 있는 한가지 방법은 이항모형에서 잔차의 이탈도(residual deviance)와 잔차의 자유도를 비교해 보는 것이다. 이 비가 1이 훨씬 넘게 되면 과산포가 있다는 것을 시사해준다. 이 모형에서 계산해보면 다음과 같다.
이 비율이 1에 가까우므로 과산포는 없다는 것을 시사해준다. 과산포가 있는지 검정하는 방법은 먼저 glm()을 통한 모형적합을 두번 하는데 먼저 family=binomial로 한번 적합시키고(예를 들어 결과를 fit로 저장) family=quasibinomial로 한번 더 적합시켜(fit.od로 저장) 다음 검정을 시행하면 된다.
pchisq(summary(fit.od)$dispersion * fit$df.residual, fit$df.residual,lower=F)
이 검정 결과 p값이 0.05이하인 경우 과산포가 있다고 할수 있다. 대장암 데이타에 이검정을 적용해보면 다음과 같다.
fit=glm(formula = status ~ rx + obstruct + adhere + nodes + extent +
surg, family = binomial, data = colon1)
fit.od=glm(formula = status ~ rx + obstruct + adhere + nodes + extent +
surg, family = quasibinomial, data = colon1)
pchisq(summary(fit.od)$dispersion*fit$df.residual,
fit$df.residual,lower=F)[1] 0.2803691검정 결과 p 값이 0.28 로 명백히 유의하지 않으므로(p > 0.05) 과산포는 없다고 확신할 수 있다.
Plot for Odds Ratios
Odds Ratios를 그림으로 직관적으로 전달하기 위해 plot.OR이라는 함수를 만들어 보았다. 이 함수는 base에있는 그래프 명령어만을 사용하여 ggplot2패키지를 사용한 것과 비슷안 플롯을 만들어준다. 사용하는 방법은 plot.OR() 함수에 glm() 함수의 결과를 전달해 주기만 하면 된다. 그림은 모두 세가지로 그릴 수 있으며(type=1,2,3, default는 1) show.OR(dafault값은 TRUE) 및 show.CI인수(dafault값은 FALSE)를 조절할수 있으므로 모두 12가지 조합으로 그릴 수 있으나 대표적으로 네가지만 소개하겠다. plot.OR함수는 내부적으로 plot() 함수를 사용하므로 plot()함수에서 사용하는 인수들을 사용할 수 있다.첫번째 예에서는 main인수를 써서 plot의 제목을 인쇄한다.
plot.OR(reduced.model, main="Plot for Odds Ratios of Reduced Model") 
plot.OR 함수는 기본적으로 odds.ratios 그래프 위에 odds.ratio와 함께 p value의 significancy 를 표시해준다. 만일 confidence interval을 같이 표시해주고 싶으면 show.CI인수를 TRUE로 지정해준다. 또한 이때에는 중복을 피하기 위해 show.OR은 FALSE로 표시해주는 것이 좋다.
plot.OR(reduced.model,type=2,show.OR=FALSE,show.CI=TRUE,
main="Plot for Odds Ratios; type=2, show.CI=TRUE") 
Odds ratios를 bar graph형식으로 표시하고 싶을 때에는 type=3을 지정해주면 된다.
plot.OR(reduced.model,type=3,main="Bar Plot for ORs with type=3")
Odds ratios를 type=3 으로 표시할때에도 show.CI인수를 TRUE로 지정하여 95%신뢰구간을 같이 표시할 수 있다.
plot.OR(reduced.model,type=3,show.CI=TRUE,main="Bar plot for ORs with 95% CI")
다음 장에서는 이 plot을 그린 plot.OR 함수에 대해 자세히 설명할 예정이다. 자신만의 plot을 만드는데 관심이 있는 독자분에게는 큰 도움이 될 것이다.