일반화 선형 모델

---

원본: PDF

---

이전 절 회귀에서, 우리는 응답 y가 연속적이고 그리고 정규 분포(normal distribution)를 가진다고, 적어도 근접하게, 일반적으로 추정하는 경우인 전통적인 설정에 주로 집중했다. 어쨌든, 많은 데이터 분석 환경의 예에서, 모델이 만들어진 데이터는 확실히 비-정규이다. 예를 들면, 우리는 이진 응답을 기록했을 것이다(e.g. 생존 또는 죽음, 또는 감염 또는 비감염, 등.) 0-1 변수는 명백히 비-정규이다. 많은 다른 예에서, 관심을 가지는 응답은 횟수(count)이다, e.g. 둥지에 알의 수 또는 나무에서 포도 덩굴 수. 포아송 분포(Poisson regression)는 데이터를 셈하는 모델에 종종 사용되어 진다 그리고 포아송 회귀는 예측에서 수량의 응답에 연관되어 사용될 수 있다. 이 장에서는 비 정규 분포를 따르는 응답에 적용하기 위해 일반화 선형 모델(generalized linear model)을 소개한다. 단어 "linear"는 여기에 나타난다 왜냐하면 응듭은 예측자의 선형 조합으로 여진히 모델이 만들어 진다 그러나 연관은 이전 장에서 처럼 직접 연관이 필요하지 않다. 대신에, 우리는 선형 예측자에 응답을 연결하는 link 함수를 소개할 필요가 있다.

일반화 선형 모델을 정의하기 전에, 우리는 특별한 경우를 소개하는 것으로 시작한다 - 로지스틱 회귀(logistic regression).

---

1 로지스틱 회귀 (Logistic Regression)

전통적인 회귀 프레임워크에서, 우리는 하나 또는 더 많은 예측 변수가 있는 함수에서 연속적인 응답 변수 y의 모델 구축에 관심을 가졌다. 대부분의 회귀 문제는 이런 형식이다. 어쨌든, 이진 응답을 제외하고 관심있는 응답이 연속적이지 않은 경우는 다양한 예가 있다. 관심이 있는 측정된 출력이 성공 또는 실패인 실험을 고려해라, 우리는 1 또는 0으로 코드화할 수 있다. 성공 또는 실패의 확률은 예측자 변수의 집합에 의존할 것이다. 이런 데이터로 모델을 생성하는 법법으로 하나의 생각은 예측자의 어떤 값이 주어진 성공의 확률을 추정할 목적을 가지는 회귀를 단순히 적합해야 한다. 어쨌든 이 접근은 작동하지 않을 것이다 왜냐하면 확률은 0 과 1 사이의 값으로 제약 된다. 전통적인 회귀에서는 연속적인 응답을 가지고 출발한다, 예측된 변수는 모든 실수에서 범위를 가질 수 있다. 그러므로 다른 모델링 기술이 필요하다.

이전 장 회귀에서, 모델에 대해 생각할 수 있는 하나의 방법은 조건부 기대(conditional expectation)의 용어에 있다: E[y|x]는 주어진 x에서 y의 조건부 기대이다. 단순 선형 회귀에서, 우리는 E[y|x] = β₀ + β₁x 를 가졌다: 이것은 조건부 기대(conditional expectation)는 x의 선형 함수이다. 만약 y가 값 0과 1을 가지는 이진이면, 그러면

E[y|x] = P(y = 1 | x).

이것은, 이진 출력에서, x에서 y의 회귀는 조건부 확률(conditional probability)이다. 만약 우리가 성공을 y = 1로 표기하면, 그러면 목표는 x에서 주어지는 성공의 확률 모델을 생성해야 한다. 여기에서 설명된 이런 문제의 접근은 로지스틱 회귀(logistic regression)로 알려졌다. 다른 접근법 또한 가능한다는 것에 주목해라, 유명한 프라빗 회귀(probit regression).

위험 평가 연구에서, dose-response 곡선은 주어진 약물 투여량 또는 독소 노출 그리고 인간 또는 동물에서 그것의 효과 사이의 관계를 설명하는데 사용된다. 우리는 예가 있는 문제로 설명해야 한다.

Example. 딱정벌레는 5 시간 동안에 다양한 농도(mf/L)에서 가스 이황화 탄소에 노출되었다.그리고 죽은 딱정벌레 수가 기록되었다. 데이터는 아래의 테이블이다:

Dose(량)	# Exposed	# Killed	Proportion
49.1	59	6	0.102
53.0	60	13	0.217
56.9	62	18	0.290
60.8	56	28	0.500
64.8	63	52	0.825
68.7	59	53	0.898
72.6	62	61	0.984
76.5	60	60	1.000

만약 우리가 CS₂의 농도를 나타낸다고 하면, 그러면 데이블에서 죽음의 확률은 x의 수준의 증가로 증가한다. 사망률(mortality rate)은 량의 증가로 증가함을 주목해라. 우리의 목적은 회귀 모델을 이용하여 량 x의 함수로 사망의 확률을 모델로 생성해야 한다.

만약 우리가 CS₂의 가장 낮은 농도 49.1 mf/L 에서 딱정벌레 노출된 순간을 고려한다면, 우리는 전체 n = 59 딱정벌레가 노출되고 그리고 그 중 6 마리가 죽었다. 이항 확률(binomial probability) 모델은 이런 타입의 결과의 모델을 생성하기 위해 사용된다. 이항 실험은 (ii)각 시행은 두개의 가능한 출력 중 하나가 결과인 경우에 (i) n번의 독립 시행으로 구성 된다 (일반적으로 각각 1 또는 0 으로 기록되는 성공 또는 실패라 불림), 그리고 (iii) 성공 확률 p는 각 시행에서 같게 나타난다. 가장 낮은 농도의 딱정벌레 예에서, 우리는 n = 59의 딱정벌레를 가졌다. 각 딱정벌레에서, 결과는 죽음(성공) 이거나 생존(실패)이다. 우리는 p를 개별 딱정벌레가 이 농도에서 5시간 후에 죽을 것인지 나타내도록 놓을 수 있다. 만약 우리가 y를 죽은 딱정벌레의 수를 나타내게 놓는다면, 그러면 y는 이항 분포(binomial distribution)를 가지는 랜덤 변수(random variable)이다. k = 0, 1, ... , 59 인 경우에 y는 값 k를 가정하는 확률을 보여주게 할 수 있는 하나는, 아래의 식으로 주어진다:

현재 설정에서, 상황은 지금 보다 더 복잡해 진다 왜냐하면 각 농도에서, 우리는 성공 확률 p가 x에 의존하는 경우에 상응하는 이항 분포를 가지기 때문이다, x는 CS₂ 농도. 그러므로, 우리의 목적은 그러면 모델 p(x)가 된다, CS₂ 농도의 노출에서 사망 확률은 x와 일치한다. 위에서 다룬 것 처럼, 모델은

p(x) = β₀ + β₁x

p(x)는 0과 1사이의 값을 가져야 하므로 동작하지 않을 것이다. 이런 상황에서 모델링의 표준 방법의 하나는 로지스틱 회귀 함수(logistic regression function)를 사용해야 한다:

p(x) = e^β₀+β₁x
(1 + e^β₀+β₁x) (2)

β₀ 그리고 β₁ 은 관심있는 회귀 변수인 경우 (단순한 선형 회귀의 교차점과 기울기 변수와 비슷). 해석에 주의해라, (2)에서 p(x)의 값은 요구되는 0과 1 사이에 놓이는 제약 사항이다.

여기에 (2)의 로지스틱 회귀 함수에 대한 몇가지 주의 사항:

1. 만약 β₁ = 0 이면 p(x)는 상수이다. 이것은, 성공 확률은 x에 의존하지 않을 것이다.
2. 만약 β₁ > 0 이면 p(x)가 증가 함수이다 그리고 β₁ < 0 이면 감소 함수이다.
3. 확률 p에서, 함수 p/(1 − p)는 교차 비(odds ratio)라 불린다.
4. 하나는 아래로 표시할 수 있다
   log[p(x) / (1 - p(x))] = β₀ + β₁x.
   함수 log[p(x) / (1 - p(x))]는 p(x)의 로짓(logit)으로 불린다:
   logit(p(x)) = log[p(x) / (1 - p(x))] = β₀ + β₁x,
   그리고 로지스틱 회귀에서 링크(link) 함수로 알려졌다. p(x)의 로짓은 일반적의 선형 회귀 표현을 만든다는 것을 주목하라.

