統計データ解析

第4講サンプルコード

Published

October 24, 2024

準備

以下で利用する共通パッケージを読み込む．

library(conflicted) # 関数名の衝突を警告
conflicts_prefer(   # 優先的に使う関数を指定
  dplyr::filter(),
  dplyr::select(),
  dplyr::lag(),
  )
library(tidyverse)  
library(broom)      # 解析結果を tibble 形式に集約
library(gt)         # 表の作成
library(gtsummary)  # 分析結果の表の作成

推定量の性質を調べる数値実験

問題

最小二乗推定量の性質を数値実験 (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 <- set_names(c(-1, 2), # 回帰係数 (切片,係数の真値)
                  c("beta0","beta1"))
sigma <-  sqrt(2) # 誤差の標準偏差(分散の平方根)を設定
mc_trial <- function(){ 
    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(coef(toy_lm)) # 推定された係数だけ返す
}

少数回の実験で確認してみる．同じ実験を複数回行うには関数 base::replicate() を利用するとよい．

mc_num <- 5 # 実験回数
replicate(mc_num, mc_trial())

                 [,1]       [,2]      [,3]       [,4]      [,5]
(Intercept) -1.326977 -0.8713116 -1.454163 -0.8811262 -1.726232
x1           1.945990  1.9140482  2.116954  1.9916924  2.057346

数値実験を行う．

mc_num <- 5000 # 実験回数
mc_data <- replicate(mc_num, mc_trial()) |> # mc_num回試行を行う
    t() |> as_tibble()                      # 得られる結果を転置してデータフレームにする
names(mc_data) <- names(beta) # 列名を変更(上書き)
mc_data # 内容を確認

# A tibble: 5,000 × 2
      beta0 beta1
      <dbl> <dbl>
 1 -2.08     2.14
 2 -0.873    1.94
 3 -1.89     2.10
 4 -2.46     2.14
 5 -1.83     2.04
 6 -0.00339  1.99
 7 -1.58     2.08
 8 -0.344    2.00
 9 -0.567    1.97
10 -0.745    1.95
# ℹ 4,990 more rows

データフレーム mc_data には mc_mun 回の実験結果が保存されている．

回帰係数の分布(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の真値 (水平線)

Tip

軸名をギリシャ文字にしたい場合は以下を加えればよい(Figure 1)．

last_plot() + # 直前のプロットを指す
    labs(x = expression(beta[0]), y = expression(beta[1]))

Tip

周辺分布のヒストグラムを追加する場合は関数 ggExtra::ggMarginal() が利用できる(Figure 2)．

ggExtra::ggMarginal(last_plot(), type = "histogram")

Tip

各回帰係数の周辺分布は以下のようにしても描ける．

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 する
}

Tip

あるいは変数名を評価して列を指定することもできる．¹

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(formula, toy_data)
summary(toy_lm)$coefficients        # 係数に関する情報はリストの要素として保管されている
summary(toy_lm)$coefficients[,2]    # 列番号での指定
summary(toy_lm)$coef[,"Std. Error"] # 列名での指定．coef と省略してもよい
tidy(toy_lm) # 関数 broom::tidy() でも同様に取得できる

広告費と売上データを用いて係数の精度を議論する．
東京の気候データを用いて係数の精度を議論する．

解答例

人工データによる標準誤差と真の誤差の比較を行うために，以下のような重回帰モデル(適宜変更して構わない)での試行を設定する．

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(){ 
    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(set_names(summary(toy_lm)$coef[,"Std. Error"],  # 標準誤差を
                     c("beta0.se","beta1.se","beta2.se"))) # 名前を付けて返す
}

これを用いて数値実験を行う．

mc_num <- 5000 # 実験回数
mc_data <- 
    replicate(mc_num, mc_trial()) |> # mc_num回の試行
    t() |> as_tibble()               # データフレームの作成

各回帰係数の標準誤差の分布を描くと以下のようになる．

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)

推定値とその標準誤差を比較する．

summary(adv_lm1)$coef[,1:2]

              Estimate  Std. Error
(Intercept) 7.03259355 0.457842940
TV          0.04753664 0.002690607

summary(adv_lm2)$coef[,1:2]

             Estimate Std. Error
(Intercept) 9.3116381 0.56290050
radio       0.2024958 0.02041131

summary(adv_lm3)$coef[,1:2]

              Estimate  Std. Error
(Intercept) 2.92109991 0.294489678
TV          0.04575482 0.001390356
radio       0.18799423 0.008039973

Tip

関数 gtsummary::tbl_regression() を利用すると上記を表として簡単にまとめられる．

