예시

• 정규분포. 평균 \(\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)라고 한다.

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}} \]

추정 절차

• 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} \]

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

  vcov.glm <- function(obj){
      so <- summary(obj, corr = F)
      so$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\)).
  • se.fit=TRUE: 표준오차 반환

• 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]} \]
    • 이탈도의 signed summand.

• 가우시안족에 대해서는 4종의 잔차가 동일

• 이탈도와 Pearson 잔차는 보통 매우 비슷

• 지렛대값의 개념과 적합에 대한 효과는 선형회귀에서와 같이 중요 (Davison and Snell, 1991)

  • GLM의 용처에서는 단순 회귀 분석이나 설계된 실험에서와 같이 높은 레버리지가 발생할 수 없는 경우가 많음

오프셋

• offset()의 인수를 평가하여 선형 예측변수에 더해줌.

  • 모형공식에서 I()로 인수를 감싸주고 계수를 1로 강제하는 것과 동등

• 예제) MASS::anorexia 자료: 거식증 환자의 세 가지 치료 방법에 대한 임상 시험의 치료 전후 체중이 포함된 자료

  ax.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt), 
              family = gaussian, data = anorexia)

• offset을 사용하지 않고 모형공식의 좌변에 Postwt - Prewt를 사용하여 모형을 적합하면 이를 이용한 예측은 실제 치료 후 체중이 아닌 체중 증가를 예측

• offset을 이용하면 적합된 모형으로 관심 있는 반응변수에 대한 직접적인 예측이 가능

살충제 자료

• 담배나방(Heliothis virescens) 유충이 살충제인 pyrenthroid trans-cypermethrin에 저항을 보이기 시작하는 용량 실험. 측정치는 각 20 마리의 실험 배치에서 죽거나 쓰러진 숫자.

• 용량은 2의 거듭제곱이므로 log2(dose)를 이용해 로지스틱 회귀 모형을 적합

• 분석 결과는 성별 간의 기울기 차이에 대한 약간의 증거를 보여준다.

• 성별의 주효과는 유의하지 않으나, sex:ldose 교호항의 주변범주이므로 후자를 포함한다면 모형에 포함되어야 한다.

• 예측

  ld <- seq(0, 5, 0.1)
  df_pred <- data.frame(ldose = rep(ld, 2), 
                        sex = factor(c(rep("M", length(ld)), 
                                       rep("F",  length(ld))), 
                                     levels = levels(sex)))
  df_pred$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 교호 작용이 크게 유의하지 않으므로, (로짓 척도에서) 공통 기울기를 갖고 성별에 따라 절편이 다른 모형을 다시 적합한다.

  budworm.lg0 <- 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

• MASS 라이브러리에는 \(\widehat{\xi}_p\), \(\widehat{\xi}_p\)의 점근 표준오차를 계산하는 함수가 제공

  MASS::dose.p

  function (obj, cf = 1:2, p = 0.5) 
  {
      eta <- family(obj)$linkfun(p)
      b <- coef(obj)[cf]
      x.p <- (eta - b[1L])/b[2L]
      names(x.p) <- paste("p = ", format(p), ":", sep = "")
      pd <- -cbind(1, x.p)/b[2L]
      SE <- sqrt(((pd %*% vcov(obj)[cf, cf]) * pd) %*% c(1, 1))
      res <- structure(x.p, SE = SE, p = p)
      class(res) <- "glm.dose"
      res
  }
  <bytecode: 0x1207bf990>
  <environment: namespace:MASS>

  • family 함수를 사용하여 연결함수를 추출한 다음 선형 예측변수의 값을 얻는 데 사용하는 방법에 주목

• 사분위수 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 \]

• 생물시료 분석에서는 프로빗 연결이 전통적으로 사용되어 왔다. 추정된 확률들이 꼬리에 집중되어 있지 않으면 프로빗과 로짓 연결은 유사한 결과가 나오는 경향이 있다.

  dose.p(update(budworm.lg0, family = binomial(link = probit)), 
         cf = c(1, 3), p = 1:3/4)  # female

                Dose        SE
  p = 0.25: 2.191229 0.2384478
  p = 0.50: 3.257703 0.2240685
  p = 0.75: 4.324177 0.2668745

