plsr.Rmd
data(biometric, package = "dmtr")| 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 |
회귀분석에서와 같이 하나의 종속변수에 영향을 주는 \(k\)개의 독립변수가 있다고 하자. 모든 변수는 평균조정되었다고 간주한다. 본 장에서 다루고자 하는 부분최소자승법(partial least squares: PLS)는 앞에서 다룬 주성분 회귀분석(PCR)과 유사하나, 도출되는 새로운 잠재변수들이 다르다.
종속변수가 하나만 존재하는 경우에는 데이터 행렬 \(\mathbf{X}\)와 종속변수벡터 \(\mathbf{y}\)가 동일한 잠재변수로 설명된다고 가정할 수 있다. \((n \times k)\) 데이터 행렬 \(\mathbf{X}\)와 종속변수벡터 \(\mathbf{y}\)에 대하여 동시에 \(A\)개의 잠재변수벡터 \(\mathbf{t}_1, \ldots, \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{2.1}\\ \mathbf{y} &=& \mathbf{t}_1 b_1 + \mathbf{t}_2 b_2 + \cdots + \mathbf{t}_A b_A + \mathbf{f} \tag{2.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{2.3}\\ \mathbf{y} &=& \mathbf{T}\mathbf{b} + \mathbf{f} \tag{2.4} \end{eqnarray}\]
위 반복 알고리즘은 함수는 nipals_plsr()로 구현되어 있다.
y <- biometric %>% select(weight) %>% scale()
X <- biometric %>% select(age, height) %>% scale()
plsr_fit <- nipals_plsr(X, y, ncomp = 2L)
plsr_fit
#> $T
#> LV1 LV2
#> [1,] -0.7858087 0.3080356
#> [2,] 0.3229184 -1.4835503
#> [3,] 0.6295298 0.3450656
#> [4,] -1.9750173 -0.5667679
#> [5,] 0.9178231 -1.5567950
#> [6,] 0.6082079 1.7210543
#> [7,] -0.2725870 -0.3890125
#> [8,] -0.6425632 0.8890449
#> [9,] 0.7325346 -0.4071198
#> [10,] 0.4649624 1.1400451
#>
#> $W
#> LV1 LV2
#> age 0.7992341 -0.6028787
#> height 0.6039375 0.7978329
#>
#> $P
#> LV1 LV2
#> age 0.8321183 -0.6028787
#> height 0.5545982 0.7978329
#>
#> $b
#> LV1 LV2
#> 0.9937095 0.0417354
#>
#> $ncomp
#> [1] 2회귀계수 b와 잠재변수 스코어 행렬 T, 그리고 종속변수에 적용된 평균 및 분산조정 수치를 이용하여, 아래와 같이 분산조정 이전의 종속변수에 대한 회귀식을 추정할 수 있다.
c(
attr(y, "scaled:center"),
plsr_fit$b * attr(y, "scaled:scale") +
rowSums(solve(t(plsr_fit$T) %*% plsr_fit$T) %*% t(plsr_fit$T) * attr(y, "scaled:center"))
)
#> weight LV1 LV2
#> 65.4000000 5.5069003 0.2312876위 추정된 회귀식은 아래와 같이 잠재변수 스코어 행렬 T를 독립변수로 하여 평균조정 이전의 종속변수를 설명하는 회귀모형의 추정 결과이다.
T_A <- as_tibble(plsr_fit[["T"]])
reg_data <- biometric %>% select(weight) %>% bind_cols(T_A)
lm_fit <- fit_linear_regression(reg_data, weight, all_of(names(T_A)))
lm_fit
#> $betas
#> (Intercept) LV1 LV2
#> 65.4000020 5.5069325 0.2312937
#>
#> $hessian
#> (Intercept) LV1 LV2
#> (Intercept) 2.000000e+01 -1.340594e-11 1.398881e-11
#> LV1 -1.340594e-11 1.494215e+01 -6.557255e-12
#> LV2 1.398881e-11 -6.557255e-12 2.105785e+01
#>
#> $mse
#> [1] 7.038464
#>
#> $df
#> [1] 7
#>
#> $se
#> (Intercept) LV1 LV2
#> 0.8389556 0.9706158 0.8176115
#>
#> $sst
#> [1] 276.4
#>
#> $rsq
#> [1] 0.8217466
#>
#> $rsqadj
#> [1] 0.770817
#>
#> $n
#> [1] 10위 과정을 함수 fit_plsr()로 아래와 같이 실행할 수 있다.
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 LV1 5.51 0.971 5.67 7.56e- 4
#> 3 LV2 0.231 0.818 0.283 7.85e- 1위 부분최소자승 회귀분석에 단계별 선택방법을 적용할 경우, 첫 번째 잠재변수만을 독립변수로 사용한 회귀모형이 최종적으로 선택되며, 이 때의 수정결정계수가 두 개의 주성분을 모두 이용한 모형보다 높음을 확인할 수 있다.
selected_lv <- select_variables_stepwise(reg_data, weight, names(T_A))
#> Inital variable: LV1 , p-value = 0.0003124394 < 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_lv))
#> $betas
#> (Intercept) LV1
#> 65.399974 5.506908
#>
#> $hessian
#> (Intercept) LV1
#> (Intercept) 2.000000e+01 1.057487e-11
#> LV1 1.057487e-11 1.494215e+01
#>
#> $mse
#> [1] 6.229061
#>
#> $df
#> [1] 8
#>
#> $se
#> (Intercept) LV1
#> 0.7892440 0.9131028
#>
#> $sst
#> [1] 276.4
#>
#> $rsq
#> [1] 0.8197088
#>
#> $rsqadj
#> [1] 0.7971724
#>
#> $n
#> [1] 10\(m\)개의 종속변수가 존재하여, 종속변수 데이터가 벡터가 아닌 (\(n \times m\)) 행렬
\[\mathbf{Y} = \left[ \mathbf{y}_1 \, \cdots \, \mathbf{y}_m \right]\]
으로 표현될 때, 각각의 종속변수에 대해 따로 잠재변수를 산출하기보다는, 여러 종속변수를 설명하는 공통의 잠재변수행렬 \(\mathbf{T}\)를 산출하는 것이 합리적이라 할 수 있다.
앞 절의 모형을 일반화하여 아래와 같은 모형을 가정한다.
\[\begin{eqnarray} \mathbf{X} &=& \mathbf{T} \mathbf{P}^\top + \mathbf{E} \tag{2.5}\\ \mathbf{Y} &=& \mathbf{U} \mathbf{Q}^\top + \mathbf{F} \tag{2.6}\\ \mathbf{U} &=& \mathbf{T} \mathbf{B} + \mathbf{H} \tag{2.7} \end{eqnarray}\]
식 (2.5)의 모형은 앞서 하나의 종속변수의 경우에서 살펴본 모형식 (2.3)와 동일하다. 식 (2.6)에서 \((n \times A)\) 행렬 \(\mathbf{U}\)는 \(\mathbf{Y}\)를 설명하는 \(A\)개의 잠재변수를 나타내는 행렬이며, \((m \times A)\) 행렬 \(\mathbf{Q}\)는 종속변수행렬 \(\mathbf{Y}\)와 잠재변수행렬 \(\mathbf{U}\)간의 선형관계를 나타내는 로딩행렬이다. 또한 식 (2.7)는 잠재변수행렬 \(\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}\)는 오차항에 해당하는 행렬이다.
다수의 종속변수가 존재하는 경우에도 NIPALS 알고리즘을 이용하여 모형을 추정한다. 이때는 각 잠재변수 \(\mathbf{t}_a\)를 추출할 때 추출한 잠재변수의 수렴 여부를 확인할 필요가 없었던 위 2.1.1절의 경우와는 달리, 각 잠재변수 \(\mathbf{t}_a\)와 \(\mathbf{u}_a\)를 추출하는 과정에서 반복적인(iterative) 기법으로 두 잠재변수 벡터들을 업데이트하며 수렴 여부를 확인하여야 한다.
위 반복 알고리즘은 함수는 nipals_plsr_n()으로 구현되어 있다.