はじめに

　飼育・栽培条件と動植物の生長の関係など、ある変数の変化が別の変数に影響を与える場合がある。このような変数間の関係をモデル化するための統計手法として回帰分析（regression analysis）が挙げられる。変数間の関係を統計的にモデル化することで、変数間に存在する因果関係について理解したり、一方の変数から他方の変数を予測したりすることができるようになる。 　今回の講義では、2つの変数間の関係を直線的な関係としてモデル化する単回帰分析（simple regression analysis）について解説する。単回帰分析では、ある変数\(Y\)を別の変数\(X\)に関連付ける。このとき、数学的には、\(Y\)は\(X\)の関数（a function of \(X\)）であると言うが、回帰分析を行う場合には、その代わりに、「\(Y\)の\(X\)の上への回帰（the regression of \(Y\) on \(X\)）」という言葉がよく用いれる。例えば、栽培試験を行った際の肥料の量と収量の関係を考える場合、収量Yは、肥料の量Xの影響をうけて変化する、と考える。このとき、回帰分析では、これら2つの変数を区別して、\(X\)を独立変数（independent variable）、\(Y\)を従属変数（dependent variable）とよぶ。\(X\)と\(Y\)の関係を回帰式として表すことで、例えば、肥料の量と収量の関係を理解し、さらには、単位施肥量あたりの収量の伸びを求めることも可能となる。なお、\(X\)と\(Y\)の関係は、必ずしも\(X\)が原因で\(Y\)が結果でなくても構わない。ある2つの変量間の関係を回帰式として表し、一方から他方を予測することが目的である場合もある。例えば、人の胴回りのサイズ\(X\)と体重\(Y\)との関係を回帰式で表し、\(X\)から\(Y\)を予測するような場合が考えられる。この場合、\(X\)は\(Y\)の変動の原因ではない。 　回帰は、様々な状況で用いられる。因果関係の理解や、予測の他にも、ある変数の影響で別の変数が影響を受けるのを補正するような利用例も考えられる。例えば、種子のサイズ\(X\)と、そこから発芽した幼植物のサイズ\(Y\)の関係もこうした例の一つである。幼植物に何らかの処理をする実験において、種子のサイズ\(X\)の違いによる幼植物のサイズ\(Y\)への影響を補正できれば、実験結果のより正確な評価が可能になるかもしれない。このように様々な場面で利用される回帰分析について、ここでは、その基礎となる単回帰分析について説明する。

単回帰分析

　施肥量と収量の関係を調べたコムギの栽培試験データ（David 1938, Snedecor and Cochran 1987）を例として用いて、単回帰分析について解説していく。

Nitrogen 40 60 80 100 120 140 160
Yield 15.9 18.8 21.6 25.2 28.7 30.4 30.7

まず、データを準備する。

nitrogen <- c(40, 60, 80, 100, 120, 140, 160)
yield <- c(15.9, 18.8, 21.6, 25.2, 28.7, 30.4, 30.7)
head(ny <- data.frame(nitrogen, yield))
##   nitrogen yield
## 1       40  15.9
## 2       60  18.8
## 3       80  21.6
## 4      100  25.2
## 5      120  28.7
## 6      140  30.4

次に、両者の関係を図示する。

plot(ny)
図1 施肥量（X）とコムギの収量（Y）の関係

図1 施肥量（X）とコムギの収量（Y）の関係

施肥量が多くなるほど収量が高くなる傾向が見てとれる。

次に、収量を施肥量によって説明する単回帰モデルを作成してみる。

model <- lm(yield ~ nitrogen, data = ny)

回帰分析の結果（推定されたモデル）は、modelに代入されている。

回帰分析の結果を表示させるには関数summaryを用いる。

summary(model)
## 
## Call:
## lm(formula = yield ~ nitrogen, data = ny)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -0.5679 -0.3357 -0.2036  0.7286  1.5607  0.5929 -1.7750 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.13214    1.19945   9.281 0.000244 ***
## nitrogen     0.13339    0.01114  11.978 7.15e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.179 on 5 degrees of freedom
## Multiple R-squared:  0.9663, Adjusted R-squared:  0.9596 
## F-statistic: 143.5 on 1 and 5 DF,  p-value: 7.153e-05

上のコマンドを実行して表示された結果について順に説明していく。

Call:
lm(formula = yield ~ nitrogen, data = ny)

これは先ほど入力したコマンドが繰り返されたものである。入力した直後にこの出力が得られても、有用な情報でないように思われるが、後で述べるように複数の回帰モデルを作って比較をする場合などには、どのようなモデルを想定して得られた結果であるかを再確認する必要がある。

なお、ここでは、\(i\)番目の肥料水準における収量を\(y_i\)、その水準での施肥量を\(x_i\)として、 \[ y_i = \mu + \beta x_i + \epsilon \] というモデルを想定して回帰分析を行っている。先述したように、\(x_i\)のことを独立変数（independent variable）または説明変数（explanatory variable）、\(y_i\)のことを従属変数（dependent variable）または応答変数（response variable）とよぶ。\(\mu\)や\(\beta\)を回帰モデルのパラメータ（parameter）、\(\epsilon_i\)を誤差（error）とよぶ。また、\(\mu\)を母切片（population intercept）、\(\beta\)を母回帰係数（population regression coefficient）とよぶ。