신생아 저체중 자료

    as.tibble(MASS::birthwt)

# 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 == other
  • smoke = 임신 중 흡연상태 (0/1)
  • ptl = 이전 조기 진통 횟수
  • ht = 고혈압 병력(0/1)
  • ui = 자궁 과민증 유무(0/1)
  • ftv = 첫 3개월 동안 의사 방문 횟수
  • bwt = 실제 출생 체중(그램)

• 실제 출생 체중을 사용할 수 있으나 여기서는 나머지 변수를 통해 출생 체중이 낮은지 예측하는 데 중점을 둔다.

• 이항(0/1) 반응으로 로지스틱 회귀를 사용

• 변수를 사용하는 방법을 신중하게 고려

  • ptl에 선형적인 반응을 기대하는 것은 무리. 1보다 큰 값을 가진 표본은 너무 작기 때문에 유무만 확인.
  • ftv도 세 가지 수준으로 줄일 수 있음
  • 비가우시안 GLM의 경우 처치 대비를 사용하는 것이 일반적

• 반응 변수가 이진 (binary)이기 때문에 모형이 정확하더라도 이탈도가 대략적으로 카이제곱 분포를 가질 것이라는 보장은 없다. 그러나 값이 자유도와 거의 일치하므로 적합에 대해 심각한 문제는 없어 보인다.

모형 선택

• 주효과 모형

  birthwt.step <- MASS::stepAIC(birthwt.glm, trace = F) 
  birthwt.step$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차 교호작용 추가

  birthwt.step2 <- MASS::stepAIC(birthwt.glm, ~ .^2 +  
                                   I(scale(age)^2) + I(scale(lwt)^2), 
                                 trace = F)
  birthwt.step2$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

• 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).

• 코펜하겐에 거주하는 1681명의 가구주를 대상으로 자신이 거주하고 있는 임대주택 유형, 다른 입주자들과의 접촉 (만남 및 연락) 정도, 아파트 관리의 영향력에 대한 생각, 주거 수준에 대한 만족도 등을 조사한 4원 분류

as_tibble(MASS::housing)

# 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를 설명인자로 하는 분석 실시

모형 선택

