Non-Linear Parameters in the Covariates (McCullagh §11.4)

Box-Tidwell 선형화 · 반복 2중 루프 · 분산 보정 · 지수합 모형의 식별성 · 약물 혼합의 \(\log(x_1 + \theta x_2)\)

GLM 의 설계행렬 안에 모수가 선형으로 들어오지 않는 상황, 즉 공변량 자체가 \(g(x;\theta)\) 꼴로 미지 모수 \(\theta\) 를 담는 경우를 다룬다. Box-Tidwell(1962) 선형화를 유도하고, 보조 공변량 \(v=\partial g/\partial \theta\) 를 더한 이중 IRLS 가 어떻게 \(\theta\) 를 끌어올리는지, 분산 행렬이 왜 최종 반복에서 \(\hat{\beta}v\) 로 재계산되어야 하는지, 지수합 \(\sum \beta_j e^{k_j x_j}\) 에서 식별성이 쉽게 붕괴되는 이유, 약물 혼합 \(\log(x_1+\theta x_2)\) 에서 profile RSS 곡선이 선형화의 대안이 되는 맥락을 직관과 수식으로 함께 풀어낸다.

Statistics
GLM
Engineering
Optimization
저자

Kwangmin Kim

공개

2026년 04월 19일

1 서론 — 공변량 안에 숨은 모수

GLM 의 선형 예측자(linear predictor)는 이름 그대로 모수 \(\beta\) 에 대해 선형이다.

\[\eta = \sum_{j=1}^{p} \beta_j x_j .\]

그런데 데이터 분석 현장에서 자주 마주치는 상황은 정작 “\(x_j\)” 자리에 들어갈 공변량 자체가 미지 모수를 품는 경우다.

  • 약리학: 두 약물이 혼합되어 있을 때 유효 농도는 \(\log(x_1 + \theta x_2)\)\(\theta\) 는 약물 2 의 상대 효능 (potency)이다.
  • 농학: 비료 반응면에서 토양에 이미 남아있는 양분 \(\alpha_i\) 를 보정한 \(u_i = 1/(x_i + \alpha_i)\) 가 진짜 공변량이다.
  • 경제학: 수익의 로그-오즈(log-odds)가 \(\log(\text{income} - \theta)\) 처럼 임계 소득 \(\theta\) 를 기준점으로 작동한다.
  • 역학: 노출의 비선형 용량-반응이 \(x^\theta\) (Box-Tidwell) 또는 \(e^{kx}\) (지수형 반응)로 모델링된다.

이 문제의 껄끄러운 점은 \(\theta\) 가 알려진 값이면 그대로 \(g(x;\theta)\) 를 설계행렬에 꽂아 평범한 GLM 으로 돌아온다는 것이다. \(\theta\)데이터로부터 추정해야 하는 순간, 모델은 비선형이 된다. 선형 예측자가 모수 \((\beta, \theta)\) 에 대해 비선형이기 때문이다.

§11.2 가 분산 함수 안의 모수, §11.3 이 링크 함수 안의 모수를 다룬다면, 이번 §11.4 는 공변량 안의 모수라는 세 번째 배출구 (outlet)를 다룬다. 기법 자체는 §11.3 과 놀랄 만큼 닮았다 — Taylor 1차 전개로 보조 공변량 하나를 추가해 문제를 GLM 안으로 되돌리는 Box-Tidwell (1962) 의 선형화가 핵심이다 (McCullagh & Nelder, 1989, §11.4).

직관적으로 보면 다음 그림이다.

  • \(\theta\) 는 공변량의 모양 을 조절하는 다이얼이다 (\(\theta\) 가 바뀌면 \(g(x;\theta)\) 전체 곡선이 움직인다).
  • 현재 추정치 \(\theta_0\) 부근에서 곡선이 어느 방향으로 얼마나 움직이는지 를 재는 것이 편미분 \(\partial g/\partial\theta\) 이다.
  • 그 “움직임 벡터” 를 보조 공변량 \(v\) 로 회귀에 끼워 넣으면, 회귀가 \(v\) 에 할당한 계수가 곧 “다이얼을 얼마나 돌려야 하는가”를 알려준다.

