Parameters in the Link Function — Box-Cox 멱 링크와 Pregibon 선형화

1 서론 — 왜 링크 모수를 추정하는가

표준 GLM 은 링크 함수 \(g(\mu) = \eta\) 를 분포에서 정해지는 고정 함수 로 본다 — 포아송에는 \(\log\), 이항에는 logit, 감마에는 \(1/\mu\). 이론적으로는 canonical link 가 우아하지만, 실무에서는 세 가지 의문이 남는다.

1. 정말 그 링크가 맞는가? 감마 오차에 \(1/\mu\) 와 \(\log \mu\) 중 어느 쪽이 데이터에 더 적합한지는 판단 근거가 필요하다
2. 두 링크 사이에 연속적 선택지가 있는가? logit 과 cloglog 사이 어딘가에 더 적합한 링크가 있을 수 있다
3. 링크 선택의 불확실성은 얼마인가? “로그 링크” 라고 단정했을 때, 데이터가 다른 링크를 얼마나 배제하는가?

McCullagh & Nelder (1989, §11.3) 의 답은 링크 함수를 모수 가족(parametric family) 으로 확장하고, 그 모수를 추정하는 것이다.

1.1 직관: 링크 가족의 ‘GPS 좌표’

\(\mu\) 에서 \(\eta\) 로 가는 “변환” 은 여러 후보가 있다. 로그, 역수, 제곱근, 항등, … 이들을 각각 독립된 선택지로 보는 대신, 한 축 위의 연속적 점 으로 배치하면 어떨까.

\[ g(\mu; \lambda) = \mu^\lambda = \begin{cases} 1/\mu & \lambda = -1 \\ \log\mu & \lambda = 0 \\ \sqrt{\mu} & \lambda = 0.5 \\ \mu & \lambda = 1 \end{cases} \]

\(\lambda\) 는 링크 공간의 “GPS 좌표”. 데이터가 \(\lambda\) 를 골라 주면, 그 좌표가 곧 최적 링크다. 95% 신뢰구간이 여러 표준 링크를 포함하면 — 그들이 모두 데이터와 양립 가능하다는 진단도 자동으로 나온다.

2 멱 가족 — Box-Cox 의 GLM 버전

2.1 정의와 연속성 보정

\[ g(\mu; \lambda) = \begin{cases} \mu^\lambda & \lambda \neq 0 \\ \log\mu & \lambda = 0 \end{cases} \]

\(\lambda \to 0\) 극한에서 \(\mu^\lambda \to 1\) 이고 \(1\) 이 아닌 \(\log\mu\) 로 가야 의미가 있다. 이 불연속을 제거한 Box-Cox 변환

\[ g(\mu; \lambda) = \frac{\mu^\lambda - 1}{\lambda} \]

은 L’Hôpital 규칙으로 \(\lambda \to 0\) 에서 \(\log\mu\) 에 수렴한다. GLM 적합에서는 절편이 자동으로 \(-1/\lambda\) 를 흡수하므로 \(\mu^\lambda\) 형과 수치적으로 동등하다.

2.2 표준 링크와의 대응

\(\lambda\)	\(g(\mu)\)	표준 이름	전형적 용도
\(-1\)	\(1/\mu\)	역수	감마 canonical
\(-0.5\)	\(1/\sqrt{\mu}\)	—	분산 안정화 (중간)
\(0\)	\(\log\mu\)	로그	포아송 canonical, 곱셈 모형
\(0.5\)	\(\sqrt{\mu}\)	제곱근	포아송 분산 안정화
\(1\)	\(\mu\)	항등	정규 선형모형

2.3 왜 이 가족이 특히 유용한가

• 대부분의 실무 링크가 포함됨 — 데이터가 \(\lambda\) 를 선택하면 어느 표준 링크에 가까운지 즉시 알 수 있다
• \(\lambda\) 에 대해 미분 가능 — Pregibon 선형화가 작동하는 전제
• 해석 가능성 — \(\lambda\) 가 “곱셈성(\(\log\)) 과 덧셈성(\(\mu\)) 사이의 어디” 인지를 말해 준다

3 Pregibon (1980) 의 선형화 — 미지 \(\lambda\) 를 보조 공변량으로

3.1 문제 설정

