必要なパッケージ

require(plotly)

重回帰分析

　回帰分析を行おうとする実験において、一つ以上の共変量が利用できる場合が少なくない。例えば、作物の収量を制御あるいは予測したいときに、収量を説明するのに、施肥量だけでなく、潅水量や日射量も用いることができれば、制御や予測がより正確になる可能性がある。あるいは、ある個体の遺伝的特性を予測するのに、1つの遺伝子だけでなく、複数の遺伝子に関する情報を用いたほうが、遺伝的特性をより正確に予測できる可能性がある。

　今、変数\(x_1,x_2,...,x_k\)（共変量covariatesとも呼ばれる）があり、これら共変量に対する応答も含めた\(n\)セットの計測データ\(x_{i1},x_{i2},...x_{ik},y_i,\:\:i=1,...,n\)があるとする。このとき、重回帰は、共変量\(x_i\)の線形結合に切片と誤差を加えて、応答を表現する。

　すなわち、 \[ y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \cdots + + \beta_kx_{ik} + \epsilon_i \] 　ここで、誤差\(\epsilon_i\)は、それぞれ独立に平均0、分散\(\sigma^2\)の正規分布に従うと仮定する。\(k\)は、共変量の数であり、パラメータ数は切片である\(\beta_0\)を含めて全部で\(p = k+1\)となる。なお、重回帰分析では、多くの場合、共変量の数\(k\)ではなく、パラメータ数\(p\)を、次元の表現や統計量の計算に用いる。

　重回帰を行う場合も、単変量の場合と同じように、回帰係数や切片の推定や検定、応答の予測やその信頼区間に興味があると考えられる。ここから、その方法について説明する。

回帰モデルの行列による表記

　\(n\)セットの観察値全て（\(i = 1,...,n\)）についての回帰式は、以下のように記述できる。 \[ y_1 = \beta_0 + \beta_1x_{11} + \beta_2x_{12} + \cdots + + \beta_kx_{1k} + \epsilon_1 \] \[ y_2 = \beta_0 + \beta_1x_{21} + \beta_2x_{22} + \cdots + + \beta_kx_{2k} + \epsilon_2 \] \[ \vdots \] \[ y_n = \beta_0 + \beta_1x_{n1} + \beta_2x_{n2} + \cdots + + \beta_kx_{nk} + \epsilon_n \]

　これら式は、それらをベクトルや行列に「束ねる」ことで、ベクトルと行列を用いたかたちで表すことができる。すなわち、 \[ \boldsymbol{y} = \boldsymbol{X\beta}+\boldsymbol{\epsilon} \]

ここで、\(\boldsymbol y\)は、 \[ \boldsymbol y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} \] また、\(\boldsymbol X\)は、 \[ \boldsymbol X = \begin{bmatrix} \boldsymbol x_1 \\ \boldsymbol x_2 \\ \vdots \\ \boldsymbol x_n \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & x_{12} & \dots & x_{1k} \\ 1 & x_{21} & x_{22} & \dots & x_{2k} \\ \vdots & & &\vdots & \\ 1 & x_{n1} & x_{n2} & \dots & x_{nk} \\ \end{bmatrix} \] \(\beta\)と\(\epsilon\)は、 \[ \boldsymbol \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{bmatrix} \]

\[ \boldsymbol \epsilon = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix} \] である。

なお、\(\boldsymbol{y, \epsilon}\)は、\(n\times1\)のベクトル、\(\boldsymbol{X}\)は、\(n\times p\)の行列、\(\beta\)は\(p \times 1\)のベクトルである。また、\(p=k+1\)である。

ここで、\(\boldsymbol \beta\)の最小二乗推定値を求めるには、次の残差平方和\(E\)を最小化することになる。 \[ E=\sum_{i=1}^n e_i^2 = \sum_{i=1}^n \left[y_i - (\beta_0+\beta_1x_{i1}+\cdots+\beta_kx_{ik})\right]^2 \]