　回帰モデルのパラメータ\(\mu\)や\(\beta\)は、標本から推定されるが、その推定値を、それぞれ、標本切片（sample intercept）および標本回帰係数（sample regression coefficient）とよぶ。標本から推定された\(\mu\)、\(\beta\)の値を、以降、それぞれ、\(m\)、\(b\)で表す。 　

Residuals:
      1       2       3       4       5       6       7 
-0.5679 -0.3357 -0.2036  0.7286  1.5607  0.5929 -1.7750

この出力は、残差を表している。データ数が多い場合には、残差の最大値、最小値、四分位点が表示される。この値をもとに、残差の大きさのチェックをすることができる。例えば、上述した単回帰モデルでは誤差の期待値（平均）は0となることを想定しているが、中央値（median）がそこから大幅にはずれていないか確認することができる。また、誤差の最大値と最小値、または、25%点と75%点がほぼ同じ値をとっているかどうかを確認することで、誤差が0を中心として左右対称の分布をしているかを確認できる。この例では、データ数が少ないために、個々の標本における残差が示されている。

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.13214    1.19945   9.281 0.000244 ***
nitrogen     0.13339    0.01114  11.978 7.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

回帰モデルのパラメータ\(\mu\)、\(\beta\)の推定値\(m\)、\(b\)と、その標準誤差、\(t\)値、\(p\)値が表示されている。

Residual standard error: 1.179 on 5 degrees of freedom
Multiple R-squared:  0.9663,    Adjusted R-squared:  0.9596 
F-statistic: 143.5 on 1 and 5 DF,  p-value: 7.153e-05

最初の行は、残差の標準偏差を表している。これは、誤差分散\(\sigma^2\)の推定値を\(s^2\)とすると、\(s\)で表される値である。

　2行目は、決定係数R^2である。また、補正\(R^2\)（Adjusted R-squared）は、自由度調整済み決定係数とよばれる統計量である。いずれも回帰が説明する変動の度合いを表しており、高いほうが回帰の説明力が高い。

　3行目は、回帰モデルの有意性を表す\(F\)検定の結果である。回帰係数\(\beta\)が0であるという帰無仮説（\(H_0: \beta = 0\)）のもとでの検定であり、この\(p\)値が非常に小さい場合には、帰無仮説\(H_0\)を棄却して対立仮説（\(H_1: \beta <, \ne, > 0\)）を採択すべきであると解釈される。

回帰分析の結果を図示して眺めてみる。まず、散布図を描き、そこに回帰直線を引く。

plot(ny$yield ~ ny$nitrogen, ylim = c(15, 35))
abline(model, col = "red")
図2. 散布図に回帰直線を加えた図

図2. 散布図に回帰直線を加えた図

次に、回帰モデルにデータをあてはめた場合のyの値を計算し、図示する。

yield.fit <- fitted(model) # calculate fitted values
plot(ny$yield ~ ny$nitrogen, ylim = c(15, 35))
abline(model, col = "red")
points(ny$nitrogen, yield.fit, pch = 3, col = "green")
図3. モデルをあてはめて計算されるyの値は全て直線上に乗る

図3. モデルをあてはめて計算されるyの値は全て直線上に乗る

観察値yは、回帰モデルで説明される部分（モデルを当てはめたときの値）と、回帰で説明されない残差の和として表される。そこで、残差を図示して、その関係を確認してみる。

plot(ny$yield ~ ny$nitrogen, ylim = c(15, 35))
abline(model, col = "red")
points(ny$nitrogen, yield.fit, pch = 3, col = "green") 
segments(ny$nitrogen, yield.fit, 
         ny$nitrogen, yield.fit + resid(model), col = "gray")
図4. yの値は、モデルをあてはめて計算されるyの値（緑の点）とモデルの残差（灰色の線）の和として表される

図4. yの値は、モデルをあてはめて計算されるyの値（緑の点）とモデルの残差（灰色の線）の和として表される

　では、観察されていない\(x（50, 70, …, 150）\)に対して、回帰モデルを用いて\(y\)を予測してみる。

plot(ny$yield ~ ny$nitrogen, ylim = c(15, 35))
abline(model, col = "red")
points(ny$nitrogen, yield.fit, pch = 3, col = "green") 
segments(ny$nitrogen, yield.fit, 
         ny$nitrogen, yield.fit + resid(model), col = "gray")
yield.pred <- predict(model, data.frame(nitrogen = seq(50, 150, 20)))
points(seq(50, 150, 20), yield.pred, pch = 2, col = "blue")
図5. 予測値（青三角）は全て回帰直線の上に乗る

図5. 予測値（青三角）は全て回帰直線の上に乗る

単回帰の統計モデル

　単回帰モデルでは、\(n\)組の観察データ\((x_1,y_1),\dots,(x_n,y_n)\)に対して、\(y_i\)を次のような\(x_i\)の線形関数と誤差からなるとモデル化する。 \[ y_i=\mu+\beta x_i+\epsilon_i \]

ここで、誤差はデータ間で独立していて、それぞれ、平均0、分散\(\sigma^2\)の正規分布にしたがっていると考える。すなわち、 \[ \epsilon_i \sim N(0,\sigma^2) \] したがって、 \[ y_i\sim N(\mu+\beta x_i,\sigma^2) \] となる。

