社会科学のためのデータ分析入門
第4章 予測
高橋 慎
2021-11-03
4.2.4 平均への回帰
ゴルトンの平均への回帰
回帰の効果ともいう．
「両親の身長の平均値が全体平均よりも低い（高い）とき，その子供は両親よりも身長が高い（低い）傾向にある」
ゴルトン卿の計算例：
両親の身長の平均が70インチのとき子供の身長が70インチ以上である割合 \[\frac{(8+5+3)}{(1+2+3+5+8+9+9+8+5+3)} \ = \ 0.30\]
両親の身長の平均が68インチのとき子供の身長が68インチ以上である割合 \[\frac{(13+10+7+3+1)}{(3+7+11+13+14+13+10+7+3+1)} \ = \ 0.41\]
平均への回帰は，予測変数の値が分布の平均から大きく離れている観察において，結果変数の値が平均に近づく傾向にあるという経験的事象を表す．
この傾向は，偶然のみによって説明できる．
4.2.5 Rにおけるデータの結合
大統領選挙において平均への回帰現象が見られるかどうか検証する．
そのために，2008年の選挙におけるオバマの得票率を用いて，彼の2012年における得票率を予測する．
2012年の選挙結果が記録されたデータセット
pres12.csvを2008年の選挙データセットに結合（merge）して用いる．
pres08 <- read.csv("pres08.csv") # 2008年のデータを読み込む
pres12 <- read.csv("pres12.csv") # 2012年のデータを読み込む
## 2つのデータセットをざっと見る
head(pres08)## state.name state Obama McCain EV
## 1 Alabama AL 39 60 9
## 2 Alaska AK 38 59 3
## 3 Arizona AZ 45 54 10
## 4 Arkansas AR 39 59 6
## 5 California CA 61 37 55
## 6 Colorado CO 54 45 9## state Obama Romney EV
## 1 AL 38 61 9
## 2 AK 41 55 3
## 3 AZ 45 54 11
## 4 AR 37 61 6
## 5 CA 60 37 55
## 6 CO 51 46 92012年の選挙データセットに含まれる変数
state：州の略称Obama：オバマの得票率（パーセンテージ）Romney：ロムニーの得票率（パーセンテージ）EV：その州の選挙人票の数
Rでは，merge()関数を用いて2つのデータセットを結合できる．データを結合するために，2つのデータセットに含まれる
state変数を用いる．
## state state.name Obama.x McCain
## Length:51 Length:51 Min. :33.00 Min. : 7.00
## Class :character Class :character 1st Qu.:43.00 1st Qu.:40.00
## Mode :character Mode :character Median :51.00 Median :47.00
## Mean :51.37 Mean :47.06
## 3rd Qu.:57.50 3rd Qu.:56.00
## Max. :92.00 Max. :66.00
## EV.x Obama.y Romney EV.y
## Min. : 3.00 Min. :25.00 Min. : 7.00 Min. : 3.00
## 1st Qu.: 4.50 1st Qu.:40.50 1st Qu.:41.00 1st Qu.: 4.50
## Median : 8.00 Median :51.00 Median :48.00 Median : 8.00
## Mean :10.55 Mean :49.06 Mean :49.04 Mean :10.55
## 3rd Qu.:11.50 3rd Qu.:56.00 3rd Qu.:58.00 3rd Qu.:11.50
## Max. :55.00 Max. :91.00 Max. :73.00 Max. :55.00データフレーム同士が同じ名前の変数を含んでいれば，結合されたデータフレームの当該変数には末尾に
.xや.yが加えられ，それぞれの変数がどちらのデータフレームのものかが示される．変数名が異なる場合は
by.xとby.yという引数を用いて，それぞれのデータフレームで結合用の変数名を指定することもできる．デフォルトでは結合されたデータフレームには引数
by.xで指定されたデータフレームxの変数名が保存される．
## 例示のため変数名を変更
names(pres12)[1] <- "state.abb"
## 異なる名前の変数を用いてデータセットを結合
pres <- merge(pres08, pres12, by.x = "state", by.y = "state.abb")
summary(pres)## state state.name Obama.x McCain
## Length:51 Length:51 Min. :33.00 Min. : 7.00
## Class :character Class :character 1st Qu.:43.00 1st Qu.:40.00
## Mode :character Mode :character Median :51.00 Median :47.00
## Mean :51.37 Mean :47.06
## 3rd Qu.:57.50 3rd Qu.:56.00
## Max. :92.00 Max. :66.00
## EV.x Obama.y Romney EV.y
## Min. : 3.00 Min. :25.00 Min. : 7.00 Min. : 3.00
## 1st Qu.: 4.50 1st Qu.:40.50 1st Qu.:41.00 1st Qu.: 4.50
## Median : 8.00 Median :51.00 Median :48.00 Median : 8.00
## Mean :10.55 Mean :49.06 Mean :49.04 Mean :10.55
## 3rd Qu.:11.50 3rd Qu.:56.00 3rd Qu.:58.00 3rd Qu.:11.50
## Max. :55.00 Max. :91.00 Max. :73.00 Max. :55.00cbind()という複数のデータフレームの列を結合する関数を用いることもできるが，2点の問題がある．データフレーム間で対応する観察が同じ行にきちんと並んでいることを前提としている．
2つの列が同じ情報を含む同じ変数を表すものであっても，データフレームに含まれるすべての列を保存してしまう．
## state.name state Obama McCain
## Length:51 Length:51 Min. :33.00 Min. : 7.00
## Class :character Class :character 1st Qu.:43.00 1st Qu.:40.00
## Mode :character Mode :character Median :51.00 Median :47.00
## Mean :51.37 Mean :47.06
## 3rd Qu.:57.50 3rd Qu.:56.00
## Max. :92.00 Max. :66.00
## EV state.abb Obama Romney
## Min. : 3.00 Length:51 Min. :25.00 Min. : 7.00
## 1st Qu.: 4.50 Class :character 1st Qu.:40.50 1st Qu.:41.00
## Median : 8.00 Mode :character Median :51.00 Median :48.00
## Mean :10.55 Mean :49.06 Mean :49.04
## 3rd Qu.:11.50 3rd Qu.:56.00 3rd Qu.:58.00
## Max. :55.00 Max. :91.00 Max. :73.00
## EV
## Min. : 3.00
## 1st Qu.: 4.50
## Median : 8.00
## Mean :10.55
## 3rd Qu.:11.50
## Max. :55.00## state.name state Obama McCain EV state.abb Obama Romney EV
## 8 D.C. DC 92 7 3 DE 59 40 3
## 9 Delaware DE 62 37 3 DC 91 7 3## state state.name Obama.x McCain EV.x Obama.y Romney EV.y
## 8 DC D.C. 92 7 3 91 7 3
## 9 DE Delaware 62 37 3 59 40 3結合されたデータフレームを用いて，アメリカ大統領選挙データで平均への回帰が実際に見られるかどうかを調べる．
アメリカ政治の分極化の進展を考慮して，選挙年ごとに各州の得票率のz得点を計算し，標準化（平均を引いて標準偏差で割ること）を行う．
こうすることで，オバマの各州における得票結果を，彼の当該選挙年における平均的な結果と比較して測定することができる．
これは
scale()関数を使って簡単に行うことができる．
- オバマの2012年の標準化得票率を，2008年の標準化得票率で回帰する．
##
## Call:
## lm(formula = Obama2012.z ~ Obama2008.z, data = pres)
##
## Coefficients:
## (Intercept) Obama2008.z
## -3.521e-17 9.834e-01両者の間には強い正の線形関係がみられ，オバマは，2008年により多くの票を得た州から2012年でもより多くの票を得ている．
結果変数と予測変数を共に標準化すると切片は０になる点に注意したい．
- 標準化すると両変数の標本平均\(\bar{Y}\)と\(\bar{X}\)は共に０になるから，切片の予測値は\(\hat{\alpha} = \bar{Y} - \hat{\beta} \bar{X} = 0\)となる．
式に
-1を加えることでモデルを切片なしで推定することもできる．
##
## Call:
## lm(formula = Obama2012.z ~ -1 + Obama2008.z, data = pres)
##
## Coefficients:
## Obama2008.z
## 0.9834- 当てはめられた回帰直線とデータの座標の両方をプロットすると，強い線形関係を観察することができる．
par(family = "HiraginoSans-W3") # 日本語が文字化けしないようにおまじない
plot(pres$Obama2008.z, pres$Obama2012.z, xlim = c(-4, 4), ylim = c(-4, 4),
xlab = "オバマの標準化得票率（2008年）",
ylab = "オバマの標準化得票率（2012年）")
abline(fit1) # 回帰直線を描く
オバマが2012年に2008年のときよりも高い標準化得票率を記録した州の割合を計算する．
平均への回帰現象が存在するのであれば，こうした州の割合は下位25パーセンタイルのほうが上位25パーセンタイルよりも多いはずである．
下位25パーセンタイルと上位25パーセンタイルの計算には
quantile()という関数を用いる．
## 下位25パーセンタイル
mean((pres$Obama2012.z >
pres$Obama2008.z)[pres$Obama2008.z
<= quantile(pres$Obama2008.z, 0.25)])## [1] 0.5714286## 上位25パーセンタイル
mean((pres$Obama2012.z >
pres$Obama2008.z)[pres$Obama2008.z
>= quantile(pres$Obama2008.z, 0.75)])## [1] 0.4615385- 結果は，平均への回帰現象を裏づけるものである．
4.2.6 モデルの当てはまり
モデルの当てはまりは，どの程度モデルがデータに当てはまっているか，つまり，どれくらい正確にモデルが観察を予測しているかの尺度である．
モデルの当てはまりは，決定係数（coefficient of determination），すなわち\(R^2\)を見ることで評価できる．
\[ R^2 \ = \ 1 - \frac{\textsf{SSR}}{\textsf{総平方和（total sum of squares：TSS）}} \ = \ 1 - \frac{\sum_{i=1}^n \hat\epsilon_i^2}{\sum_{i=1}^n (Y_i - \overline{Y})^2} \]
\(R^2\)の値は０（結果変数と予測変数の相関が０）から１（相関が１）の間をとり，線形モデルがどの程度うまくデータに当てはまっているかを示す．
決定係数はモデルの当てはまりの尺度で，予測変数によって説明された結果変数の変動の割合を表し，１から残差平方和（SSR）と総平方和（TSS）との比を引いたものと定義される．
わかりやすい例として，2000年のアメリカ大統領選挙でのフロリダ州における群ごとの開票結果を，1996年の結果から予測する．
フロリダには67の郡があり，
florida.csvにはこれら2つの選挙における候補者別得票数が含まれている．
## county Clinton96 Dole96 Perot96 Bush00 Gore00 Buchanan00
## 1 Alachua 40144 25303 8072 34124 47365 263
## 2 Baker 2273 3684 667 5610 2392 73
## 3 Bay 17020 28290 5922 38637 18850 248
## 4 Bradford 3356 4038 819 5414 3075 65
## 5 Brevard 80416 87980 25249 115185 97318 570
## 6 Broward 320736 142834 38964 177323 386561 788データファイルに含まれる変数
county：群の名前Clinton96：1996年のクリントンの得票数Dole96：1996年のドールの得票数Perot96：1996年のペローの得票数Bush00：2000年のブッシュの得票数Gore00：2000年のゴアの得票数Buchanan00：2000年のブキャナンの得票数
ここでは，共にリバタリアン（自由至上主義）で1996年に立候補したロス・ペロー（Ross Perot）と，彼と同じ政党から2000年に立候補したパット・ブキャナンに注目し，前者への投票から後者への投票を予測する．
##
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = florida)
##
## Coefficients:
## (Intercept) Perot96
## 1.34575 0.03592この回帰モデルからTSSとSSRを計算し，\(R^2\)を算出する．
resid()関数を用いて回帰分析の出力から残差のベクトルを抽出する．
## TSS（総平方和）とSSR（残差平方和）を計算する
TSS2 <- sum((florida$Buchanan00 - mean(florida$Buchanan00))^2)
SSR2 <- sum(resid(fit2)^2)
## 決定係数
(TSS2 - SSR2) / TSS2## [1] 0.5130333結果から，ブキャナンの2000年の得票の変動の51%がペローの1996年の票によって説明できることがわかる．
行った計算を関数にして，異なる回帰モデルの決定係数を簡単にできるようにする．
この関数は，様々な要素を含むリスト形式のオブジェクトである
lm()関数の出力を入力する．結果変数の値は，回帰分析の出力オブジェクトを用いて
fitted()関数から得られる当てはめ値と，各観察の残差を足し合わせることで再計算できる．
R2 <- function(fit) {
resid <- resid(fit) # 残差
y <- fitted(fit) + resid # 結果変数
TSS <- sum((y - mean(y))^2)
SSR <- sum(resid^2)
R2 <- (TSS - SSR) / TSS
return(R2)
}
R2(fit2)## [1] 0.5130333- この関数を用いる代わりに，
summary()関数の入力としてlm()関数の出力を用いることで\(R^2\)を得ることもできる．
## [1] 0.5130333- 比較のために，先に見たオバマの州レベルの回帰分析の\(R^2\)を，回帰分析結果を表す出力オブジェクトである
fit1を用いて計算する．
## [1] 0.9671579フロリダ州の（ペローとブキャナンの）回帰分析の決定係数は，州レベルの（オバマの）回帰分析のものよりもかなり低いことがわかる．
当てはまりが悪いことから，フロリダ州の回帰分析の残差をより詳細に調査する．
そのためには，残差を当てはめ値に対してプロットした残差プロット（residual plot）を作成する．
par(family = "HiraginoSans-W3") # 日本語が文字化けしないようにおまじない
plot(fitted(fit2), resid(fit2), xlim = c(0, 1500), ylim = c(-750, 2500),
xlab = "当てはめ値", ylab = "残差")
abline(h = 0)
図からは，極めて大きな残差すなわち外れ値があり，2000年の選挙でブキャナンが予想されたよりも2000票も多く得票していることがわかる．
残差の値が最大値をとる群の名前を抽出する．
## [1] "PalmBeach"- 実は，外れ値を示したパームビーチ郡では，いわゆるチョウ型投票用紙（butterfly ballot）がこの選挙で使われていた．
パームビーチ群のチョウ型投票用紙
有権者は，投票したい候補者に対応する穴を開けることになっているが，この投票用紙はとても紛らわしい．
この郡の多くのアル・ゴア支持者は，上から３番目の穴ではなく，２番目のブキャナン候補の穴を間違って開けたようである．
2000年の選挙では，フロリダ州でブッシュが僅差で勝利したために大統領に当選したが，このパームビーチ郡での疑問の残る開票結果によってゴアは大統領のイスを逃したと広く信じられている．
同じモデルをパームビーチ群を除いて推定する．
## パームビーチ群を除くデータ
florida.pb <- subset(florida, subset = (county != "PalmBeach"))
fit3 <- lm(Buchanan00 ~ Perot96, data = florida.pb)
fit3##
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = florida.pb)
##
## Coefficients:
## (Intercept) Perot96
## 45.84193 0.02435## [1] 0.8511675- パームビーチ群を除くと，決定係数は0.51から0.85へと大幅に向上した．
par(family = "HiraginoSans-W3") # 日本語が文字化けしないようにおまじない
## 残差プロット
plot(fitted(fit3), resid(fit3), xlim = c(0, 1500), ylim = c(-750, 2500),
xlab = "当てはめ値", ylab = "残差",
main = "パームビーチ群を除いた残差プロット")
abline(h = 0) # ０における水平線
モデルの当てはまりの改善は，残差プロットからも見てとれる．
パームビーチ群を含めた場合と含めなかった場合での回帰直線を比較する．
par(family = "HiraginoSans-W3") # 日本語が文字化けしないようにおまじない
## 散布図と回帰直線
plot(florida$Perot96, florida$Buchanan00, xlab = "1996年のペローの得票数",
ylab = "2000年のブキャナンの得票数") # 散布図
abline(fit2, lty = "dashed") # パームビーチ群を含んだ回帰直線
abline(fit3) # パームビーチ群を除いた回帰直線
text(30000, 3250, "パームビーチ群")
text(30000, 1500, "パームビーチ群を含んだ回帰直線")
text(30000, 400, "パームビーチ群を除いた回帰直線")
パームビーチ群を除くことで回帰直線が大きく変化していることから，回帰直線がパームビーチ群の影響を非常に強く受けていたことがわかる．
新しい回帰直線は，パームビーチ群を除いた残りの観察値によく当てはまっている．
ここで扱ったモデルの当てはまりは，サンプル内予測（in-sample prediction）に基づくものであり，サンプル外予測（out-of-sample prediction）ではない点に留意する必要がある．
つまり，決定係数などのモデルの当てはまりを示す統計量は，あるモデルがどの程度そこにある標本に当てはまるかを示すものである．
特定の標本に適合するよう過度に調整されたオーバーフィッティング（overfitting）のモデルは，別の標本に対してはそれほど正確な予測はできないかもしれない．
他のデータにも適用可能な一般的なモデルを求めるのであれば，モデルを特定の標本へオーバーフィッティングさせないよう注意を払う必要がある．