重回帰分析をする際，説明変数の中に互いに相関が高い変数が含まれる場合，通常の最小2乗法 では回帰係数の推定精度が悪くなる（多重共線性）という問題がある。
このような場合の対処法として，以下のものがある。
主成分回帰は，どのような場合にも適用できるものではないことに注意が必要である。
具体的には，以下のような場合に適用する。3 番目は重要（必然的に 2 番目も）。
重回帰分析では，回帰超平面上の予測値と実測値の差の二乗和を最小にするという考え方をとった。
主成分回帰では，回帰超平面に下ろした垂線の長さの二乗和を最小にするという考え方をする。
両者の解析例の違いを示すためのテストデータを作成する。
library(MASS) set.seed(53) m <- c(50, 80) sd <- c(10, 5) d <- mvrnorm(50, m, matrix(c(1, 0.7, 0.7, 1), 2), empirical=TRUE) d <- t(t(d)*sd+m) x <- d[,1] y <- d[,2]
以下のようなデータである（最初の 6 対）。
## x y ## [1,] 552.1985 483.2704 ## [2,] 551.9537 475.2877 ## [3,] 552.1223 485.4826 ## [4,] 561.5995 480.3897 ## [5,] 569.2573 488.7343 ## [6,] 531.6893 475.4973
Fig 1. で，赤が直線回帰，青が主成分回帰である。常に，主成分回帰直線の傾きの方が大きい。
plot(x, y, pch=19, main="Fig. 1") lm.ans <- lm(y~x) abline(lm.ans, col="red") a <- prcomp(d) r1 <- a$rotation[1,1] r2 <- a$rotation[2,1] abline(mean(y)-r2/r1*mean(x), r2/r1, col="blue")
iris データセットを使って説明する。
相関係数は以下のようになる。独立変数のうちに互いに相関の高いものがある。
round(cor(iris[,1:4]), 3)
## Sepal.Length Sepal.Width Petal.Length Petal.Width ## Sepal.Length 1.000 -0.118 0.872 0.818 ## Sepal.Width -0.118 1.000 -0.428 -0.366 ## Petal.Length 0.872 -0.428 1.000 0.963 ## Petal.Width 0.818 -0.366 0.963 1.000
重回帰分析の結果は以下のようになるが，この結果はそのまま受け入れることはできない。
lm.iris <- lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, iris) summary(lm.iris)
## ## Call: ## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, ## data = iris) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.82816 -0.21989 0.01875 0.19709 0.84570 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.85600 0.25078 7.401 9.85e-12 ## Sepal.Width 0.65084 0.06665 9.765 < 2e-16 ## Petal.Length 0.70913 0.05672 12.502 < 2e-16 ## Petal.Width -0.55648 0.12755 -4.363 2.41e-05 ## ## Residual standard error: 0.3145 on 146 degrees of freedom ## Multiple R-squared: 0.8586, Adjusted R-squared: 0.8557 ## F-statistic: 295.5 on 3 and 146 DF, p-value: < 2.2e-16
実は，この分析での独立変数のトレランスは以下のようになる。
(tol <- 1/diag(solve(cor(iris[,2:4]))))
## Sepal.Width Petal.Length Petal.Width ## 0.78689664 0.06623581 0.07025267
Petal.Length, Petal.Width のトレランスはそれぞれ 0.066，0.070 となり，一般的な基準である 0.1 より低いため多重共線性の存在が示唆される。
重回帰分析による予測結果は Fig. 2 のようになる。
lm.prediction <- predict(lm.iris) plot(lm.prediction, iris$Sepal.Length, col=rep(c(1,2,4), each=50), pch=19, main="Fig.2") abline(0, 1)
library(pls)
## ## Attaching package: 'pls' ## ## 以下のオブジェクトは 'package:stats' からマスクされています: ## ## loadings
pcr.iris <- pcr(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, ncomp=3, data=iris) pcr.prediction <- predict(pcr.iris) pcr.prediction <- cbind(pcr.prediction[,,1:3])
PCR 回帰の結果は以下のように示される。
summary(pcr.iris)
## Data: X dimension: 150 3 ## Y dimension: 150 1 ## Fit method: svdpc ## Number of components considered: 3 ## TRAINING: % variance explained ## 1 comps 2 comps 3 comps ## X 95.09 99.12 100.00 ## Sepal.Length 74.29 82.11 85.86
予測値の最初の部分を示す。
head(pcr.prediction)
## 1 comps 2 comps 3 comps ## 1 4.880974 4.987828 5.015416 ## 2 4.899464 4.717888 4.689997 ## 3 4.858013 4.789387 4.749251 ## 4 4.929821 4.808353 4.825994 ## 5 4.877276 5.041815 5.080499 ## 6 4.996961 5.360115 5.377194
第 1 主成分だけを用いて予測すると Fig. 3 のようになる。
plot(pcr.prediction[,1], iris$Sepal.Length, col=rep(c(1,2,4), each=50), pch=19, xlim=c(4, 8), main="Fig. 3") abline(0, 1)
第 1，第 2 主成分だけを用いて予測すると Fig. 4 のようになる。
plot(pcr.prediction[,2], iris$Sepal.Length, col=rep(c(1,2,4), each=50), pch=19, xlim=c(4, 8), main="Fig. 4") abline(0, 1)
第 1，第 2，第 3 主成分の全部を用いて予測すると Fig. 5 のようになる。
plot(pcr.prediction[,3], iris$Sepal.Length, col=rep(c(1,2,4), each=50), pch=19, xlim=c(4, 8), main="Fig. 5") abline(0, 1)
all.equal(lm.prediction, pcr.prediction[,3])
## [1] TRUE
まず，重回帰分析に使った 3 つの独立変数の主成分分析を行い，それぞれの主成分得点を求める。
prcomp.ans <- prcomp(iris[,2:4]) prcomp.predict <- predict(prcomp.ans) # 主成分得点（prcomp.ans$x と同じ）
主成分得点の最初の方を示す。
round(head(prcomp.predict), 3)
## PC1 PC2 PC3 ## [1,] -2.592 -0.183 0.032 ## [2,] -2.543 0.311 -0.032 ## [3,] -2.654 0.117 -0.046 ## [4,] -2.461 0.208 0.020 ## [5,] -2.602 -0.281 0.044 ## [6,] -2.280 -0.621 0.020
第 1 主成分得点だけを使って，重回帰分析を行った結果を示す。
lm1 <- lm(iris$Sepal.Length ~ prcomp.predict[,1]) summary(lm1)
## ## Call: ## lm(formula = iris$Sepal.Length ~ prcomp.predict[, 1]) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.28826 -0.29759 -0.02007 0.28029 1.06985 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 5.84333 0.03440 169.87 <2e-16 ## prcomp.predict[, 1] 0.37123 0.01795 20.68 <2e-16 ## ## Residual standard error: 0.4213 on 148 degrees of freedom ## Multiple R-squared: 0.7429, Adjusted R-squared: 0.7412 ## F-statistic: 427.6 on 1 and 148 DF, p-value: < 2.2e-16
第 1，第 2 主成分得点だけを使って，重回帰分析を行った結果を示す。
lm2 <- lm(iris$Sepal.Length ~ prcomp.predict[,1:2]) summary(lm2)
## ## Call: ## lm(formula = iris$Sepal.Length ~ prcomp.predict[, 1:2]) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.03052 -0.23189 -0.00612 0.21751 0.77520 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 5.84333 0.02879 202.936 < 2e-16 ## prcomp.predict[, 1:2]PC1 0.37123 0.01503 24.704 < 2e-16 ## prcomp.predict[, 1:2]PC2 -0.58457 0.07295 -8.014 3.16e-13 ## ## Residual standard error: 0.3527 on 147 degrees of freedom ## Multiple R-squared: 0.8211, Adjusted R-squared: 0.8186 ## F-statistic: 337.3 on 2 and 147 DF, p-value: < 2.2e-16
第 1 〜 第 3 主成分得点を使って，重回帰分析を行った結果を示す。
lm3 <- lm(iris$Sepal.Length ~ prcomp.predict[,1:3]) summary(lm3)
## ## Call: ## lm(formula = iris$Sepal.Length ~ prcomp.predict[, 1:3]) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.82816 -0.21989 0.01875 0.19709 0.84570 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 5.84333 0.02568 227.519 < 2e-16 ## prcomp.predict[, 1:3]PC1 0.37123 0.01340 27.697 < 2e-16 ## prcomp.predict[, 1:3]PC2 -0.58457 0.06506 -8.984 1.22e-15 ## prcomp.predict[, 1:3]PC3 0.86983 0.13969 6.227 4.80e-09 ## ## Residual standard error: 0.3145 on 146 degrees of freedom ## Multiple R-squared: 0.8586, Adjusted R-squared: 0.8557 ## F-statistic: 295.5 on 3 and 146 DF, p-value: < 2.2e-16
prcomp.lm <- cbind(predict(lm1), predict(lm2), predict(lm3)) head(cbind(prcomp.lm, pcr.prediction), 10)
## 1 comps 2 comps 3 comps ## 1 4.880974 4.987828 5.015416 4.880974 4.987828 5.015416 ## 2 4.899464 4.717888 4.689997 4.899464 4.717888 4.689997 ## 3 4.858013 4.789387 4.749251 4.858013 4.789387 4.749251 ## 4 4.929821 4.808353 4.825994 4.929821 4.808353 4.825994 ## 5 4.877276 5.041815 5.080499 4.877276 5.041815 5.080499 ## 6 4.996961 5.360115 5.377194 4.996961 5.360115 5.377194 ## 7 4.898979 4.957293 4.894684 4.898979 4.957293 4.894684 ## 8 4.918727 4.970316 5.021245 4.918727 4.970316 5.021245 ## 9 4.903162 4.663900 4.624913 4.903162 4.663900 4.624913 ## 10 4.915514 4.784900 4.881642 4.915514 4.784900 4.881642
左の 3 列は主成分得点（第 1 〜 3 主成分得点）を使った重回帰分析の予測値，右の 3 列は pcr 回帰（第 1 〜 3 主成分）を使った予測値である。それぞれ対応する予測値が同じであることがわかる。
つまり，主成分回帰と，主成分得点を用いる重回帰分析は同じだということである。