• 초기 모형: 모든 Type \(\times\) Cont \(\times\) Infl 조합에 대해 Sat의 조건부 분포가 동일한 모형. 즉, 만족도는 다른 설명 인자와 무관.

  house.glm0  <-  glm(Freq ~ Infl*Type*Cont +  Sat, family = poisson, 
                      data = MASS::housing)
  summary(house.glm0,  cor = F)

  
  Call:
  glm(formula = Freq ~ Infl * Type * Cont + Sat, family = poisson, 
      data = MASS::housing)
  
  Coefficients:
                                      Estimate Std. Error z value Pr(>|z|)    
  (Intercept)                        3.136e+00  1.196e-01  26.225  < 2e-16 ***
  InflMedium                         2.733e-01  1.586e-01   1.723 0.084868 .  
  InflHigh                          -2.054e-01  1.784e-01  -1.152 0.249511    
  TypeApartment                      3.666e-01  1.555e-01   2.357 0.018403 *  
  TypeAtrium                        -7.828e-01  2.134e-01  -3.668 0.000244 ***
  TypeTerrace                       -8.145e-01  2.157e-01  -3.775 0.000160 ***
  ContHigh                           1.402e-15  1.690e-01   0.000 1.000000    
  Sat.L                              1.159e-01  4.038e-02   2.871 0.004094 ** 
  Sat.Q                              2.629e-01  4.515e-02   5.824 5.76e-09 ***
  InflMedium:TypeApartment          -1.177e-01  2.086e-01  -0.564 0.572571    
  InflHigh:TypeApartment             1.753e-01  2.279e-01   0.769 0.441783    
  InflMedium:TypeAtrium             -4.068e-01  3.035e-01  -1.340 0.180118    
  InflHigh:TypeAtrium               -1.692e-01  3.294e-01  -0.514 0.607433    
  InflMedium:TypeTerrace             6.292e-03  2.860e-01   0.022 0.982450    
  InflHigh:TypeTerrace              -9.305e-02  3.280e-01  -0.284 0.776633    
  InflMedium:ContHigh               -1.398e-01  2.279e-01  -0.613 0.539715    
  InflHigh:ContHigh                 -6.091e-01  2.800e-01  -2.176 0.029585 *  
  TypeApartment:ContHigh             5.029e-01  2.109e-01   2.385 0.017083 *  
  TypeAtrium:ContHigh                6.774e-01  2.751e-01   2.462 0.013811 *  
  TypeTerrace:ContHigh               1.099e+00  2.675e-01   4.106 4.02e-05 ***
  InflMedium:TypeApartment:ContHigh  5.359e-02  2.862e-01   0.187 0.851450    
  InflHigh:TypeApartment:ContHigh    1.462e-01  3.380e-01   0.432 0.665390    
  InflMedium:TypeAtrium:ContHigh     1.555e-01  3.907e-01   0.398 0.690597    
  InflHigh:TypeAtrium:ContHigh       4.782e-01  4.441e-01   1.077 0.281619    
  InflMedium:TypeTerrace:ContHigh   -4.980e-01  3.671e-01  -1.357 0.174827    
  InflHigh:TypeTerrace:ContHigh     -4.470e-01  4.545e-01  -0.984 0.325326    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for poisson family taken to be 1)
  
      Null deviance: 833.66  on 71  degrees of freedom
  Residual deviance: 217.46  on 46  degrees of freedom
  AIC: 610.43
  
  Number of Fisher Scoring iterations: 5

  • 높은 잔차 이탈도는 이 단순 모형이 부적절하다는 것을 분명히 나타내므로 확률은 설명 인자에 따라 달라지는 것으로 보임.

• 이러한 종류의 변형을 허용하는 가장 간단한 항을 모형에 추가

  MASS::addterm(house.glm0,  ~. + Sat:(Infl+Type+Cont), test = "Chisq")

  Single term additions
  
  Model:
  Freq ~ Infl * Type * Cont + Sat
           Df Deviance    AIC     LRT   Pr(Chi)    
  <none>        217.46 610.43                      
  Infl:Sat  4   111.08 512.05 106.371 < 2.2e-16 ***
  Type:Sat  6   156.79 561.76  60.669 3.292e-11 ***
  Cont:Sat  2   212.33 609.30   5.126   0.07708 .  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  • Infl 인자는 단일 항 중 가장 크게 AIC를 감소시킨다.

• 다음 단계: Sat:Infl 항을 모형에 통합하고 재평가
  • Cox and Snell (1981)에 따르면 세 가지 항이 모두 필요하므로 이들을 모두 포함하도록 초기 모형을 업데이트

  house.glm1 <- 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

• 이탈도는 모형을 잘 적합한 것을 나타내지만 가능한 상호 교호항들을 추가하는 것을 고려하자.

  MASS::addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test = "Chisq")

  Single term additions
  
  Model:
  Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont + 
      Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont
                Df Deviance    AIC     LRT Pr(Chi)  
  <none>             38.662 455.63                  
  Infl:Type:Sat 12   16.107 457.08 22.5550 0.03175 *
  Infl:Cont:Sat  4   37.472 462.44  1.1901 0.87973  
  Type:Cont:Sat  6   28.256 457.23 10.4064 0.10855  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  • 첫 번째 항인 Type \(\times\) Infl 교호작용은 약간 유의한 것으로 보이지만 AIC가 증가하므로 단순성을 이유로 포함하지 않기로 한다.

• 선택된 주효과 모형(house.glm1)은 세 설명인자 Type, Infl, Cont가 만족도 수준의 확률에 영향을 미친다는 것을 보여준다.

모형 평가

