連続変数を組み込んだ解析

　これまでの検定や分散分析の項では、説明変数はカテゴリカル変数だけを使ってきた。今回は、説明変数が連続変数であるような解析を取り扱う。連続的に変化するデータをそのまま解釈するのは難しい。標本の平均値のように、平均的な連続的変化を1本の直線や曲線で表現して、解釈したり予測するほうが、易しいだろう。このような解析を総称して回帰分析regression analysisと呼ぶ。また、しばしば回帰することをモデルに当てはめるfittingと表現することがある。特に今回解説するような、平均的な連続的変化を直線で表す解析を線形回帰linear regressionという。説明変数が連続変数になると、いままでとどのように変わるのだろうか。イメージ図を下に示した。

図にしてみるとなんとなく検定・分散分析と回帰分析は類似していることがわかるだろうか。検定や分散分析は計算上、説明変数に0/1を代入してその平均値の変化を検出するのに対して、回帰分析は0と1の間や0以下や1以上にわたって（正しくは、説明変数の範囲に限られているが）の平均値の変化を検出するイメージだ。言い方を変えれば、線形回帰は標本を直線y = βx + αに当てはめたときのαとβを推定する試みである。これは今まで紹介してきた線形モデルの記法にならって、

y ~ βx + α

と表記できることがわかるだろう。また、二元配置分散分析でも紹介したように、この時の帰無仮説は「α = 0」「β = 0」となる。「α = 0」とは、x = 0を代入した時、すなわち切片は0であるという仮説である。そして、「β = 0」とは連続変数xの変化に対して、yはほとんど変化しない、ほぼ水平な直線であることを示している。では、具体的にαとβはどのように計算し、帰無仮説を検証するのだろうか。

係数の推定：最小二乗法

　標本の平均値が、標本の代表的な値として認識されるのは、平均値が標本の「中心」あるいは、まさにすべてのデータから平均的に近い値であるからであろう。このように、線形回帰でもデータから平均的に近い直線を求める必要がある。そもそも、直線において「データから平均的に近い」とはどのようなものだろうか。データと直線の距離を測る方法はいくつか存在するが、線形回帰ではこの距離の指標として、予測値と実測値の差 = 残差の二乗和R（あるいは平方和）を用いる。例えば、x1に対応するy1、x2に対応するy2……xnに対応するyn、とサンプルサイズnの標本があったとする。yはy = βx + αの直線から予測されると考えられるのでy1の予測値^y1は^y1 = β x1 + αである。以下同様に、予測値を求める。そして予測値と実測値の差の二乗和Rを計算する。残差二乗和Rはまだ未知数のαとβを含んでいる。したがって、αとβをいろいろ変化させて = Rをαとβの関数と考えてRを最小にするαとβを探せばよいことになる（下図）。

この方法を最小二乗法 (ordinary) least squares methodと呼ぶ。まさに、残差二乗和を最小にする方法だからである。残差二乗和（残差平方和）という言葉は分散分析の項でも出てきたと思うが、平均からのずれの二乗和という、同じ考えに基づいていることがわかるだろう。また、線形回帰において残差は平均 = 0の正規分布に従い、かつその分散は一定であることを仮定している。ここで、Rは上図のようにαとβを変数として持つ平面関数となっている。関数の最小値を求める方法として微分法がある。元の関数に対し、それを微分した関数を導関数と呼ぶが、導関数 = 0となるときを確認すれば、最小値が求まる可能性があるのだった。しかし、2変数関数の場合、どうすればよいのだろうか。詳細は省くが、ほとんど普通の微分と同じ手続きで最小値を求めることができる。今回はαとβの2変数あるが、Rをαだけを変数としてみて微分（βは定数と考える）、Rをβだけを変数としてみて微分（αは定数と考える）、この2つの導関数を求めて、導関数 = 0となるαとβを求めればよい（この方法を偏微分法と呼ぶ）。以下に計算を示す。下図において、導関数 = 0となるときの^αと^βが最終的に求めたい値で、変数の状態のαとβと区別するために表記している。このように最小二乗法で求めた^αと^βを最小二乗推定量least squares estimator (OLS estimator)と呼ぶ。あるいは、確率分布を特徴づける変数や今後登場する他の方法で推定された推定量をまとめてパラメータparameterと呼ぶこともある。

上記のように、一見複雑な式を解かなければならないようだが、最終的には比較的シンプルな形で表記でき、平均をEと表現すれば、^α = E(y) - ^β×E(x)、^β = cov(x, y)/var(x)となる。 これらの値が、いわば分散分析などで求めた平均値の差にあたる要素である。

