R에서의 통계 모델(Statistical models in R)

이 장은 통계 분석을 R로 표현하는 핵심 틀인 모델 수식(formula) 과 그 모델을 적합(fit)하는 함수들의 세계로 여러분을 안내해요. 회귀분석과 분산분석 정도는 익숙하다는 전제로 출발하고, 뒤로 갈수록 일반화선형모형과 비선형 회귀까지 살짝 언급하니까, 필요한 만큼만 챙겨 가면 됩니다.

출처: R 공식 매뉴얼

본문

이 절은 독자가 통계 방법론, 특히 회귀분석과 분산분석에 어느 정도 익숙하다고 전제합니다. 뒤로 가면서 조금 더 야심 찬 전제를 하나 얹죠. 일반화선형모형(generalized linear model)과 비선형 회귀에 대해 조금은 알고 있다는 전제입니다.

통계 모델을 적합한다는 작업은 그 요구사항이 충분히 잘 정의되어 있어서, 넓은 범위의 문제에 두루 적용되는 범용 도구를 만들어 낼 수 있을 만큼 명확해요. R은 바로 그 점을 노려서, 통계 모델 적합을 아주 간단하게 만들어 주는 시설들을 서로 맞물린 형태로 갖추고 있습니다. 서문에서도 언급했듯이 기본 출력은 최소한이고, 상세한 내용은 필요할 때 추출 함수(extractor function) 를 불러서 요구해야 해요.

11.1 통계 모델 정의하기: 수식(formulae)

통계 모델의 기본 템플릿은 독립이고 등분산(homoscedastic)인 오차를 가진 선형 회귀 모델입니다.

y_i = sum_{j=0}^p beta_j x_{ij} + e_i,     i = 1, ..., n,

여기서 e_i 는 NID(0, sigma^2) 입니다. 행렬 표기로 쓰면

y = X  beta + e

이렇게 되는데, y 는 반응 벡터(response vector) 이고 X 는 모델 행렬(model matrix) 또는 설계 행렬(design matrix) 이라고 하며, 열 x_0, x_1, ..., x_p 즉 결정 변수(determining variable) 들로 이루어집니다. 아주 흔하게 x_0 은 1로만 채워진 열, 즉 절편(intercept) 항을 정의하는 열이에요.

예시(Examples)

정식 정의를 내리기 전에, 예시 몇 개를 보면 그림이 훨씬 선명해져요.

y, x, x0, x1, x2, … 가 숫자형 변수이고 X 는 행렬, A, B, C, … 는 요인(factor)이라고 합시다. 아래 왼쪽의 수식은 오른쪽에서 설명하는 통계 모델을 지정하며, 그 수식은 다음과 같습니다.

y ~ x
y ~ 1 + x

둘 다 x 에 대한 y 의 단순 선형 회귀 모델을 뜻해요. 차이는 첫 번째가 절편을 묵시적으로 포함하고 두 번째가 명시적으로 포함한다는 점뿐입니다.

y ~ 0 + x
y ~ -1 + x
y ~ x - 1

원점을 지나는(즉 절편 항이 없는) x 에 대한 y 의 단순 선형 회귀입니다.

log(y) ~ x1 + x2

변환된 변수 log(y) 를 x1 과 x2 에 회귀시키는 다중 회귀입니다(묵시적인 절편 항 포함).

y ~ poly(x,2)
y ~ 1 + x + I(x^2)

x 에 대한 y 의 2차 다항 회귀입니다. 첫 번째 형태는 직교 다항식(orthogonal polynomial)을, 두 번째 형태는 명시적인 거듭제곱을 기저로 씁니다.

y ~ X + poly(x,2)

행렬 X 로 이루어진 모델 행렬에 x 의 2차 다항 항을 더한 y 의 다중 회귀입니다.

y ~ A

A 의 수준(level)으로 클래스를 나누는 y 의 단일 분류 분산분석(analysis of variance) 모델입니다.

y ~ A + x

A 의 수준으로 클래스를 나누고 x 를 공변량(covariate)으로 둔 y 의 단일 분류 공분산분석(analysis of covariance) 모델입니다.

y ~ A*B
y ~ A + B + A:B
y ~ B %in% A
y ~ A/B

A 와 B 의 두 요인에 대한 비가법(non-additive) 모델입니다. 처음 두 수식은 같은 교차 분류(crossed classification)를, 뒤의 두 수식은 같은 중첩 분류(nested classification)를 지정해요. 추상적으로 보면 네 수식 모두 같은 모델 부분공간을 가리킵니다.

y ~ (A + B + C)^2
y ~ A*B*C - A:B:C