では、\(\mu = 0, \beta = 1, \sigma^2 = 2\)のときに従う分布を図示してみる。

mu <- 0
beta <- 1
sigma2 <- 2
x <- 1:10
epsilon <- rnorm(10, sd = sqrt(sigma2))
y <- mu + beta * x + epsilon
plot(x, y, ylim = c(min(y) - 5, max(y) + 5), col = "blue")
abline(mu, beta)
s <- seq(-3, 3, 0.1)
t <- dnorm(s, 0, sqrt(sigma2))
for(i in 1:length(x)) {
  lines(x[i] - t, mu + beta * x[i] + s, col = "green")
}
model <- lm(y ~ x)
lines(x, fitted(model), col = "red", lty = "dashed")

　ここで、黒線が単回帰の統計モデルにおける\(x\)と\(\mu + \beta x_i\)の関係であり、緑の曲線が\(y\)が従う正規分布\(N(\mu + \beta x_i, \sigma^2)\)である。また、赤の破線は、推定されたパラメータ（母数）に基づく回帰、すなわち、推定された回帰直線を表す。次の節では、どのようにしてパラメータが推定されるかを説明する。

パラメータ（母数）の推定

　ここでは、回帰パラメータの推定法について解説する。

　先述したように単回帰のモデルは、 \[ y_i=\mu+\beta x_i+\epsilon_i \] として表現される。この式は、観察値\(y_i\)が、回帰方程式で説明される部分\(\mu+\beta x_i\)と、回帰方程式では説明されない誤差\(\epsilon_i\)から成ることを意味している。\(y_i\)は与えられているので、上式の、\(\mu\)や\(\beta\)に具体的な値\(m\)と\(b\)を代入すると、それに伴って誤差\(\epsilon_i\)の具体的な値\(e_i\)は変化する。このとき、どのようにして母切片\(\mu\)と母回帰係数\(\beta\)の推定値を求めればよいだろうか？

　何をもって「最適」とするかについては様々な基準が考えられるが、ここでは、誤差\(e_i\)の値を小さくすることを考えてみる。誤差\(e_i\)は正負両方の値をとるため、単純に和をとると互いに相殺されてしまう。そこで、\(e_i\)の2乗和（sum of squared error: SSE）を最小にすることを考える。すなわち、 \[ SSE=\sum_{i=1}^n e_i^2 =\sum_{i=1}^n \left[y_i-(m+b x_i)\right]^2 \] を最小にするような\(m\)と\(b\)を考えてみる。

\(SSE\)を最小化する\(m\)は、 \[ \frac {\partial SSE}{\partial m}=-2 \sum_{i=1}^n (y_i-m-b x_i) = 0 \\ \Leftrightarrow　\sum_{i=1}^n y_i -nm -b \sum_{i=1}^n x_i =0 \\ \Leftrightarrow　m = \frac {\sum_{i=1}^n y_i} n - b \frac {\sum_{i=1}^n x_i} n =\bar y-b\bar x \] として計算される。

また、\(SSE\)を最小化する\(b\)は、 \[ \frac {\partial SSE}{\partial b}=-2 \sum_{i=1}^n x_i(y_i-m-b x_i) = 0 \\ \Leftrightarrow \sum_{i=1}^n x_iy_i-m \sum_{i=1}^n x_i - b \sum_{i=1}^n x_i^2 =0\\ \Leftrightarrow \sum_{i=1}^n x_iy_i- n(\bar y- b \bar x)\bar x- b \sum_{i=1}^n x_i^2 =0\\ \Leftrightarrow \sum_{i=1}^n x_iy_i - n\bar x \bar y- b \left(\sum_{i=1}^n x_i^2 -n\bar x ^2\right)=0\\ \Leftrightarrow b = \frac {\sum_{i=1} n x_i y_i - n\bar x \bar y}{\sum_{i=1} n x_i^2 - n\bar x^2}= \frac {SSXY}{SSX} \] ここで、\(SSXY\)と\(SSX\)は、\(x\)と\(y\)の偏差積和と\(x\)の偏差平方和で、それぞれ、 \[ SSXY = \sum_{i=1}^n (x_i - \bar x)(y_i - \bar y) \\ = \sum_{i=1}^n x_iy_i - \bar x \sum_{i=1}^n y_i - \bar y \sum_{i=1}^n x_i + n\bar x \bar y\\ = \sum_{i=1}^n x_iy_i - n\bar x \bar y - n\bar y\bar x + n\bar x\bar y\\ = \sum_{i=1}^n x_iy_i - n\bar x \bar y \] \[ SSX = \sum_{i=1}^n (x_i - \bar x)^2 \\ = \sum_{i=1}^n x_i^2 - 2 \bar x \sum_{i=1}^n x_i + n\bar x \\ = \sum_{i=1}^n x_i^2 - 2n\bar x^2 + n\bar x^2 \\ = \sum_{i=1}^n x_i^2 - n\bar x^2 \] として計算される。

\(SSE\)を最少にする\(m\)と\(b\)をこれらパラメータの推定値とし、以降、単に\(m\)と\(b\)で表すことにする。すなわち、 \[ b=\frac {SSXY}{SSX} \\ m=\bar y- b \bar x \]