자연스럽게 문제가 지금 도출된다 - 우리는 로지스틱 회귀 모델에서 변수 β₀와 β₁ 추정하기 위해 어떻게 사용할 것인가? 대답은 최대 우도 추정(maximum likelihood estimation)의 방법을 이용해야 한다. 최대 우도 추정 뒤에 논리값은 무엇이 관측된 데이터가 가장 잘 발생하는지 β₀와 β₁의 값으로 결정해야 한다. 최대 우도 추정 방법은 로지스틱 회귀를 제외한 많은 통계적 응용에 가장 넓게 사용된다. 최대 우도 추정은 종종 데이터의 가장 효율적인 사용의 측면에서 다른 타입의 추정 프로시저 보다 더 좋은 성능을 보인다. 그러므로, 최대 우도 추정은 통계 연습에서 가정 유명한 추정 방법이다. 우리는 최대 우도 추정에 대해 간략한 소개를 지금 제공했다.

1.1 최대 우도 추정에 대한 간략한 소개
(A Brief Introduction to Maximum Likelihood Estimation)

이론을 설명하기 위해, 예제를 제공한다. 새의 어떤 종을 2개의 다른 지역에서 찾았다고 가정하라, A 그리고 B로 말함. 새의 30%는 지역 A에서 웨스트 나일 바이러스(West Nile virus)에 감염되었다 그리고 새의 50%는 지역 B에서 바이러스에 감염되었다. 지금 지역 중 하나에서 n = 100 마리의 새의 표본을 가졌다고 가정해라, 그러나 어느 지역인지 알지 못한다. n = 100 마리의 새를 검사하고 그리고 그들 중 40%가 웨스트 나일 바이러스를 지녔다고 알아 냈다. p를 새의 표본을 얻은 모집단에서 새의 비율을 나타내라. 그러면 p = 0.3 또는 p = 0.5 중 하나는 새가 각각 지역 A 또는 B의 여부에 의존. 주어진 관측 데이터 (i.e. 100 마리의 감염된 새 중 40), p의 값은 어느쪽에 더 가까운가? 최대 우도의 사상으로 관측된 데이터가 더 잘 발생하는 변수의 값을 찾아야 한다(이 경우 p)

만약 p 무작위로 선택한 새가 감염될 확률이면, 그러면 1 - p는 새가 감염되지 않을 확률이다. 이항 밀도를 사용하여 (1), 그러면 n 마리의 새 중 k 마리가 감염으로 관측되는 확률은

L(p) = ( n
k ) p^k(1 - p)^n-k

이다.

우리는 지금 p의 함수로 L을 생각한다, 성공 확률 변수 그리고 그것을 우도(likelihood)라 부른다. 최대 우도 추정의 목적은 L을 최대로 하는 p의 값을 찾는데 있다, 우도(likelihood). 그러므로 p에 단지 2개의 가능한 값이 있다, 우리는 p = 0.3 그리고 p = 0.5에서 평가한다:

L(0.3) = ( 100
40 ) 0.3⁴⁰(1 - 0.3)⁶⁰ = 0.00849.

그리고

L(0.5) = ( 100
40 ) 0.5⁴⁰(1 - 0.5)⁶⁰ = 0.01084.

그러므로, p = 0.5는 감염이 있는 100 마리 중 40마리의 새를 관측으로 주어진 더 선호하는 값이다. 그것은, p의 최대 우도 추정은 p^ = 0.5 이다.

계산 방법은 이항 분포 확률 계산하는 위키 참조: URL

이항 확률이 적용된 대부분의 최대 우도 추정의 응용(applications)에서, 성공 확률은 완전하게 알려지지 않는다. 예로, 이전 예에서 관심있는 모집단에서 n = 100 마리의 새의 표본을 단순히 가졌고 그리고 그들 중 k = 40 마리가 웨스트 나일 바이러스에 감염되었다고 가정해라. p를 감염이 있는 모집단에서 새의 비율로 나타내라. 우리는 이 수정된 표본에서 p를 추정하기 위해 최대 우도를 사용할 수 있다. 여기서, 우리 모두는 0 < p ≤ 1 임을 알고 있다. 다시 한번, 우도(likelihood)는

L(p) = ( n
k ) p^k(1 - p)^n-k = ( 100
40 ) p⁴⁰(1 - p)⁶⁰

이다.

목적은 이 표현에서 L(p)를 초대로하는 p의 값을 찾아야 한다. 동일하게, 우리는 이 표현에서 아래의 식으로 작업을 쉽게 하는 자연 로그(natural logarithm)로 최대화 할 수 있다:

l(p) = log[L(p)] = log[ ( 100
40 ) ] + 40 log[p] + 60 log[1 - p]

p에 기대되는 이 식의 유도를 가지고 그리고 그것을 0 으로 설정함으로써 우리가 p의 최대 우도 추정자(maximum likelihood estimator: mle)를 찾게 해 준다. 계산을 수행한 후에, mle는

p^ = k ⁄ n

임을 찾았다, p의 자연(로그) 추정자가 단순히 표본 비율인 경우. 이 예에서, mle는 p^ = 40/100 = 0.4 임을 알아냈다.

이 예는 최대 우도 추정에 상당히 기초적인 설명이다. 최대 우도 이론은 많은 변수를 가지는 매우 상당히 복잡한 모델에서 잘 확장되어 사용된다. 최대 우도 이론에서, 그것은 가장 표준인 상황을 따른다, 최대 우도 추정자는 상대적으로 큰 표본이 제공되는 근사한 정규 분포를 가진다. 최대 우도 추정자는 또한 최대 우도 추정을 매우 유명한 통계적 추정 과정으로 만들어 지는 다른 편향되지 않은(unbiased) 추정자 보다 더 작은 분산을 가지는 경향이 있다.

1.2 단순한 로지스틱 회귀 (Simple Logistic Regression)

지금 우리는 최대 우도 추정이 로지스틱 회귀 모델에서 어떻게 사용될 수 있는지에 대한 생각을 얻을 수 있다. 단순한 로지스틱 회귀 모델에서, 성공 확률 p는 회귀 변수 x에 의존한다, 따라서 우리는 (2)에서 정의되는 p(x)로 나타낸다. 우리의 데이터는 x₁, x₂, ... , m 의 값에서 수집되었다. 딱정벌레 사망 예에서, 가장 낮은 가스의 농도는 x₁ = 49.1 이다, y₁ = 6 의 죽음이 있는 이 농도에 노출된 n₁ = 59 마리를 가졌다. 딱정벌레가 생존해 있는지 아닌지는 다른 하나에 의존하는 가정에서, 우리는 우도 함수(likelihood function)를 아래와 같이 쓸 수 있다

L = ( n₁
y₁ ) p(x₁)^y₁(1 - p(x₁))^n₁−y₁ × ... × ( n_m
y_m ) p(x_m)^y_m(1 - p(x_m))^{n_m − y_m}

이 우도 함수 L은 실제로 β₀와 β₁의 함수임을 주목해라, 로지스틱 회귀 변수

p(x_i) = e^β₀+β₁x
(1 + e^β₀+β₁x) , i = 1, 2, ... , m 에서.

이기 때문이다.

불운하게도, 단순한 선형 회귀의 경우에서 처럼 근접한 형식으로 로지스틱 회귀의 추정치 β₀와 β₁가 주어지는 식이 존재하지 않는다. 대신에, 반복 알고리즘은 β₀와 β₁의 최대 우도 추정을 결정하기 위해 필요하다. 많은 소프트웨어 패키지는 로지스틱 회귀 모델을 적합하는 능력을 가진다. 최대 우도 추정을 찾기 위한 가장 유명한 알고리즘 2개는 Newton-Ralphson algorithm 그리고 the iterated re-weighted least squares algorithm 이다. 우리는 generalized linear model 표준인 R에서 "glm" 함수를 사용하는 로지스틱 회귀 모델의 적합을 설명할 것이다. (SAS에서 PROC LOGISTIC을 사용하는 로지스틱 회귀 모델을 적합할 수 있다.)

Example. 우리는 위에서 딱정벌레 예제로 돌아간다 지금. 로지스틱 회귀 모델을 적합하기 위한 R 코드는 (데이터는 forward 패키지에 있음):

bliss=read.table("bliss.dat", header=T)
alive=bliss[,2]-bliss[,3]
deaths=bliss[,3]
concentration=bliss[,1]
phat=deaths/(deaths+alive)
fit <- glm(cbind(deaths,alive) ~ concentration, family=binomial)
summary(fit)

이다

glm 함수에서, 명령어 "family = binomial"는 로지스틱 회귀 모델을 적합하라고 R 에게 말한다; glm은 다른 모델을 잘 적합할 수 있다, 우리가 보는 바와 같이. 위의 명령어는 아래의 출력을 만든다:

Call:
glm(formula = cbind(deaths, alive) ~ concentration, family = binomial)

Deviance Residuals:
    Min       1Q  Median      3Q    Max
-1.2746  -0.4668  0.7688  0.9544  1.2990

Coefficients:
               Estimate  Std. Error  z value     Pr(>|z|)
