まず、ビールの出荷量に関する分析を実行する(説明は省略する)。
#dir.create("data") # dataディレクトリがないなら作る
## データを持っていないならダウンロード
#download.file(url = "https://git.io/fA6Zk",
# destfile = "data/beer2010.csv")
Beer <- read_csv("data/beer2010.csv") # データを読み込む
#glimpse(Beer) # データの中身を確認
## ビールの売り上げを気温に回帰する
fit_beer <- lm(beer ~ temp, data = Beer)
summary(fit_beer) # 結果を表示
##
## Call:
## lm(formula = beer ~ temp, data = Beer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.860 -20.242 0.882 11.037 76.939
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.546 21.735 4.028 0.00241
## temp 2.104 1.165 1.805 0.10117
##
## Residual standard error: 31.98 on 10 degrees of freedom
## Multiple R-squared: 0.2458, Adjusted R-squared: 0.1704
## F-statistic: 3.259 on 1 and 10 DF, p-value: 0.1012
この分析の残差プロットを描く。
res_beer <- tibble(res = fit_beer$residuals,
fitted = fit_beer$fitted.values) %>%
ggplot(aes(x = fitted, y = res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x = "予測値", y = "残差")
print(res_beer)
はっきりとしたパタンは見られないものの、隣り合う予測値同士の残差が近いところが多いように見える。また、分散は均一ではないようにも見えるが、一つだけある外れ値を除けば、分散が均一だと言えないこともなさそうだ。標本サイズ(観測数)が12しかないので、はっきりしたことはわからない。
次に、正規QQプロットを描いてみよう。
## 残差を標準化する
z_res <- with(fit_beer, (residuals - mean(residuals)) / sd(residuals))
## 図を描く
qq_beer <- tibble(z_res) %>%
ggplot(aes(sample = z_res)) +
geom_abline(intercept = 0, slope = 1, # 45度線
linetype = "dashed", color = "gray") +
geom_qq() +
labs(x = "標準正規分布", y = "標準化した残差の分布") +
xlim(-4, 4) + ylim(-4, 4) +
coord_fixed() # 横軸と縦軸を1:1にする
print(qq_beer)
点が45度線からずれており、誤差の分布が正規分布に従っていないように見える。標本サイズが小さいので、やはりはっきりしたことはわからないが、このモデルの正しさを積極的に支持することは難しそうである。
Rds形式の衆院選データ (hr-data.Rds) を読み込む。 手元にない場合はまずダウンロードする。
# dir.create("data") # dataディレクトリがない場合は作る
#download.file(url = "https://git.io/fp00p",
# destfile = "data/hr-data.Rds")
HR <- read_rds("data/hr-data.Rds")
## Rdsファイルの読み込みがうまくいかない場合は以下を実行してCSVファイルを使う
#download.file(url = "https://git.io/fxhQU",
# destfile = "data/hr-data.csv")
#HR <- read_csv("data/hr-data.csv")
正しく読み込めたかどうか確認する。
## Rows: 8,803
## Columns: 22
## $ year <dbl> 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996,…
## $ ku <chr> "aichi", "aichi", "aichi", "aichi", "aichi", "aichi", "aich…
## $ kun <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,…
## $ status <fct> 現職, 元職, 現職, 新人, 新人, 新人, 新人, 現職, 元職, 新人,…
## $ name <chr> "KAWAMURA, TAKASHI", "IMAEDA, NORIO", "SATO, TAISUKE", "IWA…
## $ party <chr> "NFP", "LDP", "DPJ", "JCP", "others", "kokuminto", "indepen…
## $ party_code <dbl> 8, 1, 3, 2, 100, 22, 99, 8, 1, 3, 2, 10, 100, 99, 22, 8, 1,…
## $ previous <dbl> 2, 3, 2, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0,…
## $ wl <fct> 当選, 落選, 落選, 落選, 落選, 落選, 落選, 当選, 落選, 復活…
## $ voteshare <dbl> 40.0, 25.7, 20.1, 13.3, 0.4, 0.3, 0.2, 32.9, 26.4, 25.7, 12…
## $ age <dbl> 47, 72, 53, 43, 51, 51, 45, 51, 71, 30, 31, 44, 61, 47, 43,…
## $ nocand <dbl> 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7,…
## $ rank <dbl> 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5,…
## $ vote <dbl> 66876, 42969, 33503, 22209, 616, 566, 312, 56101, 44938, 43…
## $ eligible <dbl> 346774, 346774, 346774, 346774, 346774, 346774, 346774, 338…
## $ turnout <dbl> 49.2, 49.2, 49.2, 49.2, 49.2, 49.2, 49.2, 51.8, 51.8, 51.8,…
## $ exp <dbl> 9828097, 9311555, 9231284, 2177203, NA, NA, NA, 12940178, 1…
## $ expm <dbl> 9.828097, 9.311555, 9.231284, 2.177203, NA, NA, NA, 12.9401…
## $ vs <dbl> 0.400, 0.257, 0.201, 0.133, 0.004, 0.003, 0.002, 0.329, 0.2…
## $ exppv <dbl> 28.341505, 26.851941, 26.620462, 6.278449, NA, NA, NA, 38.2…
## $ smd <fct> 当選, 落選, 落選, 落選, 落選, 落選, 落選, 当選, 落選, 落選,…
## $ party_jpn <chr> "新進党", "自民党", "民主党", "共産党", "その他", "国民党",…
2012年衆院選の自民党候補だけを抜き出してデータフレームを作る。
得票率 voteshare、選挙費用 expm 、年齢 age の関係を、次の線形モデルで表す。 \[得票率_i = \beta_0 + \beta_1 選挙費用_i + \beta_2 年齢_i + \epsilon_i.\]
母数について、以下の2組の仮説を立てる。
有意水準を5% (0.05) に設定し、回帰分析でこのモデルの母数である \(\beta_0\)、\(\beta_1\)、\(\beta_2\) の推定と検定を行う。
##
## Call:
## lm(formula = voteshare ~ expm + age, data = LDP2012)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.392 -8.661 -1.616 5.744 39.938
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.95257 3.60177 7.206 5.30e-12
## expm 0.63259 0.18087 3.498 0.000545
## age 0.26858 0.06361 4.223 3.26e-05
##
## Residual standard error: 11.94 on 282 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.1071, Adjusted R-squared: 0.1007
## F-statistic: 16.91 on 2 and 282 DF, p-value: 1.164e-07
包括的な仮説から検定しよう。上で示された結果のうち、最後の行にある\(F\)統計量の\(p値\) を見ると、\(1.16 \times 10^{-7} < 0.05\) だから、有意水準5%で帰無仮説を棄却する。すなわち、\(\beta_1\)と\(\beta_2\)のうち、少なくとも一方は0ではない。言い換えると、(上で想定したモデルが正しいなら)選挙費用と年齢のうち少なくとも一方は得票率に影響を与える。
次に、係数を個別に検定しよう。まず、過去の選挙費用 (expm) の係数の\(p\)値は、\(0.000545 < 0.05\) だから、有意水準5%で帰無仮説を棄却する。すなわち、選挙費用は得票率に影響を与える。\(\beta_1\)の推定値は約\(0.63\)だから、年齢を一定に保つと、選挙費用が100万円増えるごとに、得票率は平均すると\(0.63\)パーセントポイントずつ上昇すると考えられる。この選挙における自民党候補の選挙費用の標準偏差は、3.94百万円である。選挙費用が400万円(おおよそ1標準偏差分)異なると、\(2.53\)パーセントポイント得票率に差が出るということなので、この効果は実質的にはそれほど重要ではなさそうだ。
年齢 (age) の係数の\(p\)値は、\(3.26 \times 10^{-5} < 0.05\)だから、有意水準5%で帰無仮説を棄却する。すなわち、年齢は得票率に影響を与える。\(\beta_2\)の推定値は約\(0.27\)だから、過去の当選回数を一定に保つと、年齢が1歳上がるごとに得票率が\(0.27\)パーセントポイントずつ上がることが予測される。 一回り年齢が違うと、\(3.22\)パーセントポイントだけ得票率に差が出るということなので、無視できるほど小さい差ではないが、実質的に大きな効果があると言い切ることは難しい。30歳の候補と60歳の候補を比べると、\(8.06\)パーセントポイントの差になるので、年齢差が大きければ実質的に意味がある違いが出ることがわかる(2012年の自民候補のうち、最年少候補は28歳、最年長候補はは78歳である)。
この分析結果を図にまとめる。
coefs_q12 <- coefplot(
model = fit_q12, # 図示するモデルを指定
intercept = FALSE, # 切片を表示しない
pointSize = 4, # 点の大きさを4に
## 内側の線を50%信頼区間にする。既定値は1(「点推定値±1標準誤差」)
innerCI = qt(df = summary(fit_q12)$df[2],
p = 0.25, lower.tail = FALSE),
lwdInner = 2.5, # 内側の線の太さを2.5に
## 外側の線を95%信頼区間にする。既定値は2(「点推定値±2標準誤差」)
outerCI = qt(df = summary(fit_q12)$df[2],
p = 0.025, lower.tail = FALSE),
lwdOuter = 0.5, # 外側の線の太さを0.5に
## 説明変数を日本語で表示する
newNames = list(expm = "選挙費用(100万円)", age = "年齢"),
xlab = "係数の推定値",
ylab = "説明変数",
title = "2012年衆院選における自民党候補の得票率 (%) を予測するモデル"
)
print(coefs_q12)
この図が推定結果をまとめている。縦軸にはモデルに含まれる説明変数が並んでいる。横軸は、係数の推定値の大きさを表している。図中の点は、各説明変数の係数の点推定値である。また、水平に引かれた線は、内側の太い線が50パーセント信頼区間、外側の細い線が95パーセント信頼区間である。いずれの変数についても、95パーセント信頼区間がゼロを跨いでいないので、推定された効果は5%の有意水準で統計的に有意であることがわかる。この分析の標本サイズは285、自由度調整済み決定係数は0.1である。
上のモデルの残差プロットを作る。
res_ldp <- tibble(res = fit_q12$residuals,
fitted = fit_q12$fitted.values) %>%
ggplot(aes(x = fitted, y = res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x = "予測値", y = "残差")
print(res_ldp)
特に目立ったパタンはなさそうである。残差の分散は、予測値が小さいとき(予測得票率が40%以下)にそれ以外の場合よりも小さく、誤差の分散は均一ではないかもしれない。
最後に、正規QQプロットで誤差の分布を診断する。
## 残差を標準化する
z_res <- with(fit_q12, (residuals - mean(residuals)) / sd(residuals))
## 図を描く
qq_ldp <- tibble(z_res) %>%
ggplot(aes(sample = z_res)) +
geom_abline(intercept = 0, slope = 1, # 45度線
linetype = "dashed", color = "gray") +
geom_qq(pch = 16, size = 1) +
labs(x = "標準正規分布", y = "標準化した残差の分布") +
xlim(-4, 4) + ylim(-4, 4) +
coord_fixed() # 横軸と縦軸を1:1にする
print(qq_ldp)
分布の端で残差が45度線から外れており、誤差の分布が正規分布ではない可能性が示唆される。