　また、^αの結果をもとの回帰式に入れて整理すると、

　　y - E(y) = ^β(x - E(x))

となり、回帰直線は必ず、点（E(x), E(y)）を通ることがわかる。

　共分散covarianceは本項で始めて登場した統計量だ。名前通り、共分散も分散の仲間である。標本分散は(1/n)Σ(xi - E(x))^2と計算する。その一方で、標本共分散は、2変数に関する分散ともいえる値であり、(1/n)Σ(xi - E(x))(yi - E(y))と計算する。分散がデータのばらつきの指標にであったのに対し、共分散は2変数の関係の強さ（相関）を示す値だ。関係の強さとは具体的に何かというと、xが増加した時にyも増加するか？減少するか？あるいは変化がないのか？を表す。実際、xが増加した時にyが増加すれば共分散は正に大きな値、減少すれば負に大きな値、変化がなければ0付近の値を示す。詳細は相関と因果の項で紹介しよう。

線形回帰の検定統計量

　上記で求めた係数が0から有意に異なっているかを明らかにするためには、これまでと同様に検定統計量が必要である。そのために、係数の分散を計算し、平方根（標準偏差）を求める。そして、回帰係数の推定量/(係数の標準偏差)がt分布に従うことを利用して検定を行う。ここでなんとt検定が登場した。確かに、t検定の検定統計量は(平均値の差)/(標本の標準偏差)であり、類似した検定統計量を計算していることがわかるだろう。係数の分散の導出は、本稿最後の追記することにして、ここでは避ける。残差が従う正規分布の分散をσ^2、xの残差平方和をSxxとして、最小二乗推定量^αおよび^βの分散は下記のように表現できる。

σは不偏推定量の項で話したように、我々が直接知ることはできない未知の値である。そこでσは不偏推定量^σで代用するしかない。母分散の不偏推定量が一ひねりあったように、線形回帰の残差が従う正規分布の母分散の推定も一ひねりいる。残差平方和をサンプルサイズで割るのではなく、サンプルサイズ - 回帰係数の推定量の数 = 自由度で割る。今回の回帰係数の推定量の数はαとβの二つなのでn - 2で割るのだ。不偏分散ではn - 1で割ったじゃないか、と思われるだろうが、脱線が進むので本稿の最後に追記する。とりあえず、以上の要素の計算ができれば、線形回帰の検定が可能となる。

逆問題とメカニズムの想起

　では、次のような設定で解析に使うデータの生成を行ってみよう。

------------------------------------------------------

library(plotn)

a <- 10 #切片

b <- 0 #傾き

s <- 5 #残差の母標準偏差

n <- 30 #サンプルサイズ

x <- runif(n, 0, 30)

y <- b * x + rnorm(n = n, mean = a, sd = s) #bx + aに従うデータの生成、残差は正規分布

#y ~ b * x + a + rnorm(n = n, mean = 0, sd = s)としてもよい。

fit1 <- lm(y ~ x) #線形回帰

fit2 <- summary(fit1) #fittingの要約

plotn(y ~ x) #図1の描画

ae <- coef(fit2)[1,1] #あてはめた結果得られた係数a

be <- coef(fit2)[2,1] #あてはめた結果得られた係数b

overdraw(abline(a = ae, b = be))

------------------------------------------------------

切片a = 10、傾きb = 0と指定した。すなわち、切片だけ有意に0から異なり、傾きはないことが期待される。実際、図1のようにデータのばらつきに対して、直線を当てはめると切片が10付近で、ほとんど傾きがないこ直線が推定されていることがわかる。

　さて、ここでデータの生成過程に注目すると、yにはb * x + rnorm(n = n, mean = a, sd = s)として生成したデータを格納している。このデータ生成過程は、まさに線形回帰が想定している仮定となっている。もちろん、これは線形回帰をするための仮定に合わせたデータ生成を私が選択したからだ。しかし、実際の解析ではデータの生成過程は未知である。実際の解析では、図1のデータ（xと対になるy）だけを確認し、その背後に存在するデータの生成過程（メカニズム）を想像することになる。xがyを引き起こす原因であるならば、xを入力input、yを出力outputとするような関数y = f(x)を想像するようなものである。このように入出力が与えられ、その関係を満たすような関係を解く解析を逆問題inverse problemと呼ぶ（一方、入力と関数が与えられ出力を得る問題は順問題direct probrem）。逆問題の解答はただ一つになるわけでなく、むしろたくさんの解答がありうる。図1で言えば、xの小さい順から点同士を線で結んだって、一つの答えである（図2）。