(Intercept)   -14.82300     1.28959   -11.49  <2e-16 ***
concentration   0.24942     0.02139    11.66  <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.2024 on 7 degrees of freedom
Residual deviance:   7.3849 on 6 degrees of freedom
AIC: 37.583

Number of Fisher Scoring iterations: 4

출력에서 우리는 교차점과 기울기에서 최대 우도 추정치가 β₀^ = −14.8230 그리고 β₁^ = 0.2494 를 본다, 아래의 추정된 로지스틱 회귀 모델의 결과:

p(x_i) = e^{−14.8 + 0.25x}
1 + e^{−14.8 + 0.25x}

또한, 살충제 농도에 대한 기울기 추정은 상당한 의미가 있다, 이것은 놀랍지 않다. x의 값에 연결하여, CS₂ 농도, 5시간 노출 후에 그 량에서 죽는 딱정벌레의 추정된 확률을 산출한다. 단순한 비율과 추정된 로지스틱 회귀 곡선의 그림은 Fingure 1에 그려졌다.

Figure 1: 추정된 로지스틱 회귀 곡선을 따라 딱정벌레 사망 비율 대 이황화 탄소 농도

기울기에 대한 추론 (Inference for the Slope). 단순한 로지스틱 회귀에서 관심의 주요 질문은(i.e. 오로지 하나의 예측 변수를 가지는 로지스틱 회귀) 만약 기울기 β₁이 0 이 아닌 경우이다. 만약 β₁ = 0 인 경우를 상기하라, 그러면 예측자는 성공 또는 실패 확률을 가지고 있지 않다. 딱정벌레 예에서, 그것은 CS₂이 증가함에 따라 그것은 상당히 분명하다, 게다가 사망 확률은 증가한다. β₁의 중요성을 검증하는 단순한 방법은 Wald 검증을 사용해야 한다. 큰 표본 크기에서, 로지스틱 회귀의 β₁ 같은 최대 우도 추정자는 전형적으로 근사적인 정규 분포를 따른다. 만약 우리가 가설 검증을 원한다면

H₀: β₁ = 0 대 H_α: β₁ ≠ 0,

그러면 통계량은

z = β₁^
se^(β₁^)

H₀이 참이고 그리고 표본 크기가 충분히 클 때 근사적으로 표준 정규 분포를 따를 것이다. 만약 H₀이 거짓이면, 그러면 |z|는 크게 되는 경향일 것이다 그리고 우리는 표준 정규 분포의 임계 값으로 z를 비교할 수 있다. 만약 당신이 표준 정규 랜덤 변수를 제곱한다면, 1의 자유도에서 카이-제곱 랜덤 변수를 얻는다. 그러므로, 대안으로, 우리는 1의 자유도에서 카이-제곱 분포의 임계값으로 z²를 비교 가능하다. 만약 H₀가 거짓이면 (i.e β₁ ≠ 0), z²는 커지게 되는 경향일 것이다.

딱정벌레 예에서, β₁^ = 0.2494 그리고 se^(β₁^) = 0.0214 이다. 이 표준 에러는 최대 우도 이론을 사용하여 찾아진다 (로그-우도 함수의 2번째 부분 유도를 가지도록 포함하는 이론). z 검증 통계량은 그려면

이다, 높은 중요도 (p < 2^-16) 이다.

1.3 로지스틱 회귀 계수의 해석
(Interpreting the Logistic Regression Coeficients)

딱정벌레 사망 표본에서 추정된 기울기는 β₁^ = 0.2494 이다. CS₂ 량과 사망율 간에 우리에게 말하는 것은 무엇인가? 단순 선형 모델(linear regression)을 기억하라, 기울기 β₁은 x에서 하나의 단위 변화에서 평균 응답 y의 변화를 측정이다. 로지스틱 회귀에서 기울기 β₁의 해석은 어쨌든 그렇게 간단하지 않다. 로지스틱 회귀 기울기의 이해와 해석을 돕기 위해, 우리는 첫번째로 예측자 변수 x가 이분법인 경우인 매우 단순한 예제를 볼 것이다.

전립선 암 예제(Prostate Cancer Example). n = 53 인 전립선 암 환자의 데이터 (Collet, 1991)가 수집되었다. 만약 암이 림프절 주의로 퍼졌는지 확인하기 위하여 각 환자에게 개복술(laparectomy) 수술을 했다. 목적은 만약 종양의 크기가 암이 림프절에 퍼졌는지 아닌지 예측할수 있는지 확인하기 위함이다. 지시자(표시) 회귀 변수(indicator regressor variable) x를 다음으로 정의