2 Box-Tidwell 선형화 — 수식 구성

공변량이

\[g(x; \theta)\]

꼴이고, 선형 예측자의 해당 항이

\[\beta\, g(x;\theta)\]

라고 하자. 현재 추정치 \(\theta_0\) 근방에서 \(g\) 를 Taylor 1차로 전개한다.

\[ g(x;\theta) \simeq g(x;\theta_0) + (\theta - \theta_0)\,\left[\frac{\partial g}{\partial \theta}\right]_{\theta=\theta_0}. \]

그러면 \(\beta g(x;\theta)\)

\[ \beta g(x;\theta) \simeq \beta\, g(x;\theta_0) + \beta(\theta-\theta_0)\left[\frac{\partial g}{\partial \theta}\right]_{\theta=\theta_0} = \beta u + \gamma v, \]

단,

\[u = g(x;\theta_0), \qquad v = \left[\frac{\partial g}{\partial \theta}\right]_{\theta=\theta_0}, \qquad \gamma = \beta(\theta-\theta_0).\]

핵심 관찰: 비선형 항 하나가 선형 항 두 개 로 분해된다. 둘 다 GLM 이 다룰 수 있는 “보통 공변량” 이다. \(u\)\(\theta_0\) 에서 평가한 현재의 공변량이고, \(v\)\(\theta\) 가 움직일 때 \(g\) 가 어느 방향으로 반응하는지 알려주는 방향 벡터 다.

이제 평범한 GLM 을 \(u\)\(v\) 로 적합하면 회귀 계수 \(\hat\beta\)\(\hat\gamma\) 가 나온다. Taylor 전개의 정의에 의해 \(\gamma = \beta(\theta-\theta_0)\) 이므로,

\[ \boxed{\;\theta_1 = \theta_0 + \frac{\hat\gamma}{\hat\beta}\;} \]

갱신된 추정치 다.

분자·분모의 의미. \(\hat\beta\) 는 “공변량 \(u\) 에 대한 \(y\) 의 기울기” — 반응이 모수 \(\beta\) 에 얼마나 민감한가. \(\hat\gamma\) 는 “보조 공변량 \(v\) 로 잡힌 추가 기울기” — \(\theta_0\) 에 맞춘 \(u\) 만으로는 설명되지 않고 \(\theta\) 방향으로 더 움직여야 할 신호 의 크기다. 비율 \(\hat\gamma/\hat\beta\) 로 나누는 행위는 이 추가 신호를 “기울기 단위” 로 환산하는 것 — \(\hat\beta\) 가 크면 작은 \(\hat\gamma\) 로도 충분하고, \(\hat\beta\) 가 0 근처면 \(\hat\gamma\) 가 작아도 보정이 크게 튀어 수렴이 불안정해진다. 이것이 Box-Tidwell 이터레이션의 외부 루프다. 내부 루프는 평범한 IRLS (GLM fit)이므로, 전체 구조는 GLM IRLS 를 감싸는 또 하나의 반복 이 된다.

왜 “선형화” 라고 부르나

\(\theta\)\(g\) 에 대해 비선형이지만, \(\gamma\) 는 설계행렬의 새 열 \(v\) 에 대해 선형 이다. Taylor 1차 전개로 비선형 파라미터 공간을 선형 회귀 공간에 투영하여, 매 반복마다 방향·크기 를 GLM 이 알려주도록 만든 것이다. 그래서 이 아이디어는 §11.3 의 Pregibon 링크 선형화와 정확히 같은 수학적 뼈대 를 공유한다 — 단지 Taylor 전개의 대상이 링크 \(g(\mu)\) 대신 공변량 \(g(x)\) 로 바뀌었을 뿐이다.