------------------------------------------------------

plotn(y ~ x)#図2の描画

d <- data.frame(x = x, y = y)

d <- d[order(x),]

for(i in 1:(nrow(d)-1)){

  overdraw(segments(d[i,"x"], d[i,"y"], d[i+1,"x"], d[i+1,"y"]))#xの小さい順に線で結ぶ

}

x_new <- runif(5, 0, 30)#同じ条件で5つ新たなデータを生成

y_new <- b * x_new + rnorm(n = 5, mean = a, sd = s)#同じ条件で5つ新たなデータを生成

overdraw(points(x_new, y_new, col = "red", pch = 16))#新たなデータの描画

------------------------------------------------------

しかし、それが納得できる答えかといえば、ほとんどの人はそう思わないだろう。なぜなら、もし新たなデータ点が得られた時、その点はどの線上にも乗らない可能性が高く、解答として不適切になることが容易に想像できてしまう（図2、この考え方はモデル選択にも通じる）。さらには、新たなデータが追加されるたびに答えを変える羽目になる。したがって、もっとデータの生成過程に寄り添った思考が必要になる。例えば、データ点はある直線的関係にあり、そこに誤差が加わることで観測データになっている（先行研究からそう考えるに足る根拠がある）、と考えることができれば、回帰分析のような手段を検討するだろう。ほかにも、複数の要因が観測値に影響を与えている可能性があるなら、分散分析や重回帰も候補に挙がる。直線関係でないなら、非線形回帰や一般化線形モデルの利用を考えるだろう。このように、得られた入出力からその背後にある現象を読み解こうとする逆問題的思考において統計モデリングは重要な役割を果たすのである。

線形回帰の検出力

　前置きが長くなったが、本格的にシミュレートしてみよう。次のように切片≠ 0、傾き = 0のデータを10000回生成して、線形回帰を行う。

------------------------------------------------------

a <- 10

b <- 0

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- b * x + rnorm(n = n, mean = a, sd = s)

  fit1 <- lm(y ~ x) #線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #切片aのp値

  p2 <- c(p2, coef(fit2)[2,4]) #傾きbのp値

}

plotn(y ~ x) #図3の描画

ae <- coef(fit2)[1,1]

be <- coef(fit2)[2,1]

overdraw(abline(a = ae, b = be))

histn(p1, xlab = "P value", ylab = "頻度") #図4の描画

histn(p2, xlab = "P value", ylab = "頻度") #図5の描画

sum(p1 < 0.05)

## [1] 9979

sum(p2 < 0.05)

## [1] 504

------------------------------------------------------

上記のように、切片の検出力は99%以上、傾きの危険率は5%程度となっている。

　次は、傾きだけが0ではない場合を考えよう。

------------------------------------------------------

a <- 0

b <- 1

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- b * x + rnorm(n = n, mean = a, sd = s)

  fit1 <- lm(y ~ x) #線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #切片aのp値

  p2 <- c(p2, coef(fit2)[2,4]) #傾きbのp値

}

plotn(y ~ x) #図6の描画

ae <- coef(fit2)[1,1]

be <- coef(fit2)[2,1]

overdraw(abline(a = ae, b = be))

histn(p1, xlab = "P value", ylab = "頻度") #図7の描画

histn(p2, xlab = "P value", ylab = "頻度") #図8の描画

sum(p1 < 0.05)

## [1] 488

sum(p2 < 0.05)

## [1] 10000

------------------------------------------------------

上記のように、切片の危険率は5%程度、傾きの検出力は100%程度となっている。

　次は、切片と傾きがともに0ではない場合を考えよう。

------------------------------------------------------

a <- 10

b <- 1

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- b * x + rnorm(n = n, mean = a, sd = s)

  fit1 <- lm(y ~ x) #線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #切片aのp値

  p2 <- c(p2, coef(fit2)[2,4]) #傾きbのp値

}

plotn(y ~ x) #図9の描画

ae <- coef(fit2)[1,1]

be <- coef(fit2)[2,1]

overdraw(abline(a = ae, b = be))

histn(p1, xlab = "P value", ylab = "頻度") #図10の描画

histn(p2, xlab = "P value", ylab = "頻度") #図11の描画

sum(p1 < 0.05)

## [1] 9983

sum(p2 < 0.05)

## [1] 10000

------------------------------------------------------