GLM 적합식 \(g(\mu; \lambda) = \mathbf x^T \boldsymbol\beta\) 에서 \(\lambda\) 가 미지이면 \(\boldsymbol\beta\) 의 IRLS 를 여러 \(\lambda\) 격자에서 반복해야 한다. 이를 한 번의 IRLS 로 \(\lambda\) 까지 함께 갱신 하는 묘수가 Pregibon 의 기여다.

3.2 Taylor 전개

초기값 \(\lambda_0\) 근방에서 일차 전개:

\[ g(\mu; \lambda) \simeq g(\mu; \lambda_0) + (\lambda - \lambda_0) \left.\frac{\partial g}{\partial \lambda}\right|_{\lambda_0} \]

멱 가족에서 \(\partial \mu^\lambda / \partial \lambda = \mu^\lambda \log \mu\) 이므로

\[ \mu^\lambda \simeq \mu^{\lambda_0} + (\lambda - \lambda_0) \mu^{\lambda_0} \log \mu \tag{11.1} \]

3.3 GLM 적합식의 재배치

\(g(\mu; \lambda) = \sum_j \beta_j x_j\) 를 (11.1) 로 대입:

\[ \mu^{\lambda_0} + (\lambda - \lambda_0) \mu^{\lambda_0} \log \mu \simeq \sum_j \beta_j x_j \]

\((\lambda - \lambda_0) \mu^{\lambda_0} \log \mu\) 를 우변으로 옮기고 \(\mu\) 를 현재 추정치 \(\hat\mu_0\) 로 고정:

\[ \mu^{\lambda_0} \simeq \sum_j \beta_j x_j - (\lambda - \lambda_0) \hat\mu_0^{\lambda_0} \log \hat\mu_0 \]

이것은 원래 공변량 \(\{x_j\}\) 에 추가 공변량 \(v = -\hat\mu_0^{\lambda_0} \log \hat\mu_0\) 를 더한 \(\lambda_0\)-링크 GLM 과 같다. 새 계수 \(\gamma\) 는 정의상 \(\gamma = (\lambda - \lambda_0)\).

3.4 알고리즘

초기값: lambda_0 설정 (예: 1)
반복:
  (a) lambda_0-링크 GLM 적합 → hat_beta, hat_mu_0
  (b) v_i = -hat_mu_0^lambda_0 * log(hat_mu_0) 계산
  (c) (X, v) 확장 설계행렬로 lambda_0-링크 GLM 재적합 → gamma_hat
  (d) lambda_new = lambda_0 + gamma_hat
  (e) |lambda_new - lambda_0| < tol 이면 수렴
  (f) lambda_0 ← lambda_new

3.5 직관: 왜 ’공변량 하나 추가’로 \(\lambda\) 가 추정되는가

핵심은 “\(\lambda_0\)-링크에서 남은 불일치의 방향” 을 설계행렬에 추가하는 것이다.

\(\lambda_0 = 1\) (항등 링크)로 적합한 모형이 데이터와 안 맞는다면, 그 불일치 패턴은 \(\mu \log \mu\) 형태로 나타난다 — 이는 \(\lambda\) 를 1 에서 조금 움직였을 때 링크가 바뀌는 방향 벡터. 이 방향을 공변량으로 추가하면 새 계수가 곧 “\(\lambda\) 가 1 에서 얼마나 벗어나는가” 를 측정한다.

이는 score test 의 기하학과 동일하다:

• 귀무가설 \(H_0: \lambda = \lambda_0\)
• Score 방향: \(\partial \ell / \partial \lambda |_{\lambda_0}\)
• 대안에 대한 검정력은 이 방향으로의 잔차 적합으로 얻어진다

Pregibon 의 보조 공변량 \(v\) 는 정확히 이 score 방향이다.

4 Goodness-of-Link 검정

4.1 Score 검정으로서의 해석

\(H_0: \lambda = \lambda_0\) 가 \(H_1: \lambda \neq \lambda_0\) 에 대해 검정 가능한가.

이탈도 차이:

\[ T = D(\lambda_0) - D(\hat\lambda) \;\overset{H_0}{\sim}\; \chi^2_1 \]

\(D(\lambda_0)\) 는 \(\lambda_0\)-링크로 적합한 이탈도, \(D(\hat\lambda)\) 는 \(\lambda\) 까지 최적화한 이탈도. 자유도 1 은 \(\lambda\) 가 추가된 1차원.

4.2 더 빠른 버전 — Rao score test

Pregibon 선형화 결과에서 보조 공변량 \(v\) 의 계수 \(\hat\gamma\) 에 대한 Wald 통계량

