library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
#' 日本語を用いるので macOS ではフォントの設定を行う
if(Sys.info()["sysname"] == "Darwin") { # macOS か調べて日本語フォントを指定
theme_update(text = element_text(family = "HiraginoSans-W4"))
update_geom_defaults("text", list(family = theme_get()$text$family))
update_geom_defaults("label", list(family = theme_get()$text$family))}統計データ解析I
第13講 練習問題 解答例
準備
以下で利用する共通パッケージを読み込む．
回帰分析の確率シミュレーション
問題
回帰分析におけるモデルの推定量の精度に関する 確率シミュレーションを考えなさい．
解答例
人工データによる回帰モデルの推定を行う． 以下のモデルのパラメタは適当に変更してよい．
alpha <- 2 # 切片
beta <- 3 # 回帰係数
n <- 20 # データ数
sigma <- 0.5 # 誤差の標準偏差データを生成する．
x <- runif(n, min = -1, max = 1) # 説明変数 (区間[-1,1]を想定)
epsilon <- rnorm(n, sd = sigma) # 誤差 (正規分布)
y <- alpha + beta * x + epsilonデータの視覚化を行う
toy_data <- tibble(x = x, y = y)
gg <- # 後で描き加えるため保存しておく
toy_data |>
ggplot(aes(x = x, y = y)) +
geom_point(colour = "forestgreen")
print(gg)
回帰式の推定を行う．
toy_lm <- lm(y ~ x, data = toy_data)
coef(toy_lm) # 推定された係数の取得(Intercept) x
1.979740 2.701823 推定された回帰分析の結果を表示する． これらの方法については回帰係数の有意性検定でも詳述する．
#' base R での情報の表示
toy_lm |> summary() # さまざまな情報がlist形式
Call:
lm(formula = y ~ x, data = toy_data)
Residuals:
Min 1Q Median 3Q Max
-0.8331 -0.4578 -0.1656 0.4227 1.1536
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.9797 0.1525 12.98 1.42e-10 ***
x 2.7018 0.2456 11.00 2.02e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6664 on 18 degrees of freedom
Multiple R-squared: 0.8705, Adjusted R-squared: 0.8633
F-statistic: 121 on 1 and 18 DF, p-value: 2.023e-09#' tidyverse での情報の表示 (tibble形式)
toy_lm |> broom::tidy() # 推定された係数の情報# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1.98 0.153 13.0 1.42e-10
2 x 2.70 0.246 11.0 2.02e- 9toy_lm |> broom::glance() # 推定に関するさまざまな情報# 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.870 0.863 0.666 121. 0.00000000202 1 -19.2 44.4 47.4
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>toy_lm |> broom::augment() # 推定に用いられたデータの情報# A tibble: 20 × 8
y x .fitted .resid .hat .sigma .cooksd .std.resid
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 5.76 0.970 4.60 1.15 0.215 0.609 0.524 1.95
2 -0.113 -0.863 -0.351 0.238 0.122 0.683 0.0101 0.381
3 1.92 -0.192 1.46 0.456 0.0505 0.676 0.0131 0.702
4 1.93 0.255 2.67 -0.736 0.0704 0.660 0.0497 -1.15
5 3.26 0.543 3.45 -0.182 0.112 0.684 0.00532 -0.290
6 0.372 -0.673 0.163 0.209 0.0896 0.684 0.00533 0.329
7 2.25 0.131 2.33 -0.0867 0.0595 0.685 0.000568 -0.134
8 3.00 0.686 3.83 -0.833 0.141 0.650 0.149 -1.35
9 -0.322 -0.694 0.104 -0.425 0.0928 0.677 0.0230 -0.670
10 -0.608 -0.657 0.204 -0.812 0.0873 0.654 0.0778 -1.28
11 2.82 0.462 3.23 -0.404 0.0980 0.678 0.0221 -0.638
12 -1.25 -0.919 -0.504 -0.743 0.134 0.658 0.111 -1.20
13 -0.194 -0.957 -0.605 0.412 0.142 0.677 0.0369 0.667
14 0.518 -0.463 0.728 -0.211 0.0648 0.684 0.00370 -0.327
15 1.93 -0.386 0.937 0.993 0.0587 0.639 0.0736 1.54
16 1.90 0.0694 2.17 -0.266 0.0556 0.683 0.00496 -0.410
17 3.69 0.221 2.58 1.11 0.0670 0.627 0.107 1.72
18 0.882 -0.201 1.44 -0.555 0.0506 0.672 0.0195 -0.855
19 4.06 0.825 4.21 -0.149 0.175 0.685 0.00640 -0.246
20 0.610 -0.815 -0.222 0.832 0.113 0.651 0.112 1.33 推定結果の視覚化を行う．
gg +
geom_abline(intercept = alpha,
slope = beta,
colour = "red") + # 真の回帰式
geom_abline(intercept = coef(toy_lm)[1],
slope = coef(toy_lm)[2],
colour = "blue") # 推定された回帰式
Monte-Carlo 実験を行う．
mc_trial <- function(){
epsilon <- rnorm(n, sd = sigma)
y <- alpha + beta * x + epsilon # 説明変数は固定しておく
est <- lm(y ~ x) # データフレームにせずに直接 x,y を渡す
#' データフレームにする場合は以下のようにすればよい
#' lm(y ~ x, data = tibble(x = x, y = y))
return(coef(est))
}
mc_data <- # 実験結果をデータフレームに変換
replicate(2000, mc_trial()) |> t() |> as_tibble()推定値の分布を視覚化する．
mc_data |> # 切片
ggplot(aes(x = `(Intercept)`)) +
geom_density(fill = "pink") +
geom_vline(xintercept = alpha, colour = "orange") +
labs(x = expression(hat(alpha)), title = "切片の分布")
mc_data |> # 傾き
ggplot(aes(x = x)) +
geom_density(fill = "palegreen") +
geom_vline(xintercept = beta, colour = "darkgreen") +
labs(x = expression(hat(beta)), title = "傾きの分布")
推定された回帰式のばらつきを表示する．
mc_data |>
slice_sample(n = 40) |> # Monte-Carlo 実験から40個ランダムに選択
rowid_to_column() |> # 番号列(rowid)を作成
ggplot() +
geom_abline(aes(intercept = `(Intercept)`,
slope = x,
colour = as_factor(rowid)), # 色を変える
alpha = 0.4, # alpha値を低くして色を薄くする
show.legend = FALSE) + # 凡例は表示しない
geom_abline(intercept = alpha,
slope = beta,
colour = "red",
linewidth = 1.1) + # 真の回帰式(太め)
xlim(-1,1) + ylim(alpha-beta*1.1, alpha+beta*1.1) + # 描画範囲の指定
labs(title = "推定された回帰式のばらつき")
同じ生成モデルでもデータによって推定結果がばらつくことがわかる．
回帰モデルの点推定
問題
東京の気象データを用いて， 必要であれば適当な期間を抽出し， 日射量から気温を説明する回帰モデルを構成しなさい．
解答例
データの読み込みを行う．
tw_data <- read_csv("data/tokyo_weather.csv")データの散布図 (1年分) を描く．
tw_data |>
ggplot(aes(x = solar, y = temp)) +
geom_point(colour = "olivedrab") +
labs(x = "日射量", y = "気温")
回帰式の推定を行い，その結果を表示する．
tw_lm <- lm(temp ~ solar, # 目的変数 ~ 説明変数
data = tw_data) # 気温を日射量で説明
tw_lm |> summary() # 結果の要約
Call:
lm(formula = temp ~ solar, data = tw_data)
Residuals:
Min 1Q Median 3Q Max
-15.755 -7.149 0.513 7.468 14.345
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.19919 0.89429 12.523 < 2e-16 ***
solar 0.46433 0.05749 8.077 9.82e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.707 on 364 degrees of freedom
Multiple R-squared: 0.152, Adjusted R-squared: 0.1497
F-statistic: 65.24 on 1 and 364 DF, p-value: 9.825e-15tw_lm |> broom::tidy() # 係数とその統計量# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 11.2 0.894 12.5 3.64e-30
2 solar 0.464 0.0575 8.08 9.82e-15tw_lm |> broom::glance() # その他の統計量# 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.152 0.150 7.71 65.2 9.82e-15 1 -1266. 2538. 2549.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>回帰直線を前の図に重ね描きして図示する．
last_plot() +
geom_abline(intercept = coef(tw_lm)[1], # 切片
slope = coef(tw_lm)[2], # 回帰係数
colour = "slateblue",
linewidth = 1.2) 
Note
関数 geom_smooth() を用いることもできる．
tw_data |>
ggplot(aes(x = solar, y = temp)) +
geom_point(colour = "olivedrab") +
geom_smooth(method = "lm", # 関数 lm を用いて信頼区間(平滑化方法)を計算
se = FALSE, # 信頼区間を付けない
colour = "slateblue") +
labs(x = "日射量", y = "気温")
期間を限って分析する． 夏のモデルのためのデータの散布図を描く．
tw_data |>
filter(month %in% 7:9) |> # 7月-9月を抽出
ggplot(aes(x = solar, y = temp)) +
geom_point(colour = "olivedrab") +
labs(x = "日射量", y = "気温")
回帰式の推定を行い，その結果を表示する．
tw_lm2 <- lm(formula(tw_lm), # 前の式を利用
data = tw_data,
subset = month %in% 7:9) # 7月-9月を抽出
tw_lm2 |> summary()
Call:
lm(formula = formula(tw_lm), data = tw_data, subset = month %in%
7:9)
Residuals:
Min 1Q Median 3Q Max
-6.2489 -0.8316 0.3418 1.0467 4.1866
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.43954 0.51095 45.874 < 2e-16 ***
solar 0.28037 0.02872 9.761 8.95e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.732 on 90 degrees of freedom
Multiple R-squared: 0.5143, Adjusted R-squared: 0.5089
F-statistic: 95.28 on 1 and 90 DF, p-value: 8.951e-16
Note
期間を限った分析は以下のように書くこともできる．
tw_data |>
filter(month %in% 7:9) |> # 7月-9月を抽出
lm(formula(tw_lm), data = _) # パイプ演算の内容を data に渡す
Call:
lm(formula = formula(tw_lm), data = filter(tw_data, month %in%
7:9))
Coefficients:
(Intercept) solar
23.4395 0.2804 推定結果の可視化を行う．
last_plot() +
geom_abline(intercept = coef(tw_lm2)[1],
slope = coef(tw_lm2)[2],
colour = "tomato",
linewidth = 1.2) 
Note
関数 geom_smooth() を用いてもよい．
tw_data |>
filter(month %in% 7:9) |>
ggplot(aes(x = solar, y = temp)) +
geom_point(colour = "olivedrab") +
geom_smooth(method = "lm",
se = FALSE,
colour = "tomato") +
labs(x = "日射量", y = "気温")
全データのモデルと夏のモデルを比較する．
tw_data |>
mutate(summer = ifelse(month %in% 7:9,
"Jul-Sep", # 7月-9月のラベル
"others")) |>
ggplot(aes(x = solar, y = temp, colour = summer)) +
geom_point() +
geom_abline(intercept = coef(tw_lm)[1],
slope = coef(tw_lm)[2],
colour = "slateblue",
linewidth = 1.2) +
geom_abline(intercept = coef(tw_lm2)[1],
slope = coef(tw_lm2)[2],
colour = "tomato",
linewidth = 1.2) +
labs(x = "日射量", y = "気温", colour = NULL)
Note
関数 geom_smooth() を用いる場合は以下のようにすればよい．
tw_data |>
mutate(summer = ifelse(month %in% 7:9,
"7-9月", # 7月-9月のラベル
"その他")) |>
ggplot(aes(x = solar, y = temp, colour = summer)) +
geom_point() +
geom_smooth(method = "lm",
se = FALSE,
colour = "slateblue") +
geom_smooth(data = \(x)filter(x, summer == "7-9月"), # 7-9月に限定
method = "lm",
se = FALSE,
colour = "tomato") +
labs(x = "日射量", y = "気温")
data に一変数関数を渡すと関数 geom_smooth() 内での処理を限定することができる． 上記は base R での無名関数の記述を用いているが，tidyverse では
~ filter(., summer == "7-9月")
. %>% filter(summer == "7-9月")のように書くこともできる．
回帰モデルの区間推定
問題
前問で作成した回帰モデルについて区間推定を行いなさい．
解答例
前問で作成した1年分のモデル tw_lm を用いて区間推定を行う．
confint(tw_lm) 2.5 % 97.5 %
(Intercept) 9.4405734 12.9578089
solar 0.3512825 0.5773738区間推定を視覚化する． 関数 geom_smooth() を用いると簡潔に記述できる．
tw_data |>
ggplot(aes(x = solar, y = temp)) +
geom_point(aes(colour = factor(month))) + # 月ごとに色を変える
geom_smooth(method = "lm", # 関数 lm を用いて信頼区間(平滑化方法)を計算
level = 0.95, # 既定値なので無くても良い
colour = "royalblue", # 線の色
fill = "steelblue") + # 塗り潰しの色
labs(x = "日射量", y = "気温", colour = NULL)
Note
上記の描画を愚直に行うには以下のような手続きを考えればよい．
視覚化したい信頼区間に合わせて適切な説明変数と信頼区間を作成する．
tw_conf <-
tw_lm |>
broom::augment(newdata = tibble(solar = tw_data |> pull(solar) |>
range() |> # 日射量の範囲を取得
pretty(n = 50)), # 50個程度の区間を作成
interval = "confidence",
conf.level = 0.95) # 既定値なので無くてもよいデータ点と信頼区間を描画する．
tw_data |>
ggplot(aes(x = solar, y = temp)) +
geom_point(aes(colour = factor(month))) +
geom_abline(intercept = coef(tw_lm)[1],
slope = coef(tw_lm)[2],
colour = "royalblue",
linewidth = 1.2) +
geom_ribbon(data = tw_conf, # temp 列の代わりに .fitted列がある
aes(x = solar, y = .fitted,
ymin = .lower, ymax = .upper),
fill = alpha("steelblue", 0.3)) +
labs(x = "日射量", y = "気温", colour = NULL)
データ数がある程度多い場合はデータの説明変数をそのまま用いてもよい．
tw_lm |>
broom::augment(interval = "confidence") |>
ggplot(aes(x = solar, y = temp)) +
geom_point(colour = "olivedrab") + # augmentの出力にmonthは含まれていないため
geom_line(aes(y = .fitted),
colour = "royalblue") +
geom_ribbon(aes(ymin = .lower, ymax = .upper),
fill = alpha("steelblue", 0.3)) +
labs(x = "日射量", y = "気温")
データに信頼区間を追加する方法もある．
tw_data |>
mutate(broom::augment(tw_lm, interval = "confidence")) |>
ggplot(aes(x = solar, y = temp)) +
geom_point(aes(colour = factor(month))) +
geom_line(aes(y = .fitted),
colour = "royalblue") +
geom_ribbon(aes(ymin = .lower, ymax = .upper),
fill = alpha("steelblue", 0.3)) +
labs(x = "日射量", y = "気温", colour = NULL)
tw_data |>
mutate(predict(tw_lm, # 推定したモデル
interval = "confidence", # 信頼区間の推定を指定
level = 0.95) |>
as_tibble()) |> # tibble形式に変換して追加する (列名に注意)
ggplot(aes(x = solar, y = temp)) +
geom_point(aes(colour = factor(month))) + # augmentの出力にmonthは含まれていないため
geom_line(aes(y = fit), # predict の出力は .fitted でない
colour = "royalblue") +
geom_ribbon(aes(ymin = lwr, ymax = upr), # predict の出力は .lower/.upper でない
fill = alpha("steelblue", 0.3)) +
labs(x = "日射量", y = "気温", colour = NULL)
夏のモデル tw_lm2 で同様に区間推定を行う．
confint(tw_lm2) 2.5 % 97.5 %
(Intercept) 22.4244444 24.4546376
solar 0.2233102 0.3374357視覚化を行う
model.frame(tw_lm2) |> # モデルの作成に用いたデータ
ggplot(aes(x = solar, y = temp)) +
geom_point(colour = "orange") +
geom_smooth(method = "lm", level = 0.95,
colour = "red", fill = "pink") +
labs(x = "日射量", y = "気温", colour = NULL)
1年のデータも重ねて表示する．
tw_lm2 |> broom::augment() |> # この関数でも取り出せる
ggplot(aes(x = solar, y = temp)) +
geom_point(data = tw_data, # 1年分のデータを使用
colour = alpha("grey", 0.5)) + # 薄い灰色で表示
geom_point(colour = "orange") +
geom_smooth(method = "lm", level = 0.95,
colour = "red", fill = "pink") +
labs(x = "日射量", y = "気温", colour = NULL)
1年のデータから限定して信頼区間を描いてもよい．
tw_data |>
ggplot(aes(x = solar, y = temp)) +
geom_point(aes(colour = factor(month))) +
geom_smooth(data = \(x)filter(x, month %in% 7:9),
method = "lm", level = 0.95,
colour = "red", fill = "pink") +
labs(x = "日射量", y = "気温", colour = NULL)
回帰モデルの係数の検定
問題
前問で作成した回帰モデルについて係数の検定を行いなさい．
解答例
前問で構成したモデルを用いる． 1年分のモデル tw_lm を用いた場合の base R での操作は以下のようになる．
summary(tw_lm)
Call:
lm(formula = temp ~ solar, data = tw_data)
Residuals:
Min 1Q Median 3Q Max
-15.755 -7.149 0.513 7.468 14.345
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.19919 0.89429 12.523 < 2e-16 ***
solar 0.46433 0.05749 8.077 9.82e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.707 on 364 degrees of freedom
Multiple R-squared: 0.152, Adjusted R-squared: 0.1497
F-statistic: 65.24 on 1 and 364 DF, p-value: 9.825e-15情報が多いので，整理してみる
summary(tw_lm)$coef # 名前は識別できれば途中まででも可 Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.1991911 0.89428702 12.523039 3.641899e-30
solar 0.4643281 0.05748562 8.077292 9.824715e-15summary(tw_lm)$coef["solar",c("t value","Pr(>|t|)")] t value Pr(>|t|)
8.077292e+00 9.824715e-15 summary(tw_lm)$coef[2,3:4] # 上と同じ t value Pr(>|t|)
8.077292e+00 9.824715e-15 summary(tw_lm)$fstat # F統計量 (モデルの有意性の評価) value numdf dendf
65.24265 1.00000 364.00000 anova(tw_lm)$'Pr(>F)'[1] # F統計量のp値 (summaryからは取り出し難い)[1] 9.824715e-15パイプを使う場合は以下のようにすればよい．
tw_lm |> summary() |> _$coef Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.1991911 0.89428702 12.523039 3.641899e-30
solar 0.4643281 0.05748562 8.077292 9.824715e-15tw_lm |> summary() |> _$coef["solar",c("t value","Pr(>|t|)")] t value Pr(>|t|)
8.077292e+00 9.824715e-15 tw_lm |> summary() |> _$coef[2,3:4] t value Pr(>|t|)
8.077292e+00 9.824715e-15 tidyverse での操作は以下のようになる．
tw_lm |> broom::tidy() # 係数の情報をまとめたデータフレーム# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 11.2 0.894 12.5 3.64e-30
2 solar 0.464 0.0575 8.08 9.82e-15tw_lm |> broom::tidy() |> # solarの係数のt統計量とp値を抽出
filter(term == "solar") |> select(statistic, p.value)# A tibble: 1 × 2
statistic p.value
<dbl> <dbl>
1 8.08 9.82e-15tw_lm |> broom::tidy() |> _[2, c("statistic","p.value")]# A tibble: 1 × 2
statistic p.value
<dbl> <dbl>
1 8.08 9.82e-15tw_lm |> broom::tidy() |> _[2,4:5]# A tibble: 1 × 2
statistic p.value
<dbl> <dbl>
1 8.08 9.82e-15tw_lm |> broom::glance() # モデルの情報をまとめたデータフレーム# 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.152 0.150 7.71 65.2 9.82e-15 1 -1266. 2538. 2549.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>tw_lm |> broom::glance() |> # F統計量に関する情報
select(statistic, p.value, df, df.residual)# A tibble: 1 × 4
statistic p.value df df.residual
<dbl> <dbl> <dbl> <int>
1 65.2 9.82e-15 1 364夏のモデル tw_lm2 の base R での操作は以下のようになる．
summary(tw_lm2)
Call:
lm(formula = formula(tw_lm), data = tw_data, subset = month %in%
7:9)
Residuals:
Min 1Q Median 3Q Max
-6.2489 -0.8316 0.3418 1.0467 4.1866
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.43954 0.51095 45.874 < 2e-16 ***
solar 0.28037 0.02872 9.761 8.95e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.732 on 90 degrees of freedom
Multiple R-squared: 0.5143, Adjusted R-squared: 0.5089
F-statistic: 95.28 on 1 and 90 DF, p-value: 8.951e-16coef(summary(tw_lm2)) # 関数coefでも可 Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.439541 0.51095263 45.874196 3.264630e-64
solar 0.280373 0.02872275 9.761356 8.950908e-16coef(summary(tw_lm2))["solar",c("t value","Pr(>|t|)")] t value Pr(>|t|)
9.761356e+00 8.950908e-16 coef(tw_lm2) # 推定された係数のみ取り出す場合(Intercept) solar
23.439541 0.280373 coef(summary(tw_lm2))[,"Estimate"] # 上と同じ(Intercept) solar
23.439541 0.280373 tidyverse での操作は以下のようになる．
tw_lm2 |> broom::tidy() # 係数の情報# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 23.4 0.511 45.9 3.26e-64
2 solar 0.280 0.0287 9.76 8.95e-16tw_lm2 |> broom::glance() # モデルの情報# 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.514 0.509 1.73 95.3 8.95e-16 1 -180. 366. 374.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>決定係数による回帰モデルの検討
問題
前問で作成した回帰モデルについて 決定係数を確認しなさい． また，説明変数として降水量を用いた回帰モデルについて 検討しなさい．
解答例
前問で構成した1年分のモデル tw_lm を用いる．
summary(tw_lm) # 全情報の表示
Call:
lm(formula = temp ~ solar, data = tw_data)
Residuals:
Min 1Q Median 3Q Max
-15.755 -7.149 0.513 7.468 14.345
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.19919 0.89429 12.523 < 2e-16 ***
solar 0.46433 0.05749 8.077 9.82e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.707 on 364 degrees of freedom
Multiple R-squared: 0.152, Adjusted R-squared: 0.1497
F-statistic: 65.24 on 1 and 364 DF, p-value: 9.825e-15summary(tw_lm)$r.squared[1] 0.1519948summary(tw_lm)$adj.r.squared[1] 0.1496651tw_lm |> broom::glance() # 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.152 0.150 7.71 65.2 9.82e-15 1 -1266. 2538. 2549.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>tw_lm |> broom::glance() |> # 1:2列が決定係数
select(r.squared, adj.r.squared)# A tibble: 1 × 2
r.squared adj.r.squared
<dbl> <dbl>
1 0.152 0.150夏のモデル tw_lm2 の場合は以下のようになる．
summary(tw_lm2) # 全情報の表示
Call:
lm(formula = formula(tw_lm), data = tw_data, subset = month %in%
7:9)
Residuals:
Min 1Q Median 3Q Max
-6.2489 -0.8316 0.3418 1.0467 4.1866
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.43954 0.51095 45.874 < 2e-16 ***
solar 0.28037 0.02872 9.761 8.95e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.732 on 90 degrees of freedom
Multiple R-squared: 0.5143, Adjusted R-squared: 0.5089
F-statistic: 95.28 on 1 and 90 DF, p-value: 8.951e-16summary(tw_lm2)$r.squared[1] 0.5142594summary(tw_lm2)$adj.r.squared[1] 0.5088623tw_lm2 |> broom::glance() # 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.514 0.509 1.73 95.3 8.95e-16 1 -180. 366. 374.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>tw_lm2 |> broom::glance() |>
select(r.squared, adj.r.squared)# A tibble: 1 × 2
r.squared adj.r.squared
<dbl> <dbl>
1 0.514 0.509降水量と気温の関係を調べる．
tw_lm3 <- lm(temp ~ rain, data = tw_data)
tw_lm4 <- lm(formula(tw_lm3), # 上の式を用いる
data = tw_data,
subset = month %in% 7:9) # 夏(7-9月)のモデル1年分のモデル tw_lm3 を調べる．
tw_lm3 |> summary()
Call:
lm(formula = temp ~ rain, data = tw_data)
Residuals:
Min 1Q Median 3Q Max
-16.6101 -7.9009 0.1991 6.7916 14.9991
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.30090 0.46512 37.196 <2e-16 ***
rain 0.06598 0.03136 2.104 0.0361 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.319 on 364 degrees of freedom
Multiple R-squared: 0.01201, Adjusted R-squared: 0.009298
F-statistic: 4.426 on 1 and 364 DF, p-value: 0.03609model.frame(tw_lm3) |> # 推定に用いたデータを利用
ggplot(aes(x = rain, y = temp)) +
geom_point(colour = "blue") +
geom_smooth(method = "lm",
colour = "red", fill = "pink") +
labs(x = "降水量", y = "気温")
1年分のモデル tw_lm3 に有意性はないことがわかる．
夏のモデル tw_lm4 を調べる．
tw_lm4 |> summary()
Call:
lm(formula = formula(tw_lm3), data = tw_data, subset = month %in%
7:9)
Residuals:
Min 1Q Median 3Q Max
-7.7497 -0.8979 0.5539 1.7860 4.0860
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.21401 0.28153 100.218 <2e-16 ***
rain -0.01429 0.01491 -0.958 0.34
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.472 on 90 degrees of freedom
Multiple R-squared: 0.0101, Adjusted R-squared: -0.0008977
F-statistic: 0.9184 on 1 and 90 DF, p-value: 0.3405model.frame(tw_lm4) |> # 推定に用いたデータを利用
ggplot(aes(x = rain, y = temp)) +
geom_point(colour = "orange") +
geom_smooth(method = "lm",
colour = "red", fill = "pink") +
labs(x = "降水量", y = "気温")
夏場は雨が降ると気温が下がる傾向が有意にあることが読み取れる． 決定係数が低いのはそもそも気温のばらつきが大きいことに起因すると考えられる．
補遺
ggplot2 を用いた線形回帰分析(単回帰)の例
#' パッケージの読み込み
#' 以下のほかに 'ggrepel' 'broom' 'magrittr' 'MASS' を利用
library(tidyverse)
library(ggfortify) データの読み込み (“MASS::Animals”を用いる) を行う．
data(Animals, package = "MASS")以下 Animals で参照可能となる．
データの内容を確認する．
help(Animals, package = "MASS") # 内容の詳細を表示
print(Animals) # データの表示 body brain
Mountain beaver 1.350 8.1
Cow 465.000 423.0
Grey wolf 36.330 119.5
Goat 27.660 115.0
Guinea pig 1.040 5.5
Dipliodocus 11700.000 50.0
Asian elephant 2547.000 4603.0
Donkey 187.100 419.0
Horse 521.000 655.0
Potar monkey 10.000 115.0
Cat 3.300 25.6
Giraffe 529.000 680.0
Gorilla 207.000 406.0
Human 62.000 1320.0
African elephant 6654.000 5712.0
Triceratops 9400.000 70.0
Rhesus monkey 6.800 179.0
Kangaroo 35.000 56.0
Golden hamster 0.120 1.0
Mouse 0.023 0.4
Rabbit 2.500 12.1
Sheep 55.500 175.0
Jaguar 100.000 157.0
Chimpanzee 52.160 440.0
Rat 0.280 1.9
Brachiosaurus 87000.000 154.5
Mole 0.122 3.0
Pig 192.000 180.0データのプロット (normal plot) を行う．
Animals |>
ggplot(aes(body, brain)) + # x軸，y軸に用いる列の指定
geom_point(colour = "royalblue") + # 点の追加
labs(title = "Brain and Body Weights (normal plot)",
x = "body [kg]", y = "brain [g]") # タイトルと軸名の追加
データのプロット (log plot) を行う．
Animals |>
ggplot(aes(body, brain)) +
geom_point(colour = "royalblue") +
scale_x_log10() + scale_y_log10() + # log-log plot を指定
labs(title = "Brain and Body Weights (log-log plot)",
x = "body [kg]", y="brain [g]")
データの分布から両対数変換が分析においては適切であることがわかる．
回帰分析 (単回帰) を行う．
bb_lm <- lm(log(brain) ~ log(body), # 対数変換した変数で線形回帰
data = Animals)
bb_lm |> summary() # 分析結果のまとめを表示
Call:
lm(formula = log(brain) ~ log(body), data = Animals)
Residuals:
Min 1Q Median 3Q Max
-3.2890 -0.6763 0.3316 0.8646 2.5835
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.55490 0.41314 6.184 1.53e-06 ***
log(body) 0.49599 0.07817 6.345 1.02e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.532 on 26 degrees of freedom
Multiple R-squared: 0.6076, Adjusted R-squared: 0.5925
F-statistic: 40.26 on 1 and 26 DF, p-value: 1.017e-06信頼区間付きで回帰式を表示する．
Animals |>
rownames_to_column() |>
ggplot(aes(body, brain)) +
geom_point(colour = alpha("royalblue", 0.75)) +
ggrepel::geom_text_repel(aes(label = rowname), # 各点の名前を追加
size = 3) +
geom_smooth(method = "lm", # 回帰式
colour = "dodgerblue",
fill = "dodgerblue") +
scale_x_log10() + scale_y_log10() + # log-log plot
labs(title = "Brain and Body Weights",
x = "body [kg]", y = "brain [g]")
信頼区間および予測区間付きで回帰式を表示する．
Animals |>
rownames_to_column() |>
mutate( # 予測区間の情報を追加する
broom::augment(bb_lm,
interval = "prediction", # 予測区間を指定
conf.level = 0.95) |> # 既定値なので省略可能
select(.fitted, .lower, .upper) |> exp() # 必要な列を選択して対数変換を戻す
) |>
ggplot(aes(body, brain)) +
geom_point(colour = "royalblue") +
ggrepel::geom_text_repel(aes(label = rowname), # 各点の名前を追加
size = 3) + # 文字の大きさを調整
geom_smooth(method = "lm", # 回帰式
colour = "dodgerblue",
fill = "dodgerblue") +
geom_ribbon(aes(ymin = .lower, ymax = .upper), # 予測区間の描画
fill = alpha("dodgerblue", 0.2)) +
scale_x_log10() + scale_y_log10() + # log-log plot
labs(title = "Brain and Body Weights",
x = "body [kg]", y = "brain [g]")
Note
予測区間の計算は関数 predict() を用いてもよい．
mutate( # 予測区間の情報を追加する
predict(bb_lm,
newdata = Animals, # 対数変換のため指定しないと警告が出る
interval = "prediction") |>
exp() |> # 対数変換を戻す
as_tibble() # tibble形式でfit/lwr/upr列を追加
) 診断プロットを表示する．
autoplot(bb_lm,
colour = "royalblue",
smooth.colour = "gray50",
smooth.linetype = "dashed",
ad.colour = "blue",
label.size = 3, # 以下，外れ値ラベルに関する設定
label.n = 5,
label.colour = "red")
外れ値を除いた回帰分析を行う．
idx <- c(6,16,26) # 外れ値のindex (いずれも恐竜)
bb_lm <- lm(log(brain) ~ log(body),
data = Animals,
subset = -idx) # 外れ値を除く
bb_lm |> summary()
Call:
lm(formula = log(brain) ~ log(body), data = Animals, subset = -idx)
Residuals:
Min 1Q Median 3Q Max
-0.9125 -0.4752 -0.1557 0.1940 1.9303
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.15041 0.20060 10.72 2.03e-10 ***
log(body) 0.75226 0.04572 16.45 3.24e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7258 on 23 degrees of freedom
Multiple R-squared: 0.9217, Adjusted R-squared: 0.9183
F-statistic: 270.7 on 1 and 23 DF, p-value: 3.243e-14信頼区間および予測区間付き回帰式を表示する．
Animals |>
rownames_to_column() |>
mutate(
broom::augment(bb_lm,
newdata = Animals, # 外れ値を除いたので全データで計算
interval = "prediction") |>
select(.fitted, .lower, .upper) |> exp()) |>
ggplot(aes(body, brain)) +
geom_point(colour = "royalblue") +
ggrepel::geom_text_repel(aes(label = rowname), # 各点の名前を追加
size = 3) + # 文字の大きさを調整
geom_smooth(data = \(x)slice(x, -idx), # 外れ値を除く
method = "lm", # 回帰式
colour = "dodgerblue",
fill = "dodgerblue") +
geom_ribbon(data = \(x)slice(x, -idx), # 外れ値を除く
aes(ymin = .lower, ymax = .upper), # 予測区間の描画
fill = alpha("dodgerblue", 0.2)) +
scale_x_log10() + scale_y_log10() + # log-log plot
labs(title = "Brain and Body Weights",
x = "body [kg]", y = "brain [g]")
恐竜の領域でも信頼・予測区間を描画するには，少し工夫が必要となる．
bb_new_body <- # 計算用のbody列を作成
tibble(body = Animals |> pull(body) |>
log() |> pretty(n=32) |> exp())
bb_new_data <- # 作成したbody列に対する区間の計算
tibble(bb_new_body,
broom::augment(bb_lm,
newdata = bb_new_body,
interval = "confidence") |> # 信頼区間
select(.fitted, .lower, .upper) |> exp() |>
magrittr::set_colnames(c("brain","c.lower","c.upper")), # 列名の指定
broom::augment(bb_lm,
newdata = bb_new_body,
interval = "prediction") |> # 予測区間
select(.lower, .upper) |> exp() |>
magrittr::set_colnames(c("p.lower","p.upper"))) # 列名の指定
Animals |>
rownames_to_column() |>
ggplot(aes(body, brain)) +
geom_point(colour = "royalblue") +
ggrepel::geom_text_repel(aes(label = rowname),
size = 3) +
geom_line(data = bb_new_data, # 回帰式の描画
colour = "dodgerblue") +
geom_ribbon(data = bb_new_data,
aes(ymin = c.lower, ymax = c.upper), # 信頼区間の描画
fill = alpha("dodgerblue", 0.4)) +
geom_ribbon(data = bb_new_data,
aes(ymin = p.lower, ymax = p.upper), # 予測区間の描画
fill = alpha("dodgerblue", 0.2)) +
scale_x_log10() + scale_y_log10() +
labs(title = "Brain and Body Weights",
x = "body [kg]", y = "brain [g]")
診断プロットを表示する．
autoplot(bb_lm,
colour = "royalblue",
smooth.colour = "gray50",
smooth.linetype = "dashed",
ad.colour = "blue",
label.size = 3, # 以下，外れ値ラベルに関する設定
label.n = 5,
label.colour = "red")