上記の条件の下では、切片と傾きの検出力はともに100%程度となっている。

線形回帰の危険率

　では、切片と傾きが0の場合について、危険率を明らかにする。

------------------------------------------------------

a <- 0

b <- 0

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- b * x + rnorm(n = n, mean = a, sd = s)

  fit1 <- lm(y ~ x) #線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #切片aのp値

  p2 <- c(p2, coef(fit2)[2,4]) #傾きbのp値

}

plotn(y ~ x) #図12の描画

ae <- coef(fit2)[1,1]

be <- coef(fit2)[2,1]

overdraw(abline(a = ae, b = be))

histn(p1, xlab = "P value", ylab = "頻度") #図13の描画

histn(p2, xlab = "P value", ylab = "頻度") #図14の描画

sum(p1 < 0.05)

## [1] 493

sum(p2 < 0.05)

## [1] 489

------------------------------------------------------

以上のように、切片および傾きの危険率は5%程度となった。このように、線形回帰では切片と傾きを分離して解析が可能である。

Rで行う線形回帰

　では、次のようなデータに関して線形回帰を実行する。まずはデータを図示しよう。

------------------------------------------------------

x <- c(7, 11, 27, 15, 5, 18, 9, 24, 4, 12, 

       19, 6, 8, 8, 24, 22, 11, 28, 0, 1, 

       12, 17, 10, 12, 9, 6, 18, 6, 29, 6)

y <- c(-20, -20, -53, -28, -4, -32, -11, -49, 3, -24, 

       -31, -11, -20, -12, -45, -45, -18, -48, 6, -3, 

       -22, -36, -15, -22, -26, -21, -21, -17, -53, -11)

plotn(y ~ x) #図15の描画

------------------------------------------------------

データを見ると、ばらつきは全体的に均等であり、右下がりの直線に当てはめることができそうだ。では、線形回帰を行ってみよう。線形回帰を行うときは、lm関数を使う。linear modelから名前をとられている。このとき、実際のモデルのようにy ~ b * x + aとするのではなく、説明変数を入力するだけ（y ~ x）で切片も自動的に解析してくれる。

------------------------------------------------------

fit1 <- lm(y ~ x) #線形回帰

fit2 <- summary(fit1)

fit1

## 

## Call:

## lm(formula = y ~ x)

## 

## Coefficients:

## (Intercept)            x  

##      0.2675      -1.8673

fit2

## 

## Call:

## lm(formula = y ~ x)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -10.0640  -3.8539   0.0378   3.2214  12.3430 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)   0.2675     1.8326   0.146    0.885    

## x            -1.8673     0.1216 -15.354 3.64e-15 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 5.297 on 28 degrees of freedom

## Multiple R-squared:  0.8938, Adjusted R-squared:   0.89 

## F-statistic: 235.7 on 1 and 28 DF,  p-value: 3.64e-15

------------------------------------------------------

fit1とfit2には同様の情報が格納されており、fit2のほうが見やすい形となっている。一番気になる情報はやはり、Coefficient（係数）の項だろう。計算は後にして、内容を確認する。まず(Intercept)は切片を示し、その中のEstimate列にモデルy ~ βx + αにおけるαの推定値が入っている。すなわち、α = 0.2675と推定されている。この時のP値が最右列のPr(>|t|)であり、P = 0.885となっている。つまり、5%水準において切片は有意に0と異なる、とは言えないことになる。xが0近辺の時はyは0.2675付近に分布しているが、それが0から大きく離れているとは言えなさそう、ということだ。一方、x行のEstimate列には傾きβの推定値が入っており、β = -1.8673だ。この時のP値は5%以下であり、傾きは有意に0と異なっている、と考えられる。つまり、xとyの間にはxが1増えれば、yは1.87程度減る、負の関係がある、といってもよいだろうと考える。これらの係数を使って図15に回帰直線を書き込んだものが図16である。

------------------------------------------------------

plotn(y ~ x) #図16の描画

ae <- coef(fit2)[1,1] #fit2にcoef関数をかけることで係数を取り出せる。

be <- coef(fit2)[2,1]

overdraw(abline(a = ae, b = be)) #abline関数は直線を描く関数。

------------------------------------------------------