3요인 실험이되, 모델에는 주효과와 2차 상호작용만 담습니다. 두 수식 모두 같은 모델을 지정해요.

y ~ A * x
y ~ A/x
y ~ A/(1 + x) - 1

A 의 각 수준 안에서 x 에 대한 y 의 개별 단순 선형 회귀 모델로, 부호화(coding) 방식만 다릅니다. 마지막 형태는 A 의 수준 개수만큼 서로 다른 절편과 기울기를 명시적으로 추정해 줘요.

y ~ A*B + Error(C)

두 처리 요인 A, B 와 요인 C 가 결정하는 오차 층(error stratum)을 가진 실험입니다. 예를 들어 요인 C 가 결정하는 전체구(whole plot, 따라서 하위구도 함께)를 가진 분할구(split plot) 실험이죠.

연산자 ~ 는 R에서 모델 수식(model formula) 을 정의할 때 씁니다. 보통의 선형 모델에서 그 형태는

response ~ op_1 term_1 op_2 term_2 op_3 term_3 ...

이며, 여기서

response 는 반응 변수를 정의하는 벡터나 행렬(또는 벡터·행렬로 평가되는 식)이고,
op_i 는 모델에 항을 포함하거나 제외함을 뜻하는 연산자로 + 또는 - 입니다(첫 번째 연산자는 생략 가능).
term_i 는 다음 중 하나입니다.
- 벡터·행렬 표현식 또는 1
- 요인(factor)
- 요인·벡터·행렬을 수식 연산자로 연결한 수식 표현식(formula expression)

어떤 경우든 각 항은 모델 행렬에 더하거나 뺄 열들의 모임을 정의해요. 1 은 절편 열을 뜻하며, 명시적으로 제거하지 않는 한 기본적으로 모델 행렬에 포함됩니다.

수식 연산자는 Glim이나 Genstat 같은 프로그램에서 쓰는 Wilkinson–Rogers 표기와 효과가 비슷합니다. 다만 한 가지 어쩔 수 없는 변화가 있는데, 원래의 연산자 . 가 : 로 바뀌었다는 점이에요. 왜냐하면 마침표(.)는 R에서 유효한 이름 문자이기 때문입니다.

이 표기는 다음과 같이 요약됩니다(Chambers & Hastie, 1992, p.29에 근거).

Y ~ M: Y 가 M 모델로 모델링된다.
M_1 + M_2: M_1 과 M_2 를 모두 포함한다.
M_1 - M_2: M_1 을 포함하되 M_2 의 항들은 빼낸다.
M_1 : M_2: M_1 과 M_2 의 텐서곱. 두 항이 모두 요인이라면 '하위클래스(subclasses)' 요인이 된다.
M_1 %in% M_2: M_1 : M_2 와 비슷하되 부호화 방식이 다르다.
M_1 * M_2: M_1 + M_2 + M_1:M_2 와 같다.
M_1 / M_2: M_1 + M_2 %in% M_1 과 같다.
M^n: M 의 모든 항에 '상호작용'을 n 차까지 포함한다.
I(M): M 을 격리(insulate)한다. M 안에서는 모든 연산자가 원래의 산술적 의미를 가지며, 그 항이 모델 행렬에 나타난다.

주의할 점은, 보통 함수 인자를 감싸는 괄호 안에서는 모든 연산자가 원래의 산술적 의미를 가진다는 거예요. I() 함수는 모델 수식 안에서 산술 연산자로 항을 정의할 수 있게 해 주는 항등 함수(identity function)입니다.

특히 중요한 건, 모델 수식은 모델 행렬의 열들의 스펙을 지정하며 모수(parameter)의 지정은 묵시적이라는 점입니다. 이는 다른 맥락, 예컨대 비선형 모델을 지정할 때와는 다른 이야기예요.

11.1.1 대비(Contrasts)

모델 수식이 어떻게 모델 행렬의 열들을 지정하는지, 적어도 어느 정도는 알아둘 필요가 있어요. 연속형 변수라면 간단합니다. 각 변수마다 모델 행렬의 열 하나를 제공하니까요(절편은 포함된다면 1로만 된 열을 제공하고요).

그럼 k-수준 요인 A 는 어떨까요? 답은 순서 없는(unordered) 요인과 순서 있는(ordered) 요인에서 다릅니다. 순서 없는 요인에 대해서는 k - 1 개의 열이 생성되는데, 이는 그 요인의 두 번째부터 k 번째 수준을 나타내는 지시자(indicator)입니다. (따라서 묵시적 모수화는 각 수준의 반응을 첫 번째 수준과 대비하는 방식이죠.) 순서 있는 요인에 대해서는 k - 1 개의 열이 1, ..., k 위의 직교 다항식(상수 항은 뺀)이 됩니다.

