시작하기 전에
다음의 패키지가 설치되어 있지 않으면 설치한다.
# install.packages("MASS") # install.packages("nnet") # install.packages("tidyverse") # install.packages("patchwork") # for side-by-side ggplots library (MASS)library (nnet)library (tidyverse)library (patchwork)sessionInfo ()
R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.5
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Asia/Seoul
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] patchwork_1.2.0 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
[5] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
[9] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0 nnet_7.3-19
[13] MASS_7.3-60.2
loaded via a namespace (and not attached):
[1] gtable_0.3.5 jsonlite_1.8.8 compiler_4.4.0 tidyselect_1.2.1
[5] scales_1.3.0 yaml_2.3.8 fastmap_1.1.1 R6_2.5.1
[9] generics_0.1.3 knitr_1.46 htmlwidgets_1.6.4 munsell_0.5.1
[13] pillar_1.9.0 tzdb_0.4.0 rlang_1.1.3 utf8_1.2.4
[17] stringi_1.8.3 xfun_0.43 timechange_0.3.0 cli_3.6.2
[21] withr_3.0.0 magrittr_2.0.3 digest_0.6.35 grid_4.4.0
[25] rstudioapi_0.16.0 hms_1.1.3 lifecycle_1.0.4 vctrs_0.6.5
[29] evaluate_0.23 glue_1.7.0 fansi_1.0.6 colorspace_2.1-0
[33] rmarkdown_2.26 tools_4.4.0 pkgconfig_2.0.3 htmltools_0.5.8.1
일반화 선형모형 기초
일반화 선형모형(generalized linear models, GLM)
정규분포를 따르지 않는 반응변수와 선형성을 갖도록 하는 변환을 포괄하기 위해 선형모형을 확장 (McCullagh and Nelder, 1989)
가정
반응변수 \(y\) 는 각 고정된 값의 설명변수 \(x_1, \dots, x_p\) 에 대해 독립적으로 관측된다.
선형 예측변수 \(\eta = \beta_1 x_1 + \cdots + \beta_p x_p = \boldsymbol{\beta}^T\mathbf{x}\)
\(y\) 의 분포: 밀도함수 \[
f(y_i ; \theta_i, \varphi) = \exp[A_i \{ y_i \theta_i - \gamma (\theta_i) \} / \varphi + \tau (y_i, \varphi / A_i)]
\]
\(\varphi\) : 척도 모수 (알려진 값일 수 있음)\(A_i\) : 알려진 사전 가중치모수 \(\theta_i\) 는 선형 예측변수에 의존
평균 \(\mu = \operatorname{E}[y]\) 는 선형 예측변수의 부드러운 가역 함수 \[
\mu = m(\eta), \qquad \eta = m^{-1} (\mu) = \ell (\mu)
\]
\(\ell(\cdot)\) 는 연결 함수 (link function)\(\theta = (\gamma')^{-1}(\mu)\) 의 관계
척도모수
\(\varphi\) 가 알려진 값이라면, \(y\) 의 분포는 일모수 표준 지수족이다.\(\varphi\) 가 알려지지 않은 값이라면, \(\varphi\) 는 적률법으로 장애모수로 다룬다.
예시
정규분포. 평균 \(\mu\) , 분산 \(\sigma^2\) \[
\log f(y)=\frac{1}{\varphi}\left\{y \mu-\frac{1}{2} \mu^2-\frac{1}{2} y^2\right\}-\frac{1}{2} \log (2 \pi \varphi)
\]
\(\varphi = \sigma^2\) \(\theta=\mu\) , \(\gamma(\theta)=\theta^2 / 2\)
포아송 분포. 평균 \(\mu\) \[
\log f(y) = y \log \mu - \mu - \log (y!)
\]
\(\theta = \log \mu\) , \(\varphi = 1\) , \(\gamma (\theta) = \mu = e ^ \theta\)
이항 분포. 시행 횟수 \(a\) , 성공확률 \(p\) , 반응변수 \(y = s/a\) , 성공횟수 \(s\) \[
\log f(y) = a \left[ y\log \frac{p}{1-p} + \log (1-p) \right] + \log \binom{a}{ay}
\]
\(A_i=a_i\) , \(\varphi=1\) , \(\theta = \log \frac{p}{1-p}\) , \(\gamma (\theta) = - \log (1-p) = \log (1+e^{\theta})\)
반응분포와 연결함수의 조합을 일반화선형모형족(family)라고 한다.
분포족과 연결함수
기본 연결함수는 D로 나타냄
Link
binomialGammagaussianinverse.gaussianpoisson
logitD
probit\(\bullet\)
cloglog\(\bullet\)
identity\(\bullet\) D
\(\bullet\)
inverseD
log\(\bullet\) D
1/mu^2D
sqrt\(\bullet\)
GLM에 대해 잘 알려진 결과
가능도 함수 (표본 크기 \(n\) ) \[
\ell(\theta, \varphi ; Y)=\sum_{i=1}^n \left[ A_i\left\{y_i \theta_i-\gamma\left(\theta_i\right)\right\} / \varphi+\tau\left(y_i, \varphi / A_i\right) \right]
\]
점수함수 \[
U(\theta)=A_i\left\{y_i-\gamma^{\prime}\left(\theta_i\right)\right\} / \varphi
\]
평균과 분산 \[
\operatorname{E}\left(y_i\right)= \mu_i = \gamma^{\prime}(\theta_i),
\quad
\operatorname{var}\left(y_i\right)=\frac{\varphi}{A_i} \gamma^{\prime \prime}(\theta_i)
\] 따라서 \(\theta\) 도 \(\mu\) 의 가역함수이고 \(\theta = ({\gamma^{\prime}})^{-1} (\mu) = (\ell \circ \gamma^\prime)^{-1}(\eta)\) . \[
\operatorname{E}\left(\frac{\partial^2 \ell(\theta, \varphi ; y)}{\partial \theta \partial \varphi}\right)=0
\] 즉 \(\theta\) (또는 \(\beta\) )와 \(\varphi\) 는 직교 모수 .
분산함수: \(V(\mu)=\gamma^{\prime \prime}(\theta(\mu))\)
표준연결함수와 분산함수 1
표준연결함수: \(\theta \equiv \eta\) 로 하는 연결함수 \(\ell=\left(\gamma^{\prime}\right)^{-1}\)
\(X\) 를 모형행렬이라 하면 \(\boldsymbol{\eta}=X \boldsymbol{\beta}\) 이고, \(\varphi\) 가 알려진 값일 때, \(A=\operatorname{diag} A_i\) 이다.
표준연결에 대해 \(X^T A \boldsymbol{y}\) 가 \(\boldsymbol{\beta}\) 에 대한 최소충분통계량이다.
점수함수를 이용하면 회귀 모수 \(\boldsymbol{\beta}\) 에 대한 점수 방정식은 \[
X^T A \boldsymbol{y}=X^T A \widehat{\boldsymbol{\mu}}
\]
표준연결함수와 분산함수 2
binomial\(\log (\mu /(1-\mu))\) logit\(\mu(1-\mu)\) mu(1-mu)
Gamma\(-1 / \mu\) inverse\(\mu^2\) mu^2
gaussian\(\mu\) identity1
constant
inverse.gaussian\(-2 / \mu^2\) 1 / mu^2\(\mu^3\) mu^3
poisson\(\log \mu\) log\(\mu\) mu
추정 절차
GLM의 최대 가능도 추정량의 명시적 표현식은 일반적으로 존재하지 않음. 반복법을 통해 계산.
초기 추정값 \(\widehat{\eta}_0\)
가가중치(working weight) \(W_k = \frac{A}{V}\left.\left(\frac{\mathrm{d} \mu}{\mathrm{d} \eta}\right)^2\right\vert_{\eta=\hat{\eta}_k}\)
가반응변수(working value) \(z_k = \eta+\frac{y-\mu}{\mathrm{d} \mu / \mathrm{d} \eta}\vert_{\eta=\hat{\eta}_k}\)
\(k\) 번 째 반복에서 \(\mathbf{\beta}\) 의 추정치의 새로운 근사값은 가중치 \(W_k\) 를 이용해 \(X\) 의 \(z_k\) 의 가중 선형회귀로 계산할 수 있다.
뉴턴 랩슨 방법과 비슷하나 헤시안 행렬이 그 기댓값인 피셔 정보행렬로 대체된다는 점이 다르다 (Fisher scoring).
Iterative (re-)weighted least squares (IWLS, IRLS)
\(\widehat{\boldsymbol{\beta}}\) 의 공분산행렬의 대표본 추정량은 \(\widehat{\varphi}\left(X^T \widehat{W} X\right)^{-1}\) 이다. 이는 IRLS의 부산물로 얻어진다.
\(\widehat{\boldsymbol{\beta}}\) 의 Fisher scoring에 의한 반복추정은 척도 모수 \(\varphi\) 에 의존하지 않는다. (왜 그런가?)
반복 과정은 평균함수와 분산함수를 통해서만 반응변수의 분포에 의존한다.
이탈도 분석
포화모형 (\(S\) ): 모수 \(\theta_i\) 또는 선형 예측변수 \(\eta_i\) 가 자유 모수인 모형
점수방정식으로부터 최대가능도 추정량은 \(\hat{\theta}_i = (\gamma^\prime)^{-1}(y_i)\) , \(i=1,\dotsc,n\) .
모형 \(M\) : \(p\) (< \(n\) ) 개의 회귀모수를 갖는 GLM
이탈도(deviance)
우선 척도모수 \(\varphi\) 가 1이라 가정
\(\widehat{\theta}_i\) : 모형 \(M\) 하에서 \(\theta_i\) 의 최대가능도 추정량이탈도: 포화모형 \(S\) 내에서 \(M\) 을 검정하기 위한 로그 가능비(의 2배) \[
D_M=2 \underset{i=1} { \stackrel{n}{\sum} } A_i\left[\left\{y_i \theta\left(y_i\right)-\gamma\left(\theta\left(y_i\right)\right)\right\}-\left\{y_i \widehat{\theta}_i-\gamma\left(\widehat{\theta_i}\right)\right\}\right]
\]
척도화 이탈도(잔차 이탈도): \(D_M / \varphi\)
예시: 정규분포족 (\(\varphi = \sigma^2\) )
\(D_M\) 은 잔차제곱합. 따라서 \(D_M / \varphi \sim \chi_{n-p}^2\) \(\varphi\) 의 불편추정량 \[
\widehat{\varphi}=\frac{D_M}{n-p}
\] 척도모수의 대안추정량: 표준화잔차제곱합 \[
\widetilde{\varphi}=\frac{1}{n-p} \sum_i \frac{\left(y_i-\widehat{\mu}_i\right)^2}{V\left(\widehat{\mu}_i\right) / A_i}
\] 정규분포족의 경우 \(\widetilde{\varphi} = \widehat{\varphi}\) .
정규분포족이 아닌 경우에도 척도화 이탈도의 카이제곱 근사가 유효할 수 있으며, 이 경우 \(\widehat{\varphi}\) 를 척도모수의 근사적 불편추정량으로 사용할 수 있다.
검정
\(M_0 \subset M\) 가 \(q\) 개 (\(q<p\) )의 회귀 모수들을 갖는 축소모형일 때 \(\varphi\) 가 알려진 값이면, 일반적인 가능도비 논증에 의해 \[
\frac{D_{M_0}-D_M}{\varphi} {\dot{\sim}} \chi_{p-q}^2
\] 로 \(M_0\) 를 \(M\) 내에서 검정할 수 있다.
\(\varphi\) 를 모를 경우, 관례적으로 정규분포족과 유사하게 다음의 근사 결과를 이용한다. \[
\frac{(D_{M_0}-D_M)}{\widehat{\varphi}(p-q)} {\dot{\sim}} F_{p-q, n-p}
\]
일반화 선형모형을 위한 R 함수들
glm()
glm (formula, family, data, weights, control)
family: Table 7.1의 Family Name
binomial(link=probit)와 같이 입력값을 주어야 할 수도 있다. (quasi족의 분산함수가 이렇게 쓰일 수 있다.)사용자 정의 family (MASS::negative.binomial 등)도 사용 가능
weight: 사전 가중치 \(A_i\)
control: IRLS 알고리즘 파라미터. maxit = 최대 반복 횟수. 기본값은 25; trace = 반복 절차 인쇄 여부. 기본값은 FALSE; epsilon = 반복 정지 규칙. 기본값은 1e-8. \[
|\text { deviance }^{(i)}-\text { deviance }^{(i-1)}| <\epsilon\left(\text { deviance }^{(i-1)} + 0.1 \right)
\]
glm 객체에 대한 method functions: coef, resid, print, summary, deviance
<- function (obj){<- summary (obj, corr = F)$ dispersion * so$ cov.unscaled
anova 검정: test="Chisq" = 식 (7.9)를 이용한 카이제곱 검정. test="F": 식 (7.10)을 이용한 F 검정
dispersion: 척도 모수 \(\varphi\) . summary, anova를 이용한 방법안에서와 예측을 위한 표준오차를 위해서만 쓰인다. argument dispersion을 통해 주어질 수 있다. 이항분포와 포아송족에서는 기본값 1, 다른 경우 식 (7.8).
예측 및 잔차
glm의 메소드 함수인 predict의 인자 type의 기본값은 type="link"로, 선형 예측변수 \(\eta\) 의 예측을 제공한다. 반면 type="response"는 평균 \(\mu\) 에 대한 예측을 나타낸다 (예: \(\widehat{\mu}_i\) ).
glm의 4종류 잔차 (resid 메소드의 type 인자)
response: \(y_i-\widehat{\mu}_i\) pearson: \((y_i-\widehat{\mu}_i)/\sqrt{{V}(\widehat{\mu}_i}\) working: \((y_i-\widehat{\mu}_i)/ (\mathrm{d}\mu_i/\mathrm{d}\eta_i)\) . IRLS의 마지막 단계에서 얻어짐.deviance \[
\small
\operatorname{sign}(y_i-\widehat{\mu}_i) \sqrt{A_i\left[\left\{y_i \theta\left(y_i\right)-\gamma\left(\theta\left(y_i\right)\right)\right\}-\left\{y_i \widehat{\theta}_i-\gamma\left(\widehat{\theta_i}\right)\right\}\right]}
\]
오프셋
offset()의 인수를 평가하여 선형 예측변수에 더해줌.
모형공식에서 I()로 인수를 감싸주고 계수를 1로 강제하는 것과 동등
예제) MASS::anorexia 자료: 거식증 환자의 세 가지 치료 방법에 대한 임상 시험의 치료 전후 체중이 포함된 자료
.1 <- glm (Postwt ~ Prewt + Treat + offset (Prewt), family = gaussian, data = anorexia) offset을 사용하지 않고 모형공식의 좌변에 Postwt - Prewt를 사용하여 모형을 적합하면 이를 이용한 예측은 실제 치료 후 체중이 아닌 체중 증가를 예측
offset을 이용하면 적합된 모형으로 관심 있는 반응변수에 대한 직접적인 예측이 가능
이항 자료
살충제 자료
살충제 자료 (Collett, 1991)
Sex
1
2
4
8
16
32
Male
1
4
9
13
18
20
Female
0
2
6
10
12
16
로지스틱 회귀 모형을 glm의 binomial족을 사용하여 적합하기 위해서는 시행 횟수 \(a_i\) 를 다음의 3가지 방법 중 하나로 지정해야 한다.
반응 변수가 숫자 벡터이면 데이터를 비율 형식(\(y_i = s_i/a_i\) )으로 유지하는 것으로 간주된다. 이 경우 \(a_i\) 는 인수 weights를 사용하여 가중치 벡터로 주어져야 한다. (\(a_i\) 가 모두 \(1\) 이면 기본 가중치로 충분하다.)
반응 변수가 논리형 벡터이거나 2진 인자이면 \(0/1\) 숫자 벡터로 간주되어 처리된다. 반응 변수가 다수준 인자이면 첫 번째 수준은 0 (실패)으로 처리되고 다른 모든 수준은 1 (성공)로 처리된다.
반응 변수가 2열 행렬이면 첫 번째 열에는 각 시행에 대한 성공 횟수 \(s_i\) 가 있고, 두 번째 열에는 실패 횟수 \(a_i-s_i\) 가 있다고 가정된다. 이 경우 인수 weights가 필요하지 않다. 모든 경우에 반응변수가 상대도수 \(y_i = s_i/a_i\) 이므로 평균 \(\mu_i\) 는 성공확률 \(p_i\) 이다. 따라서 적합 fitted는 이항 평균이 아닌 확률이다.
<- rep (0 : 5 , 2 )<- c (1 , 4 , 9 , 13 , 18 , 20 , 0 , 2 , 6 , 10 , 12 , 16 ) <- factor (rep (c ("M" , "F" ), c (6 , 6 )))<- cbind (numdead, numalive = 20 - numdead)<- glm (SF ~ sex* ldose, family = binomial) # option 3 summary (budworm.lg, cor = F)
Call:
glm(formula = SF ~ sex * ldose, family = binomial)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.9935 0.5527 -5.416 6.09e-08 ***
sexM 0.1750 0.7783 0.225 0.822
ldose 0.9060 0.1671 5.422 5.89e-08 ***
sexM:ldose 0.3529 0.2700 1.307 0.191
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 124.8756 on 11 degrees of freedom
Residual deviance: 4.9937 on 8 degrees of freedom
AIC: 43.104
Number of Fisher Scoring iterations: 4
분석 결과는 성별 간의 기울기 차이에 대한 약간의 증거를 보여준다.
성별의 주효과는 유의하지 않으나, sex:ldose 교호항의 주변범주이므로 후자를 포함한다면 모형에 포함되어야 한다.
예측
<- seq (0 , 5 , 0.1 )<- data.frame (ldose = rep (ld, 2 ), sex = factor (c (rep ("M" , length (ld)), rep ("F" , length (ld))), levels = levels (sex)))$ prob <- predict (budworm.lg, df_pred, type= "response" )ggplot (df_pred, aes (x = 2 ^ ldose, y = prob)) + geom_line (aes (colour = sex)) + coord_trans (x = "log2" ) + labs (x = "dose" , y = "prob" )
이탈도 분석
잔차 이탈도 4.9937은 \(\chi_8^2\) 확률변수로는 작은 값이고, 자료의 시행횟수가 적당히 크므로 모형은 매우 잘 적합되었다.
이탈도 분석을 통해 2차항이 필요한지를 확인하자.
anova (update (budworm.lg, . ~ . + sex * I (ldose^ 2 )), test = "Chisq" )
Analysis of Deviance Table
Model: binomial, link: logit
Response: SF
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 11 124.876
sex 1 6.077 10 118.799 0.0137 *
ldose 1 112.042 9 6.757 <2e-16 ***
I(ldose^2) 1 0.907 8 5.851 0.3410
sex:ldose 1 1.240 7 4.611 0.2655
sex:I(ldose^2) 1 1.439 6 3.172 0.2303
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2차항의 영향은 유의하지 않다.
LD50
\(\xi_p\) : 반응 확률이 \(p\) 인 로그 용량
\(2^{\xi_{0.5}}\) : LD50, 50% lethal dose
\(\xi_p=\frac{\ell(p)-\beta_0}{\beta_1}, \quad \frac{\partial \xi_p}{\partial \beta_0}=-\frac{1}{\beta_1}, \quad \frac{\partial \xi_p}{\partial \beta_1}=-\frac{\ell(p)-\beta_0}{\beta_1^2}=-\frac{\xi_p}{\beta_1}\)
\(\beta_0\) =절편, \(\beta_1\) =기울기 sex:ldose 교호 작용이 크게 유의하지 않으므로, (로짓 척도에서) 공통 기울기를 갖고 성별에 따라 절편이 다른 모형을 다시 적합한다.
<- glm (SF ~ sex + ldose - 1 , family = binomial) # budworm.lg0 <- glm(formula = SF ~ sex + ldose, family = binomial) summary (budworm.lg0, cor = F)$ coefficients
Estimate Std. Error z value Pr(>|z|)
sexF -3.473155 0.4685202 -7.413032 1.234445e-13
sexM -2.372412 0.3855108 -6.153945 7.557914e-10
ldose 1.064214 0.1310775 8.118971 4.701542e-16
사분위수 LD값
dose.p (budworm.lg0, cf = c (1 ,3 ), p = c (1 : 3 )/ 4 ) # female
Dose SE
p = 0.25: 2.231265 0.2499089
p = 0.50: 3.263587 0.2297539
p = 0.75: 4.295910 0.2746874
dose.p (budworm.lg0, cf = c (2 ,3 ), p = c (1 : 3 )/ 4 ) # male
Dose SE
p = 0.25: 1.196939 0.2635100
p = 0.50: 2.229262 0.2259649
p = 0.75: 3.261585 0.2549838
Probit link
\[
y = \Phi^{-1}(p) = \beta^Tx
\]
신생아 저체중 자료
# A tibble: 189 × 10
low age lwt race smoke ptl ht ui ftv bwt
<int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
1 0 19 182 2 0 0 0 1 0 2523
2 0 33 155 3 0 0 0 0 3 2551
3 0 20 105 1 1 0 0 0 1 2557
4 0 21 108 1 1 0 0 1 2 2594
5 0 18 107 1 1 0 0 1 0 2600
6 0 21 124 3 0 0 0 0 0 2622
7 0 22 118 1 0 0 0 0 1 2637
8 0 17 103 3 0 0 0 0 1 2637
9 0 29 123 1 1 0 0 0 1 2663
10 0 26 113 1 1 0 0 0 0 2665
# ℹ 179 more rows
Hosmer and Lemeshow (1989). 미국 병원에서 189명의 출생에 대한 데이터 세트를 제공. 주요 관심사: 저체중아
변수
low = 출생시 체중이 2.5kg 미만 (0/1)age = 산모 나이(세)lwt = 마지막 월경 기간에 산모의 체중(파운드)race: 1 == white/2 == black/3 == othersmoke = 임신 중 흡연상태 (0/1)ptl = 이전 조기 진통 횟수ht = 고혈압 병력(0/1)ui = 자궁 과민증 유무(0/1)ftv = 첫 3개월 동안 의사 방문 횟수bwt = 실제 출생 체중(그램)
<- birthwt %>% mutate (race = factor (race, labels = c ("white" , "black" , "other" )), ptd = factor (ptl> 0 ), ftv = fct_collapse (factor (ftv), ` 2+ ` = c ('2' , '3' , '4' , '6' )),low = factor (low),smoke = factor (smoke > 0 ),ht = factor (ht > 0 ),ui = factor (ui > 0 )%>% select (- c (bwt, ptl)) <- glm (low ~ ., family = binomial, data = bwt)summary (birthwt.glm, cor = F) # 상관행렬 생략
Call:
glm(formula = low ~ ., family = binomial, data = bwt)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.82302 1.24471 0.661 0.50848
age -0.03723 0.03870 -0.962 0.33602
lwt -0.01565 0.00708 -2.211 0.02705 *
raceblack 1.19241 0.53597 2.225 0.02609 *
raceother 0.74069 0.46174 1.604 0.10869
smokeTRUE 0.75553 0.42502 1.778 0.07546 .
htTRUE 1.91317 0.72074 2.654 0.00794 **
uiTRUE 0.68019 0.46434 1.465 0.14296
ftv1 -0.43638 0.47939 -0.910 0.36268
ftv2+ 0.17901 0.45638 0.392 0.69488
ptdTRUE 1.34376 0.48062 2.796 0.00518 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 195.48 on 178 degrees of freedom
AIC: 217.48
Number of Fisher Scoring iterations: 4
반응 변수가 이진 (binary)이기 때문에 모형이 정확하더라도 이탈도가 대략적으로 카이제곱 분포를 가질 것이라는 보장은 없다. 그러나 값이 자유도와 거의 일치하므로 적합에 대해 심각한 문제는 없어 보인다.
모형 선택
주효과 모형
<- MASS:: stepAIC (birthwt.glm, trace = F) $ anova
Stepwise Model Path
Analysis of Deviance Table
Initial Model:
low ~ age + lwt + race + smoke + ht + ui + ftv + ptd
Final Model:
low ~ lwt + race + smoke + ht + ui + ptd
Step Df Deviance Resid. Df Resid. Dev AIC
1 178 195.4755 217.4755
2 - ftv 2 1.358185 180 196.8337 214.8337
3 - age 1 1.017866 181 197.8516 213.8516
2차 교호작용 추가
<- MASS:: stepAIC (birthwt.glm, ~ .^ 2 + I (scale (age)^ 2 ) + I (scale (lwt)^ 2 ), trace = F)$ anova
Stepwise Model Path
Analysis of Deviance Table
Initial Model:
low ~ age + lwt + race + smoke + ht + ui + ftv + ptd
Final Model:
low ~ age + lwt + smoke + ht + ui + ftv + ptd + age:ftv + smoke:ui
Step Df Deviance Resid. Df Resid. Dev AIC
1 178 195.4755 217.4755
2 + age:ftv 2 12.474896 176 183.0006 209.0006
3 + smoke:ui 1 3.056805 175 179.9438 207.9438
4 - race 2 3.129586 177 183.0734 207.0734
summary (birthwt.step2)$ coef
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.58237445 1.421613247 -0.4096574 0.6820572535
age 0.07553886 0.053967335 1.3997144 0.1615988535
lwt -0.02037257 0.007497148 -2.7173759 0.0065801819
smokeTRUE 0.78004357 0.420385078 1.8555453 0.0635183923
htTRUE 2.06569557 0.748742526 2.7588864 0.0057998694
uiTRUE 1.81853023 0.667555094 2.7241650 0.0064464291
ftv1 2.92108800 2.285774273 1.2779425 0.2012697142
ftv2+ 9.24490669 2.661497465 3.4735734 0.0005135770
ptdTRUE 1.56031692 0.497001294 3.1394625 0.0016925809
age:ftv1 -0.16182450 0.096819391 -1.6714059 0.0946415389
age:ftv2+ -0.41103343 0.119144147 -3.4498835 0.0005608286
smokeTRUE:uiTRUE -1.91667524 0.973097049 -1.9696650 0.0488767714
table (bwt$ low, predict (birthwt.step2) > 0 )
FALSE TRUE
0 116 14
1 28 31
1단계에서 age와 ftv가 제거되었지만 2단계에서는 교호작용이 포함되어 있으며 세 ftv 그룹 내에서 age에 대한 기울기가 상당히 다르다.
smoke:ui 항을 추가하면 AIC가 감소하지만 해당 \(t\) 통계량은 유의수준 5%에서만 유의하다.
이항분포족 GLM에서의 주의사항
glm에서 표준오차와 \(t\) 값은 최대가능도 추정량 \(\widehat{\boldsymbol{\beta}}\) 에서 로그 가능도를 2차 테일러 전개하여 얻어진 Wald 근사에서 유도된다.
몇몇 \(\widehat{\beta}_i\) 값이 크면 \(\widehat{\boldsymbol{\beta}}\) 에서 로그 가능도의 곡률은 \(\beta_i=0\) 근처보다 훨씬 작을 수 있으므로 Wald 근사는 \(\beta_i=0\) 근방에서의 로그 가능도의 변화를 과소평가한다.
즉 \(|\widehat{\beta}_i| \rightarrow \infty\) 일 때 \(t\) 통계량이 \(0\) 에 접근하는 경향이 있고, 가능도비 검정에서 매우 유의한 계수들이 유의하지 않은 \(t\) 값 비율을 가질 수 있다 (Hauck and Donner, 1977 ).
적합된 확률이 0 또는 1에 매우 가까운 경우, Hauck-Donner 현상과 함께 수렴 문제도 발생할 수 있다.
수천 건의 사례와 약 5개의 이진 설명변수가 있는 의료 진단 문제를 고려하자. 이러한 지표 중 하나는 거의 1이 아니지만 1이면 항상 질병이 있음을 나타낸다고 하자. 그러면 해당 지표가 있는 사례의 적합된 확률은 1이 되어야 하며 이는 \(\widehat{\beta}_i = \infty\) 을 취해야만 달성될 수 있다.
glm의 결과는 경고가 되며 추정 계수는 약 \(\pm 10\) 이다. 사실 이 경우 최대 가능도 추정치가 존재하지 않는다 (Santner and Duffy, 1989 , p. 234).한 가지 해법은 작은 능형(ridge) 벌점함수를 추가하는 것이다 (Bach, 2010 ).
포아송과 다항분포 모형
포아송족
포아송족의 표준 연결함수는 log.
주요 용도: 실제로 다항분포에서 추출한 도수분포표에 포아송 로그 선형 모형을 대리로 적합
포아송 로그 선형 모형은 종종 비율(rates)에 직접 적용됨.
다원 도수분포표를 분류하는 인자를 반응인자와 설명인자로 나누는 것이 편리함.
설명인자의 주변합은 미리 (또는 추론 목적으로) 고정됨.
주요 관심사: 설명인자가 주어졌을 때 반응 인자의 조건부 확률에 있다.
알려진 사실: 합이 주어진 독립 포아송 확률 변수들의 조건부 분포 == 전체에 대한 포아송 평균의 비율로 주어진 확률을 갖는 다항분포
설명인자의 수준별 조합 내에서 다원 반응변수의 계수(counts)에 적용.
따라서 다항확률의 승법(multiplicative) 모형을 포아송 로그 선형 모형을 사용하여 적합하고 검정할 수 있음.
대리 포아송 모형에 대응되는 다항 모형
A, B, … 를 도수분포표를 분류하는 인자라 하자.최소 모형 은 모든 설명인자의 교호작용이며, 반응 인자의 총합을 고정하기 위해 분석에 포함되어야 함.반응인자와 설명인자 간의 상호 작용을 표현하기 위해 그래프 구조를 사용하는 것이 종종 도움이 된다 (Lauritzen, 1996).
코펜하겐 주거 자료
Satisfaction Low
Med.
High
Low
Med.
High
Housing Influence
Tower blocks
Low
21
21
28
14
19
37
Medium
34
22
36
17
23
40
High
10
11
36
3
5
23
Apartments
Low
61
23
17
78
46
43
Medium
43
35
40
48
45
86
High
26
18
54
15
25
62
Atrium houses
Low
13
9
10
20
23
20
Medium
8
8
12
10
22
24
High
6
7
9
7
10
21
Terraced houses
Low
18
6
7
57
23
13
Medium
15
13
13
31
21
13
High
7
5
11
5
6
13
코펜하겐에 거주하는 1681명의 가구주를 대상으로 자신이 거주하고 있는 임대주택 유형, 다른 입주자들과의 접촉 (만남 및 연락) 정도, 아파트 관리의 영향력에 대한 생각, 주거 수준에 대한 만족도 등을 조사한 4원 분류
# A tibble: 72 × 5
Sat Infl Type Cont Freq
<ord> <fct> <fct> <fct> <int>
1 Low Low Tower Low 21
2 Medium Low Tower Low 21
3 High Low Tower Low 28
4 Low Medium Tower Low 34
5 Medium Medium Tower Low 22
6 High Medium Tower Low 36
7 Low High Tower Low 10
8 Medium High Tower Low 11
9 High High Tower Low 36
10 Low Low Apartment Low 61
# ℹ 62 more rows
Cox and Snell (1981, pp. 155ff)을 따라 Sat을 반응인자로 하고 Type, Infl, Cont를 설명인자로 하는 분석 실시
모형 선택
다음 단계: Sat:Infl 항을 모형에 통합하고 재평가
Cox and Snell (1981)에 따르면 세 가지 항이 모두 필요하므로 이들을 모두 포함하도록 초기 모형을 업데이트
<- update (house.glm0, . ~ . + Sat: (Infl+ Type+ Cont)) summary (house.glm1, cor = F)
Call:
glm(formula = Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont +
Type:Cont + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont,
family = poisson, data = MASS::housing)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.135074 0.120112 26.101 < 2e-16 ***
InflMedium 0.248327 0.159979 1.552 0.120602
InflHigh -0.412645 0.184947 -2.231 0.025671 *
TypeApartment 0.292524 0.157477 1.858 0.063231 .
TypeAtrium -0.792847 0.214413 -3.698 0.000218 ***
TypeTerrace -1.018074 0.221263 -4.601 4.20e-06 ***
ContHigh -0.001407 0.169711 -0.008 0.993385
Sat.L -0.098106 0.112592 -0.871 0.383570
Sat.Q 0.285657 0.122283 2.336 0.019489 *
InflMedium:TypeApartment -0.017882 0.210496 -0.085 0.932302
InflHigh:TypeApartment 0.386869 0.233297 1.658 0.097263 .
InflMedium:TypeAtrium -0.360311 0.304979 -1.181 0.237432
InflHigh:TypeAtrium -0.036788 0.334793 -0.110 0.912503
InflMedium:TypeTerrace 0.185154 0.288892 0.641 0.521580
InflHigh:TypeTerrace 0.310749 0.334815 0.928 0.353345
InflMedium:ContHigh -0.200060 0.228748 -0.875 0.381799
InflHigh:ContHigh -0.725790 0.282352 -2.571 0.010155 *
TypeApartment:ContHigh 0.569691 0.212152 2.685 0.007247 **
TypeAtrium:ContHigh 0.702115 0.276056 2.543 0.010979 *
TypeTerrace:ContHigh 1.215930 0.269968 4.504 6.67e-06 ***
InflMedium:Sat.L 0.519627 0.096830 5.366 8.03e-08 ***
InflHigh:Sat.L 1.140302 0.118180 9.649 < 2e-16 ***
InflMedium:Sat.Q -0.064474 0.102666 -0.628 0.530004
InflHigh:Sat.Q 0.115436 0.127798 0.903 0.366380
TypeApartment:Sat.L -0.520170 0.109793 -4.738 2.16e-06 ***
TypeAtrium:Sat.L -0.288484 0.149551 -1.929 0.053730 .
TypeTerrace:Sat.L -0.998666 0.141527 -7.056 1.71e-12 ***
TypeApartment:Sat.Q 0.055418 0.118515 0.468 0.640068
TypeAtrium:Sat.Q -0.273820 0.149713 -1.829 0.067405 .
TypeTerrace:Sat.Q -0.032328 0.149251 -0.217 0.828520
ContHigh:Sat.L 0.340703 0.087778 3.881 0.000104 ***
ContHigh:Sat.Q -0.097929 0.094068 -1.041 0.297851
InflMedium:TypeApartment:ContHigh 0.046900 0.286212 0.164 0.869837
InflHigh:TypeApartment:ContHigh 0.126229 0.338208 0.373 0.708979
InflMedium:TypeAtrium:ContHigh 0.157239 0.390719 0.402 0.687364
InflHigh:TypeAtrium:ContHigh 0.478611 0.444244 1.077 0.281320
InflMedium:TypeTerrace:ContHigh -0.500162 0.367135 -1.362 0.173091
InflHigh:TypeTerrace:ContHigh -0.463099 0.454713 -1.018 0.308467
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 833.657 on 71 degrees of freedom
Residual deviance: 38.662 on 34 degrees of freedom
AIC: 455.63
Number of Fisher Scoring iterations: 4
선택된 주효과 모형(house.glm1)은 세 설명인자 Type, Infl, Cont가 만족도 수준의 확률에 영향을 미친다는 것을 보여준다.
모형 평가
모형에서 추정된 확률
<- housing %>% select (- Freq) %>% map (levels)<- predict (house.glm1, type= "response" ) # poisson means <- matrix (house.pm, ncol = 3 , byrow = T, dimnames = list (NULL , hnames[[1 ]]))<- house.pm/ drop (house.pm %*% rep (1 , 3 )) cbind (expand.grid (hnames[names (hnames) != "Sat" ]), prob = round (house.pr, 2 ))
Infl Type Cont prob.Low prob.Medium prob.High
1 Low Tower Low 0.40 0.26 0.34
2 Medium Tower Low 0.26 0.27 0.47
3 High Tower Low 0.15 0.19 0.66
4 Low Apartment Low 0.54 0.23 0.23
5 Medium Apartment Low 0.39 0.26 0.34
6 High Apartment Low 0.26 0.21 0.53
7 Low Atrium Low 0.43 0.32 0.25
8 Medium Atrium Low 0.30 0.35 0.36
9 High Atrium Low 0.19 0.27 0.54
10 Low Terrace Low 0.65 0.22 0.14
11 Medium Terrace Low 0.51 0.27 0.22
12 High Terrace Low 0.37 0.24 0.39
13 Low Tower High 0.30 0.28 0.42
14 Medium Tower High 0.18 0.27 0.54
15 High Tower High 0.10 0.19 0.71
16 Low Apartment High 0.44 0.27 0.30
17 Medium Apartment High 0.30 0.28 0.42
18 High Apartment High 0.18 0.21 0.61
19 Low Atrium High 0.33 0.36 0.31
20 Medium Atrium High 0.22 0.36 0.42
21 High Atrium High 0.13 0.27 0.60
22 Low Terrace High 0.55 0.27 0.19
23 Medium Terrace High 0.40 0.31 0.29
24 High Terrace High 0.27 0.26 0.47
확률에 가장 큰 영향을 미치는 인자는 Infl으로, 증가할수록 Sat이 낮을 확률은 감소하고 높을 확률은 증가
다음으로 중요한 인자는 Type
다른 주민과의 접촉(Cont)이 많아질수록 만족도가 낮을 확률은 낮아지고 높을 확률은 높아지는 경향이 있으나 그 효과는 상대적으로 미미하다.
비례 오즈 모형
반응변수 Sat이 순서형이므로 순서형 로지스틱 모형인 비례 오즈 모형이 자연스러움 (McCullagh, 1980 ; McCullagh and Nelder, 1989, §5.2.2)
낮은 만족도와 중간 만족도의 누적 확률의 오즈비는 설명인자에 의존하지 않음 \[
\begin{split}
\operatorname{logit}\Pr(Y \leqslant k \mid \boldsymbol{x}) &= \log\frac{\Pr(Y\le k \mid \mathbf{x})}{\Pr(Y > k \mid \mathbf{x})} =\zeta_k - \eta
,
\\
\zeta_0 &= -\infty<\zeta_1<\cdots<\zeta_K=\infty
\end{split}
\] 여기서 \(\eta\) 는 선형 예측변수.
MASS 함수 polr를 사용한 적합
<- MASS:: polr (Sat ~ Infl + Type + Cont, data = MASS:: housing, weights = Freq))
Call:
MASS::polr(formula = Sat ~ Infl + Type + Cont, data = MASS::housing,
weights = Freq)
Coefficients:
InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace
0.5663937 1.2888191 -0.5723501 -0.3661866 -1.0910149
ContHigh
0.3602841
Intercepts:
Low|Medium Medium|High
-0.4961353 0.6907083
Residual Deviance: 3479.149
AIC: 3495.149
잔차 이탈도는 이전에 적합된 다항 모형의 이탈도(3470)와 유사하며, 모수가 6개 감소한 반면 이탈도 증가량은 9.0에 불과
AIC는 크게 감소
절편의 차이인 1.1는 위에서 계산한 평균 로짓 차이인 1.22와 잘 일치
적합된 확률 계산: 다항분포족 모형과 비슷
<- predict (house.plr, expand.grid (hnames[names (hnames) != "Sat" ]), type = "probs" )cbind (expand.grid (hnames[names (hnames) != "Sat" ]), round (house.pr1, 2 ))
Infl Type Cont Low Medium High
1 Low Tower Low 0.38 0.29 0.33
2 Medium Tower Low 0.26 0.27 0.47
3 High Tower Low 0.14 0.21 0.65
4 Low Apartment Low 0.52 0.26 0.22
5 Medium Apartment Low 0.38 0.29 0.33
6 High Apartment Low 0.23 0.26 0.51
7 Low Atrium Low 0.47 0.27 0.26
8 Medium Atrium Low 0.33 0.29 0.38
9 High Atrium Low 0.19 0.25 0.56
10 Low Terrace Low 0.64 0.21 0.14
11 Medium Terrace Low 0.51 0.26 0.23
12 High Terrace Low 0.33 0.29 0.38
13 Low Tower High 0.30 0.28 0.42
14 Medium Tower High 0.19 0.25 0.56
15 High Tower High 0.10 0.17 0.72
16 Low Apartment High 0.43 0.28 0.29
17 Medium Apartment High 0.30 0.28 0.42
18 High Apartment High 0.17 0.23 0.60
19 Low Atrium High 0.38 0.29 0.33
20 Medium Atrium High 0.26 0.27 0.47
21 High Atrium High 0.14 0.21 0.64
22 Low Terrace High 0.56 0.25 0.19
23 Medium Terrace High 0.42 0.28 0.30
24 High Terrace High 0.26 0.27 0.47
음이항분포족
음이항분포와 GLM
분산이 평균보다 큰 계수 자료 분석에서 유용
관측가능한 이산변수 \(Y\) 에 대해 은닉변수 \(E\) 가 있어 \[
Y \mid E \sim \operatorname{Poi}(\mu E),
\quad
E \sim \operatorname{gamma}(\theta, 1)/\theta
\quad
\text{(shape, scale)}
\] 라 가정하면 \(Y \sim NB(\theta, \theta/(\mu + \theta))\) : \[
\begin{split}
& \operatorname{E}(Y)=\mu, \quad \operatorname{var}(Y)=\mu+\mu^2 / \theta,
\\
& f_Y(y ; \theta, \mu)= \frac{\Gamma(\theta+y)}{\Gamma(\theta) y !}
\left(\frac{\mu}{\mu + \theta}\right)^{y}\left(\frac{\theta}{\mu + \theta}\right)^{\theta}
\end{split}
\]
\(\theta\) 가 알려진 경우 이 분포는 일반적인 GLM 형식을 갖는다 (\(\theta\) 는 산포모수 \(\varphi\) 와는 다름).
Quine 자료 (6장)
계수 자료이므로 포아송 GLM을 우선 고려
glm (Days ~ .^ 4 , family = poisson, data = MASS:: quine)
Call: glm(formula = Days ~ .^4, family = poisson, data = MASS::quine)
Coefficients:
(Intercept) EthN SexM
3.0564 -0.1386 -0.4914
AgeF1 AgeF2 AgeF3
-0.6227 -2.3632 -0.3784
LrnSL EthN:SexM EthN:AgeF1
-1.9577 -0.7524 0.1029
EthN:AgeF2 EthN:AgeF3 EthN:LrnSL
-0.5546 0.0633 2.2588
SexM:AgeF1 SexM:AgeF2 SexM:AgeF3
0.4092 3.1098 1.1145
SexM:LrnSL AgeF1:LrnSL AgeF2:LrnSL
1.5900 2.6421 4.8585
AgeF3:LrnSL EthN:SexM:AgeF1 EthN:SexM:AgeF2
NA -0.3105 0.3469
EthN:SexM:AgeF3 EthN:SexM:LrnSL EthN:AgeF1:LrnSL
0.8329 -0.1639 -3.5493
EthN:AgeF2:LrnSL EthN:AgeF3:LrnSL SexM:AgeF1:LrnSL
-3.3315 NA -2.4285
SexM:AgeF2:LrnSL SexM:AgeF3:LrnSL EthN:SexM:AgeF1:LrnSL
-4.1914 NA 2.1711
EthN:SexM:AgeF2:LrnSL EthN:SexM:AgeF3:LrnSL
2.1029 NA
Degrees of Freedom: 145 Total (i.e. Null); 118 Residual
Null Deviance: 2074
Residual Deviance: 1174 AIC: 1818
평균-분산 관계 (6장)
\(\theta \approx 2\) 인 음이항 모형이 적절해 보인다.
우선 \(\theta = 2\) 로 가정. base::glm에서는 negative.binomial로 적합
<- glm (Days ~ .^ 4 , family = MASS:: negative.binomial (2 ), data = MASS:: quine) 6장의 최종 모형(에 기반한 선형 예측변수)과 비교
<- update (quine.nb, . ~ Sex/ (Age + Eth* Lrn)) anova (quine.nb0, quine.nb, test = "Chisq" )
Analysis of Deviance Table
Model 1: Days ~ Sex + Sex:Age + Sex:Eth + Sex:Lrn + Sex:Eth:Lrn
Model 2: Days ~ (Eth + Sex + Age + Lrn)^4
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 132 198.51
2 118 171.98 14 26.527 0.03433 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
이 모형(Model 1)이 고정된 \(\theta\) 가정 하에서 지나치게 단순해 보인다.
\(\theta\) 가 추정되는 경우
\(\theta\) 의 추정치와 표준오차
c (theta = quine.nb2$ theta, SE = quine.nb2$ SE)
theta SE
1.8653435 0.2580065
잔차 분석
이탈도 잔차를 이용한 진단 확인
<- data.frame (resid = resid (quine.nb2, type= "deviance" ), pred = predict (quine.nb2))<- ggplot (rs, aes (pred, resid)) + geom_point () + geom_hline (yintercept= 0 , col= "red" , linetype= "dashed" ) + xlab ("Linear predictors" ) + ylab ("Deviance residuals" )<- ggplot (rs, aes (sample = resid)) + geom_qq () + geom_qq_line () + xlab ("Theoretical Quantiles" ) + ylab ("Deviance residuals" )+ p2
특별히 후보 모형에 반대되는 증거가 보이지 않음
과대산포
산포모수 \(\varphi\)
정규분포족: \(\varphi\) 에 대한 적률 추정량은 최대가능도 추정량의 일반적인 불편(unbiased) 조정
이항분포 및 포아송족: \(\varphi=1\) (모형이 맞다면)
이항 GLM의 과대산포
모수 \(p\) 가 실제로는 평균 \(p\) 를 갖는 확률변수 \(\theta\) 인 경우
평균 \(p\) 를 정하는 선형 로지스틱 모형에 랜덤효과가 있는 경우 발생
\(\operatorname{var} \theta=\phi p(1-p)\) 라 가정하면 \[
\operatorname{E}Y=n p, \quad \operatorname{var} Y=n p(1-p)[1+(n-1) \phi]
\] 예: \(\theta \sim beta ( \alpha, \beta )\) . \(p=\mathrm{E} \theta= \alpha /(\alpha+\beta)\) 이고 \(\operatorname{var} \theta=p(1-p) /(\alpha+\beta+1)\) .
\(n_i = n\) 인 경우 \[
\operatorname{var} Y=n p(1-p)[1+(n-1) \phi]=\varphi n p(1-p)
\] \(n \equiv 1\) 에는 적용할 수 없음\(\varphi\) 에 상한이 생김: \(\phi \leqslant 1 \implies \varphi \leqslant n\) .조정된 GLM 적합은 회귀모수에 대해서도 최대 가능도 추정이 아님.
군집모형 (McCullagh and Nelder, 1989, pp. 125–6)
\(n\) 개의 자료점이 \(k\) 개의 군집에서 표집된 것으로 가정군집 내에 독립적으로 \(n/k\) 개를 이항 표집
각 군집은 랜덤 성공확률 \(\pi_i\) , \(\operatorname{E}[\pi_i] = p\) , \(\operatorname{Var}[\pi_i] = \phi\pi(1-\pi)\) 를 가짐 \[
\operatorname{E}Y=n p, \quad \operatorname{var} Y=n p(1-p)[1+(k-1) \phi]
\]
\(\varphi\) 로 척도화된 이탈도와 잔차 이탈도의 변화가 적절한 자유도를 갖는 점근 카이제곱 분포를 가짐을 시사\(\varphi\) 가 추정되어야 함 — 이탈도 분석, addterm 및 dropterm 분석에서 카이제곱 검정 대신 F 검정 사용
교환가능(exchangeable) 모형
이항 관측값을 구성하는 \(n\) 번의 시행이 교환 가능하지만 반드시 독립적이지는 않다고 가정.
임의의 시행 쌍은 상관계수 \(\delta\) 를 가짐. \[
\operatorname{var} Y=n p(1-p)[1+(n-1) \delta]=\varphi n p(1-p)
\]
\(\delta \geq 0\) 이라는 제약은 없으나 \(\delta \geqslant -1/(n-1)\) 의 제약이 있음.
관측값의 평균과 분산만 지정된다는 점에서 사실상 준가능도(quasilikelihood) 모형
R에서 이항분포 및 포아송 분포족은 \(\varphi\) 를 추정하는 것을 허용하지 않음
\(\varphi\) 가 항상 추정되는 quasibinomial과 quasipoisson족이 추가로 제공.가능도가 없고 따라서 AIC도 없으므로 모형 선택에 stepAIC를 사용할 수 없음