#' 規定の表
tbl_merge( 
    tbls = list(
        tbl_regression(adv_lm1),
        tbl_regression(adv_lm2),
        tbl_regression(adv_lm3)),
    tab_spanner = c("モデル1","モデル2","モデル3"))

規定の表では推定値・信頼区間・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 = c("モデル1","モデル2","モデル3"))

Characteristic	モデル1			モデル2			モデル3
Characteristic	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
¹ CI = Confidence Interval

Characteristic	モデル1		モデル2		モデル3
Characteristic	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
¹ SE = Standard Error

東京の気候データによる分析は以下のとおりである．

まず，データを読み込み，9月のデータを抽出して整理する．

tw_subset <-
    read_csv("data/tokyo_weather.csv") |>
    filter(month == 9)

回帰モデルの設定をする．

tw_model1 <- temp ~ solar
tw_model2 <- temp ~ solar + press
tw_model3 <- temp ~ solar + press + cloud

回帰モデルの推定を行う．

tw_lm1 <- lm(tw_model1, data = tw_subset)
tw_lm2 <- lm(tw_model2, data = tw_subset)
tw_lm3 <- lm(tw_model3, data = tw_subset)

モデルの推定値とその標準誤差は以下のとおりである．

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
Characteristic	Beta	SE ¹	Beta	SE ¹	Beta	SE ¹
(Intercept)	23	0.855	386	91.0	384	92.8
solar	0.25	0.057	0.30	0.048	0.32	0.069
press			-0.36	0.090	-0.36	0.092
cloud					0.05	0.151
¹ SE = Standard Error

関数 base::summary() を用いる場合は以下を実行すればよい．

summary(tw_lm1)$coef[,c("Estimate","Std. Error")]
summary(tw_lm2)$coef[,1:2] # 名前ではなく列番号で指定する場合
summary(tw_lm3)$coef[,1:2] # cloud の標準誤差が大きく精度が悪いことが示唆される

t 統計量の性質

問題

数値実験により t 統計量の性質を確認しなさい．

人工データを用いて t 統計量の分布を確認する．

t 統計量とそのp値は以下のようにして取り出せる．

toy_lm <- lm(formula, toy_data)
summary(toy_lm)$coef[,c("t value","Pr(>|t|)")] # 列名での指定
summary(toy_lm)$coef[,3:4]                     # 列番号での指定
tidy(toy_lm) # 関数 broom::tidy() を用いてもよい

広告費と売上データを用いて係数の有意性を議論する．
東京の気候データを用いて係数の有意性を議論する．

解答例

人工データにより 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(){ 
    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(set_names(summary(toy_lm)$coef[,"t value"], # t統計量を返す
                     c("beta0.tval","beta1.tval","beta2.tval"))) 
}

x1 の係数は 2 なので帰無仮説に従わない，x2 の係数は 0 なので帰無仮説に従うことに注意する．

数値実験を行う．

mc_num <- 5000 # 実験回数
mc_data <- 
    replicate(mc_num, mc_trial()) |> t() |> as_tibble()

各回帰係数の 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-16