사실 이 답은 이미 복잡하지만, 전부는 아닙니다. 첫째, 요인 항을 포함한 모델에서 절편을 빼면 그 첫 요인 항이 k 개의 열, 즉 모든 수준의 지시자로 부호화됩니다. 둘째, contrasts 에 대한 options 설정으로 전체 동작을 바꿀 수 있어요. R의 기본 설정은

options(contrasts = c("contr.treatment", "contr.poly"))

입니다. 이걸 굳이 언급하는 주된 이유는 R과 S가 순서 없는 요인에 대해 서로 다른 기본값을 쓰기 때문인데, S는 Helmert 대비를 씁니다. 그래서 S-PLUS 를 쓴 교과서나 논문의 결과와 비교해야 한다면

options(contrasts = c("contr.helmert", "contr.poly"))

로 설정해야 해요. 이것은 의도적인 차이로, 처리 대비(treatment contrast, R의 기본값)가 초보자가 해석하기에 더 쉽다고 여겨지기 때문입니다.

이야기는 아직 끝나지 않았어요. 사용할 대비 방식은 contrasts 와 C 함수를 써서 모델의 항별로 따로 설정할 수도 있습니다.

상호작용 항을 아직 다루지 않았는데, 이들은 성분 항을 위해 도입된 열들의 곱을 생성합니다.

세부 내용은 복잡하지만, R의 모델 수식은 한계성(marginality)이 지켜지는 한 통계 전문가가 기대하는 모델을 보통 그대로 만들어 냅니다. 예컨대 주효과 없이 상호작용만 있는 모델을 적합하면 대체로 놀라운 결과가 나오며, 이는 전문가용입니다.

11.2 선형 모델(Linear models)

보통의 다중 모델을 적합하는 기본 함수는 lm() 이며, 간결한 호출 형태는 다음과 같습니다.

> fitted.model <- lm(formula, data = data.frame)

예를 들어

> fm2 <- lm(y ~ x1 + x2, data = production)

는 y 를 x1 과 x2 에 회귀시키는 다중 회귀 모델을 적합합니다(묵시적 절편 항 포함).

중요한(기술적으로는 선택적이지만) 매개변수 data = production 은 모델을 구성하는 데 필요한 변수를 production 데이터 프레임에서 먼저 가져와야 한다고 지정해요. 이는 데이터 프레임 production 이 검색 경로에 attach 되었는지와 무관하게 항상 적용됩니다.

11.3 모델 정보를 추출하는 일반 함수(Generic functions)

lm() 의 반환값은 적합된 모델 객체, 기술적으로는 클래스 "lm" 의 결과 리스트입니다. 이 적합된 모델에 대한 정보는 클래스 "lm" 객체를 겨냥하는 일반 함수(generic function)들을 통해 표시·추출·플롯 등을 할 수 있어요. 그 함수들은

add1    deviance   formula      predict  step
alias   drop1      kappa        print    summary
anova   effects    labels       proj     vcov
coef    family     plot         residuals

이렇게 있습니다. 가장 자주 쓰이는 것 몇 가지를 간단히 설명할게요.

