library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(broom) # 解析結果を tibble 形式に集約
library(gt) # 表の作成統計データ解析
第3講 サンプルコード
準備
以下で利用する共通パッケージを読み込む．1
1 必要なパッケージは適宜追加して良い．
回帰係数の推定
問題 (広告費と売上データ)
データセット https://www.statlearning.com/s/Advertising.csv は 広告費(TV,radio,newspapers)と売上の関係を調べたもの 2 である． このデータセットを用いて以下の回帰式を推定しなさい．
2 Datasets in this presentation are taken from “An Introduction to Statistical Learning, with applications in R” (Springer, 2013) with permission from the authors: G. James, D. Witten, T. Hastie and R. Tibshirani. 参考
formula = sales ~ TV
formula = sales ~ radio
formula = sales ~ TV + radio解答例
データの読み込みを行う．
adv_data <- read_csv("https://www.statlearning.com/s/Advertising.csv")GGally パッケージを利用して散布図を描く．3
3 標準設定では図は横長で描画される．この例では縦横が同じ大きさ(7inch)になるように指定している．
GGally::ggpairs(adv_data, columns = 2:5)
“TV” の宣伝費で売上を説明する．
(adv_lm1 <- lm(sales ~ TV, data = adv_data))
Call:
lm(formula = sales ~ TV, data = adv_data)
Coefficients:
(Intercept) TV
7.03259 0.04754 adv_data |>
ggplot(aes(x = TV, y = sales)) +
geom_point(colour = "orange") +
geom_smooth(method = lm, se = FALSE)
“radio” の宣伝費で売上を説明する．
(adv_lm2 <- lm(sales ~ radio, data = adv_data))
Call:
lm(formula = sales ~ radio, data = adv_data)
Coefficients:
(Intercept) radio
9.3116 0.2025 adv_data |>
ggplot(aes(x = radio, y = sales)) +
geom_point(colour = "orange") +
geom_smooth(method = lm, se = FALSE)
両者の宣伝費で売上を説明する．
(adv_lm3 <- lm(sales ~ TV + radio, data = adv_data))
Call:
lm(formula = sales ~ TV + radio, data = adv_data)
Coefficients:
(Intercept) TV radio
2.92110 0.04575 0.18799
Tip
3次元の散布図を描くには scatterplot3d パッケージを利用することができる．
library(scatterplot3d)
s3d <- with(adv_data,
scatterplot3d(
x = TV, y = radio, z = sales,
type = "p", # plotの種類: "p"点，"l"線，"h"足付き
pch = 16, # 点の種類 (?points 参照)
angle = 45, # xy平面の見る方向 (適宜調整)
highlight.3d = TRUE # 高さ(z)ごとに色を変える
))
s3d$plane3d(adv_lm3, col = "blue", # 回帰式の定める平面の追加
draw_polygon = TRUE, # 平面の塗り潰しの設定
polygon_args = list(col = rgb(0,0,1,0.2))) 
問題 (東京の気候データ)
配布したデータセット tokyo_weather.csv は 気象庁より取得した東京の気候データを回帰分析用に整理したもの 4 である． このデータセットのうち，8月の気候データを用いて，以下の回帰式を推定しなさい．
4 参考
formula = temp ~ solar + press解答例
データの読み込みを行う．
tw_data <- read_csv("data/tokyo_weather.csv")8月の”気温”を目的変数，“日射量・気圧”を説明変数とするモデルの推定を行う．
tw_formula <- temp ~ solar + press # モデル式の定義
class(tw_formula) # formula class であることを確認[1] "formula"(tw_lm <- lm(tw_formula, # 回帰係数の推定
data = tw_data,
subset = month==9)) # 9月のデータの抽出
Call:
lm(formula = tw_formula, data = tw_data, subset = month == 9)
Coefficients:
(Intercept) solar press
386.2285 0.3020 -0.3602 (tw_df <- as_tibble(model.frame(tw_lm))) # 推定に用いたデータフレームの抽出# A tibble: 30 × 3
temp solar press
<dbl> <dbl> <dbl>
1 29.2 24.0 1012.
2 29.6 22.1 1010.
3 29.1 18.6 1011.
4 26.1 7.48 1008.
5 29.3 22.6 1005.
6 27.5 13.2 1004.
7 27 11.0 1008.
8 21.9 2.1 1008.
9 24.8 8.81 1007.
10 27.8 17.6 1009.
# ℹ 20 more rowsGGally::ggpairs(tw_df) # 散布図
散布図と回帰式の定める平面の描画の3次元プロットを行う．
s3d <- with(tw_df,
scatterplot3d(
x = solar, y = press, z = temp,
type = "p", # plotの種類: "p"点，"l"線，"h"足付き
pch = 16, # 点の種類 (?points 参照)
angle = 30, # xy平面の見る方向 (適宜調整)
highlight.3d = TRUE # 高さ(z)ごとに色を変える
))
s3d$plane3d(tw_lm, col = "blue", # 回帰式の定める平面の追加
draw_polygon = TRUE, # 平面の塗り潰しの設定
polygon_args = list(col = rgb(0,0,1,0.1))) 
最小二乗推定量の性質
問題
前問の推定結果を用いて，最小二乗推定量の以下の性質を確認しなさい．
- 推定された係数が正規方程式の解となること． \begin{equation} \boldsymbol{\hat{\beta}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} \end{equation}
- あてはめ値と残差が直交すること．
- 回帰式が標本平均を通ること．
解答例
広告費と売上データの回帰分析において， 回帰係数と正規方程式の解の一致は以下の計算で確認できる．
(beta <- coef(adv_lm3)) # 推定された回帰係数(Intercept) TV radio
2.92109991 0.04575482 0.18799423 X <- model.matrix(adv_lm3) # デザイン行列
Y <- model.frame(adv_lm3)[[1]] # 目的変数 (データフレームの1列目のベクトル)
solve(crossprod(X)) %*% crossprod(X, Y) # 正規方程式の解 [,1]
(Intercept) 2.92109991
TV 0.04575482
radio 0.18799423あてはめ値と残差の直交性はベクトルの内積を利用すればよい．
yhat <- fitted(adv_lm3) # あてはめ値
ehat <- resid(adv_lm3) # 残差
yhat %*% ehat # 直交すれば内積はO [,1]
[1,] -2.337233e-13回帰式が標本平均を通ることは説明変数の平均のあてはめ値を求めて比較すればよい．
colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値 [,1]
[1,] 14.0225mean(Y) # 目的変数の標本平均 [1] 14.0225東京の気候データにおいても同様に確認できる．
(beta <- coef(tw_lm)) # 推定された回帰係数(Intercept) solar press
386.2284847 0.3019861 -0.3602392 (X <- model.matrix(tw_lm)) # デザイン行列 (Intercept) solar press
244 1 24.01 1012.1
245 1 22.07 1010.3
246 1 18.64 1010.6
247 1 7.48 1007.5
248 1 22.58 1005.2
249 1 13.17 1003.6
250 1 11.01 1007.9
251 1 2.10 1007.8
252 1 8.81 1006.8
253 1 17.57 1009.1
254 1 17.19 1010.1
255 1 20.02 1010.0
256 1 22.00 1010.9
257 1 14.54 1009.9
258 1 9.21 1010.9
259 1 11.78 1011.5
260 1 14.84 1011.5
261 1 19.59 1011.6
262 1 19.93 1010.1
263 1 10.65 1009.3
264 1 6.65 1006.7
265 1 6.83 1008.1
266 1 4.48 1012.5
267 1 15.81 1017.2
268 1 15.49 1017.1
269 1 16.08 1012.7
270 1 11.59 1008.1
271 1 14.03 1004.7
272 1 10.11 1009.0
273 1 7.98 1007.5
attr(,"assign")
[1] 0 1 2(Y <- model.frame(tw_lm)[[1]]) # 目的変数 (データフレームの1列目のベクトル) [1] 29.2 29.6 29.1 26.1 29.3 27.5 27.0 21.9 24.8 27.8 28.1 27.7 28.0 28.2 27.4
[16] 27.9 28.7 28.9 29.0 27.2 26.7 24.8 22.1 22.2 22.4 24.6 25.3 27.4 26.3 25.6solve(crossprod(X)) %*% crossprod(X, Y) # 正規方程式の解 [,1]
(Intercept) 386.2284847
solar 0.3019861
press -0.3602392yhat <- fitted(tw_lm) # あてはめ値
ehat <- resid(tw_lm) # 残差
yhat %*% ehat # 直交すれば内積はO(に近い値)となる [,1]
[1,] 1.123776e-14colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値 [,1]
[1,] 26.69333mean(Y) # 目的変数の標本平均 [1] 26.69333
Tip
関数 stats::model.matrix() の返値は matrix class である． 関数 stats::model.frame() の返値は data.frame class であり，model.frame(...)[[k]] は data.frame ではなく vector になる． tibble(...)[[k]] は vector となるが，tibble(...)[k] は vector ではなく1列の tibble になるので注意する． 関数 base::crossprod() の引数は厳密に vector/matrix であることを要請するが， 例えば adv_data は全て数値なので行列に変換して計算できる．
crossprod(adv_data) # データフレームのままなので計算できない
crossprod(as.matrix(adv_data)) # 計算できる残差の分解
問題
前問の結果を用いて，残差の性質として以下の分解が成り立つことを確認しなさい．
\begin{equation} (\boldsymbol{y}-\bar{\boldsymbol{y}})^{\mathsf{T}} (\boldsymbol{y}-\bar{\boldsymbol{y}}) = (\boldsymbol{y}-\boldsymbol{\hat{y}})^{\mathsf{T}} (\boldsymbol{y}-\boldsymbol{\hat{y}})+ (\boldsymbol{\hat{y}}-\bar{\boldsymbol{y}})^{\mathsf{T}} (\boldsymbol{\hat{y}}-\bar{\boldsymbol{y}}) \end{equation}
\begin{equation} S_y=S+S_r \end{equation}
解答例
広告費と売上データでは以下のように確認できる．
summary(adv_lm3)
Call:
lm(formula = sales ~ TV + radio, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-8.7977 -0.8752 0.2422 1.1708 2.8328
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.92110 0.29449 9.919 <2e-16 ***
TV 0.04575 0.00139 32.909 <2e-16 ***
radio 0.18799 0.00804 23.382 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.681 on 197 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16Y <- model.frame(adv_lm3)[[1]] # 目的変数の取得
(Sy <- sum((Y-mean(Y))^2)) # 目的変数のばらつき[1] 5417.149(S <- sum(resid(adv_lm3)^2)) # 残差のばらつき[1] 556.914(Sr <- sum((fitted(adv_lm3)-mean(Y))^2)) # 回帰のばらつき[1] 4860.235S+Sr # Sy と同じになっている[1] 5417.149東京の気候データも同様に確認できる． 以下の例では目的変数を推定結果に含める lm() のオプションを利用している．
tw_formula # モデルの確認temp ~ solar + presstw_lm <- lm(tw_formula,
data = tw_data,
subset = month == 9, # 9月のデータの抽出
y = TRUE) # 目的変数をyとして返すように指定
Y <- with(tw_lm, y) # 目的変数の取得 (tw_lm$y でも可)
(Sy <- sum((Y-mean(Y))^2)) # 目的変数のばらつき[1] 151.0787(S <- sum(resid(tw_lm)^2)) # 残差のばらつき[1] 55.66168(Sr <- sum((fitted(tw_lm)-mean(Y))^2)) # 回帰のばらつき[1] 95.41699S+Sr # Sy と同じになっている[1] 151.0787
Tip
関数 base::summary() で得られる情報は， 関数 broom::tidy() および関数 broom::glance() を用いると tibble 形式のデータフレームとして取得できる．
関数 stats::resid() や stats::fitted() で得られる値は同様に 関数 broom::augumet() で取得できる．
tidy(adv_lm3) # 推定量に関する情報# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.92 0.294 9.92 4.57e-19
2 TV 0.0458 0.00139 32.9 5.44e-82
3 radio 0.188 0.00804 23.4 9.78e-59glance(adv_lm3) # 推定結果を評価する統計量# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.897 0.896 1.68 860. 4.83e-98 2 -386. 780. 794.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>augment(adv_lm3) # 残差やあてはめ値など# A tibble: 200 × 9
sales TV radio .fitted .resid .hat .sigma .cooksd .std.resid
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 22.1 230. 37.8 20.6 1.54 0.0140 1.68 0.00406 0.925
2 10.4 44.5 39.3 12.3 -1.95 0.0188 1.68 0.00871 -1.17
3 9.3 17.2 45.9 12.3 -3.04 0.0295 1.67 0.0341 -1.83
4 18.5 152. 41.3 17.6 0.883 0.0124 1.68 0.00117 0.528
5 12.9 181. 10.8 13.2 -0.324 0.00951 1.69 0.000120 -0.194
6 7.2 8.7 48.9 12.5 -5.31 0.0347 1.64 0.124 -3.22
7 11.8 57.5 32.8 11.7 0.0818 0.0129 1.69 0.0000105 0.0490
8 13.2 120. 19.6 12.1 1.09 0.00576 1.68 0.000823 0.653
9 4.8 8.6 2.1 3.71 1.09 0.0271 1.68 0.00401 0.658
10 10.6 200. 2.6 12.6 -1.95 0.0171 1.68 0.00797 -1.17
# ℹ 190 more rows決定係数によるモデルの比較
問題
決定係数を用いてモデルの比較を行いなさい．
広告費と売上データ
sales ~ TV sales ~ radio sales ~ TV + radio東京の8月の気候データ
temp ~ solar temp ~ solar + press temp ~ solar + press + cloud
解答例
広告費と売上データでのモデルの比較は以下のとおりである．
adv_formula1 <- sales ~ TV
adv_formula2 <- sales ~ radio
adv_formula3 <- sales ~ TV + radio
adv_lm1 <- lm(adv_formula1, data=adv_data, y = TRUE)
adv_lm2 <- lm(adv_formula2, data=adv_data, y = TRUE)
adv_lm3 <- lm(adv_formula3, data=adv_data, y = TRUE)
summary(adv_lm1)$adj.r.squared # 自由度調整済み決定係数[1] 0.6099148summary(adv_lm2)$adj.r.squared # (model1より減少)[1] 0.3286589summary(adv_lm3)$adj.r.squared # (model1より上昇)[1] 0.8961505予測値と実測値の比較は以下のようにすれば良い．
adv_data |>
mutate(model1 = fitted(adv_lm1), # モデルごとに予測値をデータフレームに追加
model2 = fitted(adv_lm2),
model3 = fitted(adv_lm3)) |>
pivot_longer(starts_with("model"), # モデルをラベルとして予測値をまとめる
names_to = "model", values_to = "fitted") |>
ggplot(aes(x = sales, y = fitted)) + # 実測値をx軸，予測値をy軸で表示
geom_abline(slope = 1, intercept = 0, colour = "red") + # 基準線
geom_point(aes(colour = model, shape = model)) + # 予測値をモデル別に表示
labs(y = "fitted values")
東京の気候データも同様で，モデルの比較は以下のとおりである．
tw_subset <- tw_data |> filter(month == 9) # 8月のデータの抽出
tw_formula1 <- temp ~ solar
tw_formula2 <- temp ~ solar + press
tw_formula3 <- temp ~ solar + press + cloud
tw_lm1 <- lm(tw_formula1, data=tw_subset, y = TRUE)
tw_lm2 <- lm(tw_formula2, data=tw_subset, y = TRUE)
tw_lm3 <- lm(tw_formula3, data=tw_subset, y = TRUE)
summary(tw_lm1)$adj.r.squared # 自由度調整済み決定係数[1] 0.393239summary(tw_lm2)$adj.r.squared # (model1より上昇)[1] 0.6042806summary(tw_lm3)$adj.r.squared # (model2より上昇)[1] 0.5906416比較表を作成するには，推定されたモデルの情報を取得して tibble 形式のデータフレームを作成する broom::glance() を使って，例えば次のようにすれば良い．
bind_rows( # 3つのモデルの評価指標を結合
glance(tw_lm1),
glance(tw_lm2),
glance(tw_lm3)) |>
mutate(model = map_vec(c(tw_formula1,tw_formula2,tw_formula3), deparse),
.before = 1) |> # モデル式を文字列に変換して先頭列に置く
select(1:3) |> # 決定係数と自由度調整済み決定係数を表示
set_names(c("モデル","決定係数","自由度調整済み決定係数")) |>
gt() |>
fmt_number(decimals = 3) # 表示を3桁に| モデル | 決定係数 | 自由度調整済み決定係数 |
|---|---|---|
| temp ~ solar | 0.414 | 0.393 |
| temp ~ solar + press | 0.632 | 0.604 |
| temp ~ solar + press + cloud | 0.633 | 0.591 |
予測値と実測値の比較は以下のようになる．
tw_subset |>
mutate(model1 = fitted(tw_lm1), # モデルごとに予測値をデータフレームに追加
model2 = fitted(tw_lm2),
model3 = fitted(tw_lm3)) |>
pivot_longer(starts_with("model"), # モデルをラベルとして予測値をまとめる
names_to = "model", values_to = "fitted") |>
ggplot(aes(x = temp, y = fitted)) + # 気温の実測値をx軸，予測値をy軸で表示
geom_abline(slope = 1, intercept = 0, colour = "red") + # 基準線
geom_point(aes(colour = model, shape = model)) + # 予測値をモデル別に表示
labs(y = "fitted values")