• 모형에서 추정된 확률

  hnames <- housing %>% select(-Freq) %>% map(levels)
  house.pm <- predict(house.glm1, type="response") # poisson means
  house.pm <- matrix(house.pm, ncol = 3, byrow = T, 
                     dimnames = list(NULL,  hnames[[1]]))
  house.pr <- 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)이 많아질수록 만족도가 낮을 확률은 낮아지고 높을 확률은 높아지는 경향이 있으나 그 효과는 상대적으로 미미하다.

다항분포족 모형 적합

• 라이브러리 nnet의 함수 multinom을 사용하여 대리 포아송 모델을 사용하는 대신 다항 모형을 직접 적합할 수 있다. 모형에 교호작용항이 필요하지 않다.

  house.mult <- nnet::multinom(Sat ~ Infl + Type + Cont, 
                               weights = Freq, data = housing);

  # weights:  24 (14 variable)
  initial  value 1846.767257 
  iter  10 value 1747.045232
  final  value 1735.041933 
  converged

  house.mult

  Call:
  nnet::multinom(formula = Sat ~ Infl + Type + Cont, data = housing, 
      weights = Freq)
  
  Coefficients:
         (Intercept) InflMedium  InflHigh TypeApartment TypeAtrium TypeTerrace
  Medium  -0.4192316  0.4464003 0.6649367    -0.4356851  0.1313663  -0.6665728
  High    -0.1387453  0.7348626 1.6126294    -0.7356261 -0.4079808  -1.4123333
          ContHigh
  Medium 0.3608513
  High   0.4818236
  
  Residual Deviance: 3470.084 
  AIC: 3498.084 

• 이탈도는 최소 모형의 각 셀에 대한 다항 반응이 아니라 다음의 포화 모형에 대한 것이다.

  house.mult2 <- nnet::multinom(Sat ~ Infl*Type*Cont, 
                                weights = Freq, data = housing);

  # weights:  75 (48 variable)
  initial  value 1846.767257 
  iter  10 value 1734.465581
  iter  20 value 1717.220153
  iter  30 value 1715.760679
  iter  40 value 1715.713306
  final  value 1715.710836 
  converged

  anova(house.mult,  house.mult2,  test = "none")

  Likelihood ratio tests of Multinomial Models
  
  Response: Sat
                 Model Resid. df Resid. Dev   Test    Df LR stat.
  1 Infl + Type + Cont       130   3470.084                      
  2 Infl * Type * Cont        96   3431.422 1 vs 2    34 38.66219

비례 오즈 모형

• 반응변수 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\)는 선형 예측변수.

• 비례 오즈 가정 확인: 잘 적합된 다항분포족 모형을 사용, 누적 확률의 로짓 차이 조사

  house.cpr <- house.pr %>% apply(1, cumsum)
  logit <-  function(x) log(x/(1-x))
  house.ld <- logit(house.cpr[2, ]) -  logit(house.cpr[1, ]) 
  mean(sort(house.ld))

  [1] 1.223835

  • 평균 로그 오즈비는 약 \(1.2\)이고 그 변동은 크지 않으므로 비례 오즈 모형은 유용한 단순화가 될 수 있다.

MASS 함수 polr를 사용한 적합

(house.plr <- 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와 잘 일치

• 적합된 확률 계산: 다항분포족 모형과 비슷

  house.pr1 <- 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

• 다항분포족 모형 내에서 비례 오즈 모형을 검정하기 위한 가능도비 통계량을 직접 계산

  Fr <- matrix(housing$Freq, ncol = 3,  byrow = T) 
  2 * sum(Fr * log(house.pr/house.pr1))

  [1] 9.065433

  • 위의 이탈도 차이와 일치. 자유도 6에 대해 유의하지 않음 (p=0.17).

• 비례 오즈 모형은 다항분포족 모형보다 간결할 뿐 아니라 모수의 수가 적으므로 확률에 대한 공변량의 작용을 해석하고 설명하기가 쉽다.

• stepAIC는 더 복잡한 모형 선형 예측변수를 선택한다 (왜 그런가?)

  house.plr2 <- MASS::stepAIC(house.plr, ~.^2);

  house.plr2$anova

  Stepwise Model Path 
  Analysis of Deviance Table
  
  Initial Model:
  Sat ~ Infl + Type + Cont
  
  Final Model:
  Sat ~ Infl + Type + Cont + Infl:Type + Type:Cont
  
           Step Df  Deviance Resid. Df Resid. Dev      AIC
  1                               1673   3479.149 3495.149
  2 + Infl:Type  6 22.509347      1667   3456.640 3484.640
  3 + Type:Cont  3  7.945029      1664   3448.695 3482.695

음이항분포와 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로 적합

  quine.nb <- glm(Days ~ .^4, family = MASS::negative.binomial(2), 
                  data = MASS::quine)

• 6장의 최종 모형(에 기반한 선형 예측변수)과 비교

  quine.nb0 <- 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\)가 추정되는 경우