\[ Z = \hat\gamma / \text{SE}(\hat\gamma) \;\overset{H_0}{\sim}\; N(0, 1) \]

이는 IRLS 한 번 추가 반복 으로 얻을 수 있어 여러 \(\lambda_0\) 후보를 빠르게 비교할 수 있다.

4.3 자동차 보험 예제 (Fig 11.1)

McCullagh & Nelder 의 자동차 보험 청구 데이터. 평균 청구금액을 주요 효과만으로 설명, 분산함수 \(V(\mu) \propto \mu^2\) (감마).

\(\lambda\) 를 \([-2.5, 1]\) 에서 움직이며 이탈도를 계산. 결과:

• 최솟값 \(D(\hat\lambda) = 124.51\) 이 \(\hat\lambda \approx -1\) (역수 링크) 근처
• 95% 신뢰구간: \[ \{\lambda : D(\lambda) - 124.51 < \tfrac{124.51}{108} \times 3.93\} \]
• 이 구간은 \(\lambda = 0\) (로그 링크)을 포함 — 로그 링크도 기각되지 않는다

4.4 실무적 해석

• \(\hat\lambda = -1\) 이 유일한 정답이 아님
• 역수 vs 로그는 데이터만으로는 구분 불가능한 영역
• 해석력·계산 편의·다른 유사 데이터와의 일관성 등 비-통계적 기준 으로 최종 선택

이것이 표본 크기로 링크를 식별하는 것의 현실적 한계다. 설명 계수는 쉽게 추정해도, 링크 함수의 형태는 훨씬 큰 표본이 필요하다.

5 수렴성의 함정과 대응

5.1 1차 근사의 한계

(11.1) 의 선형 근사는 \(\lambda\) 가 \(\lambda_0\) 와 가까울 때만 정확하다. \(|\lambda - \lambda_0| > 0.5\) 정도 되면 Taylor 전개의 잉여항이 커져 한 번의 갱신이 과녁을 지나칠 수 있다.

Pregibon 자신의 경고:

“이 방법은 합리적 적합을 개선하는 데 유용하지, 절망적 상황을 고칠 목적으로는 기대하지 말라.”

5.2 격자 탐색과의 결합

실무 전략:

1. 격자 탐색: \(\lambda \in \{-2, -1.5, \ldots, 1.5, 2\}\) 의 11 점에서 GLM 을 적합 (각 점에서 IRLS)
2. 대략 최솟값 위치 확인 (예: \(\lambda \approx -1\))
3. 그 근방 \(\lambda_0\) 에서 Pregibon 반복 시작 — 1차 근사가 타당한 영역

5.3 감쇠 업데이트

업데이트가 진동하면:

\[ \lambda_{\text{new}} = \alpha \lambda_0 + (1 - \alpha) (\lambda_0 + \hat\gamma), \quad \alpha \in (0, 1) \]

\(\alpha = 0.7\) 정도가 안전. 수렴은 느려지지만 발산은 막는다.

5.4 내부 반복

보조 공변량 \(v\) 의 값이 \(\hat\mu_0\) 에 의존하므로, 한 \(\lambda_0\) 에서 IRLS 를 수렴시킨 뒤 \(v\) 를 다시 계산해 추가 반복을 돌리면 더 안정적.

6 두 개 이상의 링크 모수

6.1 Shifted power family

\[ g(\mu; \alpha, \lambda) = \frac{(\mu + \alpha)^\lambda - 1}{\lambda} \]

\((\alpha, \lambda) = (0, 0)\) 이면 로그. \((\alpha, \lambda) = (0, 1)\) 이면 항등 링크 (\(\mu + 0 - 1/\lambda = \mu\), 상수 흡수).

\(\alpha\) 의 역할: \(\mu\) 를 변환하기 전에 원점을 이동. \(\mu \approx 0\) 근방에서 \(\log\) 나 멱 링크가 폭발하는 것을 방지 — 비음수 응답의 영 값 처리.

두 개의 보조 공변량:

\[ v_1 = \left.\frac{\partial g}{\partial \alpha}\right|_{\alpha_0, \lambda_0} = (\mu + \alpha_0)^{\lambda_0 - 1} \]

\[ v_2 = \left.\frac{\partial g}{\partial \lambda}\right|_{\alpha_0, \lambda_0} = \frac{(\mu + \alpha_0)^{\lambda_0} \log(\mu + \alpha_0) - [(\mu + \alpha_0)^{\lambda_0} - 1]/\lambda_0}{\lambda_0} \]