今回は練習のため、あらゆる場合において回帰直線を描いてきたが、通常、傾きが有意でない場合は図中に回帰直線を描かない。回帰線はデータの要約としてみるのに役立つ一方で、ないかもしれない関係を強調してしまうことがある。どんなに傾きが大きくてもデータのばらつきが大きいからこそ、帰無仮説を棄却できないわけで、そんな中で、あるかわからない関係性を回帰線で強調してしまうのは危険だ。論文化するときは注意したい。今回は、傾きが有意なので回帰直線を描いてもよいだろう。また、今回はスクリプトの簡便化のために実行しなかったが、回帰線は説明変数の範囲だけ描くべきだ。説明変数の範囲外はデータがないので、なんの保証もなされない。突然、傾向が変わることだってあり得る。こちらも存在しない関係性を強調しないためにも、説明変数の範囲外の描画は避けよう。

　さて、多くの場合はCoefficientを見るだろうが、この後からは各要素について計算をする。fit2の中に格納されている情報を上から計算していこう。Residualsは文字通り残差で、予測値と実測値の差=残差に関する最大、最小、四分位点の情報である。　

------------------------------------------------------

quantile(y - (ae + be*x))

#yは実測値、ae + be*xは予測値でその差（残差）をもとめて、

#quantile関数で最大最小四分位点を求める

##           0%          25%          50%          75%         100% 

## -10.06396738  -3.85390306   0.03778293   3.22140799  12.34303388

------------------------------------------------------

確かに一致した。残差の最大最小、四分位点は残差の分布を確認するために使う。仮定として残差は平均0の正規分布に従うのだから、中央値は0付近で、最大と最小の絶対値、第一四分位数と第三四分位数の絶対値はおおむね一致していることが望まれる。例えば、残差の分布は次のように確認をする。

------------------------------------------------------

es <- y - (ae + be*x)#残差を計算

histn(es, xlab = "Residuals")#図17の描画（ヒストグラム）

es <- es[order(es)]

se <- sum((y - (ae + be*x))^2)/(length(x) - 2)#残差の分散を計算

es_t <- seq(-16, 16, length = 200)

es_cum <- sapply(es_t, function(parameter, vector = es){sum(vector < parameter)})

#残差の累積分布を計算

plotn(pnorm(q = seq(-16,16, length = 200), sd = sqrt(se)) ~ 

      seq(-16,16, length = 200), type = "l",

      ylab = "CDF", xlab = "Residuals")#図18の描画（累積分布関数の比較）

overdraw(points(es_t, es_cum/length(es), type = "l", col = "red"))

------------------------------------------------------

今まで10000個データを生成して分布と見比べたときと比較し、残差はサンプルサイズ、つまり今回は30しかなく、バー（ビンbin）の広さで見た目の印象が変わるヒストグラムでは判断が難しい。そこで、累積分布関数cumulative distribution function (CDF)を使って比較をしてみよう。詳しくは本項の最後に追記するが、図18の黒線が残差が従うと予測される正規分布で、赤線が残差の実際の累積分布である。黒線と赤線が大きくかけ離れていなければ、残差は正規分布に従っていると判断してもよい。今回は極端な逸脱が目立つ場所がないので、残差の正規性はあると判断してもよいだろう。

　次にCoefficientの項だ。これは、上記で議論してきたように各係数の推定値であるEstimate項は次のように求める。

------------------------------------------------------

mean(y) - mean(x) * cov(x,y)/var(x)#切片の推定値

## [1] 0.267468

cov(x,y)/var(x)#傾きの推定値

## [1] -1.86725

------------------------------------------------------

次に各係数の分散（標準偏差）であるStd. Error項を計算しよう。

------------------------------------------------------

se <- sum((y - (ae + be*x))^2)/(length(x) - 2) #残差二乗和を自由度で割る、これが残差の不偏分散

sb <- sqrt(se/(var(x) * (length(x)-1))) #b standard error

sa <- sqrt(

  se * (1/length(x) + 

          (mean(x)^2)/(var(x) * (length(x)-1)

                       )

        )

  ) #a standard error

sa

## [1] 1.832615

sb

## [1] 0.121617

------------------------------------------------------

t値であるt value項は係数の推定値を係数の標準偏差で割るのだった。

------------------------------------------------------

ae/sa #a t-value

## [1] 0.1459489

be/sb #b t-value

## [1] -15.35352

------------------------------------------------------

t値が求まればP値が計算できる。Pr(>|t|)項は次のように計算する。

------------------------------------------------------

pt(q = ae/sa, df = length(x)-2, lower.tail = sign(ae/sa) < 0) * 2

#自由度は残差の自由度30 - 2とする。

#t値の正負によってP値の計算が変わるのでlower.tail = sign(ae/sa) < 0という引数を入れている。