anova(object_1, object_2): 부분 모델을 더 큰 모델과 비교해 분산분석표를 만들어 낸다.
coef(object): 회귀 계수(행렬)를 추출한다. 긴 형태는 coefficients(object).
deviance(object): 잔차 제곱합. 필요하면 가중된 값.
formula(object): 모델 수식을 추출한다.
plot(object): 잔차, 적합값, 몇 가지 진단을 보여 주는 네 개의 플롯을 만든다.
predict(object, newdata=data.frame): 주어진 데이터 프레임은 원래와 같은 라벨로 변수가 지정되어야 한다. 반환값은 data.frame 안의 결정 변수 값들에 대응하는 예측값의 벡터 또는 행렬.
print(object): 객체의 간결한 버전을 출력한다. 대부분 묵시적으로 쓰인다.
residuals(object): (행렬 형태의) 잔차를 추출한다. 필요하면 가중된 값. 짧은 형태는 resid(object).
step(object): 항을 더하거나 빼고 위계를 유지하면서 적절한 모델을 고른다. 단계적 탐색에서 발견된 AIC(Akaike's An Information Criterion) 값이 가장 작은 모델을 반환한다.
summary(object): 회귀분석 결과의 포괄적인 요약을 출력한다.
vcov(object): 적합된 모델 객체의 주 모수에 대한 분산-공분산 행렬을 반환한다.

11.4 분산분석과 모델 비교(Analysis of variance and model comparison)

모델 적합 함수 aov(formula, data=data.frame) 는 가장 단순한 수준에서는 lm() 과 매우 비슷하게 작동하며, 위 표에 나열된 일반 함수 대부분이 그대로 적용됩니다.

여기에 더해 aov() 는 분할구 실험이나 블록 간 정보를 회복하는 균형 불완전 블록 설계 같은 여러 오차 층을 가진 모델의 분석도 허용한다는 점에 주목하세요. 모델 수식

response ~ mean.formula + Error(strata.formula)

은 strata.formula 가 정의하는 오차 층을 가진 다층(multi-stratum) 실험을 지정합니다. 가장 단순한 경우 strata.formula 는 그냥 요인 하나인데, 이때는 그 요인의 수준 사이와 수준 안, 두 층의 실험을 정의해요.

예를 들어 모든 결정 변수가 요인일 때, 다음과 같은 모델 수식은

> fm <- aov(yield ~ v + n*p*k + Error(farms/blocks), data=farm.data)

평균 모델 v + n*p*k 와 'farms 사이', 'farms 안에서 blocks 사이', 'blocks 안' 이라는 세 개의 오차 층을 가진 실험을 기술하는 데 전형적으로 쓰입니다.

11.4.1 분산분석표(ANOVA tables)

또한 분산분석표는 일련의 적합된 모델들에 대한 것이라는 점에 주의하세요. 표시된 제곱합은 그 항을 수열의 그 위치에서 모델에 포함시킴으로써 생기는 잔차 제곱합의 감소분입니다. 따라서 포함 순서가 무의미해지는 것은 오직 직교(orthogonal) 실험에서뿐이에요.

다층 실험의 경우 절차는 먼저 반응을 오차 층 위로, 역시 순서대로 투영(project)한 다음 각 투영에 평균 모델을 적합하는 것입니다. 더 자세한 내용은 Chambers & Hastie (1992) 를 참고하세요.

기본 전체 ANOVA 표보다 더 유연한 대안은 anova() 함수로 두 개 이상의 모델을 직접 비교하는 것입니다.

> anova(fitted.model.1, fitted.model.2, ...)

그러면 순서대로 적합되었을 때 적합된 모델들 사이의 차이를 보여 주는 분산분석표가 표시됩니다. 물론 비교되는 적합 모델들은 보통 위계적인 수열이에요. 이 방식은 기본값과 다른 정보를 주는 게 아니라, 이해하고 제어하기 더 쉽게 만들어 줄 뿐입니다.

11.5 적합된 모델 갱신하기(Updating fitted models)

update() 함수는 대체로 편의 함수로, 보통 몇 개의 항만 더하거나 빼서 이전에 적합한 모델과 다른 모델을 적합하게 해 줍니다. 형태는

> new.model <- update(old.model, new.formula)

입니다. new.formula 에서 마침표 하나만으로 된 특별한 이름 . 은 '이전 모델 수식의 해당 부분' 을 대신할 수 있어요. 예를 들어

> fm05 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = production)
> fm6  <- update(fm05, . ~ . + x6)
> smf6 <- update(fm6, sqrt(.) ~ .)

는 (아마도) production 데이터 프레임의 변수들로 5변량 다중 회귀를 적합하고, 여섯 번째 회귀 변수를 포함한 추가 모델을 적합하며, 반응에 제곱근 변환을 적용한 변형 모델을 적합합니다.

특히 주의할 점은, 원래 모델 적합 함수 호출에서 data= 인자를 지정했다면 그 정보가 적합된 모델 객체를 통해 update() 와 그 포탁 함수들에 전달된다는 것입니다.

이름 . 은 다른 맥락에서도 쓰일 수 있는데, 뜻은 조금 달라요. 예를 들어

> fmfull <- lm(y ~ . , data = production)

는 반응 y 와 데이터 프레임 production 의 다른 모든 변수들을 회귀 변수로 하는 모델을 적합합니다.

증분적인 모델 수열을 탐색하는 다른 함수로는 add1(), drop1(), step() 이 있어요. 이름만 봐도 용도를 짐작할 수 있지만, 자세한 내용은 온라인 도움말을 확인하세요.

11.6 일반화선형모형(Generalized linear models)

일반화선형모형(GLM)은 선형 모델을 발전시켜, 정규분포가 아닌 반응 분포와 선형성으로의 변환까지 깔끔하고 직관적인 방식으로 수용합니다. 일반화선형모형은 다음과 같은 일련의 가정으로 기술할 수 있어요.