설계행렬에 두 열을 추가하여 GLM 적합. 두 계수가 \((\alpha, \lambda)\) 의 업데이트를 각각 제공.

6.2 이항 데이터를 위한 링크 가족

Tolerance distribution 기반 probit 분석의 일반화:

\[ g(\mu; \lambda, \delta) = \frac{\pi^{\lambda - \delta} - 1}{\lambda - \delta} - \frac{(1 - \pi)^{\lambda + \delta} - 1}{\lambda + \delta} \]

\(\pi = \mu/m\) 은 반응 비율. \((\lambda, \delta) \to (0, 0)\) 극한에서 logit 링크.

해석:

• \(\lambda\) — 링크의 “전체 스케일” (로그 멱)
• \(\delta\) — \(\pi\) 와 \(1 - \pi\) 의 비대칭성 정도

\(\delta \neq 0\) 이면 성공과 실패를 대칭적으로 다루지 않는 링크. probit, logit, cloglog 모두 포함하는 한 가족.

6.3 단일 모수 이항 링크 가족

\[ g(\mu; \lambda) = \log\!\left[\frac{(1/(1-\pi))^\lambda - 1}{\lambda}\right] \]

• \(\lambda = 1\): 로지스틱 (logit)
• \(\lambda \to 0\): complementary log-log (cloglog)

실무 사용: 생존 분석에서 “지수 붕괴 가정 (cloglog) 과 로지스틱 중 어느 쪽이 맞는가” 를 데이터로 검정. \(\hat\lambda\) 의 95% 신뢰구간이 1 만 포함하면 로지스틱, 0 만 포함하면 cloglog, 둘 다 포함하면 구분 불가.

6.4 두 모수의 식별성 문제

\((\alpha, \lambda)\) 나 \((\lambda, \delta)\) 를 동시에 추정할 때 흔한 함정:

• 고도 상관: 한 모수의 변화가 다른 모수로 거의 완전히 보상됨
• 평탄한 likelihood: 이탈도가 여러 조합에서 거의 동일
• 발산: 한 방향이 다른 방향을 먹어버려 초기값 의존성 폭발

대응: 하나씩 고정해 profile 로 확인 → 두 방향 모두 유의하면 동시 추정, 아니면 하나만 추정.

7 데이터 변환 vs 적합값 변환 — GLM 의 철학적 승리

7.1 Box & Cox (1964) 의 원래 제안

\(y\) 자체를 \(y^\lambda\) (또는 \(\log y\))로 변환한 뒤 정규-항등 선형모형 을 적용.

\[ y_i^\lambda \sim \mathcal{N}(\mathbf x_i^T \boldsymbol\beta, \sigma^2) \]

한 \(\lambda\) 로 두 목표를 동시에 달성하려는 전략:

• 체계 부분의 가법성 (선형모형 타당성)
• 무작위 부분의 등분산성·정규성

7.2 두 목표가 다른 \(\lambda\) 를 요구하는 경우

Fisher (1949) 의 결핵 데이터. Nelder & Wedderburn (1972) 재분석:

• \(\sqrt{y}\) 이면 분산이 안정 (등분산성 만족)
• \(\log y\) 이면 가법성 만족 (요인 효과 선형)

\(\lambda\) 하나로 두 목표를 모두 충족 불가. Box-Cox 는 어느 한쪽을 희생하거나 절충해야 한다.

7.3 GLM 의 분리 전략

\[ \eta_i = g(\mu_i) = \mathbf x_i^T \boldsymbol\beta, \qquad \mathrm{var}(Y_i) = \phi V(\mu_i) \]

• 가법성 → 링크 \(g\) 가 담당
• 분산 구조 → 분산함수 \(V\) 가 담당

두 역할이 수학적으로 독립적. 데이터가 가법성을 위해 로그 링크를 요구해도, 분산 구조는 별도로 \(V(\mu) \propto \mu\) (포아송적) 이나 \(V(\mu) \propto \mu^2\) (감마적) 으로 선택 가능.

7.4 Fisher 결핵 데이터의 재분석

• \(Y\) (원자료) + \(V \propto \mu\) + 로그 링크: 이탈도 \(= D_1\)
• \(\sqrt{Y}\) + \(V = \text{const}\) + 로그 링크: 이탈도 \(= D_2\)

