library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(broom) # 解析結果を tibble 形式に集約
library(ggfortify) # 診断プロットの描画
library(gt) # 表の作成
library(gtsummary) # 分析結果の表の作成第5講 回帰分析
予測と発展的なモデル
準備
以下で利用する共通パッケージを読み込む．
回帰式を用いた予測
予測は関数 broom::augment() を用いて行うことができる．
あてはめ値と予測値
#' モデルの作成
toy_train <- tibble(x1 = ..., x2 = ..., y = ...)
toy_lm <- lm(y ~ x1 + x2, data = toy_train)
toy_train_fitted <- augment(toy_lm) # あてはめ値の計算
#' 新しいデータの予測
toy_test <- tibble(x1 = ..., x2 = ...) # 予測したいデータの説明変数
toy_test_fitted <- augment(toy_lm, # 予測値の計算
newdata = toy_test)
toy_test_conf <- augment(toy_lm, newdata = toy_test, # 信頼区間
interval = "confidence", conf.level = 0.95)
toy_test_pred <- augment(toy_lm, newdata = toy_test, # 予測区間
interval = "prediction", conf.level = 0.95) 信頼区間，予測区間の水準の既定値は0.95である．
東京の気候データによる例
9,10月のデータでモデルを構築し，8,11月のデータを予測する．
#' データの整理
tw_data <- read_csv("data/tokyo_weather.csv")
tw_train <- tw_data |> # モデル推定用データ
filter(month %in% c(9,10)) # %in% は集合に含むかどうかを判定
tw_test <- tw_data |> # 予測用データ
filter(month %in% c(8,11))
#' モデルの構築
tw_model <- temp ~ solar + press # モデルの定義
tw_lm <- lm(tw_model, data = tw_train) # モデルの推定
tidy(tw_lm) # 回帰係数の評価# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 363. 112. 3.23 0.00204
2 solar 0.293 0.104 2.80 0.00689
3 press -0.341 0.111 -3.06 0.00336glance(tw_lm) # モデルの評価# 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.209 0.182 4.13 7.67 0.00111 2 -172. 351. 359.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>#' あてはめ値と予測値の計算
tw_train_fitted <- augment(tw_lm, newdata = tw_train) # あてはめ値
tw_test_fitted <- augment(tw_lm, newdata = tw_test) # 予測値グラフでの表示の例は以下のようになる．
#' 予測結果を図示
bind_rows(tw_train_fitted, tw_test_fitted) |> # 2つのデータフレームを結合
mutate(month = as_factor(month)) |> # 月を因子化して表示に利用
ggplot(aes(x = .fitted, y = temp)) +
geom_point(aes(colour = month, shape = month)) + # 月ごとに色と形を変える
geom_abline(slope = 1, intercept = 0, # 予測が完全に正しい場合のガイド線
colour = "gray") +
labs(x = "fitted", y = "observed")
関数 lubridate::month() を用いると月を文字列のラベルとすることもできる．
bind_rows(tw_train_fitted, tw_test_fitted) |>
mutate(month = month(month, label = TRUE)) |> # 文字にする場合
ggplot(aes(x = .fitted, y = temp)) +
geom_point(aes(colour = month)) + # 月ごとに色を変える
geom_abline(slope = 1, intercept = 0, # 予測が完全に正しい場合のガイド線
colour = "gray") +
labs(x = "fitted", y = "observed")
この場合は順序付きの因子になるので，‘shape’ を利用すると警告が出る．
関数 ggplot2::geom_errorbar() を用いると 区間を図示することができる．
区間表示のための関数
geom_errorbar(
mapping = NULL,
data = NULL,
stat = "identity",
position = "identity",
...,
na.rm = FALSE,
orientation = NA,
show.legend = NA,
inherit.aes = TRUE
)
#' mapping: 区間を表すために xmin,xmax または ymin,ymax を与える
#' data: データフレーム
#' ...: その他の描画オプション
#' orientation: 特別な場合に指定 (一般に向きは mapping で自動的決定)
#' 詳細は ?ggplot2::geom_errorbar関数 broom::augment() の出力の .lower/.upper 列を用いればよい．
問題
東京の気候データを用いて以下の実験を試みなさい．
- 8月のデータで回帰式を推定する
- 上記のモデルで9月のデータを予測する
特定の月のデータを取り出すには，例えば以下のようにすればよい
tw_data <- read_csv("data/tokyo_weather.csv")
tw_train <- tw_data |> filter(month == 8) # 単一の数字と比較
tw_test <- tw_data |> filter(month %in% c(9,10)) # 集合と比解答欄
解答例
東京の気候データによる分析において，信頼区間と予測区間の計算を行う．
推定用データと予測用データを準備し，推定用データを用いて回帰式の推定を行う．
tw_data <- read_csv("data/tokyo_weather.csv")
tw_train <- tw_data |> filter(month %in% 8) # 推定用データ
tw_test <- tw_data |> filter(month %in% 9) # 予測用データ
tw_model <- temp ~ solar + press + cloud # モデルの定義
tw_lm <- lm(tw_model, data = tw_train) # モデルの推定信頼区間の計算を行う．
tw_train_conf <-
augment(tw_lm, newdata = tw_train, # 推定用データに
interval = "confidence") # あてはめ値と信頼区間を付加
tw_test_conf <-
augment(tw_lm, newdata = tw_test, # 予測用データ(新規データ)に
interval = "confidence") # あてはめ値と信頼区間を付加予測区間の計算を行う．
tw_train_pred <-
augment(tw_lm, newdata = tw_train, # 推定用データに
interval = "prediction") # あてはめ値と予測区間を付加
tw_test_pred <-
augment(tw_lm, newdata = tw_test, # 予測用データ(新規データ)に
interval = "prediction") # あてはめ値と予測区間を付加8月のデータで推定したモデルで8月をあてはめた信頼区間を描画する．
tw_train_conf |>
ggplot(aes(x = day, y = temp)) +
geom_point(colour = "red", shape = 16) +
geom_point(aes(y = .fitted), colour = "blue") +
geom_errorbar(aes(ymin = .lower, ymax = .upper), colour = "royalblue") +
ylim(c(20,34)) + # 以降比較のため4つのグラフの値域を揃える
labs(x = "August", y = "Temperature", title = "Confidence Interval")
8月のデータで推定したモデルで8月をあてはめた予測区間を描画する．
tw_train_pred |>
ggplot(aes(x = day, y = temp)) +
geom_point(colour = "red", shape = 16) +
geom_point(aes(y = .fitted), colour = "blue") +
geom_errorbar(aes(ymin = .lower, ymax = .upper), colour = "steelblue") +
ylim(c(20,34)) +
labs(x = "August", y = "Temperature", title = "Prediction Interval")
8月のモデルで9月をあてはめた信頼区間を描画する．
tw_test_conf |>
ggplot(aes(x = day, y = temp)) +
geom_point(colour = "red", shape = 16) +
geom_point(aes(y = .fitted), colour = "blue") +
geom_errorbar(aes(ymin = .lower, ymax = .upper), colour = "royalblue") +
ylim(c(20,34)) +
labs(x = "September", y = "Temperature", title = "Confidence Interval")
8月のモデルで9月をあてはめた予測区間を描画する．
tw_test_pred |>
ggplot(aes(x = day, y = temp)) +
geom_point(colour = "red", shape = 16) +
geom_point(aes(y = .fitted), colour = "blue") +
geom_errorbar(aes(ymin = .lower, ymax = .upper), colour = "steelblue") +
ylim(c(20,34)) +
labs(x = "September", y = "Temperature", title = "Prediction Interval")
非線形性を含むモデルとカテゴリカル変数の扱い
線形でないモデル式を記述するためには特殊な記法がある．
- 非線形変換はそのまま関数を記述すればよい．
- 1つの変数の多項式は関数 I() を用いる．
線形でないモデル式の書き方
#' 目的変数 Y, 説明変数 X1,X2,X3
#' 交互作用を含む式 (formula) の書き方
Y ~ X1 + X1:X2 # X1 + X1*X2
Y ~ X1 * X2 # X1 + X2 + X1*X2
Y ~ (X1 + X2 + X3)^2 # X1 + X2 + X3 + X1*X2 + X2*X3 + X3*X1
#' 非線形変換を含む式 (formula) の書き方
Y ~ f(X1) # f(X1) (fは任意の関数)
Y ~ X1 + I(X2^2) # X1 + X2^2カテゴリカル変数の取り扱いには若干の注意が必要である．
- 何も宣言しなくても通常は適切に対応してくれる．
- 陽に扱う場合は関数
factor()を利用する．
因子属性の与え方の例
X <- c("A", "S", "A", "B", "D")
Y <- c(85, 100, 80, 70, 30)
toy_data1 <- tibble(X, Y)
toy_data2 <- toy_data1 |> # 因子化
mutate(X2 = factor(X)) # 関数as_factor()を用いてもよい
glimpse(toy_data2) # 作成したデータフレームの素性を見る(pillar::glimpse())Rows: 5
Columns: 3
$ X <chr> "A", "S", "A", "B", "D"
$ Y <dbl> 85, 100, 80, 70, 30
$ X2 <fct> A, S, A, B, Dtoy_data3 <- toy_data2 |> # 順序付き(levels)の因子化
mutate(X3 = factor(X, levels=c("S","A","B","C","D")))
glimpse(toy_data3) # toy_data2とはfactorの順序が異なるRows: 5
Columns: 4
$ X <chr> "A", "S", "A", "B", "D"
$ Y <dbl> 85, 100, 80, 70, 30
$ X2 <fct> A, S, A, B, D
$ X3 <fct> A, S, A, B, Dtoy_data4 <- toy_data2 |>
mutate(Y2 = factor(Y > 60)) # 条件による因子化
glimpse(toy_data4) # 条件の真偽で2値に類別されるRows: 5
Columns: 4
$ X <chr> "A", "S", "A", "B", "D"
$ Y <dbl> 85, 100, 80, 70, 30
$ X2 <fct> A, S, A, B, D
$ Y2 <fct> TRUE, TRUE, TRUE, TRUE, FALSE問題
東京の気候データ(9-11月)を用いて，気温を回帰する以下のモデルを検討しなさい．
- 日射量，気圧，湿度の線形回帰モデル
- 湿度の対数を考えた線形回帰モデル
- 最初のモデルにそれぞれの交互作用を加えたモデル
東京の気候データ(1年分)を用いて，気温を回帰する以下のモデルを検討しなさい．
- 降水の有無を表すカテゴリカル変数を用いたモデル (雨が降ると気温が変化することを検証する)
- 上記に月をカテゴリカル変数として加えたモデル (月毎の気温の差を考慮する)
解答欄
解答例
9月から11月のデータを用いた分析を行うために，データを整理する．
tw_subset <- tw_data |> filter(month %in% 9:11)日射量，気圧，湿度を用いた回帰モデルの検討を行う． まず，解析の対象となる変数の散布図を描いておく．
tw_subset |>
select(temp, solar, press, humid) |>
GGally::ggpairs(columnLabels = c("気温","日射量","気圧","湿度"))
モデルを比較する．
#' 基本となる線形回帰モデル
summary(lm(temp ~ solar + press + humid, data = tw_subset))
Call:
lm(formula = temp ~ solar + press + humid, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.7787 -2.4153 0.0586 2.5814 7.1178
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 177.64597 74.23672 2.393 0.0189 *
solar 0.78130 0.08492 9.200 1.75e-14 ***
press -0.18425 0.07274 -2.533 0.0131 *
humid 0.26864 0.03205 8.383 8.23e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.761 on 87 degrees of freedom
Multiple R-squared: 0.5976, Adjusted R-squared: 0.5838
F-statistic: 43.07 on 3 and 87 DF, p-value: < 2.2e-16#' 湿度の対数を考えた線形回帰モデル
summary(lm(temp ~ solar + press + log(humid), data = tw_subset))
Call:
lm(formula = temp ~ solar + press + log(humid), data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.868 -2.478 -0.062 2.814 6.635
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 152.12678 76.72547 1.983 0.05055 .
solar 0.73316 0.08461 8.665 2.18e-13 ***
press -0.21081 0.07366 -2.862 0.00527 **
log(humid) 16.98674 2.12197 8.005 4.83e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.837 on 87 degrees of freedom
Multiple R-squared: 0.5812, Adjusted R-squared: 0.5667
F-statistic: 40.24 on 3 and 87 DF, p-value: < 2.2e-16#' 最初のモデルにそれぞれの交互作用を加えたモデル (書き方はいろいろある)
summary(lm(temp ~ (solar + press + humid)^2, data = tw_subset))
Call:
lm(formula = temp ~ (solar + press + humid)^2, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.8008 -2.5456 -0.0834 2.5669 6.9689
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.052e+02 5.209e+02 -0.394 0.695
solar 1.333e+01 1.841e+01 0.724 0.471
press 1.979e-01 5.142e-01 0.385 0.701
humid 3.847e+00 5.796e+00 0.664 0.509
solar:press -1.271e-02 1.807e-02 -0.704 0.484
solar:humid 4.133e-03 5.894e-03 0.701 0.485
press:humid -3.586e-03 5.732e-03 -0.626 0.533
Residual standard error: 3.793 on 84 degrees of freedom
Multiple R-squared: 0.6047, Adjusted R-squared: 0.5765
F-statistic: 21.42 on 6 and 84 DF, p-value: 4.2e-15#' 更に3つの変数の積を加えたモデル
summary(lm(temp ~ solar * press * humid, data = tw_subset))
Call:
lm(formula = temp ~ solar * press * humid, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.7795 -2.6021 -0.0871 2.5201 6.9460
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.786e+01 9.223e+02 -0.052 0.959
solar -2.966e+00 8.078e+01 -0.037 0.971
press 4.228e-02 9.118e-01 0.046 0.963
humid 1.733e+00 1.175e+01 0.148 0.883
solar:press 3.413e-03 7.989e-02 0.043 0.966
solar:humid 2.331e-01 1.105e+00 0.211 0.833
press:humid -1.494e-03 1.162e-02 -0.129 0.898
solar:press:humid -2.266e-04 1.094e-03 -0.207 0.836
Residual standard error: 3.815 on 83 degrees of freedom
Multiple R-squared: 0.6049, Adjusted R-squared: 0.5716
F-statistic: 18.15 on 7 and 83 DF, p-value: 2.007e-14基本的なモデルと最後のモデルの視覚的な評価 (診断プロット) を行う．
autoplot(lm(temp ~ solar + press + humid, data = tw_subset))
autoplot(lm(temp ~ solar * press * humid, data = tw_subset))
雨の有無(カテゴリカル変数)と気温の関係を分析する． 雨の有無および月(整数値)をダミー化(因子化)1する．
1 論理値は因子として扱われるが，明示的に因子としたい場合は factor(rain > 0) としても良い． 因子の属性を与える関数としては base::as.factor() と forcats::as_factor() があるが，前者は文字列を辞書順に並べ変えてから，後者はデータフレームでの出現順に因子化する．回帰の計算においては特に問題ないが，因子化の順序はグラフなど凡例の順番に影響する．
tw_data_fact <- tw_data |>
mutate(rain = rain > 0, # 降雨の有無
month = as_factor(month)) # 月 通年でのモデルを検討する．
summary(lm(temp ~ rain, data = tw_data_fact))
Call:
lm(formula = temp ~ rain, data = tw_data_fact)
Residuals:
Min 1Q Median 3Q Max
-17.593 -6.793 -0.060 7.340 14.507
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.2932 0.4995 34.620 <2e-16 ***
rainTRUE 1.4668 0.9543 1.537 0.125
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.131 on 363 degrees of freedom
Multiple R-squared: 0.006466, Adjusted R-squared: 0.003729
F-statistic: 2.362 on 1 and 363 DF, p-value: 0.1252通年では雨と気温の関係は積極的に支持されない．
月を説明変数として加えると月毎の気温の差を考慮した回帰式が推定される．
summary(lm(temp ~ rain + month, data = tw_data_fact))
Call:
lm(formula = temp ~ rain + month, data = tw_data_fact)
Residuals:
Min 1Q Median 3Q Max
-6.794 -1.622 -0.117 1.497 8.528
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.8726 0.4460 13.168 < 2e-16 ***
rainTRUE -0.7213 0.3001 -2.404 0.0168 *
month2 1.5083 0.6409 2.354 0.0191 *
month3 7.3318 0.6246 11.738 < 2e-16 ***
month4 10.6265 0.6287 16.903 < 2e-16 ***
month5 13.3712 0.6253 21.385 < 2e-16 ***
month6 17.7021 0.6353 27.864 < 2e-16 ***
month7 22.9302 0.6241 36.743 < 2e-16 ***
month8 23.6067 0.6253 37.754 < 2e-16 ***
month9 21.0853 0.6300 33.469 < 2e-16 ***
month10 13.1445 0.6235 21.083 < 2e-16 ***
month11 8.6210 0.6288 13.710 < 2e-16 ***
month12 3.6180 0.6237 5.801 1.47e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.454 on 352 degrees of freedom
Multiple R-squared: 0.9122, Adjusted R-squared: 0.9092
F-statistic: 304.9 on 12 and 352 DF, p-value: < 2.2e-16月毎に比較した結果，雨の日の方が気温が低いことが示唆される．
診断プロットは以下のようになる．
autoplot(lm(temp ~ rain + month, data = tw_data_fact))
参考 : その他の関数・パッケージ
- 変数が増えるとモデルの比較が困難となる．
- 関数
stats::step()で自動化することができる．
モデルの探索の例
adv_data <- read_csv('https://www.statlearning.com/s/Advertising.csv')
summary(lm(sales ~ radio, data = adv_data))
Call:
lm(formula = sales ~ radio, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-15.7305 -2.1324 0.7707 2.7775 8.1810
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.31164 0.56290 16.542 <2e-16 ***
radio 0.20250 0.02041 9.921 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.275 on 198 degrees of freedom
Multiple R-squared: 0.332, Adjusted R-squared: 0.3287
F-statistic: 98.42 on 1 and 198 DF, p-value: < 2.2e-16summary(lm(sales ~ TV + radio, data = adv_data))
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-16summary(lm(sales ~ TV + radio + newspaper, data = adv_data))
Call:
lm(formula = sales ~ TV + radio + newspaper, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-8.8277 -0.8908 0.2418 1.1893 2.8292
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.938889 0.311908 9.422 <2e-16 ***
TV 0.045765 0.001395 32.809 <2e-16 ***
radio 0.188530 0.008611 21.893 <2e-16 ***
newspaper -0.001037 0.005871 -0.177 0.86
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16summary(adv_init <- lm(sales ~ TV * radio * newspaper, data = adv_data))
Call:
lm(formula = sales ~ TV * radio * newspaper, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-5.8955 -0.3883 0.1938 0.5865 1.5240
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.556e+00 4.655e-01 14.083 < 2e-16 ***
TV 1.971e-02 2.719e-03 7.250 9.95e-12 ***
radio 1.962e-02 1.639e-02 1.197 0.233
newspaper 1.311e-02 1.721e-02 0.761 0.447
TV:radio 1.162e-03 9.753e-05 11.909 < 2e-16 ***
TV:newspaper -5.545e-05 9.326e-05 -0.595 0.553
radio:newspaper 9.063e-06 4.831e-04 0.019 0.985
TV:radio:newspaper -7.610e-07 2.700e-06 -0.282 0.778
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9406 on 192 degrees of freedom
Multiple R-squared: 0.9686, Adjusted R-squared: 0.9675
F-statistic: 847.3 on 7 and 192 DF, p-value: < 2.2e-16adv_opt <- step(adv_init) # 最大のモデルから削減増加による探索Start: AIC=-16.67
sales ~ TV * radio * newspaper
Df Sum of Sq RSS AIC
- TV:radio:newspaper 1 0.070261 169.93 -18.586
<none> 169.86 -16.669
Step: AIC=-18.59
sales ~ TV + radio + newspaper + TV:radio + TV:newspaper + radio:newspaper
Df Sum of Sq RSS AIC
- radio:newspaper 1 0.19 170.12 -20.363
<none> 169.93 -18.586
- TV:newspaper 1 4.37 174.30 -15.511
- TV:radio 1 349.71 519.64 202.965
Step: AIC=-20.36
sales ~ TV + radio + newspaper + TV:radio + TV:newspaper
Df Sum of Sq RSS AIC
<none> 170.12 -20.363
- TV:newspaper 1 4.19 174.31 -17.494
- TV:radio 1 352.83 522.95 202.234summary(adv_opt) # 探索された(準)最適なモデルの確認
Call:
lm(formula = sales ~ TV + radio + newspaper + TV:radio + TV:newspaper,
data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-5.9019 -0.3818 0.1937 0.5741 1.4839
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.541e+00 2.652e-01 24.668 <2e-16 ***
TV 2.035e-02 1.605e-03 12.675 <2e-16 ***
radio 2.018e-02 9.734e-03 2.073 0.0395 *
newspaper 1.342e-02 6.377e-03 2.105 0.0366 *
TV:radio 1.136e-03 5.664e-05 20.059 <2e-16 ***
TV:newspaper -7.719e-05 3.531e-05 -2.187 0.0300 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9364 on 194 degrees of freedom
Multiple R-squared: 0.9686, Adjusted R-squared: 0.9678
F-statistic: 1197 on 5 and 194 DF, p-value: < 2.2e-16全探索ではないので最適とは限らないことに注意は必要である．
回帰分析のためのパッケージとして package::car が提供されており， 以下のことが可能である．
- 回帰モデルの評価
- 与えられたデータの再現
- 新しいデータの予測
- モデルの再構築のための視覚化
- residual plots: 説明変数・予測値と残差の関係
- marginal-model plots: 説明変数と目的変数・モデルの関係
- added-variable plots: 説明変数・目的変数をその他の変数で回帰したときの残差の関係
- component+residual plots: 説明変数とそれ以外の説明変数による残差の関係 などが用意されている
参考 : 人工データを用いた base R による予測の例
以下はあくまで例なので，これを参考に自由に数値実験を設計して下さい．
人工データの生成例
#' モデル: y = -1 + 2 * x1
set.seed(1515) # 乱数のシード値(適宜変更せよ)
n <- 50 # データ数の設定
x_obs <- tibble(x0 = 1,
x1 = runif(n), # 説明変数1
x2 = runif(n)) # 説明変数2
beta <- c(-1, 2, 0) # (切片, x1の係数, x2の係数) 実質x2は使われていない
sigma <- 1/2 # 誤差の標準偏差
epsilon <- rnorm(nrow(x_obs), sd = sigma) # 誤差項
toy_data <- x_obs |> # 目的変数の観測値を追加
mutate(y = as.vector(as.matrix(x_obs) %*% beta) + epsilon)
モデルの推定と評価の例
toy_lm1 <- lm(y ~ x1, data = toy_data) # x1のみの正しいモデル
toy_lm2 <- lm(y ~ x1 + x2, data = toy_data) # x1とx2の冗長なモデル
toy_lm3 <- lm(y ~ x2, data = toy_data) # x2のみの誤ったモデル
summary(toy_lm1)
Call:
lm(formula = y ~ x1, data = toy_data)
Residuals:
Min 1Q Median 3Q Max
-0.81627 -0.33791 -0.05144 0.29177 1.07659
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.9059 0.1171 -7.736 5.54e-10 ***
x1 1.9957 0.2155 9.259 2.96e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.466 on 48 degrees of freedom
Multiple R-squared: 0.6411, Adjusted R-squared: 0.6336
F-statistic: 85.74 on 1 and 48 DF, p-value: 2.956e-12summary(toy_lm2)
Call:
lm(formula = y ~ x1 + x2, data = toy_data)
Residuals:
Min 1Q Median 3Q Max
-0.86076 -0.34121 -0.05143 0.29323 1.05813
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.9463 0.1473 -6.422 6.16e-08 ***
x1 1.9764 0.2213 8.929 1.09e-11 ***
x2 0.1043 0.2273 0.459 0.649
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4699 on 47 degrees of freedom
Multiple R-squared: 0.6427, Adjusted R-squared: 0.6275
F-statistic: 42.27 on 2 and 47 DF, p-value: 3.139e-11summary(toy_lm3)
Call:
lm(formula = y ~ x2, data = toy_data)
Residuals:
Min 1Q Median 3Q Max
-1.77397 -0.47455 0.02159 0.51110 1.75063
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.2398 0.2020 -1.187 0.241
x2 0.4892 0.3627 1.349 0.184
Residual standard error: 0.7635 on 48 degrees of freedom
Multiple R-squared: 0.03653, Adjusted R-squared: 0.01645
F-statistic: 1.82 on 1 and 48 DF, p-value: 0.1837表としてまとめるには，例えば以下のようにすればよい．
my_gts <- function(x){
tbl_regression(x, intercept = TRUE) |>
add_significance_stars(
pattern = "{estimate} ({conf.low}, {conf.high}){stars}",
hide_ci = TRUE, hide_se = TRUE
) |>
modify_header(estimate = "**Beta (95% CI)**") |>
modify_footnote(estimate = "CI = Confidence Interval",
abbreviation = TRUE) |>
add_glance_table(include = c(r.squared,
adj.r.squared,
statistic,
p.value)) }
tbl_merge(
tbls = list(
my_gts(toy_lm1),
my_gts(toy_lm2),
my_gts(toy_lm3)),
tab_spanner = c("モデル1","モデル2","モデル3")) |>
modify_table_body( ~ .x |>
arrange(factor(variable,
levels = c("(Intercept)","x1","x2"))))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
|---|---|---|---|
| Beta (95% CI)1 | Beta (95% CI)1 | Beta (95% CI)1 | |
| (Intercept) | -0.91 (-1.1, -0.67)*** | -0.95 (-1.2, -0.65)*** | -0.24 (-0.65, 0.17) |
| x1 | 2.0 (1.6, 2.4)*** | 2.0 (1.5, 2.4)*** | |
| x2 | 0.10 (-0.35, 0.56) | 0.49 (-0.24, 1.2) | |
| R² | 0.641 | 0.643 | 0.037 |
| Adjusted R² | 0.634 | 0.627 | 0.016 |
| Statistic | 85.7 | 42.3 | 1.82 |
| p-value | <0.001 | <0.001 | 0.2 |
| 1 p<0.05; p<0.01; p<0.001 | |||
| Abbreviation: CI = Confidence Interval | |||
予測の例
新規データを作成し，予測の仕組みを確認する．
x_new <- tibble(x0 = 1,
x1 = runif(n),
x2 = runif(n,-10,10))
#' 新規データに対する目的変数の真値 (誤差なし)
y_tilde <- as.vector(as.matrix(x_new) %*% beta)
#' 各モデルでの予測
y_hat1 <- predict(toy_lm1, newdata = x_new) # x_lm1による予測値
y_hat2 <- predict(toy_lm2, newdata = x_new) # x_lm2による予測値
y_hat3 <- predict(toy_lm3, newdata = x_new) # x_lm3による予測値
可視化の例
散布図を用いて，それぞれのモデルによる予測値と真値の関係を可視化する．
tibble(obs = y_tilde,
model1 = y_hat1, model2 = y_hat2, model3 = y_hat3) |>
pivot_longer(!obs) |>
ggplot(aes(x = value, y = obs)) +
geom_point(aes(colour = name)) +
geom_abline(slope = 1, intercept = 0, colour = "gray") +
labs(x = "fitted value", y = "observed value") +
theme(legend.title = element_blank())
正しいモデルでは正確に予測が行われているが，不要な項が含まれると予測誤差が大きくなることがわかる．
数値的な評価の例
雑音のないデータ \tilde{y} と予測値の相関係数(2乗すると決定係数に相当)を用いて数値的な評価をする．
cor(y_tilde, y_hat1)^2 [1] 1cor(y_tilde, y_hat2)^2[1] 0.5110198cor(y_tilde, y_hat3)^2[1] 0.0008350369