```{r}
library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(GGally)
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月のデータを予測する．
```{r}
#' データの整理
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) # 回帰係数の評価
glance(tw_lm) # モデルの評価
#' あてはめ値と予測値の計算
tw_train_fitted <- augment(tw_lm, newdata = tw_train) # あてはめ値
tw_test_fitted <- augment(tw_lm, newdata = tw_test) # 予測値
```# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 538. 66.5 8.09 4.25e-11
2 solar 0.377 0.0652 5.78 3.18e- 7
3 press -0.513 0.0654 -7.85 1.10e-10
# 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.702 0.692 2.81 68.2 5.77e-16 2 -148. 304. 312.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>グラフでの表示の例は以下のようになる．
```{r}
#' 予測結果を図示
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() を用いると月を文字列のラベルとすることもできる．
```{r}
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)) # 集合と比解答欄
```{r}
#' ここにコードを書く
```解答例
東京の気候データによる分析において，信頼区間と予測区間の計算を行う．
推定用データと予測用データを準備し，推定用データを用いて回帰式の推定を行う．
Code
```{r}
#| code-fold: true
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) # モデルの推定
```信頼区間の計算を行う．
Code
```{r}
#| code-fold: true
tw_train_conf <-
augment(tw_lm, newdata = tw_train, # 推定用データに
interval = "confidence") # あてはめ値と信頼区間を付加
tw_test_conf <-
augment(tw_lm, newdata = tw_test, # 予測用データ(新規データ)に
interval = "confidence") # あてはめ値と信頼区間を付加
```予測区間の計算を行う．
Code
```{r}
#| code-fold: true
tw_train_pred <-
augment(tw_lm, newdata = tw_train, # 推定用データに
interval = "prediction") # あてはめ値と予測区間を付加
tw_test_pred <-
augment(tw_lm, newdata = tw_test, # 予測用データ(新規データ)に
interval = "prediction") # あてはめ値と予測区間を付加
```8月のデータで推定したモデルで8月をあてはめた信頼区間を描画する．
Code
```{r}
#| code-fold: true
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月をあてはめた予測区間を描画する．
Code
```{r}
#| code-fold: true
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月をあてはめた信頼区間を描画する．
Code
```{r}
#| code-fold: true
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月をあてはめた予測区間を描画する．
Code
```{r}
#| code-fold: true
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カテゴリカル変数の取り扱いには若干の注意が必要である．
- 何も宣言しなくても通常は適切に対応してくれる．
- 陽に扱う場合は関数 base::factor()(順序を指定可能), base::as.factor()(辞書順), forcats::asfactor()(出現順) を利用する．
因子属性の与え方の例
```{r}
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())
toy_data3 <- toy_data2 |> # 順序付き(levels)の因子化
mutate(X3 = factor(X, levels=c("S","A","B","C","D")))
glimpse(toy_data3) # toy_data2とはfactorの順序が異なる
toy_data4 <- toy_data2 |>
mutate(Y2 = factor(Y > 60)) # 条件による因子化
glimpse(toy_data4) # 条件の真偽で2値に類別される
```Rows: 5
Columns: 3
$ X <chr> "A", "S", "A", "B", "D"
$ Y <dbl> 85, 100, 80, 70, 30
$ X2 <fct> A, S, A, B, D
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, D
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年分)を用いて，気温を回帰する以下のモデルを検討しなさい．
- 降水の有無を表すカテゴリカル変数を用いたモデル (雨が降ると気温が変化することを検証する)
- 上記に月をカテゴリカル変数として加えたモデル (月毎の気温の差を考慮する)
解答欄
```{r}
#' ここにコードを書く
```解答例
9月から11月のデータを用いた分析を行うために，データを整理する．
Code
```{r}
#| code-fold: true
tw_subset <- tw_data |> filter(month %in% 9:11)
```日射量，気圧，湿度を用いた回帰モデルの検討を行う．まず，解析の対象となる変数の散布図を描いておく．
Code
```{r}
#| code-fold: true
tw_subset |>
select(temp, solar, press, humid) |>
ggpairs(columnLabels = c("気温","日射量","気圧","湿度"))
```
モデルを比較する．
Code
```{r}
#| code-fold: true
#' 基本となる線形回帰モデル
summary(lm(temp ~ solar + press + humid, data = tw_subset))
#' 湿度の対数を考えた線形回帰モデル
summary(lm(temp ~ solar + press + log(humid), data = tw_subset))
#' 最初のモデルにそれぞれの交互作用を加えたモデル (書き方はいろいろある)
summary(lm(temp ~ (solar + press + humid)^2, data = tw_subset))
#' 更に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.3585 -2.4679 -0.0807 2.8749 8.7324
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 235.32751 80.46439 2.925 0.00440 **
solar 0.74263 0.09056 8.200 1.94e-12 ***
press -0.23877 0.07796 -3.063 0.00292 **
humid 0.24709 0.03828 6.454 5.95e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.979 on 87 degrees of freedom
Multiple R-squared: 0.6048, Adjusted R-squared: 0.5912
F-statistic: 44.39 on 3 and 87 DF, p-value: < 2.2e-16
Call:
lm(formula = temp ~ solar + press + log(humid), data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.5919 -2.4348 -0.0931 2.9843 8.6557
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 188.03020 84.03399 2.238 0.02781 *
solar 0.70933 0.08885 7.984 5.34e-12 ***
press -0.24583 0.07780 -3.160 0.00217 **
log(humid) 17.03985 2.66049 6.405 7.42e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.989 on 87 degrees of freedom
Multiple R-squared: 0.6029, Adjusted R-squared: 0.5892
F-statistic: 44.02 on 3 and 87 DF, p-value: < 2.2e-16
Call:
lm(formula = temp ~ (solar + press + humid)^2, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.236 -2.528 -0.258 2.884 8.785
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.710e+02 6.057e+02 -0.613 0.542
solar 1.600e+01 1.727e+01 0.927 0.357
press 3.629e-01 5.967e-01 0.608 0.545
humid 6.392e+00 6.969e+00 0.917 0.362
solar:press -1.547e-02 1.691e-02 -0.914 0.363
solar:humid 5.139e-03 8.549e-03 0.601 0.549
press:humid -6.107e-03 6.868e-03 -0.889 0.376
Residual standard error: 4.008 on 84 degrees of freedom
Multiple R-squared: 0.613, Adjusted R-squared: 0.5854
F-statistic: 22.18 on 6 and 84 DF, p-value: 1.769e-15
Call:
lm(formula = temp ~ solar * press * humid, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-9.362 -2.231 0.115 2.455 8.913
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.909e+03 1.059e+03 -1.802 0.0752 .
solar 2.052e+02 1.089e+02 1.884 0.0630 .
press 1.882e+00 1.045e+00 1.800 0.0755 .
humid 2.593e+01 1.306e+01 1.985 0.0505 .
solar:press -2.023e-01 1.075e-01 -1.882 0.0634 .
solar:humid -2.542e+00 1.448e+00 -1.756 0.0828 .
press:humid -2.540e-02 1.290e-02 -1.970 0.0522 .
solar:press:humid 2.516e-03 1.430e-03 1.759 0.0822 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.959 on 83 degrees of freedom
Multiple R-squared: 0.6269, Adjusted R-squared: 0.5954
F-statistic: 19.92 on 7 and 83 DF, p-value: 2.028e-15基本的なモデルと最後のモデルの視覚的な評価 (診断プロット) を行う．
Code
```{r}
#| code-fold: true
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() があるが，前者は文字列を辞書順に並べ変えてから，後者はデータフレームでの出現順に因子化する．回帰の計算においては特に問題ないが，因子化の順序はグラフなど凡例の順番に影響する．
Code
```{r}
#| code-fold: true
tw_data_fact <- tw_data |>
mutate(rain = rain > 0, # 降雨の有無
month = as_factor(month)) # 月
```通年でのモデルを検討する．
Code
```{r}
#| code-fold: true
summary(lm(temp ~ rain, data = tw_data_fact))
```
Call:
lm(formula = temp ~ rain, data = tw_data_fact)
Residuals:
Min 1Q Median 3Q Max
-15.139 -8.104 -0.139 8.761 14.995
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.3045 0.5287 32.733 <2e-16 ***
rainTRUE 0.2345 1.0100 0.232 0.817
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.606 on 363 degrees of freedom
Multiple R-squared: 0.0001484, Adjusted R-squared: -0.002606
F-statistic: 0.05389 on 1 and 363 DF, p-value: 0.8166通年では雨と気温の関係は積極的に支持されない．
月を説明変数として加えると月毎の気温の差を考慮した回帰式が推定される．
Code
```{r}
#| code-fold: true
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
-9.1832 -1.6831 0.0523 1.4985 10.5365
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.8477 0.4753 14.406 < 2e-16 ***
rainTRUE -2.1266 0.3232 -6.581 1.70e-10 ***
month2 -0.2646 0.6888 -0.384 0.70105
month3 4.6423 0.6759 6.868 2.96e-11 ***
month4 9.5620 0.6819 14.022 < 2e-16 ***
month5 13.2804 0.6805 19.516 < 2e-16 ***
month6 18.7096 0.6834 27.377 < 2e-16 ***
month7 21.9961 0.6715 32.757 < 2e-16 ***
month8 23.2647 0.6721 34.617 < 2e-16 ***
month9 20.4120 0.6819 29.933 < 2e-16 ***
month10 12.4980 0.6773 18.453 < 2e-16 ***
month11 6.3367 0.6767 9.364 < 2e-16 ***
month12 2.0486 0.6721 3.048 0.00248 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.641 on 352 degrees of freedom
Multiple R-squared: 0.9087, Adjusted R-squared: 0.9056
F-statistic: 292 on 12 and 352 DF, p-value: < 2.2e-16月毎に比較した結果，雨の日の方が気温が低いことが示唆される．
診断プロットは以下のようになる．
Code
```{r}
#| code-fold: true
autoplot(lm(temp ~ rain + month, data = tw_data_fact))
```
参考 : その他の関数・パッケージ
モデルの探索の例
変数が増えるとモデルの比較が困難となるが， 関数 stats::step() で自動化することができる．
```{r}
adv_data <- read_csv('https://www.statlearning.com/s/Advertising.csv')
summary(lm(sales ~ radio, data = adv_data))
summary(lm(sales ~ TV + radio, data = adv_data))
summary(lm(sales ~ TV + radio + newspaper, data = adv_data))
summary(adv_init <- lm(sales ~ TV * radio * newspaper, data = adv_data))
adv_opt <- step(adv_init) # 最大のモデルから削減増加による探索
summary(adv_opt) # 探索された(準)最適なモデルの確認
```
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-16
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-16
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-16
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-16
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.234
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 による予測の例
人工データの生成例
以下はあくまで例なので，これを参考に自由に数値実験を設計して下さい．
```{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)
```モデルの推定と評価の例
```{r}
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)
summary(toy_lm2)
summary(toy_lm3)
```
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-12
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-11
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表としてまとめるには，例えば以下のようにすればよい．
```{r}
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 | |||
予測の例
新規データを作成し，予測の仕組みを確認する．
```{r}
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による予測値
```可視化の例
散布図を用いて，それぞれのモデルによる予測値と真値の関係を可視化する．
```{r}
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乗すると決定係数に相当)を用いて数値的な評価をする．
```{r}
cor(y_tilde, y_hat1)^2
cor(y_tilde, y_hat2)^2
cor(y_tilde, y_hat3)^2
```[1] 1
[1] 0.4615974
[1] 0.009544743