library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(broom) # 解析結果を tibble 形式に集約
library(gt) # 表の作成
library(gtsummary) # 分析結果の表の作成第4講 回帰分析
モデルの評価
準備
以下で利用する共通パッケージを読み込む．
推定量の性質を調べる数値実験
以下に人工データによる推定量の性質の確認の例を示す．
正規乱数を用いた線形単回帰モデル
set.seed(987) # 乱数のシード値を設定(再現性の担保のため)
x_obs <- tibble(x0 = 1, x1 = c(1,3,5,7)) # 説明変数の観測値
epsilon <- rnorm(nrow(x_obs), sd = 0.5) # 誤差項の生成
beta <- c(2, -3) # 真の回帰係数
toy_data <- x_obs |> # 目的変数の観測値を追加
mutate(y = as.vector(as.matrix(x_obs) %*% beta) + epsilon)
toy_formula <- y ~ x1
toy_lm <- lm(toy_formula, data = toy_data) # 回帰係数の推定
coef(toy_lm) # 回帰係数の取得(Intercept) x1
2.091690 -2.994888 summary(toy_lm) # 分析結果の概要の表示
Call:
lm(formula = toy_formula, data = toy_data)
Residuals:
1 2 3 4
-0.12531 0.02797 0.31999 -0.22265
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.09169 0.29738 7.034 0.019620 *
x1 -2.99489 0.06489 -46.151 0.000469 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2902 on 2 degrees of freedom
Multiple R-squared: 0.9991, Adjusted R-squared: 0.9986
F-statistic: 2130 on 1 and 2 DF, p-value: 0.0004692上記の例では行列とベクトルの積を用いているが， 以下のように説明変数・目的変数を個別に作成してもよい．
x_obs <- c(1,3,5,7) # 説明変数の観測値
epsilon <- rnorm(length(x_obs), sd = 0.5) # 誤差項の生成
y_obs <- beta[1] + beta[2] * x_obs + epsilon # 目的変数の観測値
toy_data <- tibble(x1 = x_obs, y = y_obs) # データフレームの作成(必要ならx0を付加しても良い)数値実験のためのコードは以下のようになる.
Monte-Carlo法の雛型
1回の試行に対応する関数を作成する．
mc_trial <- function(x) { # xは仮の引数で計算の中では使わない
## 乱数生成と推定の処理
return(tibbleクラスの返り値)} # 必要であれば関数as_tibble_row()を利用Monte-Carlo 実験を行うには関数purrr::map() を利用すればよい
mc_num <- 5000 # 実験回数を指定
mc_data <-
1:mc_num |> map(mc_trial) |> # 実験結果のリストが得られる
list_rbind() # データフレームに変換実験結果に対しては，下記の例のような統計処理・視覚化処理を行う．
mc_data |>
summarise(across(everything(), var)) # 各列の分散の計算
library(GGally)
ggpairs(mc_data) # 散布図行列の描画
tibble(x = mc_data[[k]]) |> # k列目で新しいデータフレームを作成
ggplot(aes(x = x)) + geom_histogram() # k列目のデータのヒストグラム試行の関数で仮の引数を使わない場合は関数 purrr::map() の中で無名関数を利用すれば良い．
mc_trial <- function() {
## 乱数生成と推定の処理
return(tibbleクラスの返り値)}
mc_data <-
1:mc_num |> map(\(x)mc_trial()) |> list_rbind() 関数 purrr::map() の代わりに関数 base::replicate() を利用することもできるが， 出力の整形に注意する必要がある．
mc_trial <- function() { # 1回の試行を行うプログラム
## 乱数生成と推定の処理
return(vectorクラスの返り値)}
mc_num <- 5000 # 実験回数を指定
mc_data <-
replicate(mc_num, mc_trial()) |> # Monte-Carlo実験
t() |> as_tibble() # 転置(関数t())してデータフレームに変換問題
最小二乗推定量の性質を数値実験 (Monte-Carlo法) により確認しなさい．
- 以下のモデルに従う人工データを生成する．
- 説明変数の観測データ : \begin{equation} \{1, 20, 13, 9, 5, 15, 19, 8, 3, 4\} \end{equation}
- 確率モデル : \begin{equation} y=-1+2\times x + \epsilon, \quad \epsilon\sim\mathcal{N}(0,2) \end{equation}
- 観測データから回帰係数を推定する．
- 実験を複数回繰り返し，推定値 (\hat\beta_{0},\hat\beta_{1}) の分布を調べる．
解答欄
解答例
試行を行う関数を作成する．
set.seed(2468) # 乱数のシード値 (適宜変更せよ)
x_obs <- tibble(x0 = 1, # 説明変数の観測値
x1 = c(1, 20, 13, 9, 5, 15, 19, 8, 3, 4))
beta <- c(beta0=-1, beta1=2) # 回帰係数 (切片,係数の真値)
sigma <- sqrt(2) # 誤差の標準偏差(分散の平方根)を設定
mc_trial <- function(x){
epsilon <- rnorm(nrow(x_obs), sd = sigma) # 誤差項の生成
toy_data <- x_obs |> # 目的変数の観測値を追加
mutate(y = as.vector(as.matrix(x_obs) %*% beta) + epsilon)
toy_lm <- lm(y ~ x1, data = toy_data) # 回帰係数の推定
return(as_tibble_row(coef(toy_lm))) # 推定された係数だけ返す
}少数回の実験で確認してみる．
mc_num <- 3 # 実験回数
1:mc_num |> map(mc_trial) # リストとして出力される[[1]]
# A tibble: 1 × 2
`(Intercept)` x1
<dbl> <dbl>
1 -1.33 1.95
[[2]]
# A tibble: 1 × 2
`(Intercept)` x1
<dbl> <dbl>
1 -0.871 1.91
[[3]]
# A tibble: 1 × 2
`(Intercept)` x1
<dbl> <dbl>
1 -1.45 2.12数値実験を行う．
mc_num <- 5000 # 実験回数
mc_data <-
1:mc_num |> map(mc_trial) |> # Monte-Carlo実験
list_rbind() |> # リストをtibble形式に変換
set_names(names(beta)) # 列名をbetaに合わせて変更
mc_data # 内容を確認# A tibble: 5,000 × 2
beta0 beta1
<dbl> <dbl>
1 -0.881 1.99
2 -1.73 2.06
3 -2.08 2.14
4 -0.873 1.94
5 -1.89 2.10
6 -2.46 2.14
7 -1.83 2.04
8 -0.00339 1.99
9 -1.58 2.08
10 -0.344 2.00
# ℹ 4,990 more rowsデータフレーム mc_data には mc_num 回の実験結果が保存されている．
回帰係数の分布(2次元)を調べてみる．
mc_data |>
ggplot(aes(x = beta0, y = beta1)) +
geom_point(colour = "blue", shape = 20) + # 推定値の散布図
geom_vline(xintercept = beta["beta0"], colour = "orchid") + # beta0の真値 (垂直線)
geom_hline(yintercept = beta["beta1"], colour = "orchid") # beta1の真値 (水平線)
軸名をギリシャ文字にしたい場合は以下を加えればよい(Figure 1)．
last_plot() + # 直前のプロットを指す
labs(x = expression(beta[0]), y = expression(beta[1]))
周辺分布のヒストグラムを追加する場合は関数 ggExtra::ggMarginal() が利用できる(Figure 2)．
library(ggExtra)
ggMarginal(last_plot(), type = "histogram")
各回帰係数の周辺分布は以下のようにしても描ける．
beta_cov <- sigma^2 * solve(crossprod(as.matrix(x_obs))) # 推定量の共分散行列
#' beta0 (k=1), beta1 (k=2)
for(k in 1:2) { # 同じ処理であればfor文などの利用を推奨
bar <- tibble(x = mc_data[[k]]) |>
ggplot(aes(x = x)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "blue") +
geom_vline(xintercept = beta[k], # 真の値
colour = "orchid") +
geom_function(fun = \(x) dnorm(x,
mean = beta[k],
sd = sqrt(beta_cov[k, k])),
colour = "orchid") + # 理論分布
labs(x = names(mc_data)[k])
print(bar) # for 文などの block 内でのグラフ表示は明示的に print する
}
あるいは変数名を評価して列を指定することもできる． 上記と結果は同じなので，描画はせずに以下にコードだけ記載する．
for(k in 1:2) {
foo <- sym(names(mc_data)[k]) # 文字列を扱うための手続き
bar <- mc_data |>
ggplot(aes(x = !!foo)) + # foo に入った文字列を評価する
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "blue") +
geom_vline(xintercept = beta[k], colour = "orchid") + # 真の値
geom_function(fun = \(x) dnorm(x,
mean = beta[k],
sd = sqrt(beta_cov[k, k])),
colour = "orchid") # 理論分布
print(bar)
}標準誤差の性質
問題
数値実験により標準誤差の性質を確認しなさい．
人工データを用いて標準誤差と真の誤差を比較する． 標準誤差は以下のようにして取り出せる．
toy_lm <- lm(toy_formula, toy_data) tidy(toy_lm) # 係数に関する情報をtibbleクラスとして取得 tidy(toy_lm)[2] # 列番号での要素の指定 tidy(toy_lm)["std.error"] # 列名での要素の指定広告費と売上データを用いて係数の精度を議論する．
東京の気候データを用いて係数の精度を議論する．
解答欄
解答例
人工データによる標準誤差と真の誤差の比較を行うために，以下のような重回帰モデル(適宜変更して構わない)での試行を設定する．
set.seed(1313) # 乱数のシード値 (適宜変更せよ)
x_obs <- tibble(x0 = 1,
x1 = c(1, 20, 13, 9, 5, 15, 19, 8, 3, 4), # 説明変数1
x2 = c(3, 19, 1, 4, 18, 7, 2, 10, 6, 12)) # 説明変数2
beta <- c(-1, 2, -3) # (切片，x1の係数, x2の係数)
sigma <- sqrt(2) # 誤差の標準偏差(分散の平方根)
mc_trial <- function(x){
epsilon <- rnorm(nrow(x_obs), sd = sigma) # 誤差項
toy_data <- x_obs |> # 目的変数の観測値を追加
mutate(y = as.vector(as.matrix(x_obs) %*% beta) + epsilon)
toy_lm <- lm(y ~ x1 + x2, data = toy_data) # 回帰係数の推定
return(tidy(toy_lm) |> pull(std.error) |> # 標準誤差を抽出
set_names(paste0("beta",0:2,".se")) |> # 名前を付けて返す
as_tibble_row()) # tibbleクラスに変更
}これを用いて数値実験を行う．
mc_num <- 5000 # 実験回数
mc_data <-
1:mc_num |> map(mc_trial) |> # mc_num回の試行
list_rbind() # データフレームの作成各回帰係数の標準誤差の分布を描くと以下のようになる．
beta_cov <- sigma^2*solve(crossprod(as.matrix(x_obs))) # 推定量の共分散行列
#' beta0 (k=1), beta1 (k=2), beta2 (k=3)
for(k in 1:3) {
bar <- tibble(x = mc_data[[k]]) |>
ggplot(aes(x = x)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "blue") +
geom_vline(xintercept = sqrt(beta_cov[k,k]), # 真の値
colour = "orchid") +
labs(x = names(mc_data)[k], title = "std. errors")
print(bar)
}
広告費と売上データによる分析は以下のように行うことができる． まず，データの読み込みを行う．
adv_data <-
read_csv('https://www.statlearning.com/s/Advertising.csv',
col_select = -1) # 1列目は行名と同じIDなので除去しておくモデルの推定を行う．
adv_lm1 <- lm(sales ~ TV, data = adv_data)
adv_lm2 <- lm(sales ~ radio, data = adv_data)
adv_lm3 <- lm(sales ~ TV + radio, data = adv_data)推定値とその標準誤差を比較する．
tidy(adv_lm1) |> select(1:3)# A tibble: 2 × 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 7.03 0.458
2 TV 0.0475 0.00269tidy(adv_lm2) |> select(1:3)# A tibble: 2 × 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 9.31 0.563
2 radio 0.202 0.0204tidy(adv_lm3) |> select(1:3)# A tibble: 3 × 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 2.92 0.294
2 TV 0.0458 0.00139
3 radio 0.188 0.00804関数 gtsummary::tbl_regression() を利用すると上記を表として簡単にまとめられる．
#' 規定の表
tbl_merge(
tbls = list(
tbl_regression(adv_lm1),
tbl_regression(adv_lm2),
tbl_regression(adv_lm3)),
tab_spanner = paste0("モデル",1:3))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| Beta | 95% CI | p-value | Beta | 95% CI | p-value | Beta | 95% CI | p-value | |
| TV | 0.05 | 0.04, 0.05 | <0.001 | 0.05 | 0.04, 0.05 | <0.001 | |||
| radio | 0.20 | 0.16, 0.24 | <0.001 | 0.19 | 0.17, 0.20 | <0.001 | |||
| Abbreviation: CI = Confidence Interval | |||||||||
規定の表では推定値・信頼区間・p値のみなので，必要な項目を並べるための関数を設定すると良い．
#' 表示を変更する例
my_gts <- function(x){
tbl_regression(x, intercept = TRUE) |> # 標準の表
modify_column_hide(columns = c(conf.low,p.value)) |> # 信頼区間・p値を非表示
modify_column_unhide(columns = std.error) } # 標準誤差を表示
tbl_merge(
tbls = list(
my_gts(adv_lm1),
my_gts(adv_lm2),
my_gts(adv_lm3)),
tab_spanner = paste0("モデル",1:3))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
|||
|---|---|---|---|---|---|---|
| Beta | SE | Beta | SE | Beta | SE | |
| (Intercept) | 7.0 | 0.458 | 9.3 | 0.563 | 2.9 | 0.294 |
| TV | 0.05 | 0.003 | 0.05 | 0.001 | ||
| radio | 0.20 | 0.020 | 0.19 | 0.008 | ||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||
東京の気候データによる分析は以下のとおりである．
まず，データを読み込み，9月のデータを抽出して整理する．
tw_subset <-
read_csv("data/tokyo_weather.csv") |>
filter(month == 9)回帰モデルの推定を行う．
tw_lm1 <- lm(temp ~ solar, data = tw_subset) # モデル1
tw_lm2 <- lm(temp ~ solar + press, data = tw_subset) # モデル2
tw_lm3 <- lm(temp ~ solar + press + cloud, data = tw_subset) # モデル3モデルの推定値とその標準誤差は以下のとおりである．
my_gts <- function(x){
tbl_regression(x, intercept = TRUE) |> # 標準の表
modify_column_hide(columns = c(conf.low,p.value)) |> # 信頼区間・p値を非表示
modify_column_unhide(columns = std.error) } # 標準誤差を表示
tbl_merge(
tbls = list(
my_gts(tw_lm1),
my_gts(tw_lm2),
my_gts(tw_lm3)),
tab_spanner = c("モデル1","モデル2","モデル3"))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
|||
|---|---|---|---|---|---|---|
| Beta | SE | Beta | SE | Beta | SE | |
| (Intercept) | 22 | 1.05 | 216 | 106 | 222 | 105 |
| solar | 0.35 | 0.070 | 0.33 | 0.067 | 0.22 | 0.116 |
| press | -0.19 | 0.104 | -0.19 | 0.104 | ||
| cloud | -0.29 | 0.239 | ||||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||
関数 base::summary() を用いる場合は以下を実行すればよい．
tidy(tw_lm1) |> select(term,estimate,std.error)
tidy(tw_lm2) |> select(1:3) # 名前ではなく列番号で指定する場合
tidy(tw_lm3) |> select(1:3) # cloud の標準誤差が大きく精度が悪いことが示唆されるt 統計量の性質
問題
数値実験により t 統計量の性質を確認しなさい．
人工データを用いて t 統計量の分布を確認する．
t 統計量とそのp値は以下のようにして取り出せる．
toy_lm <- lm(toy_formula, toy_data) tidy(toy_lm) |> select(term,statistic,p.value) # 列名での指定 tidy(toy_lm) |> select(1,4:5) # 列番号での指定 tidy(toy_lm) |> select(-c(2,3)) # 不要な列番号を除去広告費と売上データを用いて係数の有意性を議論する．
東京の気候データを用いて係数の有意性を議論する．
解答欄
解答例
人工データにより t 統計量の性質を調べるために，以下のような重回帰モデル(適宜変更して構わない)での試行を設定する．
set.seed(2525) # 乱数のシード値 (適宜変更せよ)
x_obs <- tibble(x0 = 1,
x1 = c(1, 20, 13, 9, 5, 15, 19, 8, 3, 4), # 説明変数1
x2 = c(3, 19, 1, 4, 18, 7, 2, 10, 6, 12)) # 説明変数2
beta <- c(-1, 2, 0) # (切片，x1の係数，x2の係数)
sigma <- sqrt(2) # 誤差の標準偏差(分散の平方根)
mc_trial <- function(x){
epsilon <- rnorm(nrow(x_obs), sd = sigma) # 誤差項
toy_data <- x_obs |> # 目的変数の観測値を追加
mutate(y = as.vector(as.matrix(x_obs) %*% beta) + epsilon)
toy_lm <- lm(y ~ x1 + x2, data = toy_data)
return(tidy(toy_lm) |> pull(statistic) |> # t統計量を返す
set_names(paste("beta",0:2,".tval")) |>
as_tibble_row())
}x1 の係数は 2 なので帰無仮説に従わない，x2 の係数は 0 なので帰無仮説に従うことに注意する．
数値実験を行う．
mc_num <- 5000
mc_data <- 1:mc_num |> map(mc_trial) |> list_rbind()各回帰係数の t 統計量の分布を図示すると以下のようになる．
n <- nrow(x_obs) # データ数 n
p <- 2 # 説明変数の次元
#' beta0 (k=1), beta1 (k=2), beta2 (k=3)
for(k in 1:3) { # 同じ処理であればfor文などの利用を推奨
bar <- tibble(x = mc_data[[k]]) |>
ggplot(aes(x = x)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "blue") +
geom_function(fun = \(x) dt(x, df = n-p-1), # 自由度 n-p-1 のt分布
colour = "orchid") +
labs(x = names(mc_data)[k])
print(bar) # for 文などの block 内でのグラフ表示は明示的に print する
}
広告費と売上データによる分析において，全てを用いたモデルと newspaper を除いたモデルを比較する．
summary(lm(sales ~ ., data = adv_data)) # "." は全て
Call:
lm(formula = sales ~ ., 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(lm(sales ~ . - newspaper, data = adv_data)) # "-" は除外
Call:
lm(formula = sales ~ . - newspaper, 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関数 gtsummary::tbl_regression() を用いて，これらの情報をまとめる．
my_gts <- function(x){
tbl_regression(x, intercept = TRUE) |>
modify_column_hide(columns = conf.low) |>
modify_column_unhide(columns = c(std.error,statistic)) |>
add_glance_table(include = c(r.squared,adj.r.squared))}
tbl_merge(
tbls = list(
my_gts(lm(sales ~ ., data = adv_data)),
my_gts(lm(sales ~ . - newspaper, data = adv_data))),
tab_spanner = c("モデル1","モデル2"))| Characteristic |
モデル1
|
モデル2
|
||||||
|---|---|---|---|---|---|---|---|---|
| Beta | SE | Statistic | p-value | Beta | SE | Statistic | p-value | |
| (Intercept) | 2.9 | 0.312 | 9.42 | <0.001 | 2.9 | 0.294 | 9.92 | <0.001 |
| TV | 0.05 | 0.001 | 32.8 | <0.001 | 0.05 | 0.001 | 32.9 | <0.001 |
| radio | 0.19 | 0.009 | 21.9 | <0.001 | 0.19 | 0.008 | 23.4 | <0.001 |
| newspaper | 0.00 | 0.006 | -0.177 | 0.9 | ||||
| R² | 0.897 | 0.897 | ||||||
| Adjusted R² | 0.896 | 0.896 | ||||||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||||
newspaperの係数のt統計量から有意性は低いと考えられる．また，自由度調整済決定係数も除いた方が若干高くなることが確認できる．
東京の気候データによる分析において，solarとpressを用いたモデルを比較する．
summary(lm(temp ~ solar + press, data = tw_subset))
Call:
lm(formula = temp ~ solar + press, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-4.2223 -0.9711 0.4653 1.6045 2.8526
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 215.70418 105.64144 2.042 0.0510 .
solar 0.33231 0.06729 4.939 3.59e-05 ***
press -0.19178 0.10443 -1.836 0.0773 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.061 on 27 degrees of freedom
Multiple R-squared: 0.5282, Adjusted R-squared: 0.4933
F-statistic: 15.12 on 2 and 27 DF, p-value: 3.937e-05summary(lm(temp ~ press, data = tw_subset))
Call:
lm(formula = temp ~ press, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-5.3175 -2.2855 0.8771 2.3231 4.1583
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 279.9332 142.0324 1.971 0.0587 .
press -0.2508 0.1406 -1.784 0.0853 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.792 on 28 degrees of freedom
Multiple R-squared: 0.1021, Adjusted R-squared: 0.07
F-statistic: 3.183 on 1 and 28 DF, p-value: 0.08527summary(lm(temp ~ solar, data = tw_subset))
Call:
lm(formula = temp ~ solar, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-5.4052 -0.9065 0.6154 1.4640 2.9743
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.72178 1.04638 20.759 < 2e-16 ***
solar 0.34644 0.06962 4.976 2.96e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.146 on 28 degrees of freedom
Multiple R-squared: 0.4693, Adjusted R-squared: 0.4504
F-statistic: 24.76 on 1 and 28 DF, p-value: 2.959e-05tbl_merge(
tbls = list(
my_gts(lm(temp ~ solar + press, data = tw_subset)),
my_gts(lm(temp ~ press, data = tw_subset)),
my_gts(lm(temp ~ solar, data = tw_subset))),
tab_spanner = c("モデル1","モデル2","モデル3"))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Beta | SE | Statistic | p-value | Beta | SE | Statistic | p-value | Beta | SE | Statistic | p-value | |
| (Intercept) | 216 | 106 | 2.04 | 0.051 | 280 | 142 | 1.97 | 0.059 | 22 | 1.05 | 20.8 | <0.001 |
| solar | 0.33 | 0.067 | 4.94 | <0.001 | 0.35 | 0.070 | 4.98 | <0.001 | ||||
| press | -0.19 | 0.104 | -1.84 | 0.077 | -0.25 | 0.141 | -1.78 | 0.085 | ||||
| R² | 0.528 | 0.102 | 0.469 | |||||||||
| Adjusted R² | 0.493 | 0.070 | 0.450 | |||||||||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||||||||
press単体では係数の推定精度も決定係数も低いが，solarと組み合わせることにより精度が上がり説明力も高くなることがわかる．
F 統計量の性質
問題
数値実験により F 統計量の性質を確認しなさい．
人工データを用いて F 統計量の分布を確認しなさい．
F 統計量とその自由度は以下のようにして取り出せる．
toy_lm <- lm(toy_formula, toy_data) glance(toy_lm)[c("statistic","p.value","df","df.residual")] # 以下全て等価 glance(toy_lm)[4:6,11] glance(toy_lm) |> select(statistic,p.value,df,df.residual) toy_lm |> glance() |> select(4:6,11)広告費と売上データのモデルの有効性を議論しなさい．
東京の気候データのモデルの有効性を議論しなさい．
解答欄
解答例
人工データにより F 統計量の性質を調べるために，前出の重回帰モデルを利用することにする．
set.seed(2525) # 乱数のシード (適宜変更せよ)
x_obs <- tibble(x0 = 1,
x1 = c(1, 20, 13, 9, 5, 15, 19, 8, 3, 4), # 説明変数1
x2 = c(3, 19, 1, 4, 18, 7, 2, 10, 6, 12)) # 説明変数2
beta <- c(-1, 0, 0) # (切片, x1の係数, x2の係数)
sigma <- sqrt(2) # 誤差の標準偏差(分散の平方根)
mc_trial <- function(x){
epsilon <- rnorm(nrow(x_obs), sd = sigma) # 誤差項
toy_data <- x_obs |> # 目的変数の観測値を追加
mutate(y = as.vector(as.matrix(x_obs) %*% beta) + epsilon)
toy_lm <- lm(y ~ x1 + x2, data = toy_data) # 回帰係数の推定
return(glance(toy_lm) |> pull(statistic) |> # F統計量を返す
as_tibble_col("fstat"))
}上記の設定では x1,x2 の係数はどちらも0なので帰無仮説が成り立つことに注意する． まず，帰無仮説が成り立つ場合の数値実験を行う．
mc_num <- 5000
mc_data <- 1:mc_num |> map(mc_trial) |> list_rbind()実験結果は1次元なので転置は不要であるが，ここでは列名の設定をしておく．
モデルの F 統計量の分布を図示すると以下のようになる．
n <- nrow(x_obs) # データ数
p <- 2 # 説明変数の次元
mc_data |>
ggplot(aes(x = fstat)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "blue") +
geom_function(fun = \(x) df(x, df1 = p, df2 = n-p-1), # 自由度p,n-p-1のF-分布
colour = "orchid") +
labs(x = "F statistic", title="null hypothesis is true")
回帰式の設定を変更して帰無仮説が成り立たない場合の数値実験を行う．
beta <- c(-1, 2, 0) # x1の係数 : 帰無仮説が成り立たない
mc_data <- 1:mc_num |> map(mc_trial) |> list_rbind()モデルの F 統計量の分布を図示すると帰無分布に従わないことがわかる．
mc_data |>
ggplot(aes(x = fstat)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "blue") +
geom_function(fun = \(x) df(x, df1 = p, df2 = n-p-1), # 自由度p,n-p-1のF-分布
colour = "orchid") +
labs(x = "F statistic", title="null hypothesis is false")
次に広告費と売上データを用いて，説明変数が1つのモデルの有効性を比較してみる．
summary(lm(sales ~ TV, data = adv_data))
Call:
lm(formula = sales ~ TV, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-8.3860 -1.9545 -0.1913 2.0671 7.2124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.032594 0.457843 15.36 <2e-16 ***
TV 0.047537 0.002691 17.67 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.259 on 198 degrees of freedom
Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099
F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16summary(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 ~ newspaper, data = adv_data))
Call:
lm(formula = sales ~ newspaper, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-11.2272 -3.3873 -0.8392 3.5059 12.7751
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.35141 0.62142 19.88 < 2e-16 ***
newspaper 0.05469 0.01658 3.30 0.00115 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.092 on 198 degrees of freedom
Multiple R-squared: 0.05212, Adjusted R-squared: 0.04733
F-statistic: 10.89 on 1 and 198 DF, p-value: 0.001148radio, newspaper は決定係数は小さく説明力は無いが，F 統計量はそれなりに小さいのでモデルの有効性が無いとは言えない．
関数 gtsummary::tbl_regression() を利用してまとめる場合は，例えば以下のようにすればよい．
my_gts <- function(x){ # 表示する項目を整理した関数を用意する
tbl_regression(x) |>
add_glance_table(include = c(adj.r.squared,statistic,p.value)) }
tbl_merge(
tbls = list(
my_gts(lm(sales ~ TV, data = adv_data)),
my_gts(lm(sales ~ radio, data = adv_data)),
my_gts(lm(sales ~ newspaper, data = adv_data))),
tab_spanner = paste0("モデル",1:3)) |>
modify_table_body(~.x |> # characteristic を並べ変える
arrange(factor(variable,
levels = c("TV","radio","newspaper"))))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| Beta | 95% CI | p-value | Beta | 95% CI | p-value | Beta | 95% CI | p-value | |
| TV | 0.05 | 0.04, 0.05 | <0.001 | ||||||
| radio | 0.20 | 0.16, 0.24 | <0.001 | ||||||
| newspaper | 0.05 | 0.02, 0.09 | 0.001 | ||||||
| Adjusted R² | 0.610 | 0.329 | 0.047 | ||||||
| Statistic | 312 | 98.4 | 10.9 | ||||||
| p-value | <0.001 | <0.001 | 0.001 | ||||||
| Abbreviation: CI = Confidence Interval | |||||||||
東京の気候データでは press, solar, rain によるモデルを検討してみる．
summary(lm(temp ~ press, data = tw_subset))
Call:
lm(formula = temp ~ press, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-5.3175 -2.2855 0.8771 2.3231 4.1583
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 279.9332 142.0324 1.971 0.0587 .
press -0.2508 0.1406 -1.784 0.0853 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.792 on 28 degrees of freedom
Multiple R-squared: 0.1021, Adjusted R-squared: 0.07
F-statistic: 3.183 on 1 and 28 DF, p-value: 0.08527summary(lm(temp ~ press + solar, data = tw_subset))
Call:
lm(formula = temp ~ press + solar, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-4.2223 -0.9711 0.4653 1.6045 2.8526
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 215.70418 105.64144 2.042 0.0510 .
press -0.19178 0.10443 -1.836 0.0773 .
solar 0.33231 0.06729 4.939 3.59e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.061 on 27 degrees of freedom
Multiple R-squared: 0.5282, Adjusted R-squared: 0.4933
F-statistic: 15.12 on 2 and 27 DF, p-value: 3.937e-05summary(lm(temp ~ press + solar + rain, data = tw_subset))
Call:
lm(formula = temp ~ press + solar + rain, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-4.2226 -0.9496 0.3838 1.6092 2.8614
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 211.106010 111.653201 1.891 0.0699 .
press -0.187276 0.110299 -1.698 0.1015
solar 0.334094 0.069502 4.807 5.59e-05 ***
rain 0.007116 0.046112 0.154 0.8785
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.099 on 26 degrees of freedom
Multiple R-squared: 0.5287, Adjusted R-squared: 0.4743
F-statistic: 9.721 on 3 and 26 DF, p-value: 0.0001777press のみではモデルの有効性があるとは言えないが， solar と組み合わせることにより有効性が確認できる． rain を加えても press の係数に変化は見られないが， solar の係数が変化し決定係数が大きくなることから solar と rain が相補的にモデルの精度を上げている可能性が示唆される．
表形式にまとめると以下のようになる．
my_gts <- function(x){ # 表示する項目を整理した関数を用意する
tbl_regression(x) |>
add_glance_table(include = c(adj.r.squared,statistic,p.value)) }
tbl_merge(
tbls = list(
my_gts(lm(temp ~ press + solar + rain, data = tw_subset)),
my_gts(lm(temp ~ press + solar, data = tw_subset)),
my_gts(lm(temp ~ press, data = tw_subset))),
tab_spanner = paste0("モデル",1:3))| Characteristic |
モデル1
|
モデル2
|
モデル3
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| Beta | 95% CI | p-value | Beta | 95% CI | p-value | Beta | 95% CI | p-value | |
| press | -0.19 | -0.41, 0.04 | 0.10 | -0.19 | -0.41, 0.02 | 0.077 | -0.25 | -0.54, 0.04 | 0.085 |
| solar | 0.33 | 0.19, 0.48 | <0.001 | 0.33 | 0.19, 0.47 | <0.001 | |||
| rain | 0.01 | -0.09, 0.10 | 0.9 | ||||||
| Adjusted R² | 0.474 | 0.493 | 0.070 | ||||||
| Statistic | 9.72 | 15.1 | 3.18 | ||||||
| p-value | <0.001 | <0.001 | 0.085 | ||||||
| Abbreviation: CI = Confidence Interval | |||||||||
参考 : 診断プロット
回帰モデルのあてはまりを視覚的に評価するものとして， 以下の関係を表示する関数が用意されている．
- Residuals vs Fitted: あてはめ値(予測値)と残差の関係
- Normal Q-Q: 残差の正規性の確認
- Scale-Location: あてはめ値と正規化した残差の関係
- Residuals vs Leverage: 正規化した残差とテコ比の関係
Tip
関数 ggfortify::autoplot() を利用すると， 上記の診断プロットを簡単に作成できる．
library(ggfortify)
#' 関数 stats::lm() による推定結果の診断プロット
tw_lm6 <- lm(temp ~ press + solar + rain, data = tw_subset)
autoplot(tw_lm6)
診断プロットは1から6まで用意されており 標準では 1,2,3,5 がまとめて表示される． 個別に表示する場合は autoplot(tw_lm6, which = 1) のように指定する． 詳細は ?ggfortify::autoplot.lm を参照せよ．