3 알고리즘 — 이중 루프 구조

θ = θ₀  (초기값)
repeat:
    u ← g(x; θ)
    v ← ∂g/∂θ 를 현재 θ 에서 평가
    (β̂, γ̂, 나머지 계수) ← GLM(y ~ u + v + 다른_공변량, family=...)
    θ ← θ + γ̂ / β̂
until |γ̂ / β̂| < tol

각 반복에서 외부 루프는 \(\theta\) 를 갱신하고, 내부 루프는 갱신된 \(\theta\) 에서 일반 IRLS 로 \(\beta\) 를 찾는다. 두 루프가 교대로 수렴에 기여 한다. 이것이 §11.3 의 Pregibon 선형화 알고리즘과 닮은 까닭이다.

3.1 수렴 보장의 한계

McCullagh & Nelder 가 명시적으로 경고한다 (p.379, 원문: “Convergence is not guaranteed for starting values arbitrarily far from the solution”). 초기값 \(\theta_0\) 가 해에서 멀리 떨어져 있으면

  • \(\hat\beta\)부호가 뒤집히거나 0 에 가까워져 \(\hat\gamma/\hat\beta\) 가 발산하고
  • \(u, v\) 의 스케일 차이가 커지면 조건수가 폭발 해 GLM 내부 IRLS 가 먼저 실패한다.

실무적 대응:

  • 초기값은 시각적 검토 (데이터 그림을 보고 \(\theta\) 가 만들 만한 곡선 형태를 어림) 또는 격자 탐색 (grid search)으로 얻는다.
  • 갱신량에 댐핑(damping) 을 건다: \(\theta_{k+1} = \theta_k + \lambda \hat\gamma/\hat\beta\), \(\lambda \in (0,1]\).
  • 수렴이 안 되면 profile likelihood (뒤에 설명)로 대체를 고려한다.

4 분산·공분산의 올바른 계산

Taylor 전개로 얻은 GLM 은 \(\beta, \gamma\) 에 대한 보통의 GLM 이다. 따라서 최종 반복에서 얻은 \((X^\top W X)^{-1}\)\(\hat\beta, \hat\gamma\) 에 대한 분산을 준다. 문제는 우리가 원하는 것이 \(\hat\theta\) 에 대한 분산이라는 점이다.

관계식

\[\gamma = \beta(\theta - \theta_0)\]

에서 델타 방법 (delta method) 을 쓸 수도 있지만, McCullagh 는 더 깔끔한 길을 제시한다. 수렴 마지막 한 번 의 반복에서 보조 공변량을 \(v\) 대신 \(\hat\beta v\) 로 바꿔 넣어 재적합한다.

이때 설계행렬의 새 열 \(\hat\beta v\) 에 대응되는 계수는

\[\gamma^* = (\theta - \theta_0)\]

가 되어, \((X^\top W X)^{-1}\) 의 해당 대각 원소가 곧바로 \(\mathrm{Var}(\hat\theta)\) 를, 해당 행이 \(\hat\theta\) 와 다른 계수들 간의 공분산 을 준다.

왜 “\(\hat\beta v\)” 로 대체해야 분산이 맞는가

직관: \(v\) 는 “\(\theta\) 가 단위만큼 움직이면 \(g\) 가 얼마 변하는가” 의 방향이지만, 실제 모델이 관측하는 민감도는 계수 \(\beta\) 가 곱해진 \(\beta v\) 다. 즉 \(y\) 의 관점에서 \(\theta\) 가 1 만큼 변하면 반응은 \(\beta \cdot (\partial g/\partial \theta) = \beta v\) 만큼 바뀐다. 이 실제 민감도 벡터 를 설계행렬에 넣어야 역정보행렬이 곧장 \(\mathrm{Var}(\hat\theta)\) 를 읽어낸다.