では、回帰係数を上述した式をもとにして計算してみよう。まずは、\(x\)と\(y\)を設定し、それらの偏差積和と偏差平方和を計算する。

x <- ny$nitrogen
y <- ny$yield
n <- length(x)          #　サンプル数をnに代入
ssxy <- sum(x * y) - n * mean(x) * mean(y)      #　偏差積和
ssx <- sum(x^2) - n * mean(x)^2         #　偏差平方和

傾き\(b\)を計算する。

(b <- ssxy / ssx)
## [1] 0.1333929

次に切片\(m\)を計算する。

(m <- mean(y) - b * mean(x))
## [1] 11.13214

計算された\(m\)と\(b\)をもとに回帰直線を描いてみる。

plot(y ~ x, ylim = c(15, 35))
abline(m, b, col = "red")

先ほど関数lmを用いて計算された回帰直線と同じものが描かれる。

なお、回帰パラメータが推定されれば、与えられた\(x_i\)に対応する\(y\)の推定値\(\hat y _i\)を計算することができるようになる。すなわち、 \[ \hat y_i=m+b x_i \] として計算できる。これにより、観察された\(x\)にモデルをあてはめたときの\(y\)の(\(\hat y\))値を計算したり、\(x\)のみが既知の場合に\(y\)を予測したりすることができる。ここでは、観察された\(x\)にモデルをあてはめたときの\(y\)の値を計算し、先ほど描いた図の上に点を散布してみる。

plot(y ~ x, ylim = c(15, 35))
abline(m, b, col = "red")
y.hat <- m + b * x
points(x, y.hat, pch = 3, col = "green")

次に、観察された\(x\)に対してモデルをあてはめたときの\(y\)の値を計算して、観察された\(y\)とモデルをあてはめたときの\(y\)間の関係を表す散布図を描いてみる。

lim <- range(c(y, y.hat))
plot(y, y.hat, xlab = "Observed", ylab = "Fitted", xlim = lim, ylim = lim)
abline(0, 1, col = "gray")
図7. 観測値とあてはめ値の間の関係

図7. 観測値とあてはめ値の間の関係

観察値とあてはめた値の一致の度合いを調べるために両者の相関係数を計算してみる。

cor(y, y.hat)
## [1] 0.9830173

実は、この相関係数の2乗が、回帰が説明する\(y\)の変動の割合（決定係数、\(R^2\)値）になっている。両者を見比べてみると、

cor(y, y.hat)^2
## [1] 0.966323

さきほど、lm関数を使って計算された\(R^2\)と同じ値となる。

Multiple R-squared:  0.9663,    Adjusted R-squared:  0.9596

回帰モデルの有意性検定

　変数間の直線的な関係が強い場合には回帰直線がよくあてはまり、両変数間の関係を回帰直線でうまくモデル化できる。しかし、変数間の直線的な関係が明瞭でない場合には、回帰直線によるモデル化がうまく行かない。ここでは、推定された回帰モデルの有効性を客観的に確認するための方法として、分散分析を用いた検定法について説明する。

まずは、再度、単回帰を行ってみる。

model <- lm(yield ~ nitrogen, data = ny)

得られた回帰モデルの有意性は、関数anovaを用いて検定できる。

anova(model)
## Analysis of Variance Table
## 
## Response: yield
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## nitrogen   1 199.289 199.289  143.47 7.153e-05 ***
## Residuals  5   6.945   1.389                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

分散分析の結果、変数nitrogenの項は高度に有意（\(p < 0.001\)）であり、施肥の量nitrogenが収量yieldに影響を与えるという回帰モデルの有効性が確認できる。

回帰モデルの分散分析では、以下に示すような計算が行われる。まず、「回帰で説明される平方和」（回帰モデルをあてはめて計算される値\(\hat y_i\)の偏差平方和）は、以下のようにして計算できる。 \[ SSR=\sum_{i=1}^n (\hat y_i-\bar y)^2 \\ =\sum_{i=1}^n \left[\mu+bx_i-(\mu+\bar bx)\right]^2 \\ =b^2 \sum_{i=1}^n (x_i-\bar x)^2 \\ =b^2\cdot SSX=b \cdot SSXY \]

また、観察値\(y\)の平均からの偏差の平方和は、回帰で説明される平方和\(SSR\)と残差平方和\(SSE\)の和として表される。すなわち、 \[ SSY=\sum_{i=1}^n (y_i-\bar y)^2 \\ =\sum_{i=1}^n (y_i-\hat y_i+\hat y_i-\bar y)^2 \\ =\sum_{i=1}^n (y_i-\hat y_i )^2 +\sum_{i=1}^n (\hat y_i-\bar y)^2 \\ =SSE+SSR \]

これは、以下の関係が成り立つことによる。 \[ 2\sum_{i=1}^n (y_i-\hat y_i)(\hat y_i-\bar y) \\ =2\sum_{i=1}^n (y_i-m-bx_i )\left[m+bx_i-(m+bx ̅ )\right] \\ =2b\sum_{i=1}^n \left[y_i-(\bar y-b\bar x)-bx_i\right](x_i-\bar x) \\ =2b\sum_{i=1}^n \left[y_i-y ̅-b(x_i-x ̅ )\right](x_i-\bar x) \\ =2b(SSXY-b\cdot SSX)=0 \]