Baker & Nelder (1978) 보고: \(D_1 \approx D_2\) — 두 분석이 거의 동등. GLM 에서 분산함수를 “\(\mu\) 비례” 로 선택하면 \(\sqrt{Y}\) 변환과 같은 효과를 얻되, 원자료 스케일에서 계수를 해석 할 수 있다는 이점.

7.5 직관: 왜 이 분리가 중요한가

\(Y\) 를 변환하면:

• 원래 단위의 해석을 잃는다 (\(\log Y\) 의 계수는 \(Y\) 에 대한 비율 해석이 가능하지만, \(\sqrt{Y}\) 의 계수는 해석 난이도 급상승)
• 예측값을 \(Y\) 스케일로 되돌릴 때 편향 발생 (\(E[\exp(\log Y)] \neq \exp(E[\log Y])\))
• 분산 구조에 대한 선택권을 잃음 (변환 후 정규 가정으로 고정)

GLM 에서 \(Y\) 는 그대로 두고 \(\mu\) 만 변환하므로 해석·예측 스케일 일관성 이 보존된다.

8 응용 분야

분야	활용	구체적 예시
손해보험	청구 분포 링크 선택	감마 오차 + 멱 링크, 역수 ↔︎ 로그 중 데이터가 결정
독성학	Tolerance 분포	logit vs cloglog 검정, 사망률 용량반응 곡선
역학	생존 확률 링크	비례 위험 가정 검증 (cloglog 근접성)
계량경제	Box-Cox 변환 검정	종속변수 변환 필요성을 이탈도로 정량화
생태학	종 발생 확률	환경 기울기에 대한 반응 곡선의 비대칭성 검정
임상시험	용량 반응	ED50 추정 모형의 링크 형태 진단
신뢰성 공학	고장 확률	Weibull ↔︎ 지수 사이 보간, cloglog 가족
마케팅	전환율	Conversion rate 의 비대칭성 (logit 타당성 확인)

9 예시 — Heart disease 연구의 logit vs cloglog

로지스틱 회귀로 심장병 발생 확률 \(\pi\) 를 나이·콜레스테롤·혈압 등으로 모형화했다. 이탈도 잔차가 패턴을 보여 링크 적절성을 의심한다.

단일 모수 이항 링크 가족 적용:

\[ g(\mu; \lambda) = \log\!\left[\frac{(1/(1-\pi))^\lambda - 1}{\lambda}\right] \]

\(\lambda\) 를 \([0, 1.5]\) 에서 격자 탐색 → 이탈도 곡선.

• \(\hat\lambda = 0.35\), 95% CI \(= [0.15, 0.60]\)
• \(\lambda = 1\) (로지스틱) 은 기각, \(\lambda = 0\) (cloglog) 는 기각 경계

해석: 로지스틱보다 비대칭 링크가 더 맞지만 완전한 cloglog 까지는 아니다. 확률이 1 에 가까워질 때 기울기가 로지스틱보다 완만해지는 패턴 — 위험 인자 누적에 따른 포화 효과.

실무 결정: 해석 편의상 cloglog (\(\lambda = 0\)) 로 진행하되, 신뢰구간이 경계에 있음을 보고서에 명시.

10 코드 예시 — Pregibon 선형화 구현

10.1 Step 1: Python — 멱 가족 \(\lambda\) 추정 (감마 오차)

import numpy as np

rng = np.random.default_rng(11)
n = 250
x = rng.uniform(0.5, 3.0, n)
# lambda_true = 0.5 (제곱근 링크)로 데이터 생성
beta_true = np.array([1.5, 0.5])
eta_true = beta_true[0] + beta_true[1] * x
mu_true = eta_true ** (1.0 / 0.5)                # eta = mu^0.5
y = rng.gamma(shape=8.0, scale=mu_true / 8.0)

def fit_power_link(x, y, lam, n_iter=100, tol=1e-9):
    """lam-멱 링크 + 감마 분산으로 GLM 적합."""
    X = np.column_stack([np.ones(len(y)), x])
    mu = np.full(len(y), y.mean())
    for _ in range(n_iter):
        eta = mu ** lam if lam != 0 else np.log(mu)
        d_eta_d_mu = lam * mu ** (lam - 1) if lam != 0 else 1.0 / mu
        V = mu ** 2                                # 감마 분산함수
        w = 1.0 / (d_eta_d_mu ** 2 * V)
        z = eta + (y - mu) * d_eta_d_mu
        WX = X * w[:, None]
        beta = np.linalg.solve(X.T @ WX, X.T @ (w * z))
        eta_new = X @ beta
        mu_new = (eta_new ** (1.0 / lam)) if lam != 0 else np.exp(eta_new)
        mu_new = np.where(mu_new > 1e-6, mu_new, mu)
        if np.max(np.abs(mu_new - mu)) < tol:
            mu = mu_new
            break
        mu = mu_new
    dev = 2 * np.sum((y - mu) / mu - np.log(y / mu))  # 감마 이탈도
    return beta, mu, dev