달리 말하면 \(\hat\gamma = \hat\beta(\hat\theta - \theta_0)\) 를 직접 재조정해서, 새 계수가 \(\hat\theta - \theta_0\) 자체 가 되도록 축을 재정규화한 것이다.

공분산 보정의 또 다른 역할: 선형화로 인해 \(\hat\theta\)\(\hat\beta\)상관 된다. 단순히 \(\theta\) 를 고정하고 \(\beta\) 만 적합하면 \(\hat\beta\) 의 표준오차가 과소평가된다. \(v\) (또는 \(\hat\beta v\)) 를 포함한 마지막 적합은 이 상관을 자동으로 조정한다 — 이는 §11.3 의 링크 선형화, §11.2 의 분산 모수 선형화와 동일한 원리다.

5 지수합 모형 — 식별성의 함정

가장 빈번하게 만나는 공변량 비선형성은 지수합 (sum of exponentials) 형태다.

\[\eta = \beta_0 + \beta_1 e^{k_1 x_1} + \beta_2 e^{k_2 x_2}.\]

약동학(compartmental pharmacokinetics), 방사선 붕괴 (radioactive decay), 경제 성장 모형 등에서 자연스럽게 등장한다. McCullagh 는 이런 모형에서 비선형 모수의 표준오차가 매우 크고, 선형 모수·서로 간에 상관이 극심 해지는 경향을 경고한다 (p.379).

5.1 직관 — 두 지수가 서로 흉내낼 수 있다

\(\beta_1 e^{k_1 x} + \beta_2 e^{k_2 x}\) 에서 \(k_1 \approx k_2\) 이면

  • \(\beta_1 e^{k_1 x} + \beta_2 e^{k_2 x} \approx (\beta_1+\beta_2) e^{k_0 x}\) (거의 하나의 지수로 축약된다)
  • 데이터는 \(\beta_1+\beta_2\) 만 정확히 말해주고, 각각의 \(\beta_1, \beta_2\)무한히 많은 조합 이 비슷한 적합도를 준다.

반대 극단으로 \(k_1\) 이 너무 크면 \(e^{k_1 x}\) 가 표본 범위 밖에서 폭발해서 수치적으로 한두 점 에 의해 좌우된다. 식별성(identifiability)이 수치적으로 깨지는 양쪽 낭떠러지 가 있다.

5.2 공변량 상관과의 공모

원래 공변량들 \(x_1, x_2\) 가 표본 내에서 상관이 있으면 (\(\mathrm{corr}(x_1, x_2)\) 가 크면) \(e^{k_1 x_1}\)\(e^{k_2 x_2}\) 도 강하게 상관되어 더 심해진다. McCullagh 는 이를 “the other covariates are themselves appreciably correlated” (p.379) 상황이라 부르며, 2개 이상의 비선형 모수 를 한 모형에 넣는 것을 경계한다.

실무 조언:

  • 모형 1: \(e^{k_1 x_1}\) 하나만 포함 — \(k_1\) 은 잘 추정됨.
  • 모형 2: \(e^{k_1 x_1} + e^{k_2 x_2}\)\(k_1, k_2\) 모두 큰 SE.
  • 모형 3: \(\beta_1 e^{k x_1} + \beta_2 e^{k x_2}\) (\(k\)공유 하도록 제약) — 자주 더 안정적이다.

제약을 통해 모수 공간을 줄이는 것 자체가 식별성에 대한 강한 무기 다.

6 약물 혼합 — \(\log(x_1 + \theta x_2)\) 의 사례

비선형 공변량이 과학적으로 의미 있게 등장하는 대표 예시는 두 약물 혼합의 joint action 모형이다 (McCullagh & Nelder §11.5.3 예고).

두 약물 \(A, B\)동일 작용기전 을 공유하면, 혼합물의 유효 용량은

\[x_1 + \theta x_2\]

