Chapter 4 부분최소자승법
회귀분석에서와 같이 하나의 종속변수에 영향을 주는 \(k\)개의 독립변수가 있다고 하자. 모든 변수는 평균조정되었다고 간주한다. 본 장에서 다루고자 하는 부분최소자승법(partial least squares: PLS)는 앞에서 다룬 주성분 회귀분석(PCR)과 유사하나, 도출되는 새로운 잠재변수들이 다르다.
독립변수와 종속변수간의 관계를 설명하기 위해, 우선 독립변수 행렬과 종속변수 행렬(또는 벡터)가 각각 서로 다른 잠재변수들에 의해 설명된다고 가정한 뒤, 두 잠재변수들간의 관계에 대한 모형을 세운다. 이 때, 본 장에서는 두 잠재변수들간의 관계가 선형인 모형(선형 PLS)만을 살펴본다.
4.1 필요 R 패키지 설치
본 장에서 필요한 R 패키지들은 아래와 같다.
| package | version |
|---|---|
| tidyverse | 1.3.1 |
| pls | 2.7-3 |
4.2 하나의 종속변수의 경우
4.2.1 기본 R 스크립트
앞 장의 주성분 회귀분석에서 사용했던 데이터에 대해 부분최소자승 회귀분석을 수행해보도록 하자.
train_df <- tribble(
~x1, ~x2, ~x3, ~y,
-3, -3, 5, -30,
-2, -3, 7, -20,
0, 0, 4, 0,
1, 2, 0, 5,
2, 2, -5, 10,
2, 2, -11, 35
)
knitr::kable(
train_df, booktabs = TRUE,
align = rep("r", ncol(train_df)),
caption = "부분최소자승 회귀분석 예제 데이터"
)| x1 | x2 | x3 | y |
|---|---|---|---|
| -3 | -3 | 5 | -30 |
| -2 | -3 | 7 | -20 |
| 0 | 0 | 4 | 0 |
| 1 | 2 | 0 | 5 |
| 2 | 2 | -5 | 10 |
| 2 | 2 | -11 | 35 |
R 패키지 pls 내의 함수 plsr()을 이용하여 PLS 모형을 아래와 같이 추정할 수 있다.
plsr_fit <- pls::plsr(y ~ x1 + x2 + x3, data = train_df, ncomp = 2)
coef(plsr_fit)## , , 2 comps
##
## y
## x1 2.475395
## x2 2.523238
## x3 -1.704636수행결과 object에 summary() 함수를 사용하여 학습된 모형의 독립변수 \(\mathbf{X}\) 및 종속변수 \(\mathbf{y}\)의 총변동에 대한 기여율을 확인할 수 있다.
summary(plsr_fit)## Data: X dimension: 6 3
## Y dimension: 6 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
## 1 comps 2 comps
## X 94.97 99.78
## y 88.28 91.46위 요약표는 하나의 잠재변수와 두 개의 잠재변수를 이용하였을 때 추정된 회귀모형들이 종속변수의 총 변량을 각각 88.2814529%와 91.4578959% 만큼을 설명함을 알려준다. 이는 앞 장에서 살펴보았던 주성분 회귀모형보다 더 높은 수치이다.
4.2.2 PLS 모형
종속변수가 하나만 존재하는 경우에는 데이터 행렬 \(\mathbf{X}\)와 종속변수벡터 \(\mathbf{y}\)가 동일한 잠재변수로 설명된다고 가정할 수 있다. (\(n \times k\)) 데이터 행렬 \(\mathbf{X}\)와 종속변수벡터 \(\mathbf{y}\)에 대하여 동시에 \(A\)개의 잠재변수벡터 \(\mathbf{t}_1, \cdots, \mathbf{t}_A\)로 설명하는 모형을 아래와 같이 기술해보자.
\[\begin{eqnarray} \mathbf{X} &=& \mathbf{t}_1 \mathbf{p}_1^\top + \mathbf{t}_2 \mathbf{p}_2^\top + \cdots + \mathbf{t}_A \mathbf{p}_A^\top + \mathbf{E} \tag{4.1}\\ \mathbf{y} &=& \mathbf{t}_1 b_1 + \mathbf{t}_2 b_2 + \cdots + \mathbf{t}_A b_A + \mathbf{f} \tag{4.2} \end{eqnarray}\]
여기서 계수벡터 \(\mathbf{p}_a\)는 \(\mathbf{X}\)에 해당하는 로딩(loading)을, 그리고 계수 \(b_a\)는 \(\mathbf{y}\)에 해당하는 로딩을 나타내며, \(\mathbf{E}\)와 \(\mathbf{f}\)는 각 모형에 해당하는 오차항(행렬 또는 벡터)이다.
위 모형을 (\(n \times A\)) 잠재변수 행렬 \(\mathbf{T} = \left[\mathbf{t}_1 \, \cdots \mathbf{t}_A \right]\)와 (\(k \times A\)) 로딩행렬 \(\mathbf{P} = \left[\mathbf{p}_1 \, \cdots \mathbf{p}_A \right]\), 그리고 로딩벡터 \(\mathbf{b} = \left[b_1 \, \cdots b_A \right]^\top\) 을 이용하여 아래와 같이 행렬식으로 나타낼 수 있다.
\[\begin{eqnarray} \mathbf{X} &=& \mathbf{T}\mathbf{P}^\top + \mathbf{E} \tag{4.3}\\ \mathbf{y} &=& \mathbf{T}\mathbf{b} + \mathbf{f} \tag{4.4} \end{eqnarray}\]
4.2.3 NIPALS 알고리즘
- [단계 0] 반복알고리즘 수행을 위한 초기화를 한다. \(a \leftarrow 1\), \(\mathbf{X}_a \leftarrow \mathbf{X}\), \(\mathbf{y}_a \leftarrow \mathbf{y}\).
- [단계 1] \(\mathbf{X}_a\)을 다중종속변수 행렬으로, \(\mathbf{y}_a\)를 독립변수 벡터로 하는 회귀모형으로부터 기울기 \(\mathbf{w}_a = [w_{a1} \, \cdots \, w_{ak}]^\top\)를 산출한다. \[\mathbf{w}_a \leftarrow \left. \mathbf{X}_a^\top \mathbf{y}_a \middle/ \mathbf{y}_a^\top \mathbf{y}_a \right. \]
- [단계 2] 기울기 벡터 \(\mathbf{w}_a\)의 크기가 1이 되도록 한다. \[\left. \mathbf{w}_a \leftarrow \mathbf{w}_a \middle/ \sqrt{\mathbf{w}_a^\top \mathbf{w}_a} \right.\]
- [단계 3] 잠재변수 \(\mathbf{t}_a\)를 행렬 \(\mathbf{X}_a\)의 각 열의 가중평균으로 구한다. 이 때, 가중치는 기울기 벡터 \(\mathbf{w}_a\)를 이용한다. \[\mathbf{t}_a \leftarrow \mathbf{X}_a \mathbf{w}_a\]
- [단계 4] 식 (4.1)와 같이 \(\mathbf{X}_a\)을 다중종속변수 행렬으로, \(\mathbf{t}_a\)를 독립변수 벡터로 하는 회귀모형으로부터 로딩벡터 \(\mathbf{p}_a\)를 구한다. \[\mathbf{p}_a \leftarrow \left. \mathbf{X}_a^\top \mathbf{t}_a \middle/ \mathbf{t}_a^\top \mathbf{t}_a \right.\]
- [단계 5] 로딩벡터 \(\mathbf{p}_a\)의 크기를 1로 조정하고, 잠재변수 벡터 \(\mathbf{t}_a\)와 기울기 벡터 \(\mathbf{w}_a\)의 크기를 그에 따라 보정한다. \[d \leftarrow \sqrt{\mathbf{p}_a^\top \mathbf{p}_a}, \, \mathbf{t}_a \leftarrow \mathbf{t}_a d, \, \mathbf{w}_a \leftarrow \mathbf{w}_a d, \, \mathbf{p}_a \leftarrow \frac{1}{d} \mathbf{p}_a \]
- [단계 6] 식 (4.2)와 같이 잠재변수 \(\mathbf{t}_a\)를 종속변수 \(\mathbf{y}_a\)에 회귀시킬 때 계수 \(b_a\)를 산출한다. \[b_a \leftarrow \left. \mathbf{y}_a^\top \mathbf{t}_a \middle/ \mathbf{t}_a^\top \mathbf{t}_a \right. \]
- [단계 7] 독립변수 행렬 \(\mathbf{X}_a\)와 종속변수벡터 \(\mathbf{y}_a\)로부터 새로 얻어진 잠재변수 벡터 \(\mathbf{t}_a\)가 설명하는 부분을 제거하고 나머지 변동만을 담은 독립변수 행렬 \(\mathbf{X}_{a + 1}\)과 종속변수벡터 \(\mathbf{y}_{a + 1}\)을 구한다. \[\mathbf{X}_{a + 1} \leftarrow \mathbf{X}_a - \mathbf{t}_a \mathbf{p}_a^\top, \, \mathbf{y}_{a + 1} \leftarrow \mathbf{y}_a - \mathbf{t}_a b_a\]
- [단계 8] \(a \leftarrow a + 1\)로 업데이트하고, [단계 1]로 돌아간다. [단계 1] - [단계 8]의 과정을 \(A\)개의 잠재변수를 얻을 때까지 반복한다.
위 NIPALS 알고리즘을 아래 nipals_plsr이라는 함수로 구현해보자. 이 때, 함수의 입력변수는 아래와 같다.
X: 평균조정된 (\(n \times k\)) 행렬y: 평균조정된 종속변수 벡터A: 잠재변수 개수
nipals_plsr <- function(X, y, A = NULL) {
if (is_empty(A) || (A > min(dim(X)))) {
A <- min(dim(X))
}
Tmat <- matrix(NA, nrow = nrow(X), ncol = A)
Wmat <- matrix(NA, nrow = ncol(X), ncol = A)
Pmat <- matrix(NA, nrow = ncol(X), ncol = A)
b <- vector("numeric", length = A)
for (a in seq_len(A)) {
# 단계 1
Wmat[, a] <- coef(lm(X ~ -1 + y))
# 단계 2
Wmat[, a] <- Wmat[, a] / sqrt(sum(Wmat[, a] ^ 2))
# 단계 3
Tmat[, a] <- X %*% Wmat[, a]
# 단계 4
Pmat[, a] <- coef(lm(X ~ -1 + Tmat[, a]))
# 단계 5
p_size <- sqrt(sum(Pmat[, a] ^ 2))
Tmat[, a] <- Tmat[, a] * p_size
Wmat[, a] <- Wmat[, a] * p_size
Pmat[, a] <- Pmat[, a] / p_size
# 단계 6
b[a] <- coef(lm(y ~ -1 + Tmat[, a]))
# 단계 7
X <- X - Tmat[, a] %*% t(Pmat[, a])
y <- y - Tmat[, a] %*% t(b[a])
}
return(list(T = Tmat, W = Wmat, P = Pmat, b = b))
}
X <- as.matrix(train_df[, c("x1", "x2", "x3")])
y <- train_df$y
nipals_fit <- nipals_plsr(X, y, A = 2)
nipals_fit## $T
## [,1] [,2]
## [1,] -6.3243200 -2.0380766
## [2,] -7.8584372 -0.5965965
## [3,] -3.6317877 1.5429616
## [4,] 0.9079469 1.9745198
## [5,] 5.7294582 0.7158171
## [6,] 11.1771398 -1.5986253
##
## $W
## [,1] [,2]
## [1,] 0.2817766 0.6579688
## [2,] 0.3130851 0.6388931
## [3,] -0.9079469 0.4245052
##
## $P
## [,1] [,2]
## [1,] 0.2537679 0.5424154
## [2,] 0.2858180 0.7279599
## [3,] -0.9240724 0.4193565
##
## $b
## [1] 2.924570 2.464658식 (4.3)과 (4.4)에서, 잠재변수 행렬 \(\mathbf{T}\)가 주어졌다 가정할 때 로딩행렬 및 벡터 \(\mathbf{P}\)와 \(\mathbf{b}\)는 아래와 같이 추정된다.
\[\begin{eqnarray} \hat{\mathbf{P}}^\top = \left(\mathbf{T}^\top \mathbf{T}\right)^{-1} \mathbf{T}^\top \mathbf{X} \tag{4.5}\\ \hat{\mathbf{b}} = \left(\mathbf{T}^\top \mathbf{T}\right)^{-1} \mathbf{T}^\top \mathbf{y} \tag{4.6} \end{eqnarray}\]
위 NIPALS 알고리즘 수행 결과에서 이를 확인해보자.
P_hat <- t(solve(t(nipals_fit$T) %*% nipals_fit$T) %*% t(nipals_fit$T) %*% X)
all(near(nipals_fit$P, P_hat))## [1] TRUEb_hat <- as.vector(t(solve(t(nipals_fit$T) %*% nipals_fit$T) %*%
t(nipals_fit$T) %*% as.matrix(y, ncol = 1)))
all(near(nipals_fit$b, b_hat))## [1] TRUE4.2.4 회귀식 변환
위 NIPALS 알고리즘 수행 결과를 원래 독립변수 \(\mathbf{X}\)와 종속변수 \(\mathbf{y}\)에 대한 식으로 변환하는 방법은 아래와 같다.
잠재변수행렬 \(\mathbf{T}\)는 아래와 같이 독립변수 행렬 \(\mathbf{X}\)와 가중치행렬 \(\mathbf{W}\), 그리고 로딩행렬 \(\mathbf{P}\)의 연산으로 표현된다.
\[\begin{equation} \mathbf{T} = \mathbf{X} \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \tag{4.7} \end{equation}\]
이를 식 (4.4)에 대입하면,
\[\begin{equation} \begin{split} \mathbf{y} &= \mathbf{X} \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \mathbf{b} + \mathbf{f}\\ &= \mathbf{X} \boldsymbol{\beta}_{PLS} + \mathbf{f} \end{split} \tag{4.8} \end{equation}\]
따라서, 원 독립변수 행렬 \(\mathbf{X}\)에 대한 회귀계수는 아래와 같이 정리된다.
\[\begin{equation} \boldsymbol{\beta}_{PLS} = \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \mathbf{b} \end{equation}\]
beta_pls <- nipals_fit$W %*%
solve(t(nipals_fit$P) %*% nipals_fit$W) %*%
as.matrix(nipals_fit$b, ncol = 1L)
beta_pls## [,1]
## [1,] 2.475395
## [2,] 2.523238
## [3,] -1.7046364.2.5 제곱합 분해
\(A\)개의 잠재변수를 사용하는 모형에 대하여 모형이 설명하는 \(\mathbf{y}\)의 변동(제곱합)을 \({SSR}\), 설명하지 못하는 변동을 \({SSE}\)라 할 때, \(\mathbf{y}\)의 전체제곱합(\({SST}\))은 다음과 같이 분해된다.
\[{SST} = {SS}(\mathbf{y}) = {SSR} + {SSE}\]
여기서 \({SS}()\)는 제곱합 함수로, 임의의 벡터 \(\mathbf{x}\)에 대해 아래와 같이 정의된다.
\[ {SS}(\mathbf{x}) = \mathbf{x}^\top \mathbf{x} \]
이 때, \({SSR}\)은 다음과 같이 산출할 수 있다.
\[\begin{equation} \begin{split} SSR &= \sum_{a = 1}^{A} SS(b_a \mathbf{t}_a)\\ &= \sum_{a = 1}^{A} b_a^2 SS(\mathbf{t}_a) \end{split} \tag{4.9} \end{equation}\]
\(a\)번째 잠재변수 \(\mathbf{t}_a\)가 \(\mathbf{y}\)를 설명하는 회귀제곱합을 \(SSR_a = b_a^2 SS(\mathbf{t}_a)\)라 할 때, \(SSR\)은 아래와 같이 분해된다.
\[ SSR = \sum_{a = 1}^{A} SSR_a \]
위 예제에서 2개의 잠재변수가 설명하는 \(\mathbf{y}\)의 총변동을 PLS 결과를 이용하여 계산하면 아래와 같다.
SSR_a <- nipals_fit$b ^ 2 * diag(t(nipals_fit$T) %*% nipals_fit$T)
SSR_a## [1] 2339.45850 84.17574이 때, 각 주성분이 설명하는 \(\mathbf{y}\) 변동의 기여율을 아래와 같이 정의한다.
\[\begin{equation} \Delta R_a^2 = \frac{SSR_a}{SST} \tag{4.10} \end{equation}\]
앞 예제에서 각 잠재변수의 \(\mathbf{y}\) 변동에 대한 기여율을 계산해보자.
SST <- sum(y ^ 2)
delta_Rsq <- SSR_a / SST
delta_Rsq## [1] 0.88281453 0.03176443잠재변수 \(A\)개를 이용한 PLS 모형이 설명하는 \(\mathbf{y}\)의 총 변동에 대한 기여도(\(SSR / SST\))은 아래와 같이 각 잠재변수의 기여도의 합으로 산출된다.
\[ R^2 = \frac{SSR}{SST} = \frac{\sum_{a = 1}^{A} SSR_a}{SST} = \sum_{a = 1}^{A} \Delta R_a^2 \]
따라서, 잠재변수 \(A\)개를 이용한 PLS 모형이 설명하는 \(\mathbf{y}\)의 총 변동에 대한 기여도(\(SSR / SST\))은 아래와 같이 각 잠재변수의 기여도의 합으로 산출된다.
앞 예제에서 잠재변수 2개를 이용한 최종모형이 설명하는 \(\mathbf{y}\)의 변동은 아래와 같다.
sum(delta_Rsq)## [1] 0.914579이는 앞 4.2.1절에서 R 패키지 pls를 이용하여 얻어진 결과와 동일함을 확인할 수 있다.
한편, 잠재변수들이 독립변수행렬 \(\mathbf{X}\)의 변동을 얼마나 설명하는지 동시에 검토할 필요가 있다. 각 잠재변수들의 제곱합 \(SS(\mathbf{t}_a)\)의 \(\mathbf{X}\)의 총변동 (\(SS(\mathbf{X})\))에 대한 비율이 그 기여율을 설명한다.
앞 예제에서 각각의 잠재변수의 \(\mathbf{X}\)에 대한 기여율은 아래와 같다.
diag(t(nipals_fit$T) %*% nipals_fit$T) / sum(diag(t(X) %*% X))## [1] 0.94972728 0.04811507잠재변수 2개를 이용한 PLS 모형의 \(\mathbf{X}\)에 대한 기여율은 아래와 같다.
sum(diag(t(nipals_fit$T) %*% nipals_fit$T)) / sum(diag(t(X) %*% X))## [1] 0.9978423잠재변수 2개가 독립변수행렬의 대부분의 변동을 설명함을 알 수 있으며, 위 결과는 역시 앞 4.2.1절에서 R 패키지 pls를 이용하여 얻어진 결과와 동일함을 확인할 수 있다.
4.2.6 독립변수의 중요도
원래의 각 독립변수가 종속변수를 설명하는 데 얼마나 영향을 미치는지는 공정분석 등에서 매우 중요하다. 식 (4.7)의 \(\mathbf{T}\)와 \(\mathbf{X}\)간의 관계식을 아래와 같이 다시 표현해보자.
\[\begin{equation} \begin{split} \mathbf{T} &= \mathbf{X} \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1}\\ &= \mathbf{X} \mathbf{W}^{*} \end{split} \end{equation}\]
이 때, \(\mathbf{W}^{*} = \left[\mathbf{w}^{*}_1 \, \cdots \, \mathbf{w}^{*}_A \right]\)를 변환가중치행렬이라 한다.
\[\begin{equation} \mathbf{W}^{*} = \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \end{equation}\]
이 때, 각 잠재변수가 설명하는 \(\mathbf{y}\)의 변동과 각 독립변수가 각 잠재변수에 기여하는 가중치를 고려하여, \(j\)번째 독립변수의 종속변수에 대한 중요도 척도로 \(VIP\)(variable importance in projection)를 다음과 같이 정의한다.
\[\begin{equation} VIP_j = \sqrt{\frac{k}{SSR} \sum_{a = 1}^{A} SSR_a \left( w^{*}_{aj} \middle/ \| \mathbf{w}^{*}_a \| \right)^2} \tag{4.11} \end{equation}\]
위 정의에 의하면 다음이 성립한다.
\[ \sum_{j = 1}^{k} VIP_j^2 = k \]
즉, 독립변수당 중요도 제곱의 평균은 1이 된다. 이에 따라, 통상 \(VIP\)가 1보다 큰 독립변수를 의미있는 변수로 간주한다.
앞 예제에 대해 각 변수의 중요도를 계산해보자.
k <- ncol(X)
Wx <- nipals_fit$W %*% solve(t(nipals_fit$P) %*% nipals_fit$W)
VIP <- sqrt(
colSums(
k / sum(SSR_a) * SSR_a *
(t(Wx ^ 2) /diag(t(Wx) %*% Wx))
)
)
VIP## [1] 0.5246196 0.5715532 1.5485804즉, \(x_3\)가 가장 영향력있는 변수라 할 수 있겠다.
4.3 다수의 종속변수의 경우
\(m\)개의 종속변수가 존재하여, 종속변수 데이터가 벡터가 아닌 (\(n \times m\)) 행렬
\[\mathbf{Y} = \left[ \mathbf{y}_1 \, \cdots \, \mathbf{y}_m \right]\]
으로 표현될 때, 각각의 종속변수에 대해 따로 잠재변수를 산출하기보다는, 여러 종속변수를 설명하는 공통의 잠재변수행렬 \(\mathbf{T}\)를 산출하는 것이 합리적이라 할 수 있다.
4.3.1 기본 R 스크립트
4개의 독립변수 및 2개의 종속변수에 대한 평균조정된 데이터가 다음과 같다.
train_df <- tribble(
~x1, ~x2, ~x3, ~x4, ~y1, ~y2,
-1, -0.5, -1, 1, 5.9, -10,
1, 1.1, -6, -6, -3.7, -2,
0, 0.3, -5, -2, 1, 11,
-3, -3.2, -9, 19, 7.7, -22,
4, 1.2, 14, -12, -7.5, 4,
-2, -2.6, -2, 9, 2.8, 1,
1, 3.7, 9, -9, -6.2, 18
)
knitr::kable(
train_df, booktabs = TRUE,
align = rep("r", ncol(train_df)),
caption = "다수의 종속변수에 대한 부분최소자승 회귀분석 예제 데이터"
)| x1 | x2 | x3 | x4 | y1 | y2 |
|---|---|---|---|---|---|
| -1 | -0.5 | -1 | 1 | 5.9 | -10 |
| 1 | 1.1 | -6 | -6 | -3.7 | -2 |
| 0 | 0.3 | -5 | -2 | 1.0 | 11 |
| -3 | -3.2 | -9 | 19 | 7.7 | -22 |
| 4 | 1.2 | 14 | -12 | -7.5 | 4 |
| -2 | -2.6 | -2 | 9 | 2.8 | 1 |
| 1 | 3.7 | 9 | -9 | -6.2 | 18 |
앞서 하나의 종속변수를 다루는 경우와 마찬가지로, R 패키지 pls 내의 함수 plsr()을 이용하여 PLS 모형을 추정할 수 있다. 이 때, formula 입력 시 종속변수의 행렬을 이용한다.
X <- as.matrix(train_df[, c("x1", "x2", "x3", "x4")])
Y <- as.matrix(train_df[, c("y1", "y2")])
plsr_multi_fit <- pls::plsr(Y ~ X, ncomp = 3)함수 수행 결과 추정된 PLS 모형으로부터 원 독립변수 \(x_1, \cdots, x_4\)와 종속변수 \(y_1, y_2\)간의 선형관계를 함수 coef()를 적용하여 아래와 같이 얻을 수 있다.
coef(plsr_multi_fit)## , , 3 comps
##
## y1 y2
## x1 -0.26509645 -2.8773570
## x2 0.08864959 2.7249194
## x3 -0.14258488 0.3377420
## x4 0.37330470 -0.7406377또한 summary() 함수를 적용하여 잠재변수들이 독립변수행렬 및 각 종속변수의 총변동에 기여하는 비율을 확인할 수 있다.
summary(plsr_multi_fit)## Data: X dimension: 7 4
## Y dimension: 7 2
## Fit method: kernelpls
## Number of components considered: 3
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps
## X 86.59 99.00 99.81
## y1 81.71 81.89 82.19
## y2 53.98 56.50 65.994.3.2 PLS 모형
앞 4.2.2절의 모형을 일반화하여 아래와 같은 모형을 가정한다.
\[\begin{eqnarray} \mathbf{X} &=& \mathbf{T} \mathbf{P}^\top + \mathbf{E} \tag{4.12}\\ \mathbf{Y} &=& \mathbf{U} \mathbf{Q}^\top + \mathbf{F} \tag{4.13}\\ \mathbf{U} &=& \mathbf{T} \mathbf{B} + \mathbf{H} \tag{4.14} \end{eqnarray}\]
식 (4.12)의 모형은 앞서 하나의 종속변수의 경우에서 살펴본 모형식 (4.3)와 동일하다. 식 (4.13)에서 (\(n \times A\)) 행렬 \(\mathbf{U}\)는 \(\mathbf{Y}\)를 설명하는 \(A\)개의 잠재변수를 나타내는 행렬이며, (\(m \times A\)) 행렬 \(\mathbf{Q}\)는 종속변수행렬 \(\mathbf{Y}\)와 잠재변수행렬 \(\mathbf{U}\)간의 선형관계를 나타내는 로딩행렬이다. 또한 식 (4.14)는 잠재변수행렬 \(\mathbf{T}\)와 \(\mathbf{U}\)간의 선형관계를 나타내는데, 특히 \(\mathbf{B}\)는 (\(A \times A\)) 대각행렬(diagonal matrix)로써, \(\mathbf{U}\)와 \(\mathbf{T}\)간에는 서로 대응하는 열 간에만 관계가 성립하며, 그 관계는 아래와 같다.
\[\begin{equation*} \mathbf{u}_a = b_a \mathbf{t}_a + \mathbf{h}_a, \, a = 1, \cdots, A \end{equation*}\]
이 때, \(b_a\)는 행렬 \(\mathbf{B}\)의 \(a\)번째 대각 원소를 나타낸다.
\[\mathbf{B} = \left[ \begin{array}{c c c c} b_{1} & 0 & \cdots & 0\\ 0 & b_{2} & & 0\\ \vdots & & \ddots & \vdots \\ 0 & 0 & \cdots & b_{A} \end{array} \right] \]
행렬 \(\mathbf{E}\), \(\mathbf{F}\), \(\mathbf{H}\)는 오차항에 해당하는 행렬이다.
4.3.3 NIPALS 알고리즘
다수의 종속변수가 존재하는 경우에도 NIPALS 알고리즘을 이용하여 모형을 추정한다. 이때는 각 잠재변수 \(\mathbf{t}_a\)를 추출할 때 추출한 잠재변수의 수렴 여부를 확인할 필요가 없었던 위 4.2.3절의 경우와는 달리, 각 잠재변수 \(\mathbf{t}_a\)와 \(\mathbf{u}_a\)를 추출하는 과정에서 반복적인(iterative) 기법으로 두 잠재변수 벡터들을 업데이트하며 수렴 여부를 확인하여야 한다.
- [단계 0] 반복알고리즘 수행을 위한 초기화를 한다. \(a \leftarrow 1\), \(\mathbf{X}_a \leftarrow \mathbf{X}\), \(\mathbf{Y}_a \leftarrow \mathbf{Y}\).
- [단계 1] 종속변수 행렬 \(\mathbf{Y}_a\)의 임의의 열 하나를 잠재변수 벡터 \(\mathbf{u}_a\)로 선정한다.
- [단계 2] \(\mathbf{X}_a\)을 다중종속변수 행렬으로, 잠재변수 \(\mathbf{u}_a\)를 독립변수 벡터로 하는 회귀모형으로부터 기울기 \(\mathbf{w}_a = [w_{a1} \, \cdots \, w_{ak}]^\top\)를 산출한다. \[\mathbf{w}_a \leftarrow \left. \mathbf{X}_a^\top \mathbf{u}_a \middle/ \mathbf{u}_a^\top \mathbf{u}_a \right. \]
- [단계 3] 기울기 벡터 \(\mathbf{w}_a\)의 크기가 1이 되도록 한다. \[\left. \mathbf{w}_a \leftarrow \mathbf{w}_a \middle/ \sqrt{\mathbf{w}_a^\top \mathbf{w}_a} \right.\]
- [단계 4] 잠재변수 \(\mathbf{t}_a\)를 행렬 \(\mathbf{X}_a\)의 각 열의 가중평균으로 구한다. 이 때, 가중치는 기울기 벡터 \(\mathbf{w}_a\)를 이용한다. \[\mathbf{t}_a \leftarrow \mathbf{X}_a \mathbf{w}_a\]
- [단계 5] \(\mathbf{Y}_a\)을 다중종속변수 행렬으로, 잠재변수 \(\mathbf{t}_a\)를 독립변수 벡터로 하는 회귀모형으로부터 기울기 (로딩벡터) \(\mathbf{q}_a = [q_{a1} \, \cdots \, q_{am}]^\top\)를 산출한다. \[\mathbf{q}_a \leftarrow \left. \mathbf{Y}_a^\top \mathbf{t}_a \middle/ \mathbf{t}_a^\top \mathbf{t}_a \right. \]
- [단계 6] 기울기 벡터 \(\mathbf{q}_a\)의 크기가 1이 되도록 한다. \[\left. \mathbf{q}_a \leftarrow \mathbf{q}_a \middle/ \sqrt{\mathbf{q}_a^\top \mathbf{q}_a} \right.\]
- [단계 7] 잠재변수 \(\mathbf{u}_a\)를 행렬 \(\mathbf{Y}_a\)의 각 열의 가중평균으로 구한다. 이 때, 가중치는 기울기 벡터 \(\mathbf{q}_a\)를 이용한다. \[\mathbf{u}_a \leftarrow \mathbf{Y}_a \mathbf{q}_a\]
- [단계 8] (수렴 확인) [단계 2]에서 [단계 7]까지의 과정을 잠재변수 벡터 \(\mathbf{u}_a\)의 모든 원소값이 수렴할 때까지 반복한다. 수렴이 확인되면 [단계 9]로 진행한다.
- [단계 9] \(\mathbf{t}_a\)를 \(\mathbf{X}_a\)에 회귀시켜, \(\mathbf{X}_a\)을 다중종속변수 행렬으로, \(\mathbf{t}_a\)를 독립변수 벡터로 하는 회귀모형으로부터 로딩벡터 \(\mathbf{p}_a\)를 구한다. \[\mathbf{p}_a \leftarrow \left. \mathbf{X}_a^\top \mathbf{t}_a \middle/ \mathbf{t}_a^\top \mathbf{t}_a \right.\]
- [단계 10] 로딩벡터 \(\mathbf{p}_a\)의 크기를 1로 조정하고, 잠재변수 벡터 \(\mathbf{t}_a\)와 기울기 벡터 \(\mathbf{w}_a\)의 크기를 그에 따라 보정한다. \[d \leftarrow \sqrt{\mathbf{p}_a^\top \mathbf{p}_a}, \, \mathbf{t}_a \leftarrow \mathbf{t}_a d, \, \mathbf{w}_a \leftarrow \mathbf{w}_a d, \, \mathbf{p}_a \leftarrow \frac{1}{d} \mathbf{p}_a \]
- [단계 11] 잠재변수벡터 \(\mathbf{u}_a\)와 \(\mathbf{t}_a\)간의 내부관계 계수 \(b_a\)를 산출한다. \[b_a \leftarrow \left. \mathbf{u}_a^\top \mathbf{t}_a \middle/ \mathbf{t}_a^\top \mathbf{t}_a \right. \]
- [단계 12] 독립변수행렬 \(\mathbf{X}_a\)와 종속변수행렬 \(\mathbf{Y}_a\)로부터 새로 얻어진 잠재변수 벡터 \(\mathbf{t}_a\)가 설명하는 부분을 제거하고 나머지 변동만을 담은 독립변수행렬 \(\mathbf{X}_{a + 1}\)과 종속변수행렬 \(\mathbf{Y}_{a + 1}\)을 구한다. \[\mathbf{X}_{a + 1} \leftarrow \mathbf{X}_a - \mathbf{t}_a \mathbf{p}_a^\top, \, \mathbf{Y}_{a + 1} \leftarrow \mathbf{Y}_a - b_a \mathbf{t}_a \mathbf{q}_a^\top \]
- [단계 13] \(a \leftarrow a + 1\)로 업데이트하고, [단계 1]로 돌아간다. [단계 1] - [단계 13]의 과정을 \(A\)개의 잠재변수를 얻을 때까지 반복한다.
위 알고리즘을 아래 nipals_plsr2이라는 함수로 구현해보자. 이 때, 함수의 입력변수는 아래와 같다.
X: 평균조정된 (\(n \times k\)) 독립변수행렬Y: 평균조정된 (\(n \times m\)) 종속변수행렬A: 잠재변수 개수
nipals_plsr2 <- function(X, Y, A = NULL) {
if (is.vector(Y)) {
Y <- as.matrix(Y, ncol = 1L)
}
if (nrow(X) != nrow(Y)) stop("X and Y must have the same numbers of observations.")
if (is_empty(A) || (A > min(dim(X)))) {
A <- min(dim(X))
}
Tmat <- matrix(NA, nrow = nrow(X), ncol = A)
Umat <- matrix(NA, nrow = nrow(X), ncol = A)
Wmat <- matrix(NA, nrow = ncol(X), ncol = A)
Pmat <- matrix(NA, nrow = ncol(X), ncol = A)
Qmat <- matrix(NA, nrow = ncol(Y), ncol = A)
Bmat <- diag(nrow = A)
for (a in seq_len(A)) {
# 단계 1
j <- sample.int(ncol(Y), size = 1L)
Umat[, a] <- Y[, j]
while (TRUE) {
# 단계 2
Wmat[, a] <- coef(lm(X ~ -1 + Umat[, a]))
# 단계 3
Wmat[, a] <- Wmat[, a] / sqrt(sum(Wmat[, a] ^ 2))
# 단계 4
Tmat[, a] <- X %*% Wmat[, a]
# 단계 5
Qmat[, a] <- coef(lm(Y ~ -1 + Tmat[, a]))
# 단계 6
Qmat[, a] <- Qmat[, a] / sqrt(sum(Qmat[, a] ^ 2))
# 단계 7
u_new <- Y %*% Qmat[, a]
# 단계 8
if (all(near(u_new, Umat[, a]))) break
Umat[, a] <- u_new
}
# 단계 9
Pmat[, a] <- coef(lm(X ~ -1 + Tmat[, a]))
# 단계 10
p_size <- sqrt(sum(Pmat[, a] ^ 2))
Tmat[, a] <- Tmat[, a] * p_size
Wmat[, a] <- Wmat[, a] * p_size
Pmat[, a] <- Pmat[, a] / p_size
# 단계 11
Bmat[a, a] <- coef(lm(Umat[, a] ~ -1 + Tmat[, a]))
# 단계 12
X <- X - Tmat[, a] %*% t(Pmat[, a])
Y <- Y - Bmat[a, a] * Tmat[, a] %*% t(Qmat[, a])
}
return(list(T = Tmat, W = Wmat, P = Pmat,
U = Umat, Q = Qmat, B = Bmat))
}
nipals_fit2 <- nipals_plsr2(X, Y, A = 3)
nipals_fit2## $T
## [,1] [,2] [,3]
## [1,] 1.5777938 -0.4117531 0.44707709
## [2,] -2.2993972 -8.0194451 -0.79803681
## [3,] 0.8145895 -5.1398023 -0.30016399
## [4,] 21.3867331 3.2065797 -0.01083136
## [5,] -17.8784038 5.7770058 -1.68928119
## [6,] 9.2701258 3.6198261 -0.09042238
## [7,] -12.8714412 0.9675888 2.44165865
##
## $W
## [,1] [,2] [,3]
## [1,] -0.1476019 0.4020251 -0.6921766
## [2,] -0.1848598 -0.4862964 0.5027723
## [3,] -0.5065102 0.8236460 0.4768770
## [4,] 0.8312518 0.4651155 0.2794808
##
## $P
## [,1] [,2] [,3]
## [1,] -0.1648502 -0.002634502 -0.5484068
## [2,] -0.1588038 -0.130002426 0.6241044
## [3,] -0.5486716 0.859489687 0.4559029
## [4,] 0.8040928 0.494337847 0.3192118
##
## $U
## [,1] [,2] [,3]
## [1,] 11.60599122 9.401433 -8.559602
## [2,] -0.03696152 3.340983 -7.252566
## [3,] -9.14720168 -11.437784 9.474719
## [4,] 22.98172912 6.021677 -4.914480
## [5,] -7.12619591 9.125035 -6.794076
## [6,] 0.47754875 -7.914057 9.259791
## [7,] -18.75490997 -8.537286 8.786214
##
## $Q
## [,1] [,2] [,3]
## [1,] 0.4832295 0.1202142 0.07948906
## [2,] -0.8754937 -0.9927480 0.99683574
##
## $B
## [,1] [,2] [,3]
## [1,] 0.8443757 0.0000000 0.000000
## [2,] 0.0000000 0.4255916 0.000000
## [3,] 0.0000000 0.0000000 3.2062994.3.4 회귀식 변환
위 NIPALS 알고리즘 수행 결과를 원래 독립변수 \(\mathbf{X}\)와 종속변수 \(\mathbf{Y}\)에 대한 식으로 변환하는 방법은 아래와 같다.
잠재변수행렬 \(\mathbf{T}\)는 하나의 종속변수일 때 살펴봤던 바와 같이 독립변수행렬 \(\mathbf{X}\)와 가중치행렬 \(\mathbf{W}\), 그리고 로딩행렬 \(\mathbf{P}\)의 연산으로 표현된다.
\[\begin{equation} \mathbf{T} = \mathbf{X} \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \end{equation}\]
이를 식 (4.14)에 대입하면,
\[\begin{equation} \mathbf{U} = \mathbf{X} \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \mathbf{B} + \mathbf{H} \end{equation}\]
이를 다시 식 (4.13)에 대입하면,
\[\begin{equation} \begin{split} \mathbf{Y} &= \mathbf{X} \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \mathbf{B} \mathbf{Q}^\top + \mathbf{H} \mathbf{Q}^\top + \mathbf{F}\\ &= \mathbf{X} \mathbf{B}_{PLS} + \mathbf{G} \end{split} \tag{4.15} \end{equation}\]
여기에서 \(\mathbf{G} = \mathbf{H} \mathbf{Q}^\top + \mathbf{F}\)는 독립변수 \(\mathbf{X}\)를 종속변수 \(\mathbf{Y}\)에 회귀시킨 뒤 남은 오차항 행렬이다. 따라서, PLS 모형을 원 독립변수 행렬 \(\mathbf{X}\)에 대한 모형으로 변환할 때의 회귀계수는 아래와 같이 정리된다.
\[\begin{equation} \mathbf{B}_{PLS} = \mathbf{W} \left(\mathbf{P}^\top \mathbf{W}\right)^{-1} \mathbf{B} \mathbf{Q}^\top \end{equation}\]
beta_pls2 <- nipals_fit2$W %*%
solve(t(nipals_fit2$P) %*% nipals_fit2$W) %*%
nipals_fit2$B %*% t(nipals_fit2$Q)
beta_pls2## [,1] [,2]
## [1,] -0.26509645 -2.8773570
## [2,] 0.08864959 2.7249194
## [3,] -0.14258488 0.3377420
## [4,] 0.37330470 -0.7406377이는 앞 4.3.1 절에서 R 패키지 pls를 통해 얻어진 결과와 동일함을 확인할 수 있다.
all(near(beta_pls2, coef(plsr_multi_fit)[, , 1]))## [1] TRUE4.3.5 제곱합 분해
\(\mathbf{Y}\)의 전체제곱합은
\[SST = SSR + SSE\]
로 분해되며, 여기서 \(SSR\)은 다음과 같이 산출된다.
\[\begin{equation} \begin{split} SSR &= \sum_{a = 1}^{A} SSR_a\\ &= \sum_{a = 1}^{A} SS(b_a \mathbf{t}_a \mathbf{q}_a^\top)\\ &= \sum_{a = 1}^{A} b_a^2 SS(\mathbf{t}_a) \end{split} \tag{4.16} \end{equation}\]
이 때, \(SSR_a\)는 잠재변수 \(\mathbf{t}_a\)에 의해 설명되는 \(\mathbf{Y}\)의 변동을 나타낸다.
SSR_a <- diag(
t(nipals_fit2$T %*% nipals_fit2$B) %*%
(nipals_fit2$T %*% nipals_fit2$B)
)
SSR_a## [1] 739.40663 26.91455 100.23860\(SSR_a\)를 전체제곱합 \(SST\)로 나누면 각 잠재변수가 \(\mathbf{Y}\)의 변동에 기여하는 비율을 산출할 수 있다.
SST <- sum(Y ^ 2)
SSR_a / SST## [1] 0.58621653 0.02133840 0.07947119또한, 잠재변수 \(\mathbf{t}_a\)가 설명하는 종속변수행렬 \(\mathbf{Y}\)의 \(j\)번째 열의 변동을 \(SSR_{aj}\)라 할 때, 이는 다음과 같이 산출된다.
\[\begin{equation} SSR_{aj} = q_{ja}^2 SSR_a \end{equation}\]
SSR_aj <- diag(SSR_a) %*% t(nipals_fit2$Q ^ 2)
SSR_aj## [,1] [,2]
## [1,] 172.6593685 566.74726
## [2,] 0.3889542 26.52560
## [3,] 0.6333587 99.60524이렇게 산출된 \(SSR_{aj}\)를 \(j\)번째 종속변수 \(\mathbf{y}_j\)의 제곱합 \(SS(\mathbf{y}_j)\)로 나누면, \(\mathbf{y}_j\)에 대한 \(\mathbf{t}_a\)의 기여도를 얻을 수 있다.
SS_j <- colSums(Y ^ 2)
SSR_aj %*% diag(1 / SS_j)## [,1] [,2]
## [1,] 0.817051715 0.53975930
## [2,] 0.001840593 0.02526248
## [3,] 0.002997154 0.09486213위의 결과에서 \(\mathbf{y}_1\)의 변동은 잠재변수 \(\mathbf{t}_1\)으로 대부분 설명되는 반면, \(\mathbf{y}_2\)의 변동은 잠재변수 \(\mathbf{t}_2\) 및 \(\mathbf{t}_3\)에 의해서도 설명됨을 볼 수 있다.