　では、上の式を用いて実際に計算してみる。まずは、回帰で説明される平方和\(SSR\)と残差平方和\(SSE\)を計算する。

(ssr <- b * ssxy)
## [1] 199.2889
(ssy <- sum(y^2) - n * mean(y)^2)
## [1] 206.2343
(sse <- ssy - ssr)
## [1] 6.945357

　次に、平方和を自由度で割った平均平方を計算する。

(msr <- ssr / 1)
## [1] 199.2889
(mse <- sse / (n - 2))
## [1] 1.389071

　最後に回帰の平均平方を誤差の平均平方で割り、\(F\)値を計算する。さらに、計算された\(F\)値に対応する\(p\)値を計算する。

(f.value <- msr / mse)
## [1] 143.4692
1 - pf(f.value, 1, n - 2)
## [1] 7.153346e-05

得られる結果は、先ほど関数anovaを用いて計算された結果と一致する。

なお、回帰の分散分析の結果は、関数summaryを用いて表示される回帰分析の結果の中にも含まれている。

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.13214    1.19945   9.281 0.000244 ***
nitrogen     0.13339    0.01114  11.978 7.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.179 on 5 degrees of freedom
Multiple R-squared:  0.9663,    Adjusted R-squared:  0.9596 
F-statistic: 143.5 on 1 and 5 DF,  p-value: 7.153e-05

「Residual standard error」は、残差の平均平方の平方根となっている。

sqrt(mse)
## [1] 1.178589

「Multiple R-squared」（\(R^2\)）は、決定係数（coefficient of determination）とよばれる値で、\(SSR\)と\(SSY\)の比である。

ssr / ssy
## [1] 0.966323

「Adjusted R-squared」（\(R_{adj}^2\)）は、自由度調整済決定係数とよばれる値で、次のように計算できる。

(ssy / (n - 1) - mse) / (ssy / (n - 1))
## [1] 0.9595876

また、「F-statistic」は、分散分析でnitrogenの効果として表されている\(F\)値とその\(p\)値に一致する。また、nitrogenの回帰係数について計算されている\(t\)値を2乗すると\(F\)値になる（9.2812 = 143.5）。

なお、\(R^2\)および\(R_{adj}^2\)は、\(SSR\)、\(SSY\)、\(SSE\)を用いて以下のように表すこともできる。 \[ R^2=\frac {SSR}{SSY}=1-\frac {SSE}{SSY} \\ R_{adj}^2=1-\frac {n-1}{n-p}\cdot\frac {SSE}{SSY} \] ここで、\(p\)はモデルに含まれるパラメータの数で、単回帰モデルでは、\(p = 2\)になります。\(R_{adj}^2\)は、モデルに含まれるパラメータの数が多ければ多いほど、調整量が大きくなる（残差平方和の小ささを低く見積もる）ことが分かる。

回帰係数の推定値が従う分布

　先述したように回帰係数\(\mu\)と\(\beta\)の推定値\(m\)と\(b\)は、標本から推定される値であり、偶然選ばれた標本に左右される確率変数である。したがって、推定値\(m\)と\(b\)は確率分布をもつ。ここでは、推定値の従う分布について考える。

まず、\(b\)について考える。\(b\)は、 \[ b=\sum_{i=1}^n \frac {(x_i-\bar x)(y_i-\bar y)}{SSX}\\ =\sum_{i=1}^n \frac{(x_i-x ̅ ) y_i}{SSX}-\bar y \sum_{i=1}^n \frac{(x_i-\bar x)}{SSX} \\ =\frac 1 {SSX} \sum_{i=1}^n y_i(x_i-\bar x) \] と表すことができる。

したがって、 推定値\(b\)の平均は、 \[ E(b)=\frac 1 {SSX} E\left(\sum_{i=1}^n y_i (x_i-\bar x)\right) \\ =\frac 1 {SSX} E\left(\sum_{i=1}^n(\mu +\beta x_i+\epsilon_i)(x_i-\bar x) \right) \\ =\frac 1 {SSX} E\left(\sum_{i=1}^n (\mu^*+\beta(x_i-\bar x)+\epsilon_i)(x_i-\bar x)\right) \\ =\frac 1 {SSX} \left[μ^* \sum_{i=1}^n (x_i-\bar x) +\beta \sum_{i=1}^n (x_i-\bar x)^2 +E\left(\sum_{i=1}^n \epsilon_i (x_i-\bar x)\right)\right] \\ =\frac 1 {SSX} [0+\beta SSX+0]= \beta \] となる。なお、回帰では、\(x_i\)は与えられた値であり、定数で確率変数ではないことに注意する。

このように、推定値\(b\)の平均は、真の値\(\beta\)に一致する。ここで、\(\mu^*\)は、\(y_i\)を\(x_i\)でなく\(x_i-\bar x\)に回帰した場合の定数項で、 \[ y_i=\mu^*+\beta (x_i-\bar x) \] と表される。