로 표현된다. 여기서 \(\theta\) 는 약물 \(B\) 를 약물 \(A\) 로 환산했을 때의 상대 효능 (relative potency) 이다. 예를 들어 \(\theta=2\) 이면 \(B\) 가 같은 양에서 \(A\) 의 두 배 효과를 낸다는 뜻이다.

로그-용량 반응 (log-dose response) 이 관행이므로, 선형 예측자는

\[\eta = \alpha + \beta \log(x_1 + \theta x_2)\]

가 된다. 이 공변량이 바로 \(g(x_1, x_2; \theta) = \log(x_1 + \theta x_2)\) 이다.

6.1 두 가지 분석 전략

전략 A: Box-Tidwell 선형화

\(\partial g/\partial \theta = x_2 / (x_1 + \theta x_2)\) 이므로, 현재 \(\theta_0\) 에서 보조 공변량은

\[u = \log(x_1 + \theta_0 x_2), \qquad v = \frac{x_2}{x_1 + \theta_0 x_2}.\]

이 둘을 넣어 GLM 을 적합하고 \(\theta_1 = \theta_0 + \hat\gamma/\hat\beta\) 로 갱신한다. 표준오차·공분산을 한 번에 얻는 장점이 있다.

전략 B: Profile deviance plot (Darby & Ellis, 1976)

\(\theta\) 를 격자 \(\{\theta^{(1)}, \dots, \theta^{(M)}\}\) 로 돌려가며 \(\theta\) 고정 상태의 MLE 을 계산한다.

\[\ell_P(\theta) = \max_{\alpha, \beta} \ell(\alpha, \beta, \theta) \quad \Rightarrow \quad \mathrm{RSS}(\theta)\]

\(\theta\) 에 대해 그린다. 곡선의 최저점\(\hat\theta\) 다. \(F\)-검정 기반의 근사 신뢰구간은

\[ \Bigl\{\theta : \mathrm{RSS}(\theta) - \mathrm{RSS}(\hat\theta) < s^2 F^*_{1, n-p, \alpha} \Bigr\}, \]

\(s^2 = D(\hat\theta)/(n-p)\) 는 잔차 평균 이탈도다 (McCullagh & Nelder, 1989, p.380).

두 전략의 트레이드오프
항목 Box-Tidwell 선형화 Profile deviance
신뢰구간 Wald (대칭, 근사) F-기반 (비대칭 허용, 정확도↑)
\(\hat\theta\) 공분산 계산 한 번에 가능 어렵다 — 추가 작업 필요
수렴 보장 초기값 민감 격자면 항상 가능
여러 비선형 모수 그대로 확장 차원의 저주
구현 난이도 IRLS + 외부 루프 적합을 수십 회 반복

McCullagh 권고: 공분산이 필요하면 선형화, 식별성·수렴에 불안감이 있으면 profile → 최저점 확인 후 선형화로 보완.

7 농업 반응면 — 토양 보정 \(1/(x_i + \alpha_i)\)

§11.5.1 이 다룰 Welch et al. (1963) coastal Bermuda grass 실험은 비선형 공변량의 또 다른 교과서적 사례다. 세 비료 \(N, P, K\) 의 역선형 반응면

\[1/\mu = \eta = \beta_0 + \sum_{i=1}^{3} \beta_i u_i, \qquad u_i = \frac{1}{x_i + \alpha_i}\]

에서 \(\alpha_i\)토양에 이미 존재하는 양분 을 모델링하는 비선형 모수다. 무처리 구(0,0,0)에서도 수확량이 0 이 아니라는 점이 \(\alpha_i > 0\) 의 필요성을 데이터가 직접 말해주는 증거다.

7.1 Box-Tidwell 적용

\(u_i = 1/(x_i+\alpha_i)\)\(\alpha_i\) 에 대해 미분하면

\[v_i = \frac{\partial u_i}{\partial \alpha_i} = -\frac{1}{(x_i+\alpha_i)^2} = -u_i^2.\]

따라서 보조 공변량은 현재 \(\alpha_i^{(0)}\) 에서