x = { 0, 작은 종양에서
1, 큰 종양에서

그리고

y = { 0, 림프절로 퍼지지 않음
1, 림프절로 퍼짐

한다.

x에서 y의 로지스틱 회귀 적합에서 최대 우도 추정자:

β₀^ = -1.4351 그리고 β₁^ = 1.6582 .

암이 림프절에 퍼져 있는지 추정 확률:

p^(1) = e^{-1.435+1.658
(1 + e-1.435+1.658)} = .5555 는 큰 종양 환자

p^(0) = e^{-1.435
(1 + e-1.435)} = .1923 는 작은 종양 환자

ODDS: 큰 종양 환자에서(i.e., x = 1), 암이 림프절로 퍼져 있을 (추정된) ODDS :

p^(1)
1 - p^(1) = .5555 ⁄ (1 - 5555 ) = 1.25

odds의 해석: 큰 종양 환자에서 암이 림프절에 펴져 있을 확률의 하나이며 그리고 암이 림프절에 퍼져 있지 않을 확률의 1.25배 이다. 작은 종양 환자에서 (i.e., x = 0), 암이 림프절에 퍼져 있을 (추정된) odds:

p^(0)
1 - p^(0) = .1923 ⁄ (1 - 1923 ) = .238

그러므로, 작은 종양 환자에서, 암이 퍼질 확률은 암이 펴지질 않을 확률의 약 1 ⁄ 4 (≈ .238) 뿐이다. 또는 반대로 생각, 우리는 작은 종양 환자에서 암이 퍼지지 않을 확률은 퍼질 확률 보다 약 4 배 더 높다.

ODDS RATIO는 x = 0의 odds 에서 x = 1 의 odds 에서의 비로 정의된다:

ODDS RATIO = p(1) ⁄ (1 - p(1))
p(0) ⁄ (1 - p(0))  

하나는 보여주기 위해 약간의 대수학으로 할 수 있다:

ODDS RATIO = e^β₁ (3)

로지스틱 회귀 기울기 β₁는 다음의 odds 비로 연관되어 진다:

LOG-ODDS = log(odds ratio) = β₁

전립선 암 예에서 우리는

β₁^ = 1.6582

찾아냈다

그러므로, odds-비 (작은 종양을 가지는 것에서 큰 종양을 가지는 비)는 다음으로 추정된다:

e^β₁ = e^1.658 = 5.25

이 예에서 odds-비의 해석은 암이 퍼지는 odds는 작은 종양을 가지고 비교된 큰 종양의 비가 약 5 배 이다.

주의 (Caution): 이러한 량은 큰 종양을 가지는 것이 작은 종양이 가지는 것과 비교해 암 확산이 5배 더 잘 발생한다는 의미가 아니다(연관 위험(relative risk)을 보라).

odds 비의 이해에 도움을 주기 위해, 작은 종양을 가지는 100명의 환자가 있다고 가정하자. 작은 종양 환장에서 암이 퍼진 odds는 거의 1/4 이다. 그러므로, 우리는 암이 퍼진 모든 20명의 환자를 기대할 수 있다, 암이 퍼지지 않는 환자는 80 퍼센트가 될 것이다. 지금 odds 비는 거의 5와 같다. 이것은, 암이 퍼저있는 odds는 약 작은 암 환자에서 보다 큰 암 환자에서 약 5 배가 더 높다:

Odds for Large Tumor Patients   = 5 × (Odds for Small Tumor patients)
= 5 × 20/80
= 100/80

해석: 큰 종양을 가지는 180명의 환자 중, 우리는 암이 퍼지지 않는 단지 80명의 환자와 비교되는 암이 퍼지는 100명의 환자를 보게 될 것이라 기대할 것이다. 만약 우리가 큰 종양을 가지는 전체 100명의 환자에서 이것을 비교하면(180명 환자 대신에) 우니는 암이 퍼지지 않은 환자 약 45명과 비교되는 암이 퍼진 약 55명의 환자를 보게 될 것이라 기대할 것이다.

독립 (Independence): 만약 odds 비가 1이면, 암이 퍼진 odds는 큰 그리고 작은 종양 환자에서 같다는 것을 주의해라. 다른 말로, 암이 퍼진 odds는 암의 크기에 의존하지 않는다. 기억해라

ODDS RATIO = e^β₁.

그러므로, 만약 odds 비가 1이면, 이것은 e⁰ = 1 이기 때문에 로지스틱 회귀 기울기 β₁ = 0 을 의미한다.

상대 위험(Relative Risk): 상대 위험의 정의:

상대 위험(Relative Risk) = p(1)
p(0)  

이 전립선 암 예제에서, 상대 위험의 추정은

p(1)^ p(0)^ = 0.5555 0.1923 = 2.89  

이것은 5.25의 odds 비 보다 상당히 작다. 상대 위험(2.89)는 큰 종양 환자에서 암이 퍼지는 확률은 작은 암 환자에서 비교하여 거의 3 배 더 높다는 것을 의미한다.

x가 연속일 때 기울기 해석(Interpreting the slope when x is continuous): 딱정벌레 사망율 예로 돌아가서, β^ = 0.2494를 상기하라. 만약 우리가 한 단위의 량을 증가하면, 그러면 x에서 변하는 한 단위에서 odds 비는 다음과 같이 보여진다 (약간의 대수학) :

p(x + 1) ⁄ (1 - p(x + 1)) p(x) ⁄ (1 - p(x)) = e^β₁  

CS₂ 농도에서 일 단위의 증가에서 추정된 odds 비는:

e^{β₁^} = e^0.2494 = 1.283,

사망 odds는 CS₂ 농도의 각 추가적인 단위 증가에서 약 1.283 배 더 크다고 나타남.

1.4 다중 로지스틱 회귀 (Multiple Logistic Regression)

전통적인 회귀의 설정에서 우리는 회귀 모델에서 하나 보다 더 많은 회귀 변수를 어떻게 포함시키는지 보았다. 다중 회귀는 게다가 로지스틱 회귀 모델로 또한 통합할 수 있다. p 회귀 변수 x₁, x₂, ... , x_p을 가진다고 상상하자. 그러면 우리는 (2) 식을 일반화 할 수 있고 그리고 다중 회귀 함수 정의는:

p(x₁, x₂, ... , x_p) = e^{β₀+β₁x₁+ ... +β_nx_n}
1 + e^{β₀+β₁x₁+ ... +β_px_p} (4)

그리고 p(x)의 로짓은

logit(p(x)) = β₀+β₁x₁+ ... +β_px_p

이다.

최대 우도는 단순한 로지스틱 회귀에서 수행한 것 처럼 β_j 그리고 다중 로지스틱 회귀 모델에서 그들의 표준 에러를 추정하기 위해 일반적으로 사용된다

다중 로지스틱 회귀의 중요성 검증 (Tests for Significance for Multiple Logistic Regression). 다만 회귀의 경우에서, 우리는 회귀 계수의 집합이 0과 다른지를 결정하기 위해 통계적 검증을 또한 수행할 수 있다. 다중 회귀 그리고 다중 로지스틱 회귀 둘 모두의 검증 절차는 같은 중요한 방법을 기반으로 한다: 전체 모델 적합 그리고 모델 축소 그리고 2개의 적합의 비교. 만약 축소된 모델이 전체 모델 만큼 비슷하게 좋다면, 그러면 축소 모델이 선호된다. 다중 로지스틱 회귀에서 검증 절차의 실제 역학은 우리가 논의한 다중 회귀와 다르다.

우도비 검증 (Likelihood Ratio Test). 로지스틱 회귀 모델은 일반화 선형 모델의 특별한 경우이다. 일반화 선형 모델에서 편차(deviance)라 불리는 통계량은 적합된 모델의 예측 값과 원시 데이터에서 실제 값과 얼마나 근접하게 일치하는지 측정되는 값을 비교한다. 최대 우도는 일반화 선형 모델에서 변수를 추정하기 위해 일반적으로 사용 된다. 우도(likelihood)는 그들의 추정치로 대체된 변수를 가지는 관측된 데이터에서 계산되는 단순한 확률 밀도(probability density) 이다. 극단적인 경우 변수의 수와 관측의 수가 같은 포화 모델(saturated model)을 적합해야 한다. 통계량의 기본적인 목적 중 하나는 가능한 작은 변수를 가지는 단순한 모델을 결정한다. 포화 모델은 가능한 관측과 같은 많은 변수를 가진다 따라서 그것은 전혀 단순하지 않게 제공된다. 어쨌든, 우리는 제안된 모델이 데이터를 얼마나 잘 적합하는지 결정하기 위해 포화 모델로 제안된 어떤 모델을 비교할 수 있다. 편차 D가 정의된다:

D = 2[log-likelihood(saturated model) − log-likelihood(proposed model)]

만약 제안된 모델이 참에 좋게 접근한다면, 그러면 편차는 상대적으로 작게 될 것이고 그리고 편차 D의 표본 분포는 근사적으로 n - p - 1 자유도에서 카이-제곱 분포(chi-squared distribution)를 따른다. 이것은 표본 크기 n이 무한대이므로 유효한 점근적인 결과(asymptotic result)의 의미이다. 어쨌든, 이 점근적인 결과(asymptotic result)는 일반적으로 그렇게 많이 사용되지 않는다. 대신에, 관심은 충첩 모델(nested models)과 비교에 있다 - 전체 모델에서 축소 모델을 비교. 중요성은 다중 회귀에서와 같다. 전체 모델

logit(p( x₁, x₂, … , x_p)) = β₀+β₁x₁+ … +β_qx_q+ … +β_q+1x_q+1+ … +β_px_p

을 고려해라 그리고 우리는 이러한 계수의 적어도 하나에서 귀무 가설(null hypothesis)

H₀: β_q+1 = ··· = β_p = 0,

대 적어도 이런 계수 중 하나는 0과 다르다는 대립 가설(alternative hypothesis) 검증을 원한다. 만약 H₀이 참이면, 그러면 회귀 변수 x_q+1 = ··· = x_p는 전체 모델에서 중복된다 그리고 제외할 수 있게 된다. 실제에서 H₀ 검증을 하기 위해, 해야할 필요가 있는 하나는 전체 모델과 축소된 모델을 적합하고 그리고 그들 각각의 편차를 비교한다. 검정 통계량:

Χ² = D_reduced - D_full

이다.

포화 모델(saturated model)에서 로그-우도의 평가 방정식 포함 그리고 우리가 편차의 차이(reduced - full)를 가졌을 때, 포화 모델의 우도를 제외하고 취소되는 둘 모두 위 방정식의 편차를 주목해라. 그러므로, 검증 통계량 Χ²은 다음으로 기술환다:

Χ² = 2[log-likelihood(full model) − log-likelihood(reduced model)] (5)

만약 H₀가 참이면, 그러면 검증 통계량 Χ²은 근사적인 카이-제곱 분포를 가진다(표본 크기가 충분히 크게 제공됨), 자유도는 전체 모델과 축소 모델 간에 변수의 수에서 차이와 같다: p - q. 만약 H₀가 거짓이면, 그러면 검증 통계량이 자유도 p - q에서 카이-제곱 분포에서 유도한다고 고려하기에 너무 크게 되는 경향이 있다. 이것은, 만약 Χ² 이 너무 크면, H₀는 기각된다. 만약 우리가 유의 수준 α에서 검증한다면, 그러면 우리는 H₀ > Χ_{α, p-q} 이면 H₀를 기각한다, p-q 자유도에서 카이-제곱 분포의 α 임계값.

(5)에서 주어진 검증 통계량은 우도 비 검정의 개념에 기반한다. 이 검정은 전체와 축소된 모델의 관측된 우도를 비교한다. 그러므로 축소 모델은 전체 모델의 제약 버전이다, 축소 모델의 우도는 전체 모델의 우도와 작거나 같게 될 것이다. 두개의 모델은 그러면 그들의 관측된 우도의 비(축소/전체) 의 조건에서 비교된다. 비의 알고리즘을 취함은 그것을 로그-우도(log-likelihoods)의 차이로 돌린다. 만약 로그-우도(log-likelihoods)에 -2를 곱하면, 그것은 귀무 가설 하에서 그것의 표본 분포는 전체와 축소 모델 간의 변수의 수의 차이와 같은 자유도를 가지는 점근적인(asymptotically) 카이-제곱이다. 검증 통계량 (5)는 때때로 로그-우도 비 통계량(log-likelihood ratio statistics) 또는 LLRS로 언급된다.

추가적으로, 유의 검정(significance tests)은 전통적인 회귀에서 부분 t-통계량(partial t-statistics)과 비슷한 Wald 통계량의 계산으로 개별적인 회귀 계수에서 수행된다 (i.e. H₀ : β_j = 0)

ω_j = β_j^
se^(β_j^)

귀무 가설(null hypothesis) β_j = 0 하에서, Wald 검증 통계량 ω_j는 근사적으로 표준 정규 분포를 따른다 (그리고 그것의 제곱은 근사적으로 자유도 1에서 카이-제곱이다). 이런 검증 아이디어를 설명하기 위해, 우리는 예제를 가지고 지금 설명한다.

Example. (open inbreeding.r) 근친교배(Inbreeding)는 식물 집단(plant populations) 내에서 자연적으로 발생한다. 연구는 식물 Yellow Monkeyflower를 이용하여 캘리포니아 나파 카운티(Napa County)의 토종 초식동물(native herbivores)에서 식물의 저항과 내성(the resistance and tolerance)에 대한 근친교배의 영향을 연구하기 위해 수행되었다. 이 연구에서 관심있는 응답 변수 y는 식물이 꽃을 피우는지 아닌지 이다 (0 없음 그리고 1 피움). 꽃 생산은 재생을 위해 필요하다. 관심있는 주요 설명 변수는 식물이 근친교배인지(1의 값) 혹은 잡종(cross-bred)인지(0의 값을 가짐) 아닌지를 가리키는 지시자 변수(indicator variable) 이다. 다른 두 공변량(covariates)을 또한 기록한다. spittlebugs 때문에 식물에서 초식 동물의 피해, 성체와 애벌레 딱정벌레, 굼벵이, 사슴 등, 이산 척도(discretized scale)를 백분률로 기록하였다 (12개 범주); 그리고 건조된 위의 식물의 토양 바이오매스(단위 그램). 데이터의 첫번째 몇 줄은 아래의 테이블에서 보여진다 (C. Ivey에 감사):

Damage	inbred	y	Biomass
0.2138	0	1	0.0525
0.2138	0	1	0.0920
0.2138	1	1	0.0239
0.2138	0	0	0.0149
0	0	0	0.0410
0.3907	0	0	0.0264
0.2138	0	1	0.1370
0.2138	1	1	0.0280
0.2138	1	1	0.0248
0.2138	1	1	0.0892
0.2138	1	0	0.0123
0	0	0	0.0105
. . .	. . .	. . .	. . .

회귀 변수 biomass 그리고 damage는 응답을 매우 다양하게 설명을 한다고 생각된다 그리고 그들은 공변량(covariates)으로 포함된다. 회귀의 주요 관심은 식물이 근친교배인지 아닌지의 지시자이다. 전체 모델에서 로짓은:

logit(p(INBRED, DAMAGE, BIOMASS) = β₀+β₁INBRED + β₂DAMAGE + β₃BIOMASS .

(이 모델은 어떤 상호작용 조건을 가지고 있지 않다는 것을 주의하라 - 의미있는 변수는 없다.) 여기에 다중 회귀 모델을 적합하는 R 코드가 있다:

dat = read.table("inbreeding.dat")
y = dat[,3]
damage=dat[,1]
inbred = dat[,2]
biomass = dat[,4]
fitfull = glm(y ~ inbred+biomass+damage, family=binomial)
summary(fitfull)

이 "전체" 모델을 실행한 R 출력은 아래에 주어진다:

Call:
glm(formula = y ~ inbred + biomass + damage, family = binomial)

Coefficients:
             Estimate  Std. Error  z value  Pr(>|z|)
(Intercept)   -1.1230      0.6795   -1.653    0.0984  .
inbred         0.3110      0.6821    0.456    0.6484  
biomass       41.3349     15.9374    2.594    0.0095  **
damage        -1.8555      2.1420   -0.866    0.3864
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 69.235 on 49 degrees of freedom
Residual deviance: 52.463 on 46 degrees of freedom
AIC: 60.463

Number of Fisher Scoring iterations: 6

이 전체 모델에서 로그-우도의 음의 2배는 잔차 편차(Residual deviance)이라 불리는 출력인 52.463 이다. inbred 그리고 damage에서 계수의 유의(수준)의 Wald 검정은 이런 회귀의 둘 모두 전체 모델에서 여유(redundant) 있게 나타나는 지시자 p-값 p = 0.6484 그리고 p = 0.3864 이 나온다. 예측자 집합의 redundancy를 위한 동시에 테스트를 진행하는 방법을 설명하기 위하여, 우리는 H₀ : β₁ = β₃ = 0을 검정하기 위해 inbred 와 damage 없는 축소된 로지스틱 회귀 모델을 적합할 수 있다

fitreduced = glm(y ~ biomass, family=binomial)
summary(fitreduced)
x2= 2*(logLik(fitfull)-logLik(fitreduced)) # log-likelihood ratio test statistic
as.numeric(x2)
pval=1-pchisq(x2,2)
as.numeric(pval)

축소 모델을 실행하고 우도 비 통계량을 계산한 출력 부분은

Coefficients:
             Estimate  Std. Error   z value  Pr(>|z|)
(Intercept)   -1.2556      0.5129    -2.448    0.0144  *
biomass       37.9157      14.7390    2.572    0.0101  *

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 69.235 on 49 degrees of freedom
Residual deviance: 53.579 on 48 degrees of freedom
AIC: 57.579

> x2= 2*(logLik(fitfull)-logLik(fitreduced)) # log-likelihood ratio test statistic
> as.numeric(x2)
[1] 1.116378
> pval=1-pchisq(x2,2)
> as.numeric(pval)
[1] 0.5722445

이 축소된 모델에서 로그-우도의 음의 2배는 53.579 이고 그리고 로그-우도 비 통계량의 값은 Χ² = 1.1164 이다. 그러므로 전체 모델과 축소된 모델은 2개의 파라메터가 다르다, 우리는 자유도 2에서 카이-제곱 분포에서 이 검정 통계량을 비교할 수 있다. 이 검정에서 p-값은 p = 0.5723 이다. 그러므로 우리는 계수 β₁ 와 β₃ 이 0 과 다르다는 증거로 불충분하다 라고 결론내린다 (귀무 가설 참). 이것은 우리가 예측자로 biomass 만 포함하는 로지스틱 회귀 모델을 사용하도록 허락한다. 주어진 biomass에서 꽃을 피우는 식물의 추정된 확률:

exp{−1.2556 + 37.9157 Biomass}
1 + exp{−1.2556 + 37.9157 Biomass}

이 분석에서, 이것은 이 식물 종은 근친 교배인지 아닌지는 그것의 생식 능력에 영향이 없다고 나타난다, 식물의 크기 그리고 초식동물에서 지속되는 손상을 고려했을 때 조차도. 이 모델에서 근친교배 (x₁)의 지시자만 사용하여 로지스틱 회귀가 실행될때 유의하는 것이 유용하다, 회귀 x₁ 은 여전히 유의하지 않다 (p = 0.9442). 적합된 로지스틱 회귀 곡선을 가지는 응답 y 대 biomass 의 그림은 Figure 2에 보여준다. 0과 하나의 응답값 간에 더 나은 구별을 위해, "jitter" 명령어가 plot 문장에 사용되었다:

plot(biomass,jitter(y,.1), main="Monkey Flower Inbreeding Data",
     ylab="Probability of Flower", xlab="Died Biomass (g)")

Figure 2: Yellow Monkey flower 식물의 로지스틱 회귀 plot. 만약 식물이 적어도 하나의 꽃을 피우면 응답은 1 그리고 꽃을 피우지 않는다면 0 이다. 응답은 식물 biomas에 대해 그려졌다. 또한 추정된 로지스틱 회귀 곡선이 보인다. 0-1 응답은 "jittered" 되었다 그래서 약간 위로 보인다.

추정된 기울기는 양수이기 때문에, 우리는 번식 확률이 식물 크기가 증가함에 따라 증가되는 것을 본다.

추정된 기울기 β₁^ = 37.9157 이고 그리고 교차 비(odds ratio)를 계산하기 위해 이것을 사용하는 것은 천문학적으로 엄청난 수를 만들어 낼 것이다. 그러므로 식물의 건조된 biomass는 0 에서 0.4 사이의 값 범위이다, 그것은 교차비를 계산하기 때문에 더 좋게 된다, biomass에서 단위 증가에서 아님, 그러나 다수의 작은 증가에서, 0.1 그래의 증가를 말함. 이 경우에, 교차비는

e^{0.1 β₁^} ≈ 44.3

이 될 것이다, 0.1 씩 추가마다 건조된 biomass 에서 증가한다, 적어도 하나의 꽃을 생산하는 교차(odds)는 인자 44.3 씩 증가한다.

2 일반화 선형 모델: 일반 설정
(Generalized Linear Models: The General Setup)

처리 방법의 비교에서 분산 모델의 일반적인 분석 뿐만 아니라 선형 회귀 모델을 포함하는 전통적인 선형 모델은 다음과 같이 표현할 수 있다:

y = β₀+β₁x₁+ ... +β_px_p + ε (6)

팩토리얼 아노바 모델(Factorial ANOVA models)은 회귀 변수가 지시자 변수일 때 결과이고 그리고 우리는 회귀가 연속적이거나 혹은 연속/지시자 변수의 혼합일 때 회귀 모델을 얻는다. (6) 식에서 응답 y 는 모델 계수의 선형 함수이다, β_j는 그런 이유로 이름이 선형 모델(linear model) 이다. (6) 식에서 에러는 종종 정규 분포를 가진다고 가정한다. 식 (6)은 에러가 비-정규인 상황을 처리하기 위해 일반화 할 수 있다. 일반화 설정에서, 응답 변수 y는 (6)에서 직접적으로 연관되지 않는 점을 제외하고 선형 조합이 β₀+β₁x₁+ ... +β_px_p 이므로 여전히 회귀에 연관된다. 대신에, 응답은 링크(link) 함수로 선형 조합 β₀+β₁x₁+ ... +β_px_p에 연결된다.

로지스틱 회귀 모델에서, 링크 함수는 로짓(logit) 이다. 로지스틱 회귀 설정에서, 응답 y는 성공 확률 p(x₁, ... , x_p)가 회귀의 집합 x = (x₁, ... , x_p)에 의존하는 0-1 베르누이 랜덤 변수이다. 식:

p(x₁, ... , x_p) = e^{β₀+β₁x₁+...+β_px_p}
1 + e^{β₀+β₁x₁+...+β_px_p} (7)

그리고 로짓 링크 함수는 다음으로 정의:

logit( p(x₁, ... , x_p)) = log[ p(x₁, ... , x_p) ⁄ (1 - p(x₁, ... , x_p)) ] = β₀+β₁x₁+...+β_px_p

전통적인 회귀 모델 식 (6)에서, 링크 함수는 단순한 식별 함수(identity function) 이다.

로지스틱 회귀 그리고 식 (6)은 일반화 선형 모델(generalized linear models)의 2가지 예이다. 일반화 선형 모델의 다른 예는 포아송 회귀(Poisson regression), (분할표(contingency tables)의 셈(count) 데이터 예에서) 로그 선형 모델(log-linear models), 생존 분석 모델 이다. 선형 모델의 중요한 3가지 기본:

1. 우리는 지수 군(exponential family)에서 독립 응답 변수의 표본 Y₁, ... , Y_n를 가진다. 기본적으로, 확률 분포의 지수 군은 밀도 함수가 지수함수로 표현한 것이다 그리고 밀도의 지지도(i.e., 밀도가 0보다 큰 포인트)는 모델의 변수에 의존하지 않는다 (정규 그리고 이항 분포는 둘 모두 지수 군이다).
2. 응답 변수 Y의 선형 모델 예측자 β₀+β₁x₁+...+β_px_p 이 있다.
3. 평균 응답은 링크 함수 g를 통해 선형 예측자에 의존한다. x_j의 값에 의존하는 이 평균 응답은 조건부 기대 μ_y|x = E[Y|x₁, ... , x_p] 로 알려져 있다 그리고
   g( μ_y|x) = β₀+β₁x₁+...+β_px_p
   식 (6)에서, 함수 g는 식별자 g(x) = x 이다 그리고 로지스틱 회귀에서, p(x)는 회귀 변수 x 값이 주어진 베르누이 응답 Y의 조건부 평균이다.

통계적 추정 그리고 일반화 선형 모델 검정은 우도 이론을 사용하여 평범하게 수행한다, 이것은 로지스틱 회귀의 경우와 같다. 일단 응답에서 분포가 결정되면, 우도 방정식은 작성할 수 있고 그리고 보통의 최대 우도 추정 절차가 그다음 따라온다. 추가적으로, 모델 검증은 분산의 차이에 기반한 우도 비 검정을 사용하여 일반적으로 수반되는 다중 로지스틱 회귀 예제를 가지고 보여 주었다. 최대 우도 추정의 사용은, 귀무 가설(null hypotheses) 하에 검정 통계량은 큰 표본 크기에서 카이-제곱 분포를 따른다.

2.1 과대산포 (Overdispersion)

전통적인 선형 모델에서(e.g. 회귀 또는 ANOVA), 에러 분산은 관측 전체에서 일정하다고 추정한다. 이 가정은 연습에서 종종 위반된다. 회귀에서, 하나는 균일하지 않은 분산의 차이를 조정하기 위해 가중 최소 제곱으로 적합할 수 있다. 일반화 선형 모델에서, 우리는 응답이 균등한 분산을 가진다고 기대하지 않았다. 예를 들어, 이 장 처음에서 딱정벌레 예에서, 이항 응답은 n_i가 살충제에 노출된 딱정벌레 수에서 n_ip(1 − p)와 같은 분산을 가진다. n_i의 다른 값은 다른 분산을 초래한다. 그럼에도 불구하고, 그것은 모델에 의해 예상되는 것을 초과하는 응답의 변화는 일반화 선형 모델에서 상당히 일반적인다. 이에 대한 용어는 과대산포(overdispersion)이다.

과대산포(overdispersion)은 몇가지 요인으로 발생할 수 있다. 예를 들어, 로지스틱 회귀에서, 관측 사이에 독립성 결여. 베르누이(i.e. 0-1) 응답 간에 양의 공분산(positive covariance)은 분산을 부풀리게 할 수 있다. 만약 과대산포가 문제이면, 과대산포 파라메터를 추정할 수 있고 그리고 로그-우도 비 검정은 대략적인 F-검정 결과인 이 추정된 과대산포(overdispersion) 분산으로 나누어 조정할 필요가 있을 것이다.

우리는 간단히 일반화 선행 모델의 주제를 간략하게 소개했다. 적합도(goodness of fit), 모델 생성 등과 같은 회귀의 많은 문제는 또한 일반화 선형 모델이 가지는 문제이다. 이 주제의 상세한 내용은 일반화 선형 모델에 초점을 맞춘 책에서 찾을 수 있을 것이다. 이 주제로 전통적인 책은 McCullagh와 Nelder의 1989년도 서적 Generalized Linear Models 이다 (McCullagh and Nelder, 1989).

우리는 일반화 선형 모델의 다른 매우 다양한 주제를 가지고 이 장을 접는다.

3 포아송 회귀 (Poisson Regression)

포아송 확률 분포(Poisson probability distribution)는 아마도 셈 데이터 모델링에서 가장 일반적으로 사용되는 이산 분포(discrete distribution) 이다. 예로, 연구는 입원치료, 응급실, 의사 방문(physician visits), 그리고 장기간 치료 방문에서 측정된 풍토성 위장염의 발생(occurrences of endemic gastroenteritis)이다. 하나의 관심은 수질 측정(water quality measures)과 연관된 위장염에 있다. 만약 우리가 Y를 특별한 시간 기간에 이 질병의 발생 수를 나타내도록 놓는다면, 그러면 포아송 분포는 Y의 분포를 모델링하는 합리적인 방법이 될 것이다. 만약 Y가 위방염 경우의 수와 같다면, 그러면 Y는 값 0, 1, 2... 으로 추정할 수 있음을 주목하라. 포아송 분포는 비율 변수 μ > 0으로 변수화 된다 그리고 확률 밀도 함수 f(y)는 다음으로 주어진다:

f(y) = P(Y = k) = e^-μ μ^k k! , k = 0, 1, .... (8)

포아송 랜덤 변수(Poisson random variable)에서 평균과 분산은 μ와 같다. 포아송 확률 모델은 아래의 3개의 조건인 상황에서 이론적으로 만들어 질 수 있다.

1. "시간" 구간에 중첩되지 않는 관심 사건의 발생은 독립적이다.
2. 작은 시간 구간에서 2개 혹은 더 많은 사건 확률은 작다, 그리고
3. 사건이 시간의 짧은 구간에서 발생하는 확률은 시간 구간의 길이에 비례한다.

하나는 포아송 분포는 확률 분포의 지수 분류에 속한다 따라서, 일반화 선형 모델 접근으로 포아송 응답을 예측자 변수와 연관하여 사용될 수 있는 것을 보여 줄 수 있다.

포아송 응답 y는 회귀 변수의 집합 x₁, ... , x_p 에 의존한다. 왜냐하면 평균은 양수이어야 한다, 예측자가 주어진 응답 y의 조건부 기대를 모델을 만드는 자연스런 방법:

μ = E[y|x₁, ... , x_p] = e^{β₀+β₁x₁+...+β_px_p} (9)

만약 우리가 양 변에 자연로그를 취한다면, 다음 식을 얻는다:

log(μ) = β₀+β₁x₁+...+β_px_p

위 식의 자연 로그는 포아송 회귀에서 링크 함수이다. 이것은 로그-선형(log-linear) 모델의 예이다.

만약 데이터 셋이 회귀 변수의 집합에서 측정치 뿐만 아니라 셈에 상응하는 응답 변수 y에서 측정치로 구성된다면, 하나는 전통적인 모델을 사용하는 이것을 단순하게 모델화 하기를 원할 것이다. 어쨌든, 만약 표준 회귀가 셈 데이터를 모델로 사용한다면, 하나는 평균이 균일하지 않은 에러 분산을 가지는 이분산성(heteroscedasticity)을 가지는 문제를 종종 접하게 된다. 표준 회귀 설정에서 가정의 하나는 회귀자 변수에 관계 없이 에러 분포의 분산이 동일하다는 것을 상기해라. 어쨌든, 위에서 강조한 것처럼, 포아송 랜덤 변수의 분산는 μ 이다. 포아송 분포는 그것이 평균과 분산이 서로 같은 점이 독특하다. 만약 포아송 평균이 식 (9)와 같이 회귀자 변수와 관계가 있다면, 그러면

var(Y) = e^{β₀+β₁x₁+...+β_px_p}

는 분명하게 분산은 게다가 회귀자 값에 의존한다고 가리킨다 따라서 동일 분산 가정은 포아송 셈 데이터(Poisson count data)를 보유하고 있지 않을 것이다.

최대 우도 검정(Maximum likelihood estimation)은 전형적으로 포아송 모델의 변수를 추정하기 위해 사용된다. 설명하기 위해, y_i를 랜덤 표본 i = 1, ... , n 에서 예측자 x_i에 의존하는 포아송 랜덤 변수로 놓아라. 그러면 우도는

L(β₀, β₁) = Π_i=1ⁿ e^−μ_i μ_i^y_i ⁄ y_i!

μ_i는 식 (9)에서 주어지는 경우. 우도에 로그를 취함여 아래의 로그-우도 함수를 얻는다:

l(β₀, β₁) = Σ_i=1ⁿ [−μ_i + y_ilog(μ_i) − log(y_i!)]

그리고 목표는 이 함수를 최대로 하는 β₀와 β₁를 찾는 것이다. 우리는 R에서 "gml" 그리고 "family=poisson"를 지정하고 사용하여 이것을 할 수 있다. 대안으로, 하나는 SAS에서 PROC GENMOD을 사용하여 할 수 있다.

따라오는 예는 포아송 회귀(Poisson regression)를 설명하기 위해 사용될 것이다.

Example. C. dubia 동물 플랑크톤(Zooplankton)을 리틀 마이에미 강(Little Miami River)에서 다른 니켈의 량(μgram⁄liter)에 노출 되었다. 하나의 량은 대조군으로 사용하는 0 이고 그리고 가장 높은 량은 34 μg/l 이다. 응답 변수는 7일간에 번식된 자손의 수 이다. 각 량에서 10 마리의 암컷이 노출되었다. 데이터 셋 zooplanktonLMR.dat는 다음 컬럼을 가진다: column 1: 량(dose), columns 2-9: 각 주어진 날짜에서 암컷이 생존했는지 죽었는지를 나타내는 지시자(0=dead, 1=alive) 그리고 column 10: 동물성 플랑크톤의 수 (data 제공: Custer, K.W. and Bur- ton, G.A.). Figure 3은 동물성 플랑크톤 대 니켈의 량에서 원시 데이터의 산점도(scatterplot)을 보여준다. 우리는 응답이 셈(count)이고 그리고 단지 보통의 최소-제곱 회귀를 적합한 것을 무시할 수 있다. Figure 3에서 명백히, 곧은 선 모델은 셈과 량 간의 관계에서 안정적으로 근사적으로 제공되지 않았다. 만약 우리가 이차 선형 모델(quadratic linear model)을 적합한다면, 우리는 Figure 4에서 보여주는 잔차 플롯을 얻는다. 이 2차 선형 모델의 적합은 명백하게 좋은 전략이 아니다, 그러나 그것은 쌍 포인트(couple points)를 설명하는데 도움이 된다. Figure 4의 잔차 플롯에서 우리는 이차 모델이 좋은 적합이 아님을 나타내는 그림에서 구조를 알아낸다. 그러므로 번식의 수는 니켈의 량에서 지수적인 감소로 나타난다, 그것은 이차 모델은 신통치 않은 성능이란 사실은 놀랍지가 않다. 추가적으로 어쨌든, 잔차 플롯(residual plot)은 포아송 타입 데이터에서 기되되는 잔차에서 매우 분명하게 비-균등 분산을 보인다.

Figure 3: 암컷이 노출된 7일 후에 동물성 플랑크톤 번식 수 대 니켈의 수준에 대한 산점 도

Figure 4: 최소-제곱 이차 모델 적합의 잔차 그림

우리는 아래의 R 코드를 이용하여 니켈 수준에서 번식 수의 회귀에 대해 단순한 포아송 회귀를 적합할 수 있다.

# lmr stands for Little Miami River
lmr <- read.table("zooplanktonLMR.dat", header=F)
y <- lmr[,10]
dose <- lmr[,1]
fit <- glm(y ~ dose, family=poisson)
summary(fit)
plot(dose, y, xlab="Nickel Dose (mu gram/liter)", pch=19,
     ylab="Number of Offspring",
     main="Zooplankton Offspring: Little Miami River Site")
nickel=seq(0,40, by=1)
yhat1=exp(coef(fit)[1]+coef(fit)[2]*nickel)
lines(nickel, yhat1, lwd=2, col=2)

이 모델을 summary 명령으로 출력:

Call:
glm(formula = y ~ dose, family = poisson)

Coefficients:
             Estimate  Std. Error  z value  Pr(>|z|)
(Intercept)   3.38150     0.04959    68.19    <2e-16 ***
dose         -0.18225     0.01092   -16.69    <2e-16 ***

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 872.94 on 59 degrees of freedom
Residual deviance: 251.32 on 58 degrees of freedom
AIC: 424.80

Number of Fisher Scoring iterations: 5

이 출력에서, 우리는 아래의 추정된 모델을 가진다:

y^ = e^{3.38 − 0.182x}

x가 니켈의 량인 경우. 량 (x)에서 계수는 0에 근접한 p-값(near-zero p-value)으로 나타나는 높게 유의(significant)하다. 원시 데이터와 적합 포화송 회귀 곡선의 그림은 Figure 5에서 보여준다.

Figure 5: 추정된 포아송 회귀 곡선에서 동물성 플랑크톤 데이터

선형 모델과 다르게, 포아송 회귀에서 기울기 계수를 해석하기 위해, 그것은 x에서 단위 증가로 (차이 대신에) 예측된 응답의 비를 보는 것이 더 잘 이해된다.

e^{β₀+β₁(x+1)}
e^β₀+β₁x = β₁

β₁^ = − 0.182를 가지는 동물성 플랑크톤 예에서, 우리는

e^{β₁^} = e^−0.182 = 0.8334

를 가진다. 그러므로, 니켈 량의 단위 증가에서, 우리는 펙터 0.8334으로 감소되는 번식의 수를 보게 되기를 기대할 수 있다.

전통적인 회귀 그리고 로지스틱 회귀 모델에서, 포아송 회귀 모델은 하나 보다 더 많은 예측자를 가지는 다중 포아송 회귀를 쉽게 일반화 할 수 있다.

3.1 과대산포 (Overdispersion)

포아송 변수의 분산은 포아송 랜덤 변수의 평균과 같다는 것을 상기해라. 포아송 회귀에서 일반적인 문제는 응답이 모델에서 기대되는 것 보다 더 잘 변화한다. 이것을 과대산포(overdisperson)라 불린다. 과대산포를 진단하는 하나의 방법은 고려되는 모델과 안정화된 모델 하에서 로그-우도에서 차이의 -2 배인 편차 통계량을 확인한다. 안정화 모델에서, 예측된 응답은 단지 실제 관측된 응답 y_i 이다. y_i^을 포아송 모델에서 예측된 응답을 나타내도록 놓아라, 하나는 편차를 보여줄 수 있다:

D = −2[log-likelihood(model) − log-likelihood(saturated)] = 2Σ[y_i log(y_i / ^y_i) − (y_i − ^y_i)] (10)

(만약 y_i = 0 이면 y_i log(y_i) = 0 인 점에 주목하라) 만약 모델이 데이터에 잘 적합한다면, 그러면 그것은 편차는 표본 크기 n 빼기(−) 추정된 계수의 수인 편차 자유도에서 대략 같게 되도록 따를 것이다: p가 예측자 수인 경우 자유도 df = n - p - 1. 점근적으로, 편차는 만약 모델이 옳바르게 명시된다면 이런 자유도에서 카이-제곱 분포를 따른다. R 에서, "glm" 함수는 이것을 "잔차 편차(Residual deviance)"라 부른다. 만약 잔차 편차가 잔차 자유도를 크기 초과한다면, 그러면 그것은 과대산포 문제를 가리킨다.

단순한 포아송 회귀에서 출력으로, 우리는 잔차 편차가 자유도 58 에서 251.32 을 확인한다. 그러므로 251.32는 58 보다 상당히 크다, 이것은 이 단순한 모델에서 과대산포 문제가 있을 수 있음을 나타낸다. 과대산포 문제가 나타났을 때, 회귀 계수의 추정된 표준 에러는 신뢰할 수 없다 그리고 그들(e.g. Wald tests)을 기반으로 한 추론은 믿을 수 없다.

추론 문제에서 하나의 해결은 실제 셈(counts)과 예측된 셈(counts) 비교에 기반한 과대산포 변수를 추정한다:

Φ = 1
(n - p - 1) ⁿΣ_i=1(y_i − ^y_i)² / ^y_i (11)

R 에서, 우리는 Φ를 계산하고 그리고 아래를 입력하여 모델을 조정할 수 있다:

> phi=sum((y-fit$fitted)^2/fit$fitted)/fit$df.residual
> summary(fit, dispersion=phi)
Coefficients:
             Estimate  Std. Error  z value  Pr(>|z|)
(Intercept)   3.38150     0.10013   33.771    <2e-16 ***
dose         -0.18225     0.02205   -8.264    <2e-16 ***

변수 추정은 윈래의 적합과 구별된다는 점에 주목하라, 그러나 표준 에러는, 과대산포의 조정 후, 더 크다 (e.g. 량 계수에서 추정된 표준 에러는 그것의 원래 값의 약 2배 이다).

과대 산포의 문제는 특성화된 모델을 유발할 수 있다. 예를 들어, 중요한 예측자 변수는 모델에서 제외될 수 있다. 과대산포의 다른 원천은 관측에서 상관관계(correlation)이다. 아마도 군집화 때문(e.g. 같은 공장 또는 계략에서 획득된 관측). 랜덤 효과는 상관관계를 이해하기 위해 이런 경우에 도입될 수 있다.

다른 이유로 과대산포는 데이터를 모델하기 위해 필요한 다른 셈(count) 분포이다. 예를 들어, 음수 이항 분포(negative binomial distribution)는 셈 데이터의 모델 생성에 필요할 수 있다, 3.2 절을 보라.

로그 변환 (Logarithm Transformation). 다른 어떤 회귀 모델를 가지고, 때때로 예측자 변수의 변환은 더 근접한 모델로 이끌 수 있다. 예를 들어, 예측자 x 에서 회귀를 통해 생성된 포아송 셈(count) y를 가정해라.

E[y|x] = e^β₀+β₁x

를 고려하는 대신에, 하나는 예측자의 로그 변환을 대신해 사용할 수 있다. 다음을 고려:

E[y|x] = e^{β₀+β₁log(x)}

그러나 이 모델은 다음 식으로 줄어든다:

E[y|x] = e^{β₀*β₁log(x)} (11)

β₀^* = log(β₀) 인 경우. 모델 (12)는 평균 응답이 예측자에 대한 지수 함수 대신에 x의 "다항(polynomial)" 함수이다.

3.2 음의 이항 모델 (Negative Binomial Model)

포아송 모델은 항상 셈(count) 응답에서 좋은 적합을 제공하지는 않는다. 대안 모델은 음의 이항 분포 이다. 성공 확률이 p 이고 그리고 우리가 성공 또는 실패의 독립 시행을 관측한 경우의 이항 실험에서 일반적인 설정을 고려해라. 우리가 r번 성공할 때까지 시행 과정을 관측한다고 가정해라. 그러면 r회 성공을 관측하기 위해 필요한 시행 W의 전체 수는 아래의 확률 분포를 가진다:

P(W = n) = ( n - 1
k - 1 ) p^r(1 - p)^n-r

응의 이항을 사용하는 회귀 적합을 위해, 포아송 대신에, 하나는 R에서 함수 "glm.nb" 함수를 사용할 수 있다. 이것은 아래의 식으로 평균 응답 μ에 연관된 추정된 계수를 생성할 것이다:

e^β₀+β₁x = μ
μ + r

r이 음의 이항 분포 함수 변수인 경우 (그리고 R에서 glm.nb의 결과에서 세타(theta)라 불림).

동물성 플랑크톤 예제를 다시 찾음, 니켈의 량에 따라 더 크게 얻는다는 것을 상기하라, 더 많은 암컷 플랑크톤은 자연스럽게 더 적은 자손을 낳고 죽는다. Y를 자손의 수를 타나내게 놓고 그리고 비율 변수 μ로 주어진 값에서, Y는 포아송 분포를 갖는다. 만약 비률 변수가 밀도 g(μ)를 가지는 감마 분포를 가지는 암컷 생존률의 우도와 관련된 랜덤 변수이면, 그러면 Y에서 결과 분포는 포아송의 무한 혼합물(infinite mixture) 이다 그리고 감마 분포는 다음으로 주어진다:

P(Y = y) = ∫₀^∞ f(y|μ)g(μ)dμ = ∫₀^∞ (e^-μμ^y / y!) g(μ)dμ

이 적분을 구함은 음의 이항 분포를 산출한다. 이러한 관점에서, 음의 이항은 포아송 회귀 모델의 훌륭한 일반화를 제공한다.

3.3 비율 모델 (Rate Models)

많은 셈(count) 회귀 모델에서, 실제 셈(count)은 관측에서 얻어지는 단위(unit)의 크기에 관계될 것이다. 예를 들어, 시설에서 찾은 알 묶음의 수는 시설에 남아있는 수와 연관된다. 동물성 플랑크톤 예에서, 그것은 7일간에 자손의 수는 암컷이 생존한 날의 수에 의존할 것이다. 예를 들어, 만약 기간이 암컷이 생존한 기간의 수와 같다면, 우리는 다음으로 기술할 수 있다:

log(y/days) = β₀ + β₁x

지수(Exponentiating), 우리는 모델

E[y|x] = e^{days + β₀ + β₁x}

을 고려할 수 있다. 이 경우에, 변수 days는 하나의 고정된 계수를 가진다. 포아송 회귀에서, 이것은 offset으로 불린다.

동물성 플랑크톤을 가지고 여기에서 표현된 예에서, 원본 데이터 집합에서, 각 암컷에서, 각 일자 (days 0-7)에서 암컷이 생존해 있는지 아닌지를 나타내기 위한 지시자 변수 (0 또는 1)을 주목하라. 우리는 생존 분석으로 알려진 통계 지류의 아이디어를 사용하여 이런 날짜 모델을 생성할 수 있다. f(y|x)를 주어진 량 x에서 암컷의 자손 수에서 확률 밀도 함수를 나타내게 놓자. d는 암컷이 생존할 날짜 수를 나타내게 놓자, 그러면 우리는 다음을 기술할 수 있다(표기의 약간의 남용):

포아송 회귀 모델은 f(y|x,d)를 추정하기 위해 사용될 수 있다 그리고 생존 분석 모델은 량 x에서 d 날 사이에 암컷 생존 확률을 모델을 생성하는데 사용할 수 있다, 위의 방정식에서 f(d|x)로 주어짐. 이 데이터 집합에서, 다수의 관측은 실험이 많은 암컷이 죽기 전에 끝을 의미하는 즉시 검열(right censored) 이다 따라서 우리는 많은 암컷에서 값 d를 알지 못한다 - 이것은 생존 분석에서 일반적으로 발생한다. 우리는 이 주제로 더 깊이 상세하게 가지 않을 것이다 그러나 관심이 있는 독자에게 생존 분석에서 몇가지 참조 서적 중 하나를 언급할 수 있다 (e.g. Hosmer and Lemeshow, 2008).

3.4 Zero-Inflated Poisson Regression

셈(count) 데이터를 가지고 작업할 때 다른 매우 일반적으로 일어나는 것은 포아송 모델과 일치하지 않는 제로 카운트의 과잉(an overabundance of zero counts) 이다. 예를 들어, 하나의 조사에서, 숲에서 나무를 관측하고 그리고 나무에서 덩쿨 수를 셈(count) 하였다. 덩쿨 수는 나무에서 나무로 0에서 30개의 덩쿨까지 다양하다. 어쨌든, zero-inflated 포아송 모델로 이끄는 덩쿨이 없는 많은 나무가 있다.

zero-inflated 포아송 모델을 생성하는 하나의 방법은 모집단이 2개의 하위 모집단으로 구성되게 명시하는 2개의 요소의 유한한 혼합을 통해서 이다: 하나의 하위 모집단은 0의 셈(counts)을 가지고 약화된다 그리고 다른 하위 모집단은 아래의 밀도를 가지는 포아송 분포를 따른다(e.g. Lambert, 1992; Min and Agresti, 2001)

p^*는 포아송 하위-모집단에서 모집단의 비율인 경우.

하나는 함수 zeroinfl를 가지는 R 패키지 pscl를 사용하여 zero-inflated 포아송 회귀 모델을 적합할 수 있다.