一般にベクトル\(\boldsymbol x = [x_1, x_2, \dots, x_n]’\)の要素の二乗和は、 \[ x_1^2 + x_2^2 + \cdots + x_n^2 = [x_1\;x_2\;\cdots x_n] \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \boldsymbol {x'x} \] と表せる。したがって、上記の残差平方和\(E\)は、 \[ E = (\boldsymbol y - \boldsymbol {X\beta})' (\boldsymbol y - \boldsymbol {X\beta}) \\ = \boldsymbol y'{\boldsymbol y} - {\boldsymbol y}'{\boldsymbol X \beta} - ({\boldsymbol X\beta})' {\boldsymbol y} + ({\boldsymbol X\beta})'{\boldsymbol X\beta} \\ = \boldsymbol {y}'{\boldsymbol y} - {\boldsymbol \beta}'({\boldsymbol X}'{\boldsymbol y}) - ({\boldsymbol \beta}'{\boldsymbol X}'){\boldsymbol y} + {\boldsymbol \beta}'{\boldsymbol X}' {\boldsymbol X\beta} \\ = \boldsymbol {y}'{\boldsymbol y} - 2 {\boldsymbol \beta}'({\boldsymbol X}'{\boldsymbol y}) + {\boldsymbol \beta}'{\boldsymbol X}' {\boldsymbol X\beta} \\ \] となることが分かる。

ここで、ベクトルの微分が、 \[ \frac {\partial}{\partial {\boldsymbol \beta}} {\boldsymbol \beta}' {\boldsymbol x} = \boldsymbol x \] \[ \frac {\partial}{\partial {\boldsymbol \beta}} {\boldsymbol \beta}' {\boldsymbol A} {\boldsymbol \beta} = {\boldsymbol A\beta} + {\boldsymbol A}'{\boldsymbol \beta} = ({\boldsymbol A} + {\boldsymbol A}'){\boldsymbol \beta} \] であることを用いると、

\[ \frac {\partial E}{\partial {\boldsymbol \beta}} = -2 \boldsymbol X' \boldsymbol y + \left[\boldsymbol X' \boldsymbol X + (\boldsymbol X' \boldsymbol X)'\right]\boldsymbol \beta \\ = - 2 \boldsymbol X' \boldsymbol y + 2 \boldsymbol X' \boldsymbol X \boldsymbol \beta = 0 \] となる。

したがって、残差平方和\(E\)を最小化する解 \[ \boldsymbol{b} = \begin{bmatrix} b_0 \\ b_1 \\ \vdots \\ b_k \end{bmatrix} \] は、 \[ \boldsymbol{X}' \boldsymbol{Xb} = \boldsymbol{X}' \boldsymbol{y} \] を満たす。

