```{r}
library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(broom) # 解析結果を tibble 形式に集約
library(gt) # 表の作成
library(gtsummary) # 分析結果の表の作成
#' 日本語表示の設定．最後の例で利用する
if (Sys.info()["sysname"] == "Darwin") { # macOSのみ設定
jp_font <- "HiraMaruProN-W4"
theme_update(text = element_text(family = jp_font))
} else {
jp_font <- NULL
}
```第4講 回帰分析
モデルの評価
準備
以下で利用する共通パッケージを読み込む．
推定量の性質
正規乱数を用いた線形単回帰モデル
以下に人工データによる推定量の性質の確認の例を示す．
```{r}
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) # 回帰係数の取得
summary(toy_lm) # 分析結果の概要の表示
```(Intercept) x1
2.091690 -2.994888
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上記の例では行列とベクトルの積を用いているが，以下のように説明変数・目的変数を個別に作成してもよい．
```{r}
#| eval: false
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}\)) の分布を調べる．
解答欄
```{r}
#' ここにコードを書く
```解答例
試行を行う関数を作成する．
Code
```{r}
#| code-fold: true
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))) # 推定された係数だけ返す
}
```少数回の実験で確認してみる．
Code
```{r}
#| code-fold: true
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数値実験を行う．
Code
```{r}
#| code-fold: true
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次元)を調べてみる．
Code
```{r}
#| code-fold: true
#| fig-width: 6
#| fig-height: 6
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)．
```{r}
#| label: fig-greek
#| fig-cap: "ギリシャ文字を利用するの例"
#| fig-width: 6
#| fig-height: 6
last_plot() + # 直前のプロットを指す
labs(x = expression(beta[0]), y = expression(beta[1]))
```
周辺分布のヒストグラムを追加する場合は関数 ggExtra::ggMarginal() が利用できる(Figure 2)．
```{r}
#| label: fig-marginal
#| fig-cap: "周辺分布を付加するの例"
#| fig-width: 6
#| fig-height: 6
library(ggExtra)
ggMarginal(last_plot(), type = "histogram")
```
各回帰係数の周辺分布は以下のようにしても描ける．
```{r}
#| label: fig-each-marginal
#| fig-cap: "周辺分布を個別に描画する例"
#| fig-show: hold
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 する
}
```
あるいは変数名を評価して列を指定することもできる． 上記と結果は同じなので，描画はせずに以下にコードだけ記載する．
```{r}
#| eval: false
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)
}
```標準誤差の性質
問題
数値実験により標準誤差の性質を確認しなさい．
人工データを用いて標準誤差と真の誤差を比較する．標準誤差は以下のようにして取り出せる．
```{r} #| ecex: no toy_lm <- lm(toy_formula, toy_data) tidy(toy_lm) # 係数に関する情報をtibbleクラスとして取得 tidy(toy_lm)[2] # 列番号での要素の指定 tidy(toy_lm)["std.error"] # 列名での要素の指定 ```# A tibble: 2 × 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 2.09 0.297 7.03 0.0196 2 x1 -2.99 0.0649 -46.2 0.000469 # A tibble: 2 × 1 estimate <dbl> 1 2.09 2 -2.99 # A tibble: 2 × 1 std.error <dbl> 1 0.297 2 0.0649広告費と売上データを用いて係数の精度を議論する．
東京の気候データを用いて係数の精度を議論する．
解答欄
```{r}
#' ここにコードを書く
```解答例
人工データによる標準誤差と真の誤差の比較を行うために，以下のような重回帰モデル(適宜変更して構わない)での試行を設定する．
Code
```{r}
#| code-fold: true
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クラスに変更
}
```これを用いて数値実験を行う．
Code
```{r}
#| code-fold: true
mc_num <- 5000 # 実験回数
mc_data <-
1:mc_num |> map(mc_trial) |> # mc_num回の試行
list_rbind() # データフレームの作成
```各回帰係数の標準誤差の分布を描くと以下のようになる．
Code
```{r}
#| code-fold: true
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)
}
```
広告費と売上データによる分析は以下のように行うことができる．まず，データの読み込みを行う．
Code
```{r}
#| code-fold: true
adv_remote <- "https://www.statlearning.com/s/Advertising.csv"
adv_file <- "data/advertising.csv"
if(file.exists(adv_file)) { # 保存してある場合
adv_data <- read_csv(file = adv_file)
} else { # 保存していない場合
adv_data <- read_csv(file = adv_remote)
write_csv(adv_data, file = adv_file) # 保存する
}
```モデルの推定を行う．
Code
```{r}
#| code-fold: true
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)
```推定値とその標準誤差を比較する．
Code
```{r}
#| code-fold: true
tidy(adv_lm1) |> select(1:3)
tidy(adv_lm2) |> select(1:3)
tidy(adv_lm3) |> 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.00269
# A tibble: 2 × 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 9.31 0.563
2 radio 0.202 0.0204
# 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() を利用すると上記を表として簡単にまとめられる．
```{r}
#| column: margin
#' 規定の表
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値のみなので，必要な項目を並べるための関数を設定すると良い．
```{r}
#| column: margin
#' 表示を変更する例
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月のデータを抽出して整理する．
Code
```{r}
#| code-fold: true
tw_subset <-
read_csv("data/tokyo_weather.csv") |>
filter(month == 9)
```回帰モデルの推定を行う．
Code
```{r}
#| code-fold: true
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
```モデルの推定値とその標準誤差は以下のとおりである．
Code
```{r}
#| code-fold: true
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) | 24 | 1.33 | 408 | 76.1 | 398 | 75.8 |
| solar | 0.14 | 0.087 | 0.14 | 0.064 | 0.04 | 0.100 |
| press | -0.38 | 0.075 | -0.37 | 0.075 | ||
| cloud | -0.26 | 0.209 | ||||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||
関数 broom::tidy() を用いる場合は以下を実行すればよい．
```{r}
#| eval: false
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値は以下のようにして取り出せる．
```{r} #| eval: false 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)) # 不要な列番号を除去 ```広告費と売上データを用いて係数の有意性を議論する．
東京の気候データを用いて係数の有意性を議論する．
解答欄
```{r}
#' ここにコードを書く
```解答例
人工データにより \(t\) 統計量の性質を調べるために，以下のような重回帰モデル(適宜変更して構わない)での試行を設定する．
Code
```{r}
#| code-fold: true
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 なので帰無仮説に従うことに注意する．
数値実験を行う．
Code
```{r}
#| code-fold: true
mc_num <- 5000
mc_data <- 1:mc_num |> map(mc_trial) |> list_rbind()
```各回帰係数の \(t\) 統計量の分布を図示すると以下のようになる．
Code
```{r}
#| code-fold: true
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 を除いたモデルを比較する．
Code
```{r}
#| code-fold: true
summary(lm(sales ~ ., data = adv_data)) # "." は全て
summary(lm(sales ~ . - newspaper, data = adv_data)) # "-" は除外
```
Call:
lm(formula = sales ~ ., data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-8.8105 -0.9008 0.2641 1.1783 2.8336
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.0052094 0.3942082 7.623 1.06e-12 ***
...1 -0.0005798 0.0020992 -0.276 0.783
TV 0.0457759 0.0013988 32.725 < 2e-16 ***
radio 0.1883832 0.0086480 21.784 < 2e-16 ***
newspaper -0.0012433 0.0059319 -0.210 0.834
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.689 on 195 degrees of freedom
Multiple R-squared: 0.8973, Adjusted R-squared: 0.8951
F-statistic: 425.7 on 4 and 195 DF, p-value: < 2.2e-16
Call:
lm(formula = sales ~ . - newspaper, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-8.7767 -0.8931 0.2663 1.1968 2.8375
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.9779046 0.3711527 8.023 9.22e-14 ***
...1 -0.0005246 0.0020775 -0.252 0.801
TV 0.0457633 0.0013941 32.827 < 2e-16 ***
radio 0.1877652 0.0081100 23.152 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.685 on 196 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8957
F-statistic: 570.4 on 3 and 196 DF, p-value: < 2.2e-16関数 gtsummary::tbl_regression() を用いて，これらの情報をまとめる．
Code
```{r}
#| code-fold: true
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) | 3.0 | 0.394 | 7.62 | <0.001 | 3.0 | 0.371 | 8.02 | <0.001 |
| ...1 | 0.00 | 0.002 | -0.276 | 0.8 | 0.00 | 0.002 | -0.252 | 0.8 |
| TV | 0.05 | 0.001 | 32.7 | <0.001 | 0.05 | 0.001 | 32.8 | <0.001 |
| radio | 0.19 | 0.009 | 21.8 | <0.001 | 0.19 | 0.008 | 23.2 | <0.001 |
| newspaper | 0.00 | 0.006 | -0.210 | 0.8 | ||||
| R² | 0.897 | 0.897 | ||||||
| Adjusted R² | 0.895 | 0.896 | ||||||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||||
newspaper の係数の\(t\)統計量から有意性は低いと考えられる．また，自由度調整済決定係数も除いた方が若干高くなることが確認できる．
東京の気候データによる分析において， solar と press を用いたモデルを比較する．
Code
```{r}
#| code-fold: true
summary(lm(temp ~ solar + press, data = tw_subset))
summary(lm(temp ~ press, data = tw_subset))
summary(lm(temp ~ solar, data = tw_subset))
```
Call:
lm(formula = temp ~ solar + press, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-2.850 -1.674 -0.122 1.589 2.789
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 408.15579 76.13755 5.361 1.16e-05 ***
solar 0.14137 0.06357 2.224 0.0347 *
press -0.37975 0.07535 -5.040 2.73e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.874 on 27 degrees of freedom
Multiple R-squared: 0.5294, Adjusted R-squared: 0.4946
F-statistic: 15.19 on 2 and 27 DF, p-value: 3.805e-05
Call:
lm(formula = temp ~ press, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-4.179 -1.055 0.017 1.514 3.309
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 410.40575 81.31819 5.047 2.44e-05 ***
press -0.37997 0.08048 -4.721 5.94e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.002 on 28 degrees of freedom
Multiple R-squared: 0.4432, Adjusted R-squared: 0.4234
F-statistic: 22.29 on 1 and 28 DF, p-value: 5.941e-05
Call:
lm(formula = temp ~ solar, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-5.0223 -1.7788 0.3453 1.9400 3.4880
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.44529 1.33286 18.34 <2e-16 ***
solar 0.14179 0.08696 1.63 0.114
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.564 on 28 degrees of freedom
Multiple R-squared: 0.08671, Adjusted R-squared: 0.05409
F-statistic: 2.658 on 1 and 28 DF, p-value: 0.1142Code
```{r}
#| code-fold: true
tbl_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) | 408 | 76.1 | 5.36 | <0.001 | 410 | 81.3 | 5.05 | <0.001 | 24 | 1.33 | 18.3 | <0.001 |
| solar | 0.14 | 0.064 | 2.22 | 0.035 | 0.14 | 0.087 | 1.63 | 0.11 | ||||
| press | -0.38 | 0.075 | -5.04 | <0.001 | -0.38 | 0.080 | -4.72 | <0.001 | ||||
| R² | 0.529 | 0.443 | 0.087 | |||||||||
| Adjusted R² | 0.495 | 0.423 | 0.054 | |||||||||
| Abbreviations: CI = Confidence Interval, SE = Standard Error | ||||||||||||
press 単体では係数の推定精度も決定係数も低いが， solar と組み合わせることにより精度が上がり説明力も高くなることがわかる．
\(F\) 統計量の性質
問題
数値実験により \(F\) 統計量の性質を確認しなさい．
人工データを用いて \(F\) 統計量の分布を確認しなさい．
\(F\) 統計量とその自由度は以下のようにして取り出せる．
```{r} #| eval: false 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) ```広告費と売上データのモデルの有効性を議論しなさい．
東京の気候データのモデルの有効性を議論しなさい．
解答欄
```{r}
#' ここにコードを書く
```解答例
人工データにより \(F\) 統計量の性質を調べるために，前出の重回帰モデルを利用することにする．
Code
```{r}
#| code-fold: true
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なので帰無仮説が成り立つことに注意する．まず，帰無仮説が成り立つ場合の数値実験を行う．
Code
```{r}
#| code-fold: true
mc_num <- 5000
mc_data <- 1:mc_num |> map(mc_trial) |> list_rbind()
```実験結果は1次元なので転置は不要であるが，ここでは列名の設定をしておく．
モデルの \(F\) 統計量の分布を図示すると以下のようになる．
Code
```{r}
#| code-fold: true
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")
```
回帰式の設定を変更して帰無仮説が成り立たない場合の数値実験を行う．
Code
```{r}
#| code-fold: true
beta <- c(-1, 2, 0) # x1の係数 : 帰無仮説が成り立たない
mc_data <- 1:mc_num |> map(mc_trial) |> list_rbind()
```モデルの \(F\) 統計量の分布を図示すると帰無分布に従わないことがわかる．
Code
```{r}
#| code-fold: true
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")
```
原点付近を拡大すると，帰無分布の密度が高いところには統計量は観測されないことがわかる．
Code
```{r}
#| code-fold: true
last_plot() + xlim(0,40) # 横軸の表示範囲を限定する
```
次に広告費と売上データを用いて，説明変数が1つのモデルの有効性を比較してみる．
Code
```{r}
#| code-fold: true
summary(lm(sales ~ TV, data = adv_data))
summary(lm(sales ~ radio, data = adv_data))
summary(lm(sales ~ newspaper, 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-16
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 ~ 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() を利用してまとめる場合は，例えば以下のようにすればよい．
Code
```{r}
#| code-fold: true
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( # characteristic を並べ変える
~.x |> 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 によるモデルを検討してみる．
Code
```{r}
#| code-fold: true
summary(lm(temp ~ press, data = tw_subset))
summary(lm(temp ~ press + solar, data = tw_subset))
summary(lm(temp ~ press + solar + rain, data = tw_subset))
```
Call:
lm(formula = temp ~ press, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-4.179 -1.055 0.017 1.514 3.309
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 410.40575 81.31819 5.047 2.44e-05 ***
press -0.37997 0.08048 -4.721 5.94e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.002 on 28 degrees of freedom
Multiple R-squared: 0.4432, Adjusted R-squared: 0.4234
F-statistic: 22.29 on 1 and 28 DF, p-value: 5.941e-05
Call:
lm(formula = temp ~ press + solar, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-2.850 -1.674 -0.122 1.589 2.789
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 408.15579 76.13755 5.361 1.16e-05 ***
press -0.37975 0.07535 -5.040 2.73e-05 ***
solar 0.14137 0.06357 2.224 0.0347 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.874 on 27 degrees of freedom
Multiple R-squared: 0.5294, Adjusted R-squared: 0.4946
F-statistic: 15.19 on 2 and 27 DF, p-value: 3.805e-05
Call:
lm(formula = temp ~ press + solar + rain, data = tw_subset)
Residuals:
Min 1Q Median 3Q Max
-3.3580 -1.2159 -0.2178 1.6558 2.4655
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 427.84269 75.59858 5.659 5.95e-06 ***
press -0.39827 0.07470 -5.331 1.40e-05 ***
solar 0.08765 0.07177 1.221 0.233
rain -0.03049 0.02037 -1.497 0.146
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.832 on 26 degrees of freedom
Multiple R-squared: 0.5668, Adjusted R-squared: 0.5168
F-statistic: 11.34 on 3 and 26 DF, p-value: 6.124e-05press のみではモデルの有効性があるとは言えないが， solar と組み合わせることにより有効性が確認できる． rain を加えても press の係数に変化は見られないが， solar の係数が変化し決定係数が大きくなることから solar と rain が相補的にモデルの精度を上げている可能性が示唆される．
表形式にまとめると以下のようになる．
Code
```{r}
#| code-fold: true
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.40 | -0.55, -0.24 | <0.001 | -0.38 | -0.53, -0.23 | <0.001 | -0.38 | -0.54, -0.22 | <0.001 |
| solar | 0.09 | -0.06, 0.24 | 0.2 | 0.14 | 0.01, 0.27 | 0.035 | |||
| rain | -0.03 | -0.07, 0.01 | 0.15 | ||||||
| Adjusted R² | 0.517 | 0.495 | 0.423 | ||||||
| Statistic | 11.3 | 15.2 | 22.3 | ||||||
| p-value | <0.001 | <0.001 | <0.001 | ||||||
| Abbreviation: CI = Confidence Interval | |||||||||
補遺 : 診断プロット
回帰モデルのあてはまりを視覚的に評価するものとして，以下の関係を表示する関数が用意されている．
- Residuals vs Fitted: あてはめ値(予測値)と残差の関係
- Normal Q-Q: 残差の正規性の確認
- Scale-Location: あてはめ値と正規化した残差の関係
- Residuals vs Leverage: 正規化した残差とテコ比の関係
関数 ggfortify::autoplot() を利用すると，上記の診断プロットを簡単に作成できる．
```{r}
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 を参照せよ．