• MASS에서는 \(\theta\)의 최대가능도 추정을 포함하여 glm()을 수정한 glm.nb()를 제공

  quine.nb <- MASS::glm.nb(Days ~ .^4, data = MASS::quine)
  quine.nb2 <- MASS::stepAIC(quine.nb) 

  quine.nb2$anova 

  Stepwise Model Path 
  Analysis of Deviance Table
  
  Initial Model:
  Days ~ (Eth + Sex + Age + Lrn)^4
  
  Final Model:
  Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + 
      Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + 
      Sex:Age:Lrn
  
                 Step Df   Deviance Resid. Df Resid. Dev      AIC
  1                                       118   167.4535 1095.324
  2 - Eth:Sex:Age:Lrn  2 0.09746244       120   167.5509 1092.728
  3     - Eth:Sex:Age  3 0.11060087       123   167.4403 1089.409

• 선택된 모형은 Lrn/(Age + Eth + Sex)^2로 쓸 수 있음

• AIC는 과적합되는 경향이 있으므로 3차 교호작용항을 검정

  MASS::dropterm(quine.nb2, test = "Chisq") 

  Single term deletions
  
  Model:
  Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + 
      Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + 
      Sex:Age:Lrn
              Df    AIC    LRT Pr(Chi)  
  <none>         1089.4                 
  Eth:Sex:Lrn  1 1092.5 5.0384 0.02479 *
  Eth:Age:Lrn  2 1091.2 5.7684 0.05590 .
  Sex:Age:Lrn  2 1091.1 5.7019 0.05779 .
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Sex:Age:Lrn, Eth:Sex:Lrn 교호작용 제거 여부 검토

  quine.nb3 <- update(quine.nb2, . ~ . - Eth:Age:Lrn - Sex:Age:Lrn)
  anova(quine.nb2, quine.nb3)

  Likelihood ratio tests of Negative Binomial Models
  
  Response: Days
                                                                                                                          Model
  1                             Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn
  2 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Eth:Lrn + Sex:Age + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + Sex:Age:Lrn
       theta Resid. df    2 x log-lik.   Test    df LR stat.    Pr(Chi)
  1 1.724987       127       -1053.431                                 
  2 1.865343       123       -1043.409 1 vs 2     4 10.02203 0.04005825

  • 두 항 중 하나만 제거하고 둘 다는 제거하지 않는 것을 시사

• \(\theta\)의 추정치와 표준오차

  c(theta = quine.nb2$theta, SE  = quine.nb2$SE)

      theta        SE 
  1.8653435 0.2580065 

잔차 분석

• 이탈도 잔차를 이용한 진단 확인

  rs <- data.frame(resid = resid(quine.nb2, type="deviance"), 
                   pred = predict(quine.nb2))
  p1 <- ggplot(rs, aes(pred, resid)) + geom_point() + 
          geom_hline(yintercept=0, col="red", linetype="dashed") + 
          xlab("Linear predictors") + ylab("Deviance residuals")
  p2 <- ggplot(rs, aes(sample = resid)) + geom_qq() + geom_qq_line() + 
          xlab("Theoretical Quantiles") + ylab("Deviance residuals")
  p1 + 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를 사용할 수 없음