## [1] 0.8850074

pt(q = be/sb, df = length(x)-2, lower.tail = sign(be/sb) < 0) * 2

## [1] 3.640119e-15

------------------------------------------------------

確かにすべて一致している。

　そして、最後の段落だが、Residual standard errorは残差の標準偏差なので、残差の不偏分散を格納したseの平方根をとって、

------------------------------------------------------

sqrt(se) #Residual standard error

## [1] 5.296698

------------------------------------------------------

Multiple R-squaredは今回、初めて登場した統計量で、R二乗値、寄与率、決定係数などと呼ばれ、0～1の値をとる。この値は、データ全体のばらつきのうちモデルによってあらわされるばらつきの割合を示している。つまり、この値が高ければモデルが良く当てはまっている、と考えていいだろう。Multiple R-squared = 1のとき、残差が0であることを意味し、すべてのデータが推定された直線上にのっていることを意味する。この値は分散分析における、要因の偏差平方和/全体の偏差平方和 = 要因の偏差平方和/(要因の偏差平方和 + 残差の偏差平方和)と同じである。次のように計算する。

------------------------------------------------------

var(ae + be*x)/var(y) #Multiple R-squared

## [1] 0.8938311

st <- sum((y - mean(y))^2) #全体の偏差平方和

#要因の偏差平方和/全体の偏差平方和 = 要因の偏差平方和/(要因の偏差平方和 + 残差の偏差平方和)

# = 1 - 残差の偏差平方和/全体の偏差平方和と考えて、下記のように計算もできる。

#残差不偏分散seは自由度で割られているので、偏差平方和に直すには自由度を掛ける。

1 - (se* (length(x) - 2))/st #Multiple R-squared

## [1] 0.8938311

------------------------------------------------------

F valueは分散分析の検定統計量として登場したモデルと残差の不偏分散の比である。線形回帰ではt値によって傾きを検定したのだが、分散分析の考え方と同様に要因 = bxが説明する平均平方和（不偏分散）と残差の平均平方和の比率がF分布に従うことを利用して、モデルが実測データに影響を与えているかを検定できる。DFは自由度で、説明変数の自由度が2 - 1 = 1、残差の自由度が28であることを示す。あるいは全体の自由度 = サンプルサイズ - 1 = 29なので、説明変数の自由度 = 全体の自由度 - 残差の自由度 = 29 - 28と考えてもよい。同時にF検定のP値も計算する。

------------------------------------------------------

sr <- sum((ae + be*x - mean(ae + be*x))^2)/1

#モデルによる予測値の平均と予測値の偏差平方和を自由度で割る。

#今回は、2つの係数を推定しているので、2 - 1 = 1が説明変数の自由度

#st - se * (length(x) - 2)と計算してもよい

sr/se #F-value

## [1] 235.7307

pf(df1 = 1, df2 = length(x) - 2, q = sr/se, lower.tail = F)

#説明変数の自由度は1、残差の自由度は28

## [1] 3.640119e-15

------------------------------------------------------

このように計算できた。さらに注目するべきは、傾きのP値とF検定によるP値が一致している。たびたび登場していることからも明らかなように分散分析、t検定、回帰分析は表裏一体なのである。重回帰の場合は、複数の説明変数を組み込んだ総合的なF検定が行われるので、係数のP値とF検定のP値は一致しなくなる。

　以上のように線形回帰における様々な統計量を計算できた。線形回帰を行えば、連続変数xとyの関係を定量的に評価可能である。特にこのような連続的な変化をモデルに落とし込むことは予測という観点から重要となる。例えば、理想的には図15のように説明変数の全体わたって万遍なくデータが取れていることが望ましい。しかし、実際の野外データでは観測ができず、説明変数の連続性が途切れてしまうこともあるし、データはとったものの説明変数/被説明変数が何かの形で失われてしまい、欠損値として扱うこともある（図19）。

------------------------------------------------------

x <- c(7, 11, 27, 15, 5, 18, 9, 24, 4, 12, 

       19, 6, 8, 8, 24, 22, 11, 28, 0, 1, 

       12, 17, 10, 12, 9, 6, 18, 6, 29, 6)

y <- c(-20, -20, -53, -28, -4, -32, -11, -49, 3, -24, 

       -31, -11, -20, -12, -45, -45, -18, -48, 6, -3, 

       -22, -36, -15, -22, -26, -21, -21, -17, -53, -11)

d <- data.frame(x = x, y = y)

NAs <- d[d$x > 6 & d$x < 16,] #欠測にするデータを格納