관심 대상인 반응 y 가 있고, 그 값들이 반응의 분포에 영향을 주는 자극 변수 x_1, x_2, … 가 있다.
자극 변수들은 오직 하나의 선형 함수를 통해서만 y 의 분포에 영향을 준다. 이 선형 함수를 선형 예측자(linear predictor) 라고 하며, 보통 다음과 같이 쓴다.

eta = beta_1 x_1 + beta_2 x_2 + ... + beta_p x_p,

따라서 x_i 는 beta_i 가 0일 때에만(그리고 그럴 때에만) y 의 분포에 영향이 없습니다.

y 의 분포는 다음과 같은 형태다.

f_Y(y; mu, phi)
  = exp((A/phi) * (y lambda(mu) - gamma(lambda(mu))) + tau(y, phi))

여기서 phi 는 척도 모수(scale parameter)(어쩌면 알려진 값)로 모든 관측에 대해 상수이고, A 는 사전 가중치(prior weight) 로 알려져 있다고 가정하되 관측마다 달라질 수 있으며, mu 는 y 의 평균입니다. 즉 y 의 분포는 평균과, 어쩌면 척도 모수에 의해 결정된다고 가정하는 거죠.

평균 mu 는 선형 예측자의 매끄럽고 역이 존재하는 함수다.

mu = m(eta),    eta = m^{-1}(mu) = ell(mu)

이 역함수 ell() 를 연결 함수(link function) 라고 합니다.

이 가정들은 실제 통계 실무에서 쓰이는 넓은 모델 군을 담을 만큼 느슨하면서도, 적어도 근사적으로 추정과 추론의 통일된 방법론을 개발할 수 있을 만큼 빡빡합니다. 자세한 내용은 McCullagh & Nelder (1989) 나 Dobson (1990) 같은 현재의 참고 문헌을 찾아보세요.

11.6.1 패밀리(Families)

R에 제공되는 시설로 다루는 일반화선형모형의 클래스에는 gaussian, binomial, poisson, inverse gaussian, gamma 반응 분포와, 반응 분포가 명시되지 않은 준우도(quasi-likelihood) 모델이 포함됩니다. 후자의 경우 분산 함수가 평균의 함수로 지정되어야 하지만, 그 외의 경우에는 이 함수가 반응 분포에 의해 함의됩니다.

각 반응 분포는 평균과 선형 예측자를 연결하는 다양한 연결 함수를 허용합니다. 자동으로 사용 가능한 것들은 다음 표와 같습니다.

반응 분포, 연결 함수, 그리고 모델링 작업을 수행하는 데 필요한 각종 정보의 조합을 일반화선형모형의 패밀리(family) 라고 합니다.

11.6.2 glm() 함수

반응의 분포가 오직 하나의 선형 함수를 통해 자극 변수에 의존하므로, 선형 모델에서 쓰던 것과 같은 메커니즘으로 일반화 모델의 선형 부분을 지정할 수 있어요. 차이는 패밀리를 다른 방식으로 지정해야 한다는 데 있습니다.

일반화선형모형을 적합하는 R 함수는 glm() 이며, 형태는

> fitted.model <- glm(formula, family=family.generator, data=data.frame)

입니다. 새로 등장한 것은 family.generator 뿐인데, 이것이 패밀리를 기술하는 도구예요. 모델과 추정 과정을 함께 정의하고 제어하는 함수와 표현식의 리스트를 생성하는 함수의 이름입니다. 처음 보면 조금 복잡해 보일 수 있지만, 사용법은 아주 간단해요.

표준 공급 패밀리 생성기의 이름은 Families 표의 "Family Name" 아래에 있습니다. 연결 함수에 선택지가 있을 때는 패밀리 이름에 괄호로 매개변수처럼 연결 함수의 이름을 붙일 수도 있어요. quasi 패밀리의 경우 분산 함수도 이 방식으로 지정할 수 있습니다.

몇 가지 예가 과정을 명확하게 만들어 줄 거예요.

gaussian 패밀리

다음과 같은 호출은

> fm <- glm(y ~ x1 + x2, family = gaussian, data = sales)

다음과 같은 결과를 냅니다.

> fm <- lm(y ~ x1+x2, data=sales)

다만 훨씬 비효율적이죠. gaussian 패밀리는 자동으로 연결 함수 선택지가 제공되지 않으므로 매개변수를 쓸 수 없다는 점에 주의하세요. 표준이 아닌 연결 함수가 필요한 gaussian 문제는, 나중에 보게 될 quasi 패밀리를 통해 보통 처리할 수 있어요.

binomial 패밀리

Silvey (1970) 에서 가져온 작고 인위적인 예를 들어볼게요.

