In youngroklee-ml/dmtr: Functions and Data for "Data Mining Techniques"

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(dmtr)
library(dplyr)

데이터 {#pca-data}

data(biometric, package = "dmtr")

knitr::kable(
  biometric, 
  booktabs = TRUE,
  align = c('r', 'r', 'r'),
  caption = '나이, 키, 몸무게 데이터'
)

변수의 변동과 제곱합 {#pca-ss}

총 $k$개의 독립변수가 있고 각 독립변수에 대하여 $n$개의 관측치가 있다고 하자. 이 때, $x_{ij}$를 $j$번째 독립변수에 대한 $i$번째 관측치라 하자. 즉, 관측데이터는 아래와 같은 행렬로 표현할 수 있다.

\begin{equation} \mathbf{X} = \left[ \begin{array}{c c c c} x_{11} & x_{12} & \cdots & x_{1k}\ x_{21} & x_{22} & \cdots & x_{2k}\ \vdots & \vdots & \ddots & \vdots \ x_{n1} & x_{n2} & \cdots & x_{nk} \end{array} \right] \end{equation}

주성분분석에서는 통상 원데이터를 그대로 사용하지 않고 적절한 변환을 취하는데, 주로 평균조정(mean-centered) 데이터를 이용한다. 이는 아래와 같이 독립변수에 대하여 표본평균을 뺌으로써 조정된 변수의 평균이 0이 되도록 하는 것이다.

\begin{equation} \tilde{x}{ij} \leftarrow x{ij} - \frac{1}{n} \sum_{l = 1}^{n} x_{lj} (#eq:pca-mean-centering) \end{equation}

이후에 별도의 언급이 없는 한, 행렬 $\mathbf{X}$ 및 변수값 $x_{ij}$는 식 \@ref(eq:pca-mean-centering)을 이용하여 평균조정된 것으로 가정한다.

이 밖에도 다른 변환이 사용되는 경우가 있는데, 특히 단위 등이 서로 상이할 경우에는 평균조정 이후 추가로 각 변수의 분산이 1이 되도록 분산조정을 한다.

\begin{equation} z_{ij} \leftarrow \frac{x_{ij}}{\sqrt{\frac{1}{n - 1} \sum_{l =1}^{n} x_{lj}^2}} (#eq:pca-scaling) \end{equation}

이 때, 식 \@ref(eq:pca-scaling)에서 분모 부분은 변수의 표본 표준편차로 $s_j$로 표현된다.

\begin{equation} s_{j} = \sqrt{\frac{1}{n - 1} \sum_{l =1}^{n} x_{lj}^2} \end{equation}

이후 분산조정을 이용하는 경우 행렬 $\mathbf{Z}$ 및 변수값 $z_{ij}$로 표현한다.

\begin{equation} \mathbf{Z} = \left[ \begin{array}{c c c c} z_{11} & z_{12} & \cdots & z_{1k}\ z_{21} & z_{22} & \cdots & z_{2k}\ \vdots & \vdots & \ddots & \vdots \ z_{n1} & z_{n2} & \cdots & z_{nk} \end{array} \right] \end{equation}

변수벡터 $\mathbf{x}j = [x{1j} \, x_{2j} \, \cdots \, x_{nj}]^\top$에 대한 제곱합의 정의는 아래와 같다.

\begin{equation} SS(\mathbf{x}j) = \mathbf{x}_j^\top \mathbf{x}_j = \sum{i = 1}^{n} x_{ij}^2 \end{equation}

NIPALS 알고리즘 {#pca-nipals}

NIPALS(Nonlinear Iterative Paritial Least Squares) 알고리즘은 반복적(iterative) 알고리즘을 이용하여 변동 기여율이 가장 큰 주성분부터 가장 작은 주성분까지 순차적으로 고유벡터와 주성분 스코어를 구하는 방법이다.

우선, 특이치 분해에서 사용한 식을 단순화하여, 분산조정된 행렬 $\mathbf{Z}$가 아래와 같이 주성분 스코어 행렬 $\mathbf{T}$와 고유벡터 행렬 $\mathbf{V}$로 분해된다고 하자. (분산조정 대신 평균조정만을 원할 경우 $\mathbf{Z}$ 대신 $\mathbf{X}$를 사용)

[ \mathbf{Z} = \mathbf{U} \mathbf{D} \mathbf{V}^\top = \mathbf{T} \mathbf{V}^\top ]

즉, 주성분 스코어 행렬 $\mathbf{T}$는 아래와 같다.

[ \mathbf{T} = \mathbf{Z} \mathbf{V} ]

NIPALS 알고리즘은 아래와 같이 주성분 스코어 행렬 $\mathbf{T}$의 열과 고유벡터행렬 $\mathbf{V}$의 열을 동시에 구한다.

• [단계 0] 반복알고리즘 수행을 위한 초기화를 한다. $h \leftarrow 1$, $\mathbf{Z}_h \leftarrow \mathbf{Z}$.
• [단계 1] 데이터 행렬 $\mathbf{Z}_h$의 임의의 열 하나를 주성분 스코어 벡터 $\mathbf{t}_h$로 선정한다.
• [단계 2] 로딩벡터를 구한다. $\mathbf{v}_h \leftarrow \mathbf{Z}_h \mathbf{t}_h \left/ \sqrt{\mathbf{t}_h^\top \mathbf{t}_h} \right.$
• [단계 3] 로딩벡터의 크기가 1이 되도록 한다. $\mathbf{v}_h \leftarrow \mathbf{v}_h \left/ \sqrt{\mathbf{v}_h^\top \mathbf{v}_h} \right.$
• [단계 4] 주성분 스코어 벡터를 로딩벡터에 기반하여 계산한다. $\mathbf{t}_h \leftarrow \mathbf{Z}_h \mathbf{v}_h$
• [단계 5] 주성분 스코어 벡터 $\mathbf{t}_h$가 수렴하였으면 [단계 6]으로 진행하고, 그렇지 않으면 [단계 2]로 돌아간다.
• [단계 6] 데이터 행렬 $\mathbf{Z}h$로부터 새로 얻어진 주성분 벡터 $\mathbf{t}_h$와 고유벡터 $\mathbf{v}_h$가 설명하는 부분을 제거하고 나머지 변동만을 담은 새로운 데이터 행렬 $\mathbf{Z}{h + 1}$을 구한다. [ \mathbf{Z}{h + 1} \leftarrow \mathbf{Z}{h} - \mathbf{t}_h \mathbf{v}_h^\top ]
• [단계 7] $h \leftarrow h + 1$로 업데이트하고, [단계 1]로 돌아간다. [단계 1] - [단계 7]의 과정을 $\mathbf{Z}$의 rank 수만큼의 주성분을 얻을 때까지 반복한다.

위 반복 알고리즘을 수행하는 함수를 아래와 같이 구성해보자. 아래 함수에서 입력변수 X는 데이터 행렬으로, 평균조정된 행렬 $\mathbf{X}$나 분산조정된 $\mathbf{Z}$ 모두 사용 가능하다. 입력변수 .pc은 추출하고자 하는 주성분의 개수이다.

print(nipals_pca)

위 함수를 이용하여 주성분 분해를 수행하는 함수 fit_pca()를 실행해보자.

fit_pca(biometric, c(age, height, weight), .pc = 2L)

주성분 회귀분석 {#pcr}

본 절에서 독립변수 행렬 $\mathbf{X}$을 평균조정 이전의 원래 독립변수 행렬이라고 하자. 종속변수 관측치 벡터 $\mathbf{y}$를 독립변수 데이터 행렬 $\mathbf{X}$로 설명하는 회귀 모형을 추정하고자 할 때, 독립변수 행렬 $\mathbf{X}$의 열벡터 간 다중공선성(multicollinearity)이 높으면 최소자승법에 의한 $\boldsymbol{\beta}$의 추정치의 분산이 커지는 문제가 있으며, $\mathbf{X}$ 행렬의 관측수보다 변수 수가 많을 때는 $\boldsymbol{\beta}$ 추정치를 구할 수 없다. 이 문제를 해결하기 위해 주성분 회귀분석(principal component regression; PCR)에서는 $\mathbf{X}$ 변동 대부분을 설명하는 $A$개 ($A \leq rank(\mathbf{X})$)의 주성분 스코어를 독립변수로 사용한다.

모형 추정

위 NIPALS 알고리즘을 통해 얻어진 첫 $A$개의 주성분으로 이루어진 주성분 스코어 행렬을 $\mathbf{T}_A$라 하자.

[ \mathbf{T}_A = \left[ \mathbf{t}_1 \, \mathbf{t}_2 \, \cdots \, \mathbf{t}_A \right] ]

\begin{equation} \mathbf{y} = q_0 + q_1 \mathbf{t}_1 + q_2 \mathbf{t}_2 + \cdots + q_A \mathbf{t}_A + \mathbf{f} (#eq:pcr-model) \end{equation}

여기서 $\mathbf{f}$는 오차항을 나타내는 벡터이며, $q_0$는 intercept, $q_1, \cdots, q_A$는 각 주성분 스코어에 해당하는 회귀계수들이다. 위 모형은 다중회귀모형으로 볼 수 있으므로, 다중회귀모형에 대한 모든 이론이 적용될 수 있다. 또한 위 모형에서 각 주성분 스코어 벡터 $\mathbf{t}_1, \ldots, \mathbf{t}_A$는 서로 선형 독립적(linearly independent)이므로, 회귀성 검정이 용이한 측면이 있다.

위 예제 데이터에서 weight를 종속변수로 하고 age와 height를 독립변수로 사용하여 주성분 회귀분석을 수행해보자. 이 때, 두 개의 주성분을 모두 회귀모형의 독립변수로 사용해보자.

y <- biometric %>% select(weight)
pc_fit <- fit_pca(biometric, c(age, height), .pc = 2L)
T_A <- as_tibble(pc_fit[["score"]])
reg_data <- y %>% bind_cols(T_A)
lm_fit <- fit_linear_regression(reg_data, weight, all_of(names(T_A)))
lm_fit

위 과정을 함수 fit_pcr()로 수행할 수 있다. 함수 fit_pcr()의 결과값은 주성분 분석 함수 fit_pca()의 결과와 회귀분석 함수 fit_linear_regression()의 결과를 동시에 포함한다.

fit <- fit_pcr(biometric, weight, c(age, height), .pc = 2L)

ttest_linear_regression(fit)

위 두 개의 주성분을 모두 이용한 회귀모형의 경우, 원래 독립변수 age와 height를 그대로 이용하여 추정한 회귀모형과 정확히 일치하는 종속변수에 대한 설명력을 지닌다. 다만, 위 주성분 회귀분석 모형의 경우 주성분 스코어가 직교하므로 회귀계수벡터 추정량의 분산-공분산행렬에서 공분산이 0임을 볼 수 있다.

fit$mse * solve(fit$hessian / 2)

위 주성분 회귀분석에 단계별 선택방법을 적용할 경우, 두 번째 주성분 스코어만을 독립변수로 사용한 회귀모형이 최종적으로 선택되며, 이 때의 수정결정계수가 두 개의 주성분을 모두 이용한 모형보다 높음을 확인할 수 있다.

selected_pc <- select_variables_stepwise(reg_data, weight, names(T_A))
fit_linear_regression(reg_data, weight, all_of(selected_pc))

반면, 첫 번째 주성분 스코어만을 이용한 모형의 경우, 종속변수를 거의 설명하지 못하며, 수정결정계수는 0보다도 작음을 아래에서 확인할 수 있다. 이는 독립변수의 분산을 더 많이 설명하는 잠재변수가 반드시 종속변수의 분산 또한 더 많이 설명하지는 않는다는 것을 보여준다.

fit_linear_regression(reg_data, weight, PC1)

회귀모형에 적합한 잠재변수를 추출하는 방법으로 부분최소자승 회귀분석이 있으며, 이는 별도의 문서에서 설명하도록 한다.

회귀계수 변환

원래 변수를 독립변수로 대한 회귀계수를 주성분 회귀모형으로부터 계산할 수 있다. 이는 원래 변수와 종속변수 간의 관계를 살펴보는 데 유용하다.

주성분 분해 시 각 변수에 적용된 평균 조정값을 $m_j$, 분산 조정값을 $s_j$라 할 때, 추정된 주성분 회귀모형은 아래와 같이 원래 변수를 독립변수로 하는 회귀모형으로 선형변환할 수 있다. 여기서 $v_{ja}$는 $a$번째 주성분에 대한 $j$번째 변수의 가중치를 나타낸다.

\begin{eqnarray} \hat{y}i &=& \hat{q}_0 + \sum{a = 1}^{A} \hat{q}a t{ia}\ &=& \hat{q}0 + \sum{a = 1}^{A} \hat{q}a \sum{j = 1}^{k} v_{ja} \frac{x_{ij} - m_j}{s_j}\ &=& \hat{q}0 - \sum{a = 1}^{A} \hat{q}a \sum{j = 1}^{k} v_{ja} \frac{m_j}{s_j} + \sum_{a = 1}^{A} \sum_{j = 1}^{k} \hat{q}a \frac{v{ja}}{s_j} x_{ij}\ &=& \hat{q}0 - \sum{a = 1}^{A} \sum_{j = 1}^{k} \hat{q}a v{ja} \frac{m_j}{s_j} + \sum_{j = 1}^{k} \sum_{a = 1}^{A} \hat{q}a \frac{v{ja}}{s_j} x_{ij}\ &=& \hat{\beta}0 + \sum{j = 1}^{k} \hat{\beta}j x{ij} \end{eqnarray}

여기에서, 아래와 같은 관계식이 성립한다.

[ \hat{\beta}0 = \hat{q}_0 - \sum{a = 1}^{A} \sum_{j = 1}^{k} \hat{q}a v{ja} \frac{m_j}{s_j} ]

[ \hat{\beta}j = \sum{a = 1}^{A} \hat{q}a \frac{v{ja}}{s_j} ]

위에서 수행한 함수 fit_pcr()의 결과값 fit은 주성분 분해 시 적용된 평균조정과 분산조정값 및 가중치 $v_{ja}$를 아래와 같이 center, scale, loadings원소에 각각 저장하고 있다.

fit$center
fit$scale
fit$loadings

따라서, 위 원소값들과 추정된 주성분 회귀계수 betas를 사용하여 아래와 같이 원래 변수를 독립변수로 하는 회귀모형의 회귀계수 추정값을 계산할 수 있다.

intercept <- fit$betas["(Intercept)"] -
  drop((t(fit$center / fit$scale) %*% fit$loadings) %*%
     fit$betas[colnames(fit$loadings)])
intercept

betas <- (1 / fit$scale) * 
  drop(fit$loadings %*% fit$betas[colnames(fit$loadings)])
betas

위 계산은 함수 fit_pcr() 내에서 수행되어 결과값의 org_betas 원소에 저장된다.

fit$org_betas

반응치의 예측

새로 관측된 독립변수에 대한 종속변수 예측을 수행하고자 할 때, 새 관측치에 대한 주성분 스코어를 구한 뒤 기존에 추정된 주성분 회귀모형에 대입하는 방식으로 수행한다.

위 나이/키/몸무게 데이터로 이루어진 회귀모형 예제에서, 나이가 40세이고 키가 170cm이며, 몸무게는 측정이 되지 않은 새로운 관측치가 있다고 하자.

new_obs <- tibble(age = 40, height = 170)

새 관측치에 대한 주성분 스코어를 계산은, 주성분 분해 시 적용된 평균조정 및 분산조정을 그대로 적용한 뒤 가중치 행렬을 적용하는 방식으로 수행한다.

new_obs %>%
  scale(center = fit$center, scale = fit$scale) %*%
  fit$loadings

위 과정은 함수 predict_pca()로 구현되어 있다.

new_score <- predict_pca(fit, new_obs)
new_score

이 주성분 스코어를 함수 predict_linear_regression()의 입력 데이터로 사용할 때, 새 관측치에 대한 종속변수의 예측값을 얻을 수 있다.

predict_linear_regression(fit, new_score, starts_with("PC"))

위 일련의 예측 수행 과정이 함수 predict_pcr()로 구현되어, 아래와 같이 사용할 수 있다.

predict_pcr(
  fit,
  .new_data = new_obs,
  .ci_interval = 0.95,
  .pi_interval = 0.95
)

youngroklee-ml/dmtr documentation built on June 12, 2022, 6:24 p.m.