def pregibon_step(x, y, lam0):
    """Pregibon 선형화 한 스텝 — 보조 공변량 추가 후 lambda 갱신."""
    beta, mu, _ = fit_power_link(x, y, lam0)
    v = -(mu ** lam0) * np.log(np.maximum(mu, 1e-12))
    X_aug = np.column_stack([np.ones(len(y)), x, v])
    # 확장 설계행렬로 lambda0-링크 GLM 적합 (한 번의 IRLS 반복 추가)
    eta = mu ** lam0
    d_eta_d_mu = lam0 * mu ** (lam0 - 1)
    w = 1.0 / (d_eta_d_mu ** 2 * mu ** 2)
    z = eta + (y - mu) * d_eta_d_mu
    WX = X_aug * w[:, None]
    coef = np.linalg.solve(X_aug.T @ WX, X_aug.T @ (w * z))
    gamma_hat = coef[-1]
    return lam0 + gamma_hat, coef

# Step A: 격자 탐색으로 대략 위치 확인
lams = np.linspace(-1.5, 1.5, 31)
devs = np.array([fit_power_link(x, y, l)[2] for l in lams])
lam_init = lams[np.argmin(devs)]
print(f"격자 최솟값 위치: lambda = {lam_init:.2f}")

# Step B: 그 근방에서 Pregibon 반복
lam = lam_init
for it in range(25):
    lam_new, coef = pregibon_step(x, y, lam)
    lam = 0.6 * lam + 0.4 * lam_new                # 감쇠
    if abs(lam_new - lam) < 1e-5:
        break
print(f"Pregibon 수렴: lambda_hat = {lam:.4f}  (true = 0.5)")

격자 탐색으로 전역 최솟값 영역을 먼저 확보한 뒤 Pregibon 으로 정밀화하는 것이 안정적 실무.

10.2 Step 2: Python — Profile deviance 곡선과 95% CI

import matplotlib.pyplot as plt
from scipy.stats import chi2

# 격자에서 이탈도 프로파일
lams_fine = np.linspace(-1.5, 1.5, 61)
devs_fine = np.array([fit_power_link(x, y, l)[2] for l in lams_fine])

dev_min = np.min(devs_fine)
lam_hat = lams_fine[np.argmin(devs_fine)]
cutoff = dev_min + chi2.ppf(0.95, 1)

# 95% 신뢰구간을 스캔으로 추출
in_ci = devs_fine < cutoff
ci_low = lams_fine[in_ci].min()
ci_high = lams_fine[in_ci].max()

plt.plot(lams_fine, devs_fine)
plt.axhline(cutoff, color="red", linestyle="--", label="95% cutoff")
plt.axvline(lam_hat, color="green", linestyle=":", label=f"lambda_hat = {lam_hat:.2f}")
plt.xlabel("lambda")
plt.ylabel("Deviance")
plt.legend()
print(f"95% CI for lambda: [{ci_low:.2f}, {ci_high:.2f}]")

이 곡선이 링크 선택의 불확실성을 한눈에 보여준다. CI 안에 0 (로그) 과 0.5 (제곱근) 이 모두 있으면 둘 다 허용 가능.

10.3 Step 3: R — glm 에 custom power link 주입

# R 에서 power link 를 make.link 로 정의
power_link <- function(lambda) {
  if (lambda == 0) return(make.link("log"))
  structure(list(
    linkfun  = function(mu)   mu^lambda,
    linkinv  = function(eta)  pmax(eta, 1e-8)^(1/lambda),
    mu.eta   = function(eta)  (1/lambda) * pmax(eta, 1e-8)^(1/lambda - 1),
    valideta = function(eta)  all(is.finite(eta) & eta > 0),
    name     = sprintf("power(%g)", lambda)
  ), class = "link-glm")
}

set.seed(11)
n <- 250
x <- runif(n, 0.5, 3.0)
mu_true <- (1.5 + 0.5 * x)^2
y <- rgamma(n, shape = 8, scale = mu_true/8)