\[u_i = \frac{1}{x_i + \alpha_i^{(0)}}, \qquad v_i = -u_i^2\]

가 되고, GLM (감마 오차, 역수 링크) 을 \(u_1, u_2, u_3, v_1, v_2, v_3\) 로 적합하면 각 \(v_i\) 의 계수 \(\hat c_i\)\(u_i\) 의 계수 \(\hat b_i\) 로부터

\[\delta \alpha_i = \hat c_i / \hat b_i\]

를 얻어 \(\alpha_i^{(1)} = \alpha_i^{(0)} + \delta \alpha_i\) 로 갱신한다 (McCullagh & Nelder, 1989, p.381).

7.2 왜 감마/역수 링크인가

“수확량의 비율적 표준편차 (proportional standard deviation) 가 일정” 이라는 실무 가정이 \(V(\mu) = \mu^2\) 를 함축하고, 이는 감마 분포의 분산 함수다. 역수 링크는 inverse polynomial 반응면이 \(1/\mu\) 에서 자연스럽게 선형화되기 때문에 선택된다. 즉 공변량의 비선형성은 \(\alpha_i\) 에 있지, 링크·분산에 있는 것이 아니다 — 이것이 §11.4 가 §11.2·11.3 과 독립적으로 서술되는 이유다.

8 다른 대표적 공변량 비선형 형태

형태 \(g(x;\theta)\) \(\partial g/\partial \theta\) 응용
Box-Tidwell 멱 \(x^\theta\) \(x^\theta \log x\) 용량-반응, 지질학, 심리측정
지수 감쇠 \(e^{-\theta x}\) \(-x e^{-\theta x}\) 약동학, 방사선
로지스틱 용량 \(\{1 + e^{-\theta(x-x_0)}\}^{-1}\) 체인룰 생리학, 임상
힐 (Hill) 함수 \(x^n / (K^\theta + x^n)\) 체인룰 효소 동역학
임계점 (threshold) \((x - \theta)_+\) \(-\mathbb{1}[x > \theta]\) 경제학, 환경
약물 혼합 \(\log(x_1 + \theta x_2)\) \(x_2/(x_1+\theta x_2)\) 약리학

공통 직관: \(g\) 의 모양을 \(\theta\) 로 바꿀 수 있는 한, 모든 형태는 같은 Box-Tidwell 틀로 들어온다. 체인룰을 써서 \(v\) 만 계산해 주면 나머지는 IRLS 에 위임한다.

9 Box-Tidwell vs Pregibon — 두 선형화의 비교

§11.3 (Pregibon) 과 §11.4 (Box-Tidwell) 는 수학적 구조가 거의 동일하지만 Taylor 전개의 좌표축 이 다르다.

차원 Pregibon §11.3 Box-Tidwell §11.4
미지 모수의 위치 링크 함수 \(g(\mu;\lambda)\) 공변량 \(g(x;\theta)\)
Taylor 전개 대상 \(g(\mu;\lambda)\)\(\lambda\) 에 대해 \(g(x;\theta)\)\(\theta\) 에 대해
보조 공변량 \(v = \partial g/\partial \lambda\), 현재 적합 \(\hat\mu\) 에서 평가 \(v = \partial g/\partial \theta\), 데이터 \(x\) 에서 평가
반복 구조 외부(링크 갱신) + 내부(IRLS) 외부(\(\theta\) 갱신) + 내부(IRLS)
분산 보정 마지막 반복에 \(\hat\beta v\) 투입 마지막 반복에 \(\hat\beta v\) 투입
대표 예시 Box-Cox 멱 링크 \(\mu^\lambda\) Box-Tidwell 멱 공변량 \(x^\theta\)

핵심 교훈: GLM 확장의 세 배출구 (분산 §11.2, 링크 §11.3, 공변량 §11.4) 중 뒤의 두 개는 동일한 선형화 엔진 으로 굴러간다. 이것이 McCullagh 가 §11.1 에서 “세 확장 모두 GLM 의 IRLS 를 재활용한다” 고 강조한 이유다.