推定値\(b\)の分散は、 \[ V(b)=\frac 1 {SSX^2} V\left(\sum_{i=1}^n y_i (x_i-\bar x)\right) \\ =\frac 1 {SSX^2} \sum_{i=1}^n (x_i-\bar x)^2 V(y_i) \\ =\frac {\sum_{i=1}^n (x_i-\bar x)^2} {SSX^2} \sigma^2=\frac {\sigma^2}{SSX} \] となる。なお、ここで\(\sigma^2\)は、残差分散\(\sigma^2=V(y_i)=V(e_i)\)である。

なお、推定値\(b\)は、 \[ a_i=\frac {x_i-\bar x}{SSX} \] とおくと、\(y_i\)の線形結合 \[ b=\sum_{i=1}^n a_i y_i \] である。ここで、\(y_i\)は正規分布に従うため、その線形結合である\(b\)もまた正規分布に従う。すなわち、 \[ b \sim N\left(\beta,\frac {\sigma^2}{SSX}\right) \]

いっぽう、推定値\(m\)は、\(m=\bar y-b\bar x\)と表せる。

したがって、平均は、 \[ E(m)=E(\bar y-b \bar x)=E(\bar y)- E(b)\bar x=\mu+ \beta \bar x- \beta \bar x=μ \] となる。推定値\(m\)の平均も、やはり、真の値\(\mu\)に一致する。

次に分散は、 \[ m = \bar y - b\bar x = \sum_{i=1}^n \frac{y_i}{n} - \sum_{i=1}^n \frac{y_i(x_i-\bar x)}{SSX} \bar x \\ = \sum_{i=1}^ny_i\left[\frac{1}{n}-\frac{(x_i-\bar x)\bar x}{SSX}\right] \] したがって、 \[ V(m) = V(y) \sum_{i=1}^n \left[\frac{1}{n}-\frac{(x_i-\bar x)\bar x}{SSX}\right]^2\\ \] となる。ここで、 \[ \sum_{i=1}^n \left[\frac{1}{n}-\frac{(x_i-\bar x)\bar x}{SSX}\right]^2\\ = \sum_{i = 1}^n \left[\frac{1}{n^2}-\frac{2(x_i-\bar x)\bar x}{n SSX} + \frac{(x_i-\bar x)^2\bar x^2}{SSX^2}\right] = \frac{1}{n} + 0 + \frac{\bar x^2}{SSX} \] と整理できるので、 \[ V(m) = \sigma^2 \left(\frac{1}{n} + \frac{\bar x^2}{SSX}\right) \] となる。

\(y_i\), \(b\)は正規分布に従うため、\(m=\bar y-b \bar x\)と表される\(m\)もまた正規分布に従う。すなわち、 \[ m \sim N\left(\mu ,\sigma^2 \left(\frac 1 n + \frac {\bar x^2}{SSX}\right)\right) \]

なお、誤差分散\(\sigma^2\)の真の値は未知だが、これを残差分散\(s^2\)で置き換えることができる。すなわち、 \[ s^2 = \frac {SSE}{n-2} \] この値は分散分析の際に計算した残差の平均平方である。

このとき、\(b\)に関する統計量 \[ t = \frac {b-b_0}{s⁄\sqrt{SSX}} \] は、帰無仮説 \[ H_0: \beta = b_0 \] のもとで、自由度\(n – 2\)の\(t\)分布に従う。

このとき、母回帰係数\(\beta\)が\(1-\alpha\)の確率で含まれる区間、すなわち、\((1-\alpha)100\)%信頼区間は以下のように計算される。 \[ \left[b-t_{n-2,1-\alpha/2} \frac s {\sqrt{SSX}},\: b+t_{n-2,1-\alpha/2} \frac s {\sqrt{SSX}}\right] \] ここで、\(t_{n-2,1-\alpha/2}\)は自由度\(n-2\)における両側\(\alpha\)水準の棄却限界値である。

また、\(m\)についても統計量 \[ t = \frac {m-m_0}{s\sqrt{\frac 1 n +\frac {\bar x^2} {SSX}}} \] は、帰無仮説 \[ H_0: \mu=m_0 \] のもとで、自由度\(n – 2\)の\(t\)分布に従う。

\(\mu\)が\(1-\alpha\)の確率で含まれる区間、すなわち、\((1-\alpha)100\)%信頼区間は以下のように計算される。 \[ \left[m-t_{n-2,1-\alpha/2} \;s \sqrt {\frac 1 n + \frac {\bar x^2}{SSX}},\: m+t_{n-2,1-\alpha/2} \;s \sqrt {\frac 1 n + \frac {\bar x^2}{SSX}}\right] \] では、ここまでに推定値を求めてきた\(\beta\)と\(\mu\)について、検定と信頼区間の計算を行ってみよう。

まず母回帰係数\(\beta\)について、帰無仮説\(H_0: \beta=0\)について検定してみる。

(t <- (b - 0) / (sqrt(mse / ssx)))
## [1] 11.97786
2 * pt(-abs(t), n - 2)   # two tailed test
## [1] 7.153346e-05

この検定の結果は、回帰分析結果として既に表示されていたものである。

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.13214    1.19945   9.281 0.000244 ***
nitrogen     0.13339    0.01114  11.978 7.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

なお、自由度\(n-2\)の\(t\)分布について97.5%点は

qt(0.975, n - 2)
## [1] 2.570582

であり、上で計算した\(t\)値（11.98）のほうが大きいことが分かる。

次に、95%信頼区間は、

half.int <- qt(0.975, n - 2) * sqrt(mse / ssx)
c(b - half.int, b + half.int)
## [1] 0.1047653 0.1620204

と計算される。この範囲に、母回帰係数\(\beta\)が95%の確率で含まれることを意味する。

上で行った仮説検定は、任意の\(\beta_0\)について行うことができます。例えば、帰無仮説\(H_0: \beta=0.1\)について検定するには、以下のようにする。

(t <- (b - 0.1) / (sqrt(mse / ssx)))
## [1] 2.998474
2 * pt(-abs(t), n - 2) # two tailed test
## [1] 0.03015207

結果は、5%水準で有意である。これは、先ほどの95%信頼区間に0.1が「含まれていない」ことと同じことを意味する。

では、次に母切片\(\mu\)についても検定と信頼区間の計算を行う。まず、帰無仮説\(H_0: \mu=0\)の検定を行う。

# test mu = 0
(t <- (m - 0) / sqrt(mse * (1/n + mean(x)^2 / ssx)))
## [1] 9.281037
2 * pt(-abs(t), n - 2)   # two tailed test
## [1] 0.0002442102

この結果もやはり、既に計算されていた。

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.13214    1.19945   9.281 0.000244 ***
nitrogen     0.13339    0.01114  11.978 7.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

次に95%信頼区間を計算してみよう。

half.int <- qt(0.975, n - 2) * sqrt(mse * (1/n + mean(x)^2 / ssx))
(int <- c(m - half.int, m + half.int))
## [1]  8.048858 14.215428

母切片\(\mu\)が95%の確率でこの区間に含まれるということを意味する。

最後に、帰無仮説\(H_0: \mu=10\)について検定してみよう。

(t <- (m - 10) / sqrt(mse * (1/n + mean(x)^2 / ssx)))
## [1] 0.9438847
2 * pt(-abs(t), n - 2)   # two tails test
## [1] 0.388568

結果は、5%水準でも有意ではない。これは、95%信頼区間に10が含まれることと同じことを意味する。

なお、母回帰係数と母切片の信頼区間は、次のようにするとRの専用関数を用いて計算することができます。

model <- lm(yield ~ nitrogen, data = ny)
confint(model) # 95%
##                 2.5 %     97.5 %
## (Intercept) 8.0488576 14.2154281
## nitrogen    0.1047653  0.1620204
confint(model, level = 0.99) # 99%
##                  0.5 %     99.5 %
## (Intercept) 6.29578758 15.9684981
## nitrogen    0.08848843  0.1782973

推定値、予測値の信頼区間

関数predictには様々な機能がある。まずは回帰モデルを単純に引数として関数を使ってみる。するとモデルをあてはめたときのyの値が計算される。その値は関数fittedで計算されるものと全く同じである。

predict(model) # predictをそのまま使うと、モデルを当てはめたときのyの値を表示
##        1        2        3        4        5        6        7 
## 16.46786 19.13571 21.80357 24.47143 27.13929 29.80714 32.47500
fitted(model) # モデルを当てはめたときのyの値を表示
##        1        2        3        4        5        6        7 
## 16.46786 19.13571 21.80357 24.47143 27.13929 29.80714 32.47500

オプションintervalとlevelを設定するとモデルをあてはめたときの\(y\)の値の信頼区間を計算できます。

predict(model, interval = "confidence", level = 0.95)
##        fit      lwr      upr
## 1 16.46786 14.40349 18.53222
## 2 19.13571 17.51629 20.75514
## 3 21.80357 20.52331 23.08384
## 4 24.47143 23.32633 25.61653
## 5 27.13929 25.85902 28.41955
## 6 29.80714 28.18772 31.42656
## 7 32.47500 30.41064 34.53936

関数predictを用いて\(y\)の信頼区間を図示してみる。

new.x <- data.frame(nitrogen = 40:160)
conf.int <- predict(model, interval = "confidence", level = 0.95, newdata = new.x)
head(conf.int)
##        fit      lwr      upr
## 1 16.46786 14.40349 18.53222
## 2 16.60125 14.56064 18.64186
## 3 16.73464 14.71767 18.75162
## 4 16.86804 14.87456 18.86151
## 5 17.00143 15.03132 18.97154
## 6 17.13482 15.18794 19.08171
plot(ny$yield ~ ny$nitrogen)
matlines(new.x$nitrogen, conf.int, lty = c(1, 2, 2), col = "red")
図8. yの信頼区間。xの平均付近は狭く、そこから離れるほど広くなる

図8. yの信頼区間。xの平均付近は狭く、そこから離れるほど広くなる

\(y\)の信頼区間は次のように計算できる。まず、\(x^*\)を与えたときの\(y\)、すなわち、\(y_m=E(y│x=x^* )=\mu+\beta x^*\)を推定することを考える。標本から推定された回帰係数を\(b\)とすると、\(y_m\)の推定値は、\(\hat y_m = m+bx^*\)となる。ここで、\(m\)も\(b\)も確率変数であるために、\(\hat y _m\)もまた確率変数となります。\(\hat y_m\)は \[ \hat y_m=m+bx^*=\bar y+b(x^*-\bar x) \] と表され、その分散は以下のように計算される。 \[ V(\hat y_m)=V(\bar y)+(x^*-\bar x)^2 V(b) \\ = \frac{\sigma^2} n +\frac{(x^*-\bar x)^2 \sigma^2} {SSX} \\ = \sigma^2 \left[\frac 1 n + \frac {(x^*-\bar x)^2}{SSX}\right] \] 残差分散\(\sigma^2\)について標本から計算される残差分散\(s^2\)と置き換えると、先ほどと同様に、統計量 \[ t = \frac{\hat y_m-y_m}{s\sqrt{\frac 1 n +\frac {(x^*-\bar x)^2}{SSX}}} \] は、自由度\(n – 2\)の\(t\)分布に従う。

したがって、真の値\(y_m=\mu+\beta x^*\)が\(1-\alpha\)の確率で含まれる区間、すなわち、\((1-\alpha)100\)%信頼区間は以下のように計算される。 \[ \left[\hat y_m-t_{n-2,1-\alpha /2} \; s \sqrt{\frac 1 n + \frac {(x^*-\bar x)^2}{SSX}}, \hat y_m+t_{n-2,1-\alpha /2} \; s \sqrt{\frac 1 n + \frac {(x^*-\bar x)^2}{SSX}}\right] \]

では、上式にしたがって、\(x^*\)をあたえたときの\(y\)の推定値の信頼区間を図示してみる。

x <- 40:160
y.hat <- m + b * x
mean.x <- mean(ny$nitrogen)
half.int <- qt(0.975, n - 2) * sqrt(mse) * sqrt(1/n + (x - mean.x)^2 / ssx)
y.hat.lower <- y.hat - half.int
y.hat.upper <- y.hat + half.int
plot(ny$yield ~ ny$nitrogen)
matlines(x, cbind(y.hat, y.hat.lower, y.hat.upper), lty = c(1, 2, 2), col = "red")

図8と同じ図が描かれる。

　なお、上述したのは\(x^*\)をあたえたときのyの推定値についての信頼区間であり、新たな\(y\)を予測した場合の信頼区間ではない。今、\(y\)の予測値を\(y_{pred}\)とすると、\(y_{pred}\)の変動の要因は2つあり、そのうちの1つは、誤差分散\(\sigma^2\)によるもの、もう一つは、\(m+bx^*\)の分散\(\sigma^2/n+((x^*-\bar x)^2 \sigma^2)/SSX\)によるものである。したがって、予測値\(\hat y_{pred}\)の分散は、 \[ V(\hat y_{pred})=σ^2 \left[1+\frac 1 n +\frac {(x^*-\bar x)^2}{SSX}\right] \] で、\(\hat y_{pred}\)は正規分布 \[ \hat y_{pred} \sim N(\mu+ \beta x^*,\sigma^2 \left[1+\frac 1 n + \frac {(x^*-\bar x )^2}{SSX}\right] \] に従う。

残差分散\(\sigma^2\)について標本から計算される残差分散\(s^2\)と置き換えると、先ほどと同様に、統計量 \[ t = \frac {\hat y_{pred}-y_{pred}}{s\sqrt {1+\frac 1 n +\frac {(x^*-x ̅ )^2}{SSX}}} \] は、自由度\(n – 2\)の\(t\)分布に従う

したがって、予測値\(y_{pred}\)が\(1-\alpha\)の確率で含まれる区間、すなわち、\((1-\alpha)100\)%予測区間は以下のように計算される。 \[ \left[\hat y_{pred}-t_{n-2,1-\alpha/2}\;s \sqrt{1+\frac 1 n + \frac{(x^*-\bar x)^2} {SSX}},\: \hat y_{pred}+t_{n-2,1-\alpha/2}\;s \sqrt{1+ \frac 1 n + \frac{(x^*-\bar x)^2}{SSX}}\right] \]

では、上式にしたがって、\(x^*\)をあたえたときの\(y\)の予測値（\(y\)の観察値が取るであろう値）と予測区間を図示してみよう。

x <- 40:160
y.hat <- m + b * x
mean.x <- mean(ny$nitrogen)
half.int <- qt(0.975, n - 2) * sqrt(mse) * sqrt(1 + 1/n + (x - mean.x)^2 / ssx)
y.hat.lower <- y.hat - half.int
y.hat.upper <- y.hat + half.int
plot(ny$yield ~ ny$nitrogen)
matlines(x, cbind(y.hat, y.hat.lower, y.hat.upper), lty = c(1, 2, 2), col = "green")
図9. yの信頼区間（赤）のyの予測区間（緑)

図9. yの信頼区間（赤）のyの予測区間（緑)

これは、関数predictを用いて以下のように計算できる。

conf.int <- predict(model, interval = "predict", level = 0.95, newdata = new.x)
head(conf.int)
##        fit      lwr      upr
## 1 16.46786 12.80174 20.13398
## 2 16.60125 12.94846 20.25404
## 3 16.73464 13.09500 20.37429
## 4 16.86804 13.24136 20.49471
## 5 17.00143 13.38754 20.61531
## 6 17.13482 13.53354 20.73610
plot(ny$yield ~ ny$nitrogen)
matlines(new.x$nitrogen, conf.int, lty = c(1, 2, 2), col = "green")

参考文献

以下は、本テキストを作成するにあたり参考にした書籍である。より詳しく学ぶためにもお薦めです。