에게 해의 칼리토스(Kalythos) 섬에서 남성 거주자들은 선천성 안질환을 앓고 있는데, 그 영향은 나이가 들수록 더 뚜렷해집니다. 여러 연령대의 섬 남성 표본을 뽑아 실명 여부를 검사해 기록했습니다. 데이터는 아래와 같아요.

| | Age: | 20 | 35 | 45 | 55 | 70 | | | No. tested: | 50 | 50 | 50 | 50 | 50 | | | No. blind: | 6 | 17 | 26 | 37 | 44 |

우리가 풀 문제는 이 데이터에 로짓(logistic) 모델과 프로빗(probit) 모델을 모두 적합하고, 각 모델의 LD50, 즉 남성 거주자가 실명할 확률이 50%가 되는 나이를 추정하는 것입니다.

x 연령에서 실명자 수를 y, 검사 수를 n 이라고 하면 두 모델 모두 y ~ B(n, F(beta_0 + beta_1 x)) 형태입니다. 프로빗의 경우 F(z) = Phi(z) 는 표준 정규 분포 함수이고, 로짓(기본값)의 경우 F(z) = e^z/(1+e^z) 입니다. 두 경우 모두 LD50 = - beta_0/beta_1 로, 즉 분포 함수의 인자가 0이 되는 지점입니다.

첫 단계는 데이터를 데이터 프레임으로 구성하는 것입니다.

> kalythos <- data.frame(x = c(20,35,45,55,70), n = rep(50,5),
                         y = c(6,17,26,37,44))

glm() 으로 이항 모델을 적합할 때 반응에는 세 가지 가능성이 있습니다.

반응이 벡터라면 이진(binary) 데이터를 담고 있다고 가정하므로 0/1 벡터여야 한다.
반응이 두 열 행렬이라면 첫 열은 시행의 성공 수, 둘째 열은 실패 수를 담는다고 가정한다.
반응이 요인(factor) 이라면 첫 번째 수준은 실패(0), 나머지 수준은 '성공'(1)으로 취한다.

여기서는 두 번째 관례가 필요하므로, 데이터 프레임에 행렬을 하나 더합니다.

> kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)

모델을 적합하려면

> fmp <- glm(Ymat ~ x, family = binomial(link=probit), data = kalythos)
> fml <- glm(Ymat ~ x, family = binomial, data = kalythos)

를 씁니다. 로짓 연결이 기본값이므로 두 번째 호출에서는 매개변수를 생략할 수 있어요. 각 적합의 결과를 보려면

> summary(fmp)
> summary(fml)

을 쓰면 됩니다. 두 모델 모두 (너무나) 잘 적합됩니다. LD50 추정값을 찾으려면 간단한 함수를 쓰면 돼요.

> ld50 <- function(b) -b[1]/b[2]
> ldp <- ld50(coef(fmp)); ldl <- ld50(coef(fml)); c(ldp, ldl)

이 데이터에서의 실제 추정값은 각각 43.663 년과 43.601 년입니다.

포아송 모델(Poisson models)

포아송 패밀리의 기본 연결 함수는 log 이며, 실제로 이 패밀리의 주요 용도는 실제 분포가 종종 다항(multinomial)인 빈도 데이터에 대리 포아송 로그-선형 모델을 적합하는 것입니다. 이는 크고 중요한 주제라 여기서 더 다루지는 않을게요. 사실 비가우시안 일반화 모델 사용의 상당 부분을 차지할 만큼 큽니다.

실제로 진짜 포아송 데이터가 등장하는 경우도 있는데, 과거에는 로그나 제곱근 변환 후 가우시안 데이터로 분석하곤 했습니다. 후자에 대한 우아한 대안으로, 포아송 일반화선형모형을 다음 예처럼 적합할 수 있어요.

> fmod <- glm(y ~ A + B + x, family = poisson(link=sqrt),
              data = worm.counts)

준우도 모델(Quasi-likelihood models)

모든 패밀리에서 반응의 분산은 평균에 의존하며 척도 모수를 승수로 가집니다. 분산이 평균에 의존하는 형태는 반응 분포의 특성인데, 예컨대 포아송 분포에서는 Var(y) = mu 입니다.

준우도 추정과 추론에서는 정확한 반응 분포를 지정하지 않고, 대신 연결 함수 하나와 평균에 의존하는 분산 함수의 형태만 지정합니다. 준우도 추정이 gaussian 분포와 형식적으로 동일한 기법을 쓰므로, 이 패밀리는 덕분에 표준이 아닌 연결 함수나 분산 함수로 gaussian 모델을 적합하는 방법이 되기도 해요.

예를 들어 비선형 회귀 y = theta_1 z_1 / (z_2 - theta_2) + e 를 생각해 봅시다. 이는 x_1 = z_2/z_1, x_2 = -1/z_1, beta_1 = 1/theta_1, beta_2 = theta_2/theta_1 일 때 y = 1 / (beta_1 x_1 + beta_2 x_2) + e 로 다시 쓸 수 있어요. 적절한 데이터 프레임을 만들었다고 하면, 이 비선형 회귀를 다음과 같이 적합할 수 있습니다.

> nlfit <- glm(y ~ x1 + x2 - 1,
               family = quasi(link=inverse, variance=constant),
               data = biochem)

더 자세한 정보가 필요하면 매뉴얼과 도움말 문서를 참고하세요.

11.7 비선형 최소제곱과 최대가능도 모델(Nonlinear least squares and maximum likelihood)

일부 비선형 모델 형태는 일반화선형모형(glm())으로 적합할 수 있어요. 하지만 대부분의 경우 비선형 곡선 적합 문제를 비선형 최적화 문제로 접근해야 합니다. R의 비선형 최적화 루틴은 optim(), nlm(), nlminb() 이며, 이들은 부적합도(lack-of-fit) 지표를 최소화하는 모수 값을 찾되, 다양한 모수 값을 반복적으로 시도하는 방식으로 작동합니다. 예컨대 선형 회귀와 달리 절차가 만족스러운 추정값으로 수렴한다는 보장은 없습니다. 모든 방법은 어떤 모수 값을 시도할지에 대한 초기 추정값이 필요하며, 수렴은 시작값의 품질에 결정적으로 의존할 수 있어요.

11.7.1 최소제곱(Least squares)

비선형 모델을 적합하는 한 방법은 오차 제곱합(SSE) 또는 잔차 제곱합을 최소화하는 것입니다. 이 방법은 관측 오차가 정규 분포에서 그럴듯하게 나왔을 때 타당해요.

Bates & Watts (1988) 51쪽의 예를 들어볼게요. 데이터는 다음과 같습니다.

> x <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.56, 0.56,
         1.10, 1.10)
> y <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200)

최소화할 적합 기준은

> fn <- function(p) sum((y - (p[1] * x)/(p[2] + x))^2)

입니다. 적합을 하려면 모수의 초기 추정값이 필요한데, 타당한 시작값을 찾는 한 방법은 데이터를 플롯하고 모수 값을 추측한 뒤 그 값으로 모델 곡선을 겹쳐 보는 것입니다.

> plot(x, y)
> xfit <- seq(.02, 1.1, .05)
> yfit <- 200 * xfit/(0.1 + xfit)
> lines(spline(xfit, yfit))

더 잘할 수도 있지만, 200과 0.1이라는 시작값이 충분해 보여요. 이제 적합을 수행합니다.

> out <- nlm(fn, p = c(200, 0.1), hessian = TRUE)

적합 후 out$minimum 은 SSE이고, out$estimate 는 모수의 최소제곱 추정값입니다. 추정값의 근사 표준오차(SE)를 얻으려면

> sqrt(diag(2*out$minimum/(length(y) - 2) * solve(out$hessian)))

를 합니다. 위 줄에서 뺀 2 는 모수의 개수입니다. 95% 신뢰구간은 모수 추정값 ± 1.96 SE가 되겠죠. 새 플롯에 최소제곱 적합을 겹칠 수 있습니다.

> plot(x, y)
> xfit <- seq(.02, 1.1, .05)
> yfit <- 212.68384222 * xfit/(0.06412146 + xfit)
> lines(spline(xfit, yfit))

표준 패키지 stats 는 최소제곱으로 비선형 모델을 적합하는 훨씬 광범위한 시설을 제공합니다. 방금 적합한 모델은 Michaelis-Menten 모델이므로

> df <- data.frame(x=x, y=y)
> fit <- nls(y ~ SSmicmen(x, Vm, K), df)
> fit
Nonlinear regression model
  model:  y ~ SSmicmen(x, Vm, K)
   data:  df
          Vm            K
212.68370711   0.06412123
 residual sum-of-squares:  1195.449
> summary(fit)

Formula: y ~ SSmicmen(x, Vm, K)

Parameters:
    Estimate Std. Error t value Pr(>|t|)
Vm 2.127e+02  6.947e+00  30.615 3.24e-11
K  6.412e-02  8.281e-03   7.743 1.57e-05

Residual standard error: 10.93 on 10 degrees of freedom

Correlation of Parameter Estimates:
      Vm
K 0.7651

를 쓸 수 있어요.

11.7.2 최대가능도(Maximum likelihood)

최대가능도는 오차가 정규가 아니어도 적용되는 비선형 모델 적합 방법입니다. 로그가능도를 최대화하는, 아니면 동등하게 음의 로그가능도를 최소화하는 모수 값을 찾는 방법이에요. Dobson (1990) 108–111쪽의 예를 들어볼게요. 이 예는 용량-반응 데이터에 로짓 모델을 적합하는데, 분명히 glm() 으로도 적합할 수 있습니다. 데이터는

> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113,
         1.8369, 1.8610, 1.8839)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)

이고, 최소화할 음의 로그가능도는

> fn <- function(p)
   sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))
           + log(choose(n, y)) ))

입니다. 타당한 시작값을 고르고 적합을 수행합니다.

> out <- nlm(fn, p = c(-50,20), hessian = TRUE)

적합 후 out$minimum 은 음의 로그가능도이고, out$estimate 는 모수의 최대가능도 추정값입니다. 추정값의 근사 SE 를 얻으려면

> sqrt(diag(solve(out$hessian)))

를 합니다. 95% 신뢰구간은 모수 추정값 ± 1.96 SE가 되겠죠.

11.8 표준이 아닌 모델 몇 가지(Some non-standard models)

이 장을 마무리하며, 특수한 회귀 및 데이터 분석 문제에 R이 제공하는 다른 시설 몇 가지를 간단히 소개합니다.

혼합 모델(Mixed models). 권장 패키지 nlme 는 선형 및 비선형 혼합효과 모델, 즉 일부 계수가 랜덤 효과에 해당하는 선형·비선형 회귀를 위한 lme() 와 nlme() 함수를 제공합니다. 이 함수들은 모델을 지정하는 데 수식을 많이 사용해요.
국소 근사 회귀(Local approximating regressions). loess() 함수는 국소 가중 회귀를 이용해 비모수 회귀를 적합합니다. 이런 회귀는 지저분한 데이터에서 추세를 두드러지게 하거나, 큰 데이터셋을 요약해 통찰을 얻는 데이터 축소에 유용합니다. loess 는 투영 추적 회귀 코드와 함께 표준 패키지 stats 에 들어 있어요.
강건 회귀(Robust regression). 데이터의 극단적 이상치 영향에 저항하도록 회귀 모델을 적합하는 함수가 여럿 있습니다. 권장 패키지 MASS 의 lqs 함수는 고저항 적합을 위한 최첨단 알고리즘을 제공하고, 저항은 덜하지만 통계적으로 더 효율적인 방법은 예컨대 MASS 패키지의 rlm 함수처럼 여러 패키지에서 이용 가능합니다.
가법 모델(Additive models). 이 기법은 결정 변수들 각각(보통 변수마다 하나)의 매끄러운 가법 함수들로부터 회귀 함수를 구성하는 것을 목표로 해요. acepack 패키지의 avas 와 ace, mda 패키지의 bruto 와 mars 함수는 사용자 기여 패키지에서 이 기법의 예를 보여 줍니다. 이의 확장인 일반화 가법 모델(Generalized Additive Models) 은 사용자 기여 패키지 gam 과 mgcv 에 구현되어 있어요.
나무 기반 모델(Tree-based models). 예측이나 해석을 위해 명시적인 전역 선형 모델을 찾는 대신, 나무 기반 모델은 결정 변수의 결정적 지점에서 데이터를 재귀적으로 이분하여, 최종적으로는 안은 가능한 한 동질적이고 사이는 가능한 한 이질적인 집단으로 데이터를 나누는 것을 추구합니다. 그 결과는 다른 데이터 분석 방법이 잘 주지 못하는 통찰을 종종 이끌어 냅니다.

모델은 다시 보통의 선형 모델 형태로 지정됩니다. 모델 적합 함수는 tree() 이지만, plot() 이나 text() 같은 많은 다른 일반 함수가 나무 기반 모델 적합 결과를 그래픽으로 잘 표현하도록 잘 적응되어 있어요.

나무 모델은 사용자 기여 패키지 rpart 와 tree 를 통해 R에서 이용 가능합니다.

더 알아보기

이 장에서 다룬 모델 수식과 적합 함수의 문법은 R에서 빠질 수 없는 뼈대예요. 실전에선 lm(), glm(), aov(), nls() 의 도움말을 먼저 열어보고, 안내 예제를 그대로 따라쳐 보는 걸 권해요.
일반화선형모형 이론이 궁금하다면 본문에서 인용한 McCullagh & Nelder (1989), Dobson (1990) 같은 참고 문헌이 좋은 출발점이에요.
혼합 모델·강건 회귀·가법 모델·나무 모델은 각각 nlme, MASS, acepack/mda/gam/mgcv, rpart/tree 패키지로 확장됩니다.