10 구현 — Python (Box-Tidwell via statsmodels)

임의의 GLM 에서 Box-Tidwell 을 일반화하는 스켈레톤을 제공한다.

코드
import numpy as np
import statsmodels.api as sm

def fit_boxtidwell_power(y, x, family=sm.families.Gaussian(),
                          theta0=1.0, tol=1e-5, max_iter=50, damp=0.7):
    """
    공변량 x 가 x^theta 형태로 들어갈 때 theta 와 beta 를 동시 추정.
    y = alpha + beta * x^theta + eps (GLM family 로 일반화)

    Parameters
    ----------
    y : (n,) 반응
    x : (n,) 양수 공변량
    family : sm.families.Family
    theta0 : 초기 theta
    damp : (0,1] 갱신 댐핑
    """
    x = np.asarray(x, dtype=float)
    assert np.all(x > 0), "x^theta 는 x > 0 을 요구한다"

    theta = theta0
    history = []
    for k in range(max_iter):
        # 1) 보조 공변량 구성
        u = x ** theta              # g(x;theta)
        v = (x ** theta) * np.log(x)  # dg/dtheta = x^theta * ln(x)

        X = sm.add_constant(np.column_stack([u, v]))

        # 2) GLM 적합
        model = sm.GLM(y, X, family=family).fit()
        alpha_hat, beta_hat, gamma_hat = model.params

        # 3) theta 갱신 (damping)
        if abs(beta_hat) < 1e-8:
            raise RuntimeError("beta 가 0 에 가까워 갱신 실패")

        delta = gamma_hat / beta_hat
        theta_new = theta + damp * delta
        history.append((k, theta, theta_new, beta_hat, gamma_hat))

        if abs(theta_new - theta) < tol:
            theta = theta_new
            break
        theta = theta_new

    # 4) 마지막 반복 — 분산 보정: v 대신 beta_hat * v 를 투입
    u = x ** theta
    v = (x ** theta) * np.log(x)
    X_final = sm.add_constant(np.column_stack([u, beta_hat * v]))
    model_final = sm.GLM(y, X_final, family=family).fit()
    # model_final.cov_params() 에서 3번째 대각 원소가 Var(theta_hat) 근사

    return {
        "theta_hat": theta,
        "beta_hat": beta_hat,
        "alpha_hat": alpha_hat,
        "se_theta": np.sqrt(model_final.cov_params()[2, 2]),
        "iterations": len(history),
        "history": history,
        "final_model": model_final,
    }

# 합성 예시
rng = np.random.default_rng(0)
n = 200
x = rng.uniform(0.5, 5.0, n)
true_theta, true_alpha, true_beta = 0.6, 1.0, 2.0
mu = true_alpha + true_beta * x ** true_theta
y = mu + rng.normal(0, 0.2, n)

res = fit_boxtidwell_power(y, x, theta0=1.5)
print(f"theta_hat = {res['theta_hat']:.4f} (true 0.60)")
print(f"  SE      = {res['se_theta']:.4f}")
print(f"beta_hat  = {res['beta_hat']:.4f}")

10.1 Profile deviance 플롯 (대안 추정법)

코드
import matplotlib.pyplot as plt

thetas = np.linspace(0.2, 1.2, 31)
rss = []
for th in thetas:
    u = x ** th
    X = sm.add_constant(u)
    fit = sm.OLS(y, X).fit()
    rss.append(fit.ssr)
rss = np.array(rss)

# F-기반 근사 신뢰구간 (95%)
from scipy.stats import f
n, p = len(y), 3  # alpha, beta, theta
s2 = rss.min() / (n - p)
F95 = f.ppf(0.95, 1, n - p)
cutoff = rss.min() + s2 * F95