lams <- seq(-1.5, 1.5, by = 0.05)
devs <- sapply(lams, function(l) {
  fit <- try(glm(y ~ x, family = Gamma(link = power_link(l))), silent = TRUE)
  if (inherits(fit, "try-error")) NA else fit$deviance
})
plot(lams, devs, type = "l", xlab = expression(lambda), ylab = "Deviance")
abline(h = min(devs, na.rm = TRUE) + qchisq(0.95, 1), col = "red", lty = 2)

10.4 Step 4: 이항 링크 가족 — logit vs cloglog 검정

import numpy as np
from scipy.optimize import minimize_scalar

rng = np.random.default_rng(7)
n = 500
x = rng.normal(size=n)
# cloglog 로 데이터 생성 (lambda* -> 0)
eta_true = -0.5 + 0.8 * x
pi_true = 1 - np.exp(-np.exp(eta_true))
y = rng.binomial(1, pi_true)

def inv_link(eta, lam, eps=1e-10):
    """lam-링크의 역함수.
    lam=1: logit (pi = expit(eta)),
    lam->0 한계: cloglog (pi = 1 - exp(-exp(eta))).
    """
    if abs(lam) < 1e-8:
        return 1.0 - np.exp(-np.exp(eta))
    # log[((1/(1-pi))^lam - 1)/lam] = eta
    # → (1/(1-pi))^lam = 1 + lam * exp(eta)
    base = np.maximum(1.0 + lam * np.exp(eta), eps)
    return 1.0 - base ** (-1.0 / lam)

def neg_loglik_lambda(lam, x, y, n_iter=40, tol=1e-8):
    """lam 고정 → beta 를 IRLS 로 적합 → 이항 로그우도의 음수 반환."""
    X = np.column_stack([np.ones(len(y)), x])
    beta = np.zeros(2)
    for _ in range(n_iter):
        eta = X @ beta
        pi = np.clip(inv_link(eta, lam), 1e-10, 1 - 1e-10)
        # d pi / d eta: 수치 미분 (custom 링크 일반화)
        h = 1e-5
        dpi_deta = (inv_link(eta + h, lam) - inv_link(eta - h, lam)) / (2 * h)
        dpi_deta = np.maximum(dpi_deta, 1e-12)
        V = pi * (1 - pi)
        w = dpi_deta ** 2 / V
        z = eta + (y - pi) / dpi_deta
        WX = X * w[:, None]
        beta_new = np.linalg.solve(X.T @ WX, X.T @ (w * z))
        if np.max(np.abs(beta_new - beta)) < tol:
            beta = beta_new
            break
        beta = beta_new
    pi = np.clip(inv_link(X @ beta, lam), 1e-10, 1 - 1e-10)
    return -np.sum(y * np.log(pi) + (1 - y) * np.log(1 - pi))

# 격자 탐색
lams = np.linspace(0.01, 1.5, 30)
nlls = np.array([neg_loglik_lambda(l, x, y) for l in lams])
lam_hat_grid = lams[np.argmin(nlls)]

# Brent 법으로 정밀화
res = minimize_scalar(neg_loglik_lambda, bounds=(0.01, 1.5),
                      args=(x, y), method='bounded',
                      options={'xatol': 1e-4})
lam_hat = res.x

print(f"lambda_hat (grid)  = {lam_hat_grid:.3f}")
print(f"lambda_hat (Brent) = {lam_hat:.3f}   (true -> 0: cloglog)")

예상 출력: lambda_hat 이 0 근처로 수렴해 cloglog 가 데이터에 맞음을 보여준다. 로지스틱 (\(\lambda = 1\)) 과의 이탈도 차이를 \(\chi^2_1\) 에 대입하면 goodness-of-link 검정이 완성된다. 이 골격은 단일 모수 이항 링크 가족 전반에 적용 가능하며, 수치 미분 부분을 해석 미분으로 바꾸면 속도가 크게 향상된다.

11 진단과 실무 가이드

진단	어떤 문제를 잡는가
\(D(\lambda)\) 곡선	여러 극소점, 평탄한 영역 (식별성 부족)
95% CI 의 폭	링크 선택의 데이터 증거 강도
Pregibon 반복의 진동	1차 근사 영역 벗어남 → 격자 재시작
잔차 vs 예측값 플롯	잘못된 링크는 잔차에 체계적 패턴 (bow-shape)
\(D(\lambda_0) - D(\hat\lambda)\)	Goodness-of-link \(\chi^2_1\) 검정

실무 체크리스트:

1. 표준 링크 (canonical 또는 통례) 로 먼저 적합하고 잔차 플롯 검사
2. 잔차가 체계적 패턴을 보이면 링크 의심 — 멱 가족 격자 탐색
3. \(\hat\lambda\) 가 0 이나 -1 같은 정수에 충분히 가까우면 그 표준 링크 채택 (해석·구현 간편)
4. \(\hat\lambda\) 가 표준 링크에서 유의하게 떨어지면 비표준 링크 사용 고려 — 단, 예측보다 해석에서 손실 감수
5. 두 개 이상 링크 모수는 대부분 과적합 — 강한 이론적 이유 없이는 피한다
6. 여러 유사 데이터셋에서 일관되게 \(\hat\lambda\) 가 같은 방향으로 벗어나면 분야별 관례 로 정착시킬 근거

12 Bibliographic Notes

• Box & Cox (1964): 원자료 변환의 멱 가족. GLM 의 아이디어를 역사적으로 예비함
• Pregibon (1980): Goodness-of-link test. 링크 모수의 보조 공변량 전개
• Nelder & Wedderburn (1972): GLM 의 원전. 데이터 변환 vs 적합값 변환의 분리
• Baker & Nelder (1978): GLIM 매뉴얼. 비표준 링크의 구현 기술
• Stukel (1988): 이항 데이터용 2-모수 링크 가족 (현대 대안)
• Nelder & Pregibon (1987): 링크·분산함수의 결합 추정 체계
• Aranda-Ordaz (1981): 대칭·비대칭 이항 링크 가족

13 정리

연결함수 속 모수 \(\lambda\) 는 “표준 링크의 선택” 이라는 이산 문제를 연속 가족 위의 추정 문제 로 바꾼다. Pregibon 의 선형화는 \(\lambda\) 를 보조 공변량의 계수로 변환하여, 표준 IRLS 프레임워크 안에서 모두 처리한다.

1. 멱 가족 \(g(\mu; \lambda) = \mu^\lambda\) 은 역수·로그·제곱근·항등의 표준 링크를 한 축으로 이음
2. Pregibon 선형화: Taylor 전개 → 보조 공변량 \(v = -\hat\mu^{\lambda_0} \log \hat\mu\) 추가 → 새 계수 \(\gamma = \lambda - \lambda_0\)
3. Goodness-of-link 검정은 이탈도 차이 \(\chi^2_1\) 또는 보조 공변량의 Wald 통계량
4. 다중 모수 가족 (shifted power, logit-cloglog 보간) 은 원리적으로 같으나 식별성 이슈에 주의
5. GLM 의 철학: 링크(가법성) 와 분산함수(이분산성) 의 분리 — Box-Cox 변환의 한 \(\lambda\) 로 두 목표 충족 불가능 문제 해결

다음 절(§11.4)은 공변량 자체의 비선형 모수 — Box-Tidwell 기법, drug mixture, Mitscherlich — 를 다룬다.

14 관련 주제

선행 지식

• Models with Additional Non-Linear Parameters — 개관 (McCullagh Ch.11) — 세 출구 전체 지도
• Parameters in the Variance Function — NB의 \(k\)·반올림 분산 (McCullagh §11.2) — 바로 이전 절
• GLM 적합 알고리즘 — IRLS 완전 유도 (McCullagh §2.5) — 선형화의 안쪽 루프
• The Components of a GLM — Random·Systematic·Link (McCullagh §2.2) — 링크 역할의 기초

후속 주제

• Non-linear Parameters in the Covariates — Box-Tidwell·drug mixture (McCullagh §11.4)
• Examples: Bermuda grass 비료·insecticide-synergist·insulin assay (McCullagh §11.5)
• Model Checking — score test 기하학 확장 (McCullagh Ch.12)

관련 개념

• Quasi-likelihood Functions 개관 (McCullagh Ch.9) — 링크와 분산함수의 분리 철학
• 이항 반응 모형 — 링크 함수·모수 해석 (McCullagh §4.3) — logit·probit·cloglog 비교
• Gamma GLM — 분산 함수·이탈도·산포 추정 (McCullagh §8.3) — 감마의 링크 선택 실무
• Joint Modelling of Mean and Dispersion (McCullagh Ch.10) — 이중 GLM 구조 전례