summary(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
Characteristic	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
¹ 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 
-2.3676 -1.1604  0.1108  1.0012  2.5560 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 386.22848   90.95731   4.246 0.000230 ***
solar         0.30199    0.04759   6.345 8.57e-07 ***
press        -0.36024    0.09025  -3.992 0.000453 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.436 on 27 degrees of freedom
Multiple R-squared:  0.6316,    Adjusted R-squared:  0.6043 
F-statistic: 23.14 on 2 and 27 DF,  p-value: 1.399e-06

summary(lm(temp ~ press, data = tw_subset))


Call:
lm(formula = temp ~ press, data = tw_subset)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.1953 -1.5311  0.1761  1.5901  3.0402 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 242.9579   136.5601   1.779   0.0861 .
press        -0.2142     0.1353  -1.584   0.1245  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.225 on 28 degrees of freedom
Multiple R-squared:  0.08221,   Adjusted R-squared:  0.04943 
F-statistic: 2.508 on 1 and 28 DF,  p-value: 0.1245

summary(lm(temp ~ solar, data = tw_subset))


Call:
lm(formula = temp ~ solar, data = tw_subset)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9840 -0.5948  0.5636  1.0317  1.8893 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 23.17563    0.85469  27.116  < 2e-16 ***
solar        0.25353    0.05699   4.449 0.000125 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.778 on 28 degrees of freedom
Multiple R-squared:  0.4142,    Adjusted R-squared:  0.3932 
F-statistic: 19.79 on 1 and 28 DF,  p-value: 0.0001248

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
Characteristic	Beta	SE ¹	Statistic	p-value	Beta	SE ¹	Statistic	p-value	Beta	SE ¹	Statistic	p-value
(Intercept)	386	91.0	4.25	<0.001	243	137	1.78	0.086	23	0.855	27.1	<0.001
solar	0.30	0.048	6.35	<0.001					0.25	0.057	4.45	<0.001
press	-0.36	0.090	-3.99	<0.001	-0.21	0.135	-1.58	0.12
R²	0.632				0.082				0.414
Adjusted R²	0.604				0.049				0.393
¹ SE = Standard Error

press単体では係数の推定精度も決定係数も低いが，solarと組み合わせることにより精度が上がり説明力も高くなることがわかる．

F 統計量の性質

問題

数値実験により F 統計量の性質を確認しなさい．

人工データを用いて F 統計量の分布を確認しなさい．

F 統計量とその自由度は以下のようにして取り出せる．

toy_lm <- lm(formula, toy_data)
summary(toy_lm)$fstat
summary(toy_lm)$fstatistic # 省略しない場合
glance(toy_lm) # 関数 broom::glance() を用いてもよい

広告費と売上データのモデルの有効性を議論しなさい．
東京の気候データのモデルの有効性を議論しなさい．

解答例

人工データにより 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(){ 
    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(summary(toy_lm)$fstat[1]) # F統計量を返す
}

上記の設定では x1,x2 の係数はどちらも0なので帰無仮説が成り立つことに注意する．まず，帰無仮説が成り立つ場合の数値実験を行う．

mc_num <- 5000 # 実験回数
mc_data <- 
    replicate(mc_num, mc_trial()) |> as_tibble_col("fstat")

実験結果は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 <-
    replicate(mc_num, mc_trial()) |> as_tibble_col("fstat")

モデルの 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-16

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-16

summary(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.001148

radio, 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 = c("モデル1","モデル2","モデル3")) |>
    modify_table_body( ~ .x |> # characteristic を並べ変える
                           arrange(factor(variable,
                                          levels = c("TV","radio","newspaper"))))

Characteristic	モデル1			モデル2			モデル3
Characteristic	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
¹ 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.1953 -1.5311  0.1761  1.5901  3.0402 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 242.9579   136.5601   1.779   0.0861 .
press        -0.2142     0.1353  -1.584   0.1245  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.225 on 28 degrees of freedom
Multiple R-squared:  0.08221,   Adjusted R-squared:  0.04943 
F-statistic: 2.508 on 1 and 28 DF,  p-value: 0.1245

summary(lm(temp ~ press + solar, data = tw_subset))


Call:
lm(formula = temp ~ press + solar, data = tw_subset)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.3676 -1.1604  0.1108  1.0012  2.5560 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 386.22848   90.95731   4.246 0.000230 ***
press        -0.36024    0.09025  -3.992 0.000453 ***
solar         0.30199    0.04759   6.345 8.57e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.436 on 27 degrees of freedom
Multiple R-squared:  0.6316,    Adjusted R-squared:  0.6043 
F-statistic: 23.14 on 2 and 27 DF,  p-value: 1.399e-06

summary(lm(temp ~ press + solar + rain, data = tw_subset))


Call:
lm(formula = temp ~ press + solar + rain, data = tw_subset)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.3862 -1.2316  0.2577  0.7963  2.4566 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 392.38652   90.53009   4.334 0.000195 ***
press        -0.36577    0.08980  -4.073 0.000386 ***
solar         0.26928    0.05504   4.893 4.46e-05 ***
rain         -0.01618    0.01393  -1.162 0.255967    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.427 on 26 degrees of freedom
Multiple R-squared:  0.6497,    Adjusted R-squared:  0.6093 
F-statistic: 16.08 on 3 and 26 DF,  p-value: 4.103e-06

press のみではモデルの有効性があるとは言えないが， 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 = c("モデル1","モデル2","モデル3"))

Characteristic	モデル1			モデル2			モデル3
Characteristic	Beta	95% CI ¹	p-value	Beta	95% CI ¹	p-value	Beta	95% CI ¹	p-value
press	-0.37	-0.55, -0.18	<0.001	-0.36	-0.55, -0.18	<0.001	-0.21	-0.49, 0.06	0.12
solar	0.27	0.16, 0.38	<0.001	0.30	0.20, 0.40	<0.001
rain	-0.02	-0.04, 0.01	0.3
Adjusted R²	0.609			0.604			0.049
Statistic	16.1			23.1			2.51
p-value	<0.001			<0.001			0.12
¹ CI = Confidence Interval