plt.plot(thetas, rss)
plt.axhline(cutoff, color="red", linestyle="--", label="95% CI cutoff")
plt.xlabel(r"$\theta$")
plt.ylabel("RSS")
plt.title("Profile RSS — Box-Tidwell power")
plt.legend(); plt.show()

곡선 최저점의 \(\theta\) 와 cutoff 선 아래 영역이 근사 95% 신뢰구간을 준다.

11 구현 — R (nls + Bermuda grass 스타일)

코드
library(stats)

# 약물 혼합 예시: y = alpha + beta * log(x1 + theta * x2) + eps
# nls 는 내부적으로 Gauss-Newton 을 사용 — Box-Tidwell 선형화의 한 구현.

set.seed(42)
n <- 150
x1 <- runif(n, 0.5, 5)
x2 <- runif(n, 0.5, 5)
true_theta <- 2.5
y <- 1.0 + 0.5 * log(x1 + true_theta * x2) + rnorm(n, 0, 0.1)

fit <- nls(y ~ alpha + beta * log(x1 + theta * x2),
           start = list(alpha = 0, beta = 1, theta = 1),
           trace = TRUE)
summary(fit)

# profile deviance
theta_grid <- seq(0.5, 5, length.out = 40)
rss <- sapply(theta_grid, function(th) {
    z <- log(x1 + th * x2)
    fit_lm <- lm(y ~ z)
    sum(residuals(fit_lm)^2)
})
plot(theta_grid, rss, type = "l", xlab = expression(theta), ylab = "RSS")
abline(v = theta_grid[which.min(rss)], col = "red", lty = 2)

nls() 는 Gauss-Newton 기반 비선형 최소제곱으로, Box-Tidwell 과 동일한 선형화 원리를 일반화한 것이다. GLM 과 결합이 필요하면 gnm (generalized nonlinear models) 패키지의 gnm(... , formula = ~ Mult(1, log(x1 + Mult(1, x2))) ...) 를 사용한다.

12 서지 노트

  • Box, G.E.P. & Tidwell, P.W. (1962). “Transformation of the Independents Variables.” Technometrics — 원전.
  • Darby, S.C. & Ellis, M.J. (1976). 약물 혼합의 profile deviance 분석.
  • Welch, L.F. et al. (1963). coastal Bermuda grass factorial experiment — McCullagh §11.5.1 데이터 출처.
  • Pregibon, D. (1980) — 링크 선형화와 연계 (§11.3 참조).
  • Seber, G.A.F. & Wild, C.J. (2003). Nonlinear Regression — 포괄적 참조.

13 정리

  • 세 번째 배출구: §11.4 는 비선형 모수가 공변량 에 있는 경우를 다룬다. \(\beta g(x;\theta)\) 에서 \(\theta\) 를 추정한다.
  • 선형화: Taylor 1차로 \(g\) 를 풀어 \(\beta u + \gamma v\) 로 분해, 보조 공변량 \(v = \partial g/\partial \theta\) 를 추가.
  • 갱신식: \(\theta_1 = \theta_0 + \hat\gamma/\hat\beta\).
  • 분산 보정: 마지막 반복에서 \(v\) 대신 \(\hat\beta v\) 를 투입해 \((X^\top W X)^{-1}\)\(\mathrm{Var}(\hat\theta)\) 를 직접 내놓도록 만든다.
  • 식별성: 지수합 모형처럼 비선형 모수가 여러 개면 상관이 심해져 SE 가 폭발 — 제약을 통해 차원을 줄이는 것이 실무 해법.
  • 대안: profile deviance plot (Darby-Ellis) — \(\theta\) 격자로 RSS 곡선을 그려 최저점을 찾고 F-기반 신뢰구간을 구성.
  • 엔진 재사용: §11.3 링크 선형화와 수학적 뼈대가 동일 — GLM IRLS 를 감싸는 외부 루프.

14 관련 주제

선행

후속

관련 개념

Subscribe

Enjoy this blog? Get notified of new posts by email: