従属変数と有意な相関を持つ独立変数を採用すれば,全ての独立変数の偏回帰係数が有意になると考えるのは間違いだ。
library(MASS) set.seed(123) n <- 44 r <- matrix(c(1.0, 0.3, 0.4, 0.35, 0.3, 1.0, 0.2, 0.4, 0.4, 0.2, 1.0, 0.5, 0.35, 0.4, 0.5, 1.0), 4) d <- data.frame(matrix(mvrnorm(n, mu=rep(0, 4), Sigma=r, empirical=TRUE), n)) colnames(d) <- c("y", "x1", "x2", "x3")
このデータで,変数間の相関係数は,以下のようになる。
y x1 x2 x3 y 1.00 0.3 0.4 0.35 x1 0.30 1.0 0.2 0.40 x2 0.40 0.2 1.0 0.50 x3 0.35 0.4 0.5 1.00
従属変数 \(y\) と独立変数 \(x_1\),\(x_2\),\(x_3\) の相関係数はそれぞれ 0.3,0.4,0.35 であり,有意水準 5% で "母相関係数 = 0" の検定を行うといずれも帰無仮説は棄却される(\(p\) 値はそれぞれ 0.048,0.007,0.020 である)。ちなみに,このデータは変数ごとに標準化されている(平均値 = 0,標準偏差 = 1)ので,偏回帰係数(Estimate)は標準化偏回帰係数でもあるが,それぞれの変数を 1 個だけ用いた単回帰分析の標準化偏回帰係数は従属変数と独立変数の相関係数と同じであり,したがって "偏回帰係数 = 0" の検定の \(p\) 値と母相関係数の検定結果の \(p\) 値は同じである。
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0 0.146 0.000 1.000 x1 0.3 0.147 2.038 0.048
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0 0.140 0.000 1.000 x2 0.4 0.141 2.828 0.007
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.00 0.143 0.000 1.00 x3 0.35 0.145 2.421 0.02
次に,3 つの変数 \(x_1\),\(x_2\),\(x_3\) 全てを独立変数として重回帰分析を行う。偏回帰係数の検定結果は,予想に反して,全ての独立変数で帰無仮説は棄却されない。
また,\(x_1\) と \(y\) との相関は,\(x_3\) と \(y\) の相関より低いが,偏回帰係数の \(p\) 値より小さい。また,\(x_3\) より \(x_1\) のほうが \(y\) に及ぼす影響が大きい(それぞれ 0.190,0.124 である)。
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.000 0.138 0.000 1.000 x1 0.190 0.152 1.251 0.218 x2 0.300 0.161 1.861 0.070 x3 0.124 0.172 0.718 0.477
このシミュレーション例でいえることは「単回帰分析の結果を積み上げても重回帰分析の結果にはならない」ということである。重回帰分析は従属変数と独立変数の関係だけではなく,独立変数相互間の関係も考慮されるので,従属変数と独立変数 1 個ずつを取り上げる単回帰分析の結果を寄せ集めたものとは異なるのである。
独立変数の候補変数がたくさんあるとき,その中から有効な変数を取り上げて意味のある重回帰モデルを作ろうとするとき,単相関係数(単回帰分析)の結果から従属変数と相関の高い変数を取り出して重回帰分析を行うという方法は誤った方法である。変数選択を行いたいなら,ステップワイズ変数選択などを行うべきなのである。
検定についての誤解の一つに,「\(p\) 値が小さいことは,主張したいことが正しいという証拠になる」というのがある。しかし,あらゆる検定において,\(p\) 値は検定統計量が帰無仮説のもとで起こりにくいことを示すが,サンプルサイズにも依存するということである。回帰の分散分析の帰無仮説は「回帰により説明される分散は誤差分散なみである」ということであり,\(p\) 値が小さいということは「回帰により説明される分散は誤差分散に比べて十分に大きい」ということである。しかし,重回帰式が使えるかどうかは決定係数を見るべきである。サンプルサイズがちょっと大きければ,決定係数が小さくても回帰の分散分析の帰無仮説は棄却される。
library(MASS) set.seed(123) n <- 310 r <- matrix(c(1.0, 0.3, 0.4, 0.35, 0.3, 1.0, 0.2, 0.4, 0.4, 0.2, 1.0, 0.5, 0.35, 0.4, 0.5, 1.0), 4) d2 <- data.frame(matrix(mvrnorm(n, mu=rep(0, 4), Sigma=r, empirical=TRUE), n)) colnames(d2) <- c("y", "x1", "x2", "x3")
このデータの相関係数は前節のデータのものと同じである。
y x1 x2 x3 y 1.00 0.3 0.4 0.35 x1 0.30 1.0 0.2 0.40 x2 0.40 0.2 1.0 0.50 x3 0.35 0.4 0.5 1.00
唯一違うのは,前節のデータはサンプルサイズが 44 であり,このデータのサンプルサイズは 310 ということである。
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.000 0.050 0.000 1.000 x1 0.190 0.055 3.459 0.001 x2 0.300 0.058 5.148 0.000 x3 0.124 0.062 1.987 0.048
独立変数 \(x_1\),\(x_2\),\(x_3\) の偏回帰係数(標準化偏回帰係数)は前節の結果と同じであるが,偏回帰係数の検定の \(p\) 値は前節の結果とは違い,全て 0.05 より小さく,帰無仮説は棄却される。更に回帰の分散分析の結果は \(F\)(3,306)=28.849 であり,\(p\) 値は 0.000 で,帰無仮説は棄却される。
しかし,決定係数 \(R^2\) は 0.220 であり,とても有効な予測ができるようなものではないことがわかる。
偏回帰係数の検定において全ての \(p\) 値が 0.05 より大きければ,回帰の分散分析の \(p\) 値も 0.05 より大きいというのが誤っているということは,第 1 節のデータの分析結果を見れば明らかである。
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.000 0.138 0.000 1.000 x1 0.190 0.152 1.251 0.218 x2 0.300 0.161 1.861 0.070 x3 0.124 0.172 0.718 0.477
回帰の分散分析結果は \(F\)(3,40)=3.771 であり,\(p\) 値は 0.018 である。
実際に計算してみなければわからない。
library(MASS) set.seed(123) n <- 44 r <- matrix(runif(16, 0.8, 0.9), 4, 4) r <- (r+t(r))/2 diag(r) <- 1 d3 <- data.frame(matrix(mvrnorm(n, mu=rep(0, 4), Sigma=r, empirical=TRUE), n)) colnames(d3) <- c("y", "x1", "x2", "x3")
このデータで,変数間の相関係数は,以下のようになる。0.8 以上のものばかりである。
y x1 x2 x3 y 1.000 0.886 0.848 0.878 x1 0.886 1.000 0.849 0.873 x2 0.848 0.849 1.000 0.828 x3 0.878 0.873 0.828 1.000
重回帰分析の結果は,以下のようになる。回帰係数の符号,大きさも特に問題はない。
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.000 0.062 0.000 1.000 x1 0.381 0.144 2.638 0.012 x2 0.232 0.125 1.854 0.071 x3 0.353 0.136 2.604 0.013
トレランス(VIF)も問題のあるレベルではない。
tolerance VIF x1 0.187 5.356 x2 0.247 4.041 x3 0.211 4.745
library(MASS) set.seed(1423) n <- 44 r <- matrix(runif(16, 0.7, 0.79), 4, 4) r <- (r+t(r))/2 r[2,3] <- r[3,2] <- 0.94 diag(r) <- 1 d4 <- data.frame(matrix(mvrnorm(n, mu=rep(0, 4), Sigma=r, empirical=TRUE), n)) colnames(d4) <- c("y", "x1", "x2", "x3")
このデータで,変数間の相関係数は,以下のようになる。0.8 以上のものは 1 個しかない。
y x1 x2 x3 y 1.000 0.707 0.756 0.754 x1 0.707 1.000 0.940 0.720 x2 0.756 0.940 1.000 0.785 x3 0.754 0.720 0.785 1.000
重回帰分析の結果は,以下のようになる。
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.000 0.094 0.000 1.000 x1 0.037 0.280 0.133 0.895 x2 0.390 0.313 1.248 0.219 x3 0.421 0.154 2.737 0.009
トレランス(VIF)からは,独立変数 \(x_2\) に問題がありそうである。
tolerance VIF x1 0.116 8.649 x2 0.092 10.827 x3 0.382 2.619