주성분분석

library(dmtr)
library(dplyr)

1 데이터

data(biometric, package = "dmtr")

Table 1.1: 나이, 키, 몸무게 데이터

age	height	weight
21	170	60
47	167	65
36	173	67
15	165	54
54	168	73
25	177	71
32	169	68
18	172	62
43	171	66
28	175	68

2 변수의 변동과 제곱합

총 \(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} \tag{2.1} \end{equation}\]

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

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

\[\begin{equation*} z_{ij} \leftarrow \frac{x_{ij}}{\sqrt{\frac{1}{n - 1} \sum_{l =1}^{n} x_{lj}^2}} \tag{2.2} \end{equation*}\]

이 때, 식 (2.2)에서 분모 부분은 변수의 표본 표준편차로 \(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}\]

3 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은 추출하고자 하는 주성분의 개수이다.

function (X, .pc = NULL) 
{
    if (rlang::is_empty(.pc) || (.pc > min(dim(X)))) {
        .pc <- min(dim(X))
    }
    Th <- matrix(NA, nrow = nrow(X), ncol = .pc)
    Vh <- matrix(NA, nrow = ncol(X), ncol = .pc)
    for (h in seq_len(.pc)) {
        j <- sample(ncol(X), 1L)
        Th[, h] <- X[, j]
        while (TRUE) {
            Vh[, h] <- t(t(Th[, h]) %*% X/(norm(Th[, h], "2")^2))
            Vh[, h] <- Vh[, h]/norm(Vh[, h], "2")
            th <- X %*% Vh[, h]
            if (all(dplyr::near(Th[, h], th))) 
                break
            Th[, h] <- th
        }
        X <- X - Th[, h] %*% t(Vh[, h])
    }
    return(list(T = Th, V = Vh))
}
<bytecode: 0x7f8b6d5deb70>
<environment: namespace:dmtr>

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

fit_pca(biometric, c(age, height, weight), .pc = 2L)
#> $eig
#> [1] 1.743468 1.172119
#> 
#> $score
#>              PC1        PC2
#>  [1,] -1.2672472 -0.2796131
#>  [2,]  0.2485837  1.4996458
#>  [3,]  0.6166785 -0.3344830
#>  [4,] -2.8174941  0.6097017
#>  [5,]  1.7218381  1.5057495
#>  [6,]  1.0580976 -1.7554054
#>  [7,]  0.1870461  0.3451973
#>  [8,] -0.9370744 -0.8763285
#>  [9,]  0.5953602  0.4324721
#> [10,]  0.5942115 -1.1469365
#> 
#> $loadings
#>              PC1         PC2
#> age    0.5684413  0.59068025
#> height 0.3574679 -0.80427252
#> weight 0.7410069 -0.06513491
#> 
#> $ncomp
#> [1] 2
#> 
#> $R2
#> [1] 0.5811558 0.3907064
#> 
#> $center
#>    age height weight 
#>   31.9  170.7   65.4 
#> 
#> $scale
#>       age    height    weight 
#> 12.982467  3.683296  5.541761

4 주성분 회귀분석

본 절에서 독립변수 행렬 \(\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})\))의 주성분 스코어를 독립변수로 사용한다.

4.1 모형 추정

위 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} \tag{4.1} \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
#> $betas
#> (Intercept)         PC1         PC2 
#>  65.3999999  -0.5332898   5.5093680 
#> 
#> $hessian
#>               (Intercept)          PC1           PC2
#> (Intercept)  2.000000e+01 7.216450e-12 -6.980527e-12
#> PC1          7.216450e-12 2.123307e+01  5.301315e-12
#> PC2         -6.980527e-12 5.301315e-12  1.476693e+01
#> 
#> $mse
#> [1] 7.038464
#> 
#> $df
#> [1] 7
#> 
#> $se
#> (Intercept)         PC1         PC2 
#>   0.8389556   0.8142308   0.9763576 
#> 
#> $sst
#> [1] 276.4
#> 
#> $rsq
#> [1] 0.8217466
#> 
#> $rsqadj
#> [1] 0.770817
#> 
#> $n
#> [1] 10

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

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

ttest_linear_regression(fit)
#> # A tibble: 3 × 5
#>   term        estimate std_error t_statistic  p_value
#>   <chr>          <dbl>     <dbl>       <dbl>    <dbl>
#> 1 (Intercept)   65.4       0.839      78.0   1.50e-11
#> 2 PC1            0.533     0.814       0.655 5.33e- 1
#> 3 PC2            5.51      0.976       5.64  7.80e- 4

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

fit$mse * solve(fit$hessian / 2)
#>               (Intercept)           PC1           PC2
#> (Intercept)  7.038464e-01  2.070131e-14 -1.785960e-13
#> PC1          2.070131e-14  6.629718e-01 -6.747667e-13
#> PC2         -1.785960e-13 -6.747667e-13  9.532742e-01

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

selected_pc <- select_variables_stepwise(reg_data, weight, names(T_A))
#> Inital variable:  PC2 , p-value =  0.0003803436  <  0.05 
#> Iteration  1 : forward selection - no additional variable gives statistically significant improvement of the fit.
fit_linear_regression(reg_data, weight, all_of(selected_pc))
#> $betas
#> (Intercept)         PC2 
#>   65.399987    5.509376 
#> 
#> $hessian
#>              (Intercept)          PC2
#> (Intercept) 2.000000e+01 7.244205e-12
#> PC2         7.244205e-12 1.476693e+01
#> 
#> $mse
#> [1] 6.536071
#> 
#> $df
#> [1] 8
#> 
#> $se
#> (Intercept)         PC2 
#>   0.8084597   0.9408672 
#> 
#> $sst
#> [1] 276.4
#> 
#> $rsq
#> [1] 0.8108228
#> 
#> $rsqadj
#> [1] 0.7871757
#> 
#> $n
#> [1] 10

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

fit_linear_regression(reg_data, weight, PC1)
#> $betas
#> (Intercept)         PC1 
#>  65.4000002  -0.5332894 
#> 
#> $hessian
#>               (Intercept)           PC1
#> (Intercept)  2.000000e+01 -7.771561e-12
#> PC1         -7.771561e-12  2.123307e+01
#> 
#> $mse
#> [1] 34.17259
#> 
#> $df
#> [1] 8
#> 
#> $se
#> (Intercept)         PC1 
#>    1.848583    1.794103 
#> 
#> $sst
#> [1] 276.4
#> 
#> $rsq
#> [1] 0.01092373
#> 
#> $rsqadj
#> [1] -0.1127108
#> 
#> $n
#> [1] 10

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

4.2 회귀계수 변환

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

주성분 분해 시 각 변수에 적용된 평균 조정값을 \(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
#>    age height 
#>   31.9  170.7
fit$scale
#>       age    height 
#> 12.982467  3.683296
fit$loadings
#>               PC1       PC2
#> age     0.7071068 0.7071068
#> height -0.7071068 0.7071068

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

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

betas <- (1 / fit$scale) * 
  drop(fit$loadings %*% fit$betas[colnames(fit$loadings)])
betas
#>       age    height 
#> 0.3291211 0.9552908

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

fit$org_betas
#>  (Intercept)          age       height 
#> -108.1671116    0.3291211    0.9552908

4.3 반응치의 예측

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

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

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

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

new_obs %>%
  scale(center = fit$center, scale = fit$scale) %*%
  fit$loadings
#>            PC1       PC2
#> [1,] 0.5755606 0.3067933

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

new_score <- predict_pca(fit, new_obs)
new_score
#> # A tibble: 1 × 2
#>     PC1   PC2
#>   <dbl> <dbl>
#> 1 0.576 0.307

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

predict_linear_regression(fit, new_score, starts_with("PC"))
#> # A tibble: 1 × 1
#>   .pred
#>   <dbl>
#> 1  67.4

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

predict_pcr(
  fit,
  .new_data = new_obs,
  .ci_interval = 0.95,
  .pi_interval = 0.95
)
#> # A tibble: 1 × 6
#>   .pred   .se .ci_lower .ci_upper .pi_lower .pi_upper
#>   <dbl> <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
#> 1  67.4  1.01      65.0      69.8      60.7      74.1