　これより、\(\boldsymbol{\beta}\)の最小二乗推定量は、 \[ \boldsymbol{b}=(\boldsymbol{X'X})^{-1}\boldsymbol{X'y} \] として計算される。

　また、推定値は、 \[ \boldsymbol{\hat y} = \boldsymbol{Xb} = \boldsymbol{X}(\boldsymbol{X'X})^{-1}\boldsymbol{X'y} \] となる。

また残差は、 \[ \boldsymbol{e} = \boldsymbol{y}-\boldsymbol{\hat y}=\boldsymbol{y}- \boldsymbol{Xb}=\boldsymbol{y}- \boldsymbol{X}(\boldsymbol{X'X})^{-1}\boldsymbol{X'y}=(\boldsymbol{I} - \boldsymbol{X}(\boldsymbol{X'X})^{-1}\boldsymbol{X'})\boldsymbol{y} \] と表せる。\(\boldsymbol{I}\)は、\(n\times n\)の単位行列である。

　ここで、上の2つの式をよく眺めてみると、残差と推定値の式の中に同じ行列が表れていることが分かる。残差と推定値の式の中にあらわれる行列 \[ \boldsymbol{H}=\boldsymbol{X}(\boldsymbol{X'X})^{-1}\boldsymbol{X'} \] は重要な行列であり、ハット（hat）行列とよばれる。

　ハット行列\(\boldsymbol{H}\)を用いて表すと、推定値と残差は、それぞれ \[ \boldsymbol{\hat y} = \boldsymbol{Hy} \] \[ \boldsymbol{e} = (\boldsymbol{I}-\boldsymbol{H})\boldsymbol{y} \] とコンパクトなかたちで書き表すことができる。

　なお、行列\(\boldsymbol{H}\)および\(\boldsymbol{I-H}\)は直交射影行列であり、\(n\)次元ベクトル\(\boldsymbol{y}\)は、\(\boldsymbol{H}\)によって\(n\)次元の推定値ベクトル\(\boldsymbol{\hat y}\)に射影される。また、\(n\)次元ベクトル\(\boldsymbol{y}\)は、\(\boldsymbol{I-H}\)によって、\(n\)次元の残差ベクトル\(\boldsymbol{e}\)に射影される。

補足：射影行列の定義
正方行列\(\boldsymbol{P}\)が、 \[ \boldsymbol{P}^2=\boldsymbol{P} \] を満たすとき、\(\boldsymbol{P}\)を実ベクトル空間上の射影行列という。さらに、 \[ \boldsymbol{P}^T=\boldsymbol{P} \] を満たすとき、直交射影行列とよぶ。

　いかなる射影行列\(\boldsymbol{A}\)もべき等（idempotent）で、\(\boldsymbol{A}^2=\boldsymbol{A}\)が成り立つ。別の言い方をすると、射影の射影は、もとの射影と同じことを意味する（一旦影となる空間へ射影したら、影から影へさらに射影しても変化がないことを意味する）。なお、幾何学的には、\(\boldsymbol{\hat y}\)と \(\boldsymbol{e}\)は、それらの射影行列の積が0になることから、互いに直交することが分かる。実際、\(\boldsymbol{H}\)はべき等であるため、\(\boldsymbol{H}(\boldsymbol{I-H})=\boldsymbol{H}-\boldsymbol{H}^2=\boldsymbol{H}-\boldsymbol{H}=\boldsymbol{0}\)が成り立つ。

　なお、残差\(\epsilon_i\)は、\(i,i'\)間で互いに独立であり、ベクトル\(\boldsymbol{\epsilon}\)の分散は\(\sigma^2\boldsymbol{I}\)となる（つまり、\(i, i'\)間の共分散は0である）。

例　ブルーベリーの糖度と酸度と嗜好性

　ブルーベリー品種Primadonnaで計測された果実の糖度（可溶性固形分、SSC）\(X_1\)と酸度（pH）\(X_2\)と官能試験から得られた嗜好性\(Y\)の関係（Gilbert et al., 2015, PLOS ONE 10:e0138494）をもとに、上述した計算を行ってみる。

x1 <- c(11.3, 11, 11.4, 10.6, 12.5, 13.9, 10.8, 12.6, 10.7, 12.4, 11.6, 11.6, 11.2, 11.1, 11.2, 12.3, 10.3, 11.1, 11, 11.3, 12.3, 13, 13, 12.6, 11.7)
x2 <- c(3.95, 4, 3.99, 3.18, 3.34, 3.42, 3.51, 3.62, 3.46, 4.11, 4.41, 3.94, 3.47, 3.55, 3.52, 3.67, 3.84, 4.06, 4.18, 3.58, 3.52, 3.76, 4.22, 3.8, 3.72)
y <- c(24, 21.6, 20.9, 27.2, 25.5, 27.5, 25, 29.7, 22.1, 22.8, 20, 15.6, 27.2, 23.8, 25.1, 24.8, 20, 24.6, 23.2, 28.5, 27.8, 26.6, 28.2, 25.2, 25.5)

　まずは、データを眺めてみる。

df <- data.frame(y, x1, x2)
pairs(df)

あまり明瞭なパターンではないが、\(y\)と\(x1\)には、一方が増えるともう一方も増えるような傾向がある。\(y\)と\(x2\)の関係は、その逆に見える。

　では、ここから、上述した回帰係数の計算を行っていく。まず、行列\(\boldsymbol X\)を準備する。最初の1列（切片に対応する列）は全て1からなる列になる。

X <- cbind(rep(1, length(x1)), x1, x2)
dim(X)
## [1] 25  3
head(X)
##          x1   x2
## [1,] 1 11.3 3.95
## [2,] 1 11.0 4.00
## [3,] 1 11.4 3.99
## [4,] 1 10.6 3.18
## [5,] 1 12.5 3.34
## [6,] 1 13.9 3.42

　回帰係数の推定値を計算してみる。solveは、逆行列を求める関数である。

b <- solve(t(X) %*% X) %*% t(X) %*% y
b
##         [,1]
##    25.371650
## x1  1.541424
## x2 -5.038987

　切片\(\beta_0\)の推定値が25.37、糖度\(X_1\)の回帰係数\(\beta_1\)の推定値が1.54、酸度\(X_2\)の回帰係数\(\beta_2\)の推定値が-5.04であることが分かる。

次に、ハット行列を計算してみる。

H <- X %*% solve(t(X) %*% X) %*% t(X)

得られた行列がべき等（\(\boldsymbol H^2 = \boldsymbol H\)）かどうか、念のため確認してみる。\(\boldsymbol H^2 - \boldsymbol H\)を計算して、全要素の最大値と最小値を調べる。

max(H %*% H - H)
## [1] 1.149081e-14
min(H %*% H - H)
## [1] -1.286256e-14

全ての要素が0であることが分かる。

では、\(\boldsymbol y\)の推定値\(\boldsymbol {\hat y}\)を計算してみよう。これは、ハット行列を使って簡単に計算できる。

yHat <- H %*% y

また、誤差\(\boldsymbol e\)を計算して、誤差\(\boldsymbol e\)と推定値\(\boldsymbol {\hat y}\)が直交していることも確認してみる。

e <-  y - yHat
head(e)
##            [,1]
## [1,]  1.1142579
## [2,] -0.5713656
## [3,] -1.9383251
## [4,]  1.5132346
## [5,] -2.3092331
## [6,] -2.0641077
t(e) %*% yHat 
##              [,1]
## [1,] 2.575485e-11

なお、誤差は、\((\boldsymbol{I}-\boldsymbol{H})\boldsymbol{y}\)として計算もできる。

e <- (diag(1, nrow(H)) - H) %*% y
head(e)
##            [,1]
## [1,]  1.1142579
## [2,] -0.5713656
## [3,] -1.9383251
## [4,]  1.5132346
## [5,] -2.3092331
## [6,] -2.0641077

\(\boldsymbol y\)と推定値\(\boldsymbol {\hat y}\)を図示する。

plot(y, yHat)
abline(0, 1)

\(y\)と\(\hat y\)相関係数も計算してみる。

cor(y, yHat)
##           [,1]
## [1,] 0.6385623

この相関係数のは、単回帰のときと同様、決定係数\(R^2\)に一致する。

cor(y, yHat)^2
##           [,1]
## [1,] 0.4077619

　なお、単回帰と同様にlm関数を用いて重回帰を行うことができる。

model <- lm(y ~ x1 + x2)
summary(model)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7986 -1.0381  0.1364  1.5132  4.0544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  25.3716     9.3883   2.702  0.01301 * 
## x1            1.5414     0.5917   2.605  0.01617 * 
## x2           -5.0390     1.7161  -2.936  0.00764 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.599 on 22 degrees of freedom
## Multiple R-squared:  0.4078, Adjusted R-squared:  0.3539 
## F-statistic: 7.574 on 2 and 22 DF,  p-value: 0.003144

　あるいは、先ほど準備しておいたdfを使って、

model <- lm(y ~ ., data = df)
summary(model)
## 
## Call:
## lm(formula = y ~ ., data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7986 -1.0381  0.1364  1.5132  4.0544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  25.3716     9.3883   2.702  0.01301 * 
## x1            1.5414     0.5917   2.605  0.01617 * 
## x2           -5.0390     1.7161  -2.936  0.00764 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.599 on 22 degrees of freedom
## Multiple R-squared:  0.4078, Adjusted R-squared:  0.3539 
## F-statistic: 7.574 on 2 and 22 DF,  p-value: 0.003144

と計算することもできる。y ~ .（\(y\)を\(y\)以外の変数でモデル化する)という表現は、 共変料が多い場合に便利なのでおぼえておくとよい。

平方和と分散分析表

　重回帰の平方和SSY, SSR, SSEは、行列表記を用いると簡単に表せる。

　まず、\(y\)の平方和であるSSYについて考える。SSYは、単回帰と同様に \[ SSY = \sum_i^n y_i^2 - n \bar y^2 \] として表すことができるが、これを、それぞれの要素が1の\(n\times n\)行列\(\boldsymbol{J}\)を導入して、以下のようにベクトルを用いて表す（計算する）こともできる。 \[ SSY = \boldsymbol{y'y} - \frac 1 n \boldsymbol{y'Jy}=\boldsymbol{y'}\biggl(\boldsymbol{I}-\frac 1 n \boldsymbol{J}\biggr)\boldsymbol{y} \]

同様に、残差\(e\)の平方和SSEは、 \[ SSE = \boldsymbol{e'e}=\boldsymbol{y'(I-H)'(I-H)y}=\boldsymbol{y'(I-H)y} \] として計算できる。これは、\(\boldsymbol{I-H}\)が対称な射影行列であるため、\(\boldsymbol{(I-H)'(I-H) = (I-H)(I-H) = I-H}\)が成り立つためである。

　SSYとSSEの差が回帰が説明する平方和となるが、それは、 \[ SSR=SSY-SSE=\boldsymbol{y'}\biggl(\boldsymbol{H}-\frac 1 n \boldsymbol{J}\biggr)\boldsymbol{y} \] と計算できる。このように、いずれの平方和も、行列計算を用いて容易に計算できる。

　なお、SSY, SSR, SSEの自由度は、\(n-1, p-1, n-p\)となる。したがって、重回帰の分散分析表は以下のようになる。

要因 自由度(df) 平方和(SS) 平均平方(MS) \(F\) \(p\)値
回帰 \(p-1\) \(SSR\) \(MSR = \frac {SSR}{p-1}\) \(F = \frac {MSR}{MSE}\) \(p(F_{p-1, n-p} > F)\)
誤差 \(n-p\) \(SSE\) \(MSE = \frac {SSE}{n-p}\)
全体 \(n-1\) \(SSY\)

となる。

　回帰の平均平方\(MSR\)が誤差の平均平方\(MSE\)に比較して大きく、大きな\(F\)値が得られたとき、帰無仮説\(H_0\)：共変量\(x_1,...,x_k\)が応答\(y\)に影響しない、が棄却される。より正確には、帰無仮説は、 \[ H_0: \beta_1 = \beta_2 = \cdots = \beta_k = 0 \] と書くことができ、対立仮説は、少なくとも1つの\(\beta_i,\:\:i=1,...,k\)が0でない、という仮説になる。

　では、上に示した式を用いて、重回帰の分散分析を行ってみる。

まずは、計算のための準備をする。

n <- length(y)
p <- ncol(X)
J <- matrix(1, nrow = n, ncol = n)
I <- diag(1, n)

次に、上に示した式を用いてSSR, SSE, SSYを計算する。

# SSR, SSE, SSY
SSE <- t(y) %*% (I - H) %*% y
SSY <- t(y) %*% (I - 1/n * J) %*% y
SSR <- SSY - SSE
#または SSR <- t(y) %*% (H - 1/n * J) %*% y

次に、SSR, SSEを自由度で割り、回帰と残差の平均平方を求める。

## devided by df
MSR <- SSR / (p-1)
MSR
##          [,1]
## [1,] 51.16792
MSE <- SSE / (n-p)
MSE
##         [,1]
## [1,] 6.75608

祭儀に、平均平方の比をとって\(F\)値を計算し、検定を行う。

(F <- MSR / MSE)
##         [,1]
## [1,] 7.57361
1 - pf(F, p-1, n-p)
##             [,1]
## [1,] 0.003143837

なお、単回帰のときと同様、\(R^2\) \[ R^2 = \frac {SSR}{SSY} = 1-\frac{SSE}{SSY} \] は決定係数とよばれ、回帰の当てはまりの良さを表す。

　ところで、より多くの変数を回帰に加えれば、その変数が被説明変数\(y\)に対して影響をしていない場合でも、確率的な摂動により\(y\)の変動の一部を説明してしまう。このことにより、新たな変数を加えると、\(R^2\)はいつも大きくなる。

　このように考えると、2つのモデルが同じような\(R^2\)をもつ場合には、オッカムの剃刀（Ockham’s razor）の指針（ある事柄を説明するためには、必要以上に多くを仮定するべきでない）に従って、より単純なモデルが選ばれるべきである。逆に言うと、回帰モデルに新たな変数を加える場合には、変数を新たに加えるだけでも\(R^2\)が増加することに対して何らかの考慮が必要となる。

　これを実現するための方法の一つが、単回帰のときにも出てきた自由度調整済みの決定係数 \[ R_{adj}^2 = 1-\frac{n-1}{n-p} \frac{SSE}{SST} \] である。この決定係数は、パラメータ数\(p\)が大きくなることにペナルティを課す。

　では、上で推定した回帰モデルについて、これらの係数を計算してみる。

(R2 <- 1 - SSE / SSY)
##           [,1]
## [1,] 0.4077619
(R2adj <- 1 - (n-1) / (n-p) * SSE / SSY)
##          [,1]
## [1,] 0.353922

　なお、単回帰のときと同様に、誤差分散\(\sigma^2\)の推定量はMSEであり、誤差標準偏差の推定量は、その平方根となる。

s <- sqrt(MSE)
s
##          [,1]
## [1,] 2.599246

　これらの係数や推定値は、lmを用いた回帰でも計算されていることを改めて確認しておく。

model <- lm(y ~ ., data = df)
summary(model)
## 
## Call:
## lm(formula = y ~ ., data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7986 -1.0381  0.1364  1.5132  4.0544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  25.3716     9.3883   2.702  0.01301 * 
## x1            1.5414     0.5917   2.605  0.01617 * 
## x2           -5.0390     1.7161  -2.936  0.00764 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.599 on 22 degrees of freedom
## Multiple R-squared:  0.4078, Adjusted R-squared:  0.3539 
## F-statistic: 7.574 on 2 and 22 DF,  p-value: 0.003144

得られた回帰モデルを、\((y, x_1, x_2)\)の3次元の世界で眺めてみる。 以下のコードは複雑だが、コードそのものではなく、描かれる図に注目してほしい。

df <- data.frame(y, x1, x2)
model <- lm(y ~ x1 + x2, data = df)
ax1 <- seq(min(x1), max(x1), 0.5)
ax2 <- seq(min(x2), max(x2), 0.1)
surface <- expand.grid(x1 = ax1, x2 = ax2)
surface$y <- predict(model, newdata = surface)
z <- matrix(surface$y, ncol= length(ax1), byrow = T)
plot_ly(x = ax1, y = ax2, z = z) %>% add_surface() %>% 
  add_trace(x = x1, y = x2, z = y, mode = "markers", type = "scatter3d") %>%
  add_trace(x = x1, y = x2, z = predict(model), mode = "markers", type = "scatter3d")

オレンジ色の点が観察された\((y, x_1, x_2)\)に対応する点、緑色の点が回帰モデルを当てはめたときの値、すなわち、\((\hat y, x_1, x_2)\)に対応する点である。

　2変量の回帰モデルは、\(a x_1 + b x_2 +c = 0\)というかたちをしているので平面となる。つまり、回帰モデルは、回帰平面として描かれる。回帰平面と観察された点の関係を見ると、観察された点から\(y\)軸にそって平面に降ろされた線が回帰平面と交わる点が\(y\)の推定値\(\hat y\)となっていることが分かる。回帰直線が、回帰平面になっただけで、基本的な仕組みは単回帰と同じである。

回帰パラメータと応答の推定

　回帰係数\(\boldsymbol{b}\)の推定量の共分散行列は、 \[ \boldsymbol{s_b}^2 = MSE \times (\boldsymbol{X'X})^{-1} \] と表される。この\(i\)番目の対角要素は、\(b_i\)の分散の推定量であり、位置\((i,j)\)にある非対角要素は、\(b_i\)と\(b_j\)の共分散の推定量である。

　信頼区間を求める、あるいは、ある特定の\(\beta_i\)に対して検定するには、単回帰のときと同様に推定値\(b_i\)と標本分散\(s_{b_i}\)を用いる。単回帰と異なる点は、検定統計量\(t\)が、\(n-2\)ではなく、\(n-p\)の自由度をもつ点である（なお、単回帰では、\(p=2\)であることから、単回帰も重回帰も同じルールに従っていることに注意する）。

(s2 <- MSE[1,1] * solve(t(X) %*% X))
##                        x1           x2
##     88.139822 -4.03102250 -10.84699831
## x1  -4.031022  0.35016054  -0.01754843
## x2 -10.846998 -0.01754843   2.94508500
(sb <- sqrt(diag(s2)))
##                  x1        x2 
## 9.3882811 0.5917436 1.7161250
(tstats <- b / sb)
##         [,1]
##     2.702481
## x1  2.604885
## x2 -2.936259
(pvals <- 2 * pt(-abs(tstats), n - p))
##           [,1]
##    0.013005983
## x1 0.016169568
## x2 0.007641036

　上で行った計算をまとめて表示してみる。

cbind(b, tstats, pvals)
##         [,1]      [,2]        [,3]
##    25.371650  2.702481 0.013005983
## x1  1.541424  2.604885 0.016169568
## x2 -5.038987 -2.936259 0.007641036

　なお、\(\boldsymbol {x_m}=(1, x_1, x_2, ..., x_k)'\)に対する回帰の応答、すなわち、\(y\)の推定値\(\hat y_m = \boldsymbol{x_m' b}\)は、以下の標本分散をもつ。 \[ s_{y_m}^2 = \boldsymbol{x_m}' \boldsymbol{s_b}^2 \boldsymbol x_m \]

　\(y\)の予測値\(\tilde y\)の分散（新たに観察される\(y\)の分散）は、単回帰の場合と同様に、誤差による分散が加わり、\(s_{y_p}^2 = s_{y_m}^2 + MSE\)となる。また、回帰の応答の推論においては、自由度\(n-p\)の\(t\)統計量が用いられる。

　例として、新たな果実が\(x_1 = 12.4\), \(x_2 = 3.52\)のとき、すなわち、共変量\(\boldsymbol{x}_m =(1, 12.4, 3.52)\)をもつ果実について、嗜好性の予測をすることを考える。まずは、この果実に対する推定値\(\hat y_m\)を計算する。

# predicted value
xm <- c(1, 12.4, 3.52)
ym <- xm %*% b
ym
##          [,1]
## [1,] 26.74807

データに含まれる\(Y\)に対して、ヒストグラムで推定された\(\hat y_m\)の大きさ比較してみましょう。

hist(y)
abline(v = ym)

このとき、推定値\(\hat y_m\)、および、予測値\(\tilde y_p\)の分散は、以下のように計算される。

# the variances of mean and individual responses
s2m <- t(xm) %*% s2 %*% xm
s2p <- s2m + MSE[1,1]
(sm <- sqrt(s2m))
##        [,1]
## [1,] 0.7792
(sp <- sqrt(s2p))
##          [,1]
## [1,] 2.713528

これにより、推定値、および、予測値の信頼区間を計算することができる。

# 95 CI's on the mean and individual responses
(cim <- c(ym - qt(0.975, n-p) * sm, ym + qt(0.975, n-p) * sm))
## [1] 25.13211 28.36403
(cip <- c(ym - qt(0.975, n-p) * sp, ym + qt(0.975, n-p) * sp))
## [1] 21.12056 32.37559

実測されている\(Y\)に対して、ヒストグラムで結果を確認してみる。

hist(y, xlim = c(15, 33))
abline(v = ym)
abline(v = cim, col = "orange") # 平均
abline(v = cip, col = "green") #　個々の応答

　以上の一連の解析は関数lmを用いて以下のように実行できる。

df <- data.frame(y, x1, x2)
model <- lm(y ~ ., data = df)
summary(model)
## 
## Call:
## lm(formula = y ~ ., data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7986 -1.0381  0.1364  1.5132  4.0544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  25.3716     9.3883   2.702  0.01301 * 
## x1            1.5414     0.5917   2.605  0.01617 * 
## x2           -5.0390     1.7161  -2.936  0.00764 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.599 on 22 degrees of freedom
## Multiple R-squared:  0.4078, Adjusted R-squared:  0.3539 
## F-statistic: 7.574 on 2 and 22 DF,  p-value: 0.003144

回帰係数を表示する。

coef(model)
## (Intercept)          x1          x2 
##   25.371650    1.541424   -5.038987

回帰モデルの分散分析を行う。

anova(model)
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## x1         1  44.088  44.088  6.5256 0.018072 * 
## x2         1  58.248  58.248  8.6216 0.007641 **
## Residuals 22 148.634   6.756                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

決定係数を表示する。

summary(model)$r.squared
## [1] 0.4077619
summary(model)$adj.r.squared
## [1] 0.353922

回帰係数の信頼区間を表示する。

confint(model)
##                  2.5 %    97.5 %
## (Intercept)  5.9015464 44.841753
## x1           0.3142228  2.768625
## x2          -8.5980125 -1.479962

次に、x1 = 12.4, x2 = 3.52に対するyの推定値\(\hat y\)とその信頼区間、および、予測値\(\tilde y\)の信頼区間の計算を行う。

new.data <- data.frame(x1 = 12.4, x2 = 3.52)
predict(model, new.data, interval = "confidence")
##        fit      lwr      upr
## 1 26.74807 25.13211 28.36403
predict(model, new.data, interval = "prediction")
##        fit      lwr      upr
## 1 26.74807 21.12056 32.37559

参考文献

以下は、本テキストを作成するにあたり参考にした書籍である。