d[d$x > 6 & d$x < 16, "y"] <- NA #データを欠測にする

plotn(y ~ x, data = d) #図19の描画

------------------------------------------------------

こんな時、不連続部分や欠損値に対してもモデルが妥当であれば予測を可能にしてくれるのが、統計の良いところなのだ。欠測値があるデータに関して、線形回帰を行う。

------------------------------------------------------

fit3 <- lm(y ~ x, data = d) #回帰分析

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -10.9545  -4.0227  -0.7727   4.0000  11.9545 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)   1.4091     2.4651   0.572    0.576    

## x            -1.9091     0.1414 -13.505 3.65e-10 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 5.86 on 16 degrees of freedom

##   ( 12 個の観測値が欠損のため削除されました )

## Multiple R-squared:  0.9194, Adjusted R-squared:  0.9143 

## F-statistic: 182.4 on 1 and 16 DF,  p-value: 3.647e-10

------------------------------------------------------

すると、求まった係数、特に傾きはほとんど同じである。もちろん、欠測がある場合、その欠測によってデータがゆがめられ、実際の関係性が崩れていることがあるので注意が必要だが（事実、元の解析とぴったり一致するわけではない）、今回はとりあえずうまく予測ができたとしよう。この分析結果があれば、欠測値がある領域にどんなデータがありえたかを予測ができる。例えば、次のように図示しよう。

------------------------------------------------------

plotn(y ~ x, data = d) #図20の描画

ae <- coef(fit4)[1,1]

be <- coef(fit4)[2,1]

overdraw(abline(a = ae, b = be)) #回帰直線

n_new <- 100000

x_new <- runif(n_new, 6, 16) #欠測がある領域、6 < x < 16のデータを10万回生成する

y_new <- be * x_new + rnorm(n = n_new, mean = ae, sd = fit4$sigma) 

#欠測がある領域、6 < x < 16のデータを10万回生成する

#fit4$sigmaとすれば、残差の標準偏差を得られる

overdraw(points(x_new, y_new, pch = 16, col = "#00000002"))#生成されたデータの描画

#半透明色を使うことでたくさん生成された部分は濃く、めったに生成されない部分は薄く表示される。

#濃い色の領域のデータが生成されやすいことを示す。

overdraw(points(NAs$x, NAs$y))#実際の欠測値

------------------------------------------------------

欠測部分の濃い色の部分がデータが生成されやすい = この範囲に現実データが収まっているだろうことを予測している。実際、欠測値を描画すると、濃い領域に多くが収まっていることがわかるだろう。このように、統計は差を検出あるいは重要な要因の抽出だけでなく、予測に使える。予測が可能になるということは、新たな知見が増えることでもある。自分のデータと向き合うときは、上記のような視点から眺めてみることで、さらなる洞察を得られるだろう。

追記1：推定量の分散

　推定量の分散は以下のような考え方で算出する。

追記2：自由度

　いままで自由度の詳細を述べずに来たが、ここでさっと紹介する。自由度degree of freedomとは、「互いに制約のある変数のうち自由な（独立な）値をとれる変数の個数」である。例えば、あるサンプルサイズnの標本Xを考える。この状態は一切の制約がないので、n個の変数は自由な値をとれ、自由度nとなる。ここで、標本Xの平均E(X)が定まっているとする。すると、変数間に制約が生まれる。平均E(X)を変えることなく、自由に選べる変数はn - 1個となる。n - 1個までは自由に選んでもよいが、n個目まで自由に選んでしまうとE(x)が変わってしまうだろう。このような時を自由度n - 1とする。

　一般に、モデルの自由度は推定が必要なパラメータの数(p) - 1、残差の自由度はサンプルサイズ(n) - pとされる。そして、全体の自由度は、モデルの自由度と残差の自由度を足してn - 1とされる。あるいは、全体の自由度は平均の関係式による制約によって、サンプルサイズに対して自由度が1小さい、と考えてもよい。例えば、1説明変数の線形回帰では、

y ~ bx + a

とすれば、aとbの2つが推定したいパラメータである。なのでこのモデルの自由度はp - 1 = 1である。残差の自由度はn - p = n - 2、全体の自由度はn - 1である。 自由度の定義から考えて、パラメータ間は1つの関係式による制約があり（^a = E(y) - ^b E(x)）、ゆえにモデルは1つのパラメータが定まればもう一方のパラメータが定まるので、自由度1と考えられる。残差は2つのパラメータの関係式があるため、サンプルサイズから2を引いた値になる。

　2水準（g1: A、B）および3水準（g2: X、Y、Z）の2説明変数とそれらの交互作用を考慮した分散分析の自由度は、モデル式に基づいて、

Xの場合：y ~ b0 + b1 × g1

Yの場合：y ~ b0 + b1 × g1 + b2_y + b3_y × g1

Zの場合：y ~ b0 + b1 × g1 + b2_z + b3_z × g1

となり、求めるべきパラメータの数は6個なので、モデルの自由度は5となる。さらに要因ごとに自由度を考えられて、2水準のg1は、g1の平均の関係式から一方の平均が求まれば、もう一方の平均も求まるので2 - 1 = 1。3水準のg2は、g2の平均の関係式から2水準の平均が求まれば、最後の1水準の平均も求まるので3 - 1 = 2。交互作用の自由度は、g1とg2の自由度を掛ければ求まり、1×2 = 2。これらのすべての要因の自由度を足し合わせれば5である。そして、各水準に含まれるサンプルサイズを足し合わせた総サンプルサイズをnとして、残差の自由度はn - 6となる。

　では、ここであえて、切片だけの線形モデルを考えてみよう。

y ~ a

このとき、求めるべきパラメータは1個なので、モデルの自由度は0である。残差の自由度はn - 1である。ではこの時の、推定値aはどんな値だろうか。定義にならって最小二乗法を試してみよう。

上記のように、^a = E(y)となった。すなわち、平均を求める作業とは別の視点からみれば、切片だけの線形モデルを仮定し、切片の最小二乗推定値を求める作業といえるのである。そしてこの時の残差平方和とは、定義からまさに偏差平方和のことである。この不偏分散とは偏差平方和をサンプルサイズ(n) - 1で除したものだったから、今回の残差平方和を残差の自由度で除したものと一致している。

　理論的には、標本の母分散の不偏推定量は、標本の偏差平方和を標本の自由度で除したものとなる。それゆえ、分散分析の各要因の平均平方和は偏差平方和を要因の自由度で除することで各要因の母分散を推定しているし、分散分析や線形回帰の残差平方和を残差の自由度で除することで残差の母分散を推定している。したがって、1標本における不偏分散はn - 1で除し、2パラメータを持つ線形回帰における残差平方和はn - 2で除するのである。

追記3：累積分布関数

　累積分布関数とは、確率変数Xの確率分布に対して、x以下となる確率の関数である（下図）。最小値は0、最大値は1となる。

過去のシミュレーションのように極めてサンプルサイズが大きいときは、ヒストグラムと確率分布の比較で十分、分布の類似性が確認できていた。しかし、サンプルサイズが小さいときはその比較は不適切になりうる。というのも、ヒストグラムはサンプルサイズが小さいときにはビンによって見た目が大きく変わるため、比較が難しくなる。そこでヒストグラムのような離散的変化の図ではなく、連続的な変化を示す図で確認するほうがより正確である。そこで用いられるのが累積分布関数である。確率密度関数が分布によって異なっているように、累積密度関数も分布によって特徴がある。よく知られた正規分布はもちろんのこと、自分の手元の標本の累積密度関数も簡単に描くことができる。例えば、標本が5 ～ 18の間に存在する連続値であるなら、5.01、5.02……、17.98、17.99という「刻み」を用意し、5.01以下のデータの個数、5.02以下のデータの個数……というように計算する。最後に各データの個数をサンプルサイズで割れば、累積分布関数になる。上記の残差の分布を見たときも、これにならって計算をした。

------------------------------------------------------

es <- y - (ae + be*x)#残差を計算

es <- es[order(es)]

se <- sum((y - (ae + be*x))^2)/(length(x) - 2)#残差の分散を計算

es_t <- seq(-16, 16, length = 200)#残差の最大最小をもとに「刻み」を用意

es_cum <- sapply(es_t, function(parameter, vector = es){sum(vector < parameter)})

#sum(vector < parameter)はparameter以下のvector内の要素の個数を数えている。

#parameter = es_t、vector = esなので、es_t[1]より小さなes内の要素の個数、

#es_t[2]よりも……となっている

plotn(pnorm(q = seq(-16,16, length = 200), sd = sqrt(se)) ~ 

      seq(-16,16, length = 200), type = "l",

      ylab = "CDF", xlab = "Residuals")#図18の描画（累積分布関数の比較）

overdraw(points(es_t, es_cum/length(es), type = "l", col = "red"))

------------------------------------------------------