```{r}
library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(broom) # 解析結果を tibble 形式に集約
library(gt) # 表の作成
#' 日本語表示の設定．最後の例で利用する
if (Sys.info()["sysname"] == "Darwin") { # macOSのみ設定
jp_font <- "HiraMaruProN-W4"
theme_update(text = element_text(family = jp_font))
} else {
jp_font <- NULL
}
```第3講 回帰分析
回帰モデルの考え方と推定
準備
以下で利用する共通パッケージを読み込む．1
1 : 必要なパッケージは適宜追加して良い．
回帰係数の推定
線形回帰のための関数
回帰係数は関数 stats::lm() を用いて推定することができる．
lm(formula, data, subset, weights, na.action,
method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, offset, ...)
#' formula: 目的変数名 ~ 説明変数名(複数ある場合は + で並べる)
#' data: 目的変数，説明変数を含むデータフレーム
#' subset: 推定に用いるデータフレームの部分集合を指定(指定しなければ全て)
#' na.action: 欠損の扱いを指定(既定値はoption("na.action")で設定された処理)
#' model,x,y,qr: 返値にmodel.frame,model.matrix,目的変数,QR分解を含むか指定実行例
配布データ wine.csv を用いた ボルドーワインの価格と気候の関係に関する回帰分析は以下で実行できる．
```{r}
bw_data <-
read_csv(file = "data/wine.csv") # データの読み込み
bw_data # データの表示(一部)
```# A tibble: 38 × 6
VINT LPRICE2 WRAIN DEGREES HRAIN TIME_SV
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1952 -0.999 600 17.1 160 31
2 1953 -0.454 690 16.7 80 30
3 1954 NA 430 15.4 180 29
4 1955 -0.808 502 17.2 130 28
5 1956 NA 440 15.6 140 27
6 1957 -1.51 420 16.1 110 26
7 1958 -1.72 582 16.4 187 25
8 1959 -0.418 485 17.5 187 24
9 1960 -1.97 763 16.4 290 23
10 1961 0 830 17.3 38 22
# ℹ 28 more rows```{r}
bw_lm1 <-
lm(formula = LPRICE2 ~ WRAIN + HRAIN, # 冬と収穫期の雨で価格を予測
data = bw_data)
bw_lm1 # 推定結果の簡単な表示
```
Call:
lm(formula = LPRICE2 ~ WRAIN + HRAIN, data = bw_data)
Coefficients:
(Intercept) WRAIN HRAIN
-8.102e-01 -5.447e-06 -4.408e-03 ```{r}
bw_lm2 <-
lm(formula = LPRICE2 ~ . - VINT, # VINTを除く全てで価格を予測
data = bw_data)
bw_lm2 # 推定結果の簡単な表示
```
Call:
lm(formula = LPRICE2 ~ . - VINT, data = bw_data)
Coefficients:
(Intercept) WRAIN DEGREES HRAIN TIME_SV
-12.145334 0.001167 0.616392 -0.003861 0.023847 パイプ (|>) を利用して， “データフレーム |> 操作” という形で統一的に書きたいのであれば， _ (プレースフォルダ; placeholder) を利用して以下のように書くことができる．
bw_fit <-
bw_data |>
lm(formula = LPRICE2 ~ . - VINT, data = _)問題 (広告費と売上データ)
データセット https://www.statlearning.com/s/Advertising.csv 2 は広告費(TV,radio,newspapers)と売上の関係を調べたものである． このデータセットを用いて以下の回帰式を推定しなさい．
2 : Datasets in this presentation are taken from “An Introduction to Statistical Learning, with applications in R” (Springer, 2013) with permission from the authors: G. James, D. Witten, T. Hastie and R. Tibshirani. 参考
formula = sales ~ TV
formula = sales ~ radio
formula = sales ~ TV + radio解答欄
```{r}
#' ここにコードを書く
```解答例
データの読み込みを行う．
adv_data <-
read_csv("https://www.statlearning.com/s/Advertising.csv")インターネット上にあるデータをローカルに保存しておくには， 例えば以下のようにすればよい． 2回目以降はローカルファイルを参照することになる．
```{r}
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) # ローカルファイルに保存する
}
```3 : 標準設定では図は横長で描画される．この例では縦横が同じ大きさ(7inch)になるように指定している．
GGally パッケージを利用して散布図を描く．3
Code
```{r}
#| code-fold: true
#| fig-width: 6
#| fig-height: 6
library(GGally)
ggpairs(adv_data, columns = 2:5)
```
“TV” の宣伝費で売上を説明する．
Code
```{r}
#| code-fold: true
adv_lm1 <- lm(sales ~ TV, data = adv_data)
adv_lm1
#' ()はprint()と等価で．以下のようにしても良い
#' (adv_lm1 <- lm(sales ~ TV, data = adv_data))
adv_data |>
ggplot(aes(x = TV, y = sales)) +
geom_point(colour = "orange") +
geom_smooth(method = lm, se = FALSE)
```
Call:
lm(formula = sales ~ TV, data = adv_data)
Coefficients:
(Intercept) TV
7.03259 0.04754 
“radio” の宣伝費で売上を説明する．
Code
```{r}
#| code-fold: true
adv_lm2 <- lm(sales ~ radio, data = adv_data)
adv_lm2
adv_data |>
ggplot(aes(x = radio, y = sales)) +
geom_point(colour = "orange") +
geom_smooth(method = lm, se = FALSE)
```
Call:
lm(formula = sales ~ radio, data = adv_data)
Coefficients:
(Intercept) radio
9.3116 0.2025 
両者の宣伝費で売上を説明する．
Code
```{r}
#| code-fold: true
adv_lm3 <- lm(sales ~ TV + radio, data = adv_data)
adv_lm3
```
Call:
lm(formula = sales ~ TV + radio, data = adv_data)
Coefficients:
(Intercept) TV radio
2.92110 0.04575 0.18799 3次元の散布図を描くには scatterplot3d パッケージを利用することができる．
```{r}
#| label: fig-scatterplot3d
#| fig-cap: "scatterplot3d による3次元散布図の例"
#| column: margin
library(scatterplot3d)
s3d <- with(adv_data,
scatterplot3d(
x = TV, y = radio, z = sales,
type = "p", # plotの種類: "p"点，"l"線，"h"足付き
pch = 16, # 点の種類 (?points 参照)
angle = 45, # xy平面の見る方向 (適宜調整)
highlight.3d = TRUE # 高さ(z)ごとに色を変える
))
s3d$plane3d(adv_lm3, col = "blue", # 回帰式の定める平面の追加
draw_polygon = TRUE, # 平面の塗り潰しの設定
polygon_args = list(col = rgb(0,0,1,0.2)))
```
問題 (東京の気候データ)
配布したデータセット tokyo_weather.csv は 気象庁 4 より取得した東京の気候データを回帰分析用に整理したものである． このデータセットのうち，9月の気候データを用いて，以下の回帰式を推定しなさい．
4 : 気象庁のサイト
formula = temp ~ solar + press解答欄
```{r}
#' ここにコードを書く
```解答例
データの読み込みを行う．
Code
```{r}
#| code-fold: true
tw_data <- read_csv("data/tokyo_weather.csv")
```9月の”気温”を目的変数，“日射量・気圧”を説明変数とするモデルの推定を行う．
Code
```{r}
#| code-fold: true
tw_formula <- temp ~ solar + press # モデル式の定義
class(tw_formula) # formula クラスであることを確認
tw_lm <- lm(tw_formula, # 回帰係数の推定
data = tw_data,
subset = month == 9) # 9月のデータの抽出
tw_lm
```[1] "formula"
Call:
lm(formula = tw_formula, data = tw_data, subset = month == 9)
Coefficients:
(Intercept) solar press
408.1558 0.1414 -0.3798 推定に用いた変数の関係を図示する．
Code
```{r}
#| code-fold: true
#| fig-width: 6
#| fig-height: 6
tw_df <- model.frame(tw_lm) |> # 推定に用いたデータフレームの抽出(詳しくは後述)
as_tibble() # tibble クラスに変換
tw_df
ggpairs(tw_df) # 散布図
```# A tibble: 30 × 3
temp solar press
<dbl> <dbl> <dbl>
1 30.4 20.8 1009.
2 31.2 23.0 1007.
3 30.3 17.9 1008.
4 26.8 5.78 1011.
5 24.7 1.24 1004.
6 27.1 23.5 1007
7 27.4 19.0 1006.
8 29.6 21.5 1006
9 29.3 14 1008.
10 28.3 5.85 1010.
# ℹ 20 more rows
散布図と回帰式の定める平面の描画の3次元プロットを行う．
```{r}
s3d <- with(tw_df,
scatterplot3d(
x = solar, y = press, z = temp,
type = "p", # plotの種類: "p"点，"l"線，"h"足付き
pch = 16, # 点の種類 (?points 参照)
angle = 30, # xy平面の見る方向 (適宜調整)
highlight.3d = TRUE # 高さ(z)ごとに色を変える
))
s3d$plane3d(tw_lm, col = "blue", # 回帰式の定める平面の追加
draw_polygon = TRUE, # 平面の塗り潰しの設定
polygon_args = list(col = rgb(0,0,1,0.1)))
```
最小二乗推定量の性質
情報の取得
関数 stats::lm() の出力には様々な情報が含まれる．
base::summary(lmの出力) # 推定結果のまとめ
stats::coef(lmの出力) # 推定された回帰係数
stats::fitted(lmの出力) # あてはめ値
stats::resid(lmの出力) # 残差
stats::model.frame(lmの出力) # modelに必要な変数の抽出 (データフレーム)
stats::model.matrix(lmの出力) # デザイン行列関数 base::summary() の tidyverse 版として package::broom が用意されている．関数 stats::lm() の結果を tibble クラスで表示することができる．
library(broom) # 解析結果を tibble クラスに集約
broom::tidy(lmの出力) # 推定結果のまとめ．coef(summary(lmの出力)) と同様
broom::glance(lmの出力) # 評価指標(統計量)のまとめ．決定係数やF値などを整理
broom::augment(lmの出力) # 入出力データのまとめ．あてはめ値・残差などを整理ベクトルと行列の計算
データフレーム以外の重要なデータ構造として以下の2つがある．
- ベクトル (
vector) : 1次元の配列 - 行列 (
matrix) : 2次元の同じデータ型の配列
データフレーム(の一部)は 必要であればベクトルと行列に明示的に変換できる．
as.vector(データフレーム[列名]) # base::as.vector()
as_vector(データフレーム[列名]) # purrr::as_vector()
as.matrix(データフレーム) # base::as.matrix()
#' データフレームを行列に変換する場合は全て同じデータ型でなくてはならない代表的な行列とベクトルの計算は以下のようにする．
行列 \(A,B\) の積 \(AB\) および ベクトル \(a,b\) の内積 \(a\cdot b\)
A %*% B # 行列の大きさは適切である必要がある a %*% b # ベクトルは同じ長さである必要がある正方行列 \(A\) の逆行列 \(A^{-1}\)
solve(A) # 他にもいくつか関数はある\(X^{\mathsf{T}}Y\) および \(X^{\mathsf{T}}X\) の計算
crossprod(X, Y) # cross product の略 #' X: 行列 (またはベクトル) #' Y: 行列 (またはベクトル) crossprod(X) # 同じものを掛ける場合は引数は1つで良い行列の転置を計算する関数
t()と積を組み合わせてもよい．crossprod(X,Y) = t(X) %*% Y #' XY^T を計算する関数 tcrossprod() もある tcrossprod(X,Y) = X %*% t(Y)
問題
前問の推定結果を用いて，最小二乗推定量の以下の性質を確認しなさい．
推定された係数が正規方程式の解となること．
\[\begin{equation} \boldsymbol{\hat{\beta}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} \end{equation}\]
あてはめ値と残差が直交すること．
回帰式が標本平均を通ること．
解答欄
```{r}
#' ここにコードを書く
```解答例
広告費と売上データの回帰分析において， 回帰係数と正規方程式の解の一致は以下の計算で確認できる．
Code
```{r}
#| code-fold: true
beta <- coef(adv_lm3) # 推定された回帰係数
beta
X <- model.matrix(adv_lm3) # デザイン行列
Y <- model.frame(adv_lm3)[[1]] # 目的変数 (データフレームの1列目のベクトル)
solve(crossprod(X)) %*% crossprod(X, Y) # 正規方程式の解
```(Intercept) TV radio
2.92109991 0.04575482 0.18799423
[,1]
(Intercept) 2.92109991
TV 0.04575482
radio 0.18799423あてはめ値と残差の直交性はベクトルの内積を利用すればよい．
Code
```{r}
#| code-fold: true
yhat <- fitted(adv_lm3) # あてはめ値
ehat <- resid(adv_lm3) # 残差
yhat %*% ehat # 直交すれば内積はO(数値誤差は考慮する必要がある)
``` [,1]
[1,] -2.337233e-13回帰式が標本平均を通ることは説明変数の平均のあてはめ値を求めて比較すればよい．
Code
```{r}
#| code-fold: true
colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値
mean(Y) # 目的変数の標本平均
``` [,1]
[1,] 14.0225
[1] 14.0225東京の気候データにおいても同様に確認できる．
Code
```{r}
#| code-fold: true
beta <- coef(tw_lm) # 推定された回帰係数
beta
X <- model.matrix(tw_lm) # デザイン行列
X |> head(n = 10) # 先頭10行だけ表示
Y <- model.frame(tw_lm)[[1]] # 目的変数 (データフレームの1列目のベクトル)
Y |> head(n = 10) # 先頭10成分だけ表示
solve(crossprod(X)) %*% crossprod(X, Y) # 正規方程式の解(回帰係数と一致する)
```(Intercept) solar press
408.1557895 0.1413665 -0.3797500
(Intercept) solar press
244 1 20.80 1008.7
245 1 23.04 1006.7
246 1 17.92 1008.3
247 1 5.78 1010.8
248 1 1.24 1004.1
249 1 23.47 1007.0
250 1 18.98 1006.2
251 1 21.50 1006.0
252 1 14.00 1008.3
253 1 5.85 1009.8
[1] 30.4 31.2 30.3 26.8 24.7 27.1 27.4 29.6 29.3 28.3
[,1]
(Intercept) 408.1557895
solar 0.1413665
press -0.3797500Code
```{r}
#| code-fold: true
yhat <- fitted(tw_lm) # あてはめ値
ehat <- resid(tw_lm) # 残差
yhat %*% ehat # 直交すれば内積はO(数値誤差は考慮する必要がある)
``` [,1]
[1,] 4.090332e-14Code
```{r}
#| code-fold: true
colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値
mean(Y) # 目的変数の標本平均
``` [,1]
[1,] 26.48
[1] 26.48関数 stats::model.matrix() の返値は matrix クラスである．関数 stats::model.frame() の返値は data.frame クラスであり， model.frame(...)[[k]] は data.frame ではなく vector になる． tibble(...)[[k]] は vector となるが， tibble(...)[k] は vector ではなく1列の tibble になるので注意する． 関数 base::crossprod() の引数は厳密に vector/matrix であることを要請するが，例えば adv_data は全て数値なので行列に変換して計算できる．
crossprod(adv_data) # データフレームのままなので計算できない
crossprod(as.matrix(adv_data)) # 計算できる残差の分解
問題
前問の結果を用いて，残差の性質として以下の分解が成り立つことを確認しなさい．
\[\begin{equation} (\boldsymbol{y}-\bar{\boldsymbol{y}})^{\mathsf{T}} (\boldsymbol{y}-\bar{\boldsymbol{y}}) = (\boldsymbol{y}-\boldsymbol{\hat{y}})^{\mathsf{T}} (\boldsymbol{y}-\boldsymbol{\hat{y}})+ (\boldsymbol{\hat{y}}-\bar{\boldsymbol{y}})^{\mathsf{T}} (\boldsymbol{\hat{y}}-\bar{\boldsymbol{y}}) \end{equation}\] \[\begin{equation} S_y=S+S_r \end{equation}\]
解答欄
```{r}
#' ここにコードを書く
```解答例
広告費と売上データでは以下のように確認できる．
Code
```{r}
#| code-fold: true
summary(adv_lm3)
Y <- model.frame(adv_lm3)[[1]] # 目的変数の取得
(Sy <- sum((Y-mean(Y))^2)) # 目的変数のばらつき
(S <- sum(resid(adv_lm3)^2)) # 残差のばらつき
(Sr <- sum((fitted(adv_lm3)-mean(Y))^2)) # 回帰のばらつき
S+Sr # Sy と同じになっている
```
Call:
lm(formula = sales ~ TV + radio, data = adv_data)
Residuals:
Min 1Q Median 3Q Max
-8.7977 -0.8752 0.2422 1.1708 2.8328
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.92110 0.29449 9.919 <2e-16 ***
TV 0.04575 0.00139 32.909 <2e-16 ***
radio 0.18799 0.00804 23.382 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.681 on 197 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16
[1] 5417.149
[1] 556.914
[1] 4860.235
[1] 5417.149東京の気候データも同様に確認できる． 以下の例では目的変数を推定結果に含める lm() のオプションを利用している．
Code
```{r}
#| code-fold: true
tw_formula # モデルの確認
tw_lm <- lm(tw_formula,
data = tw_data,
subset = month == 9, # 9月のデータの抽出
y = TRUE) # 目的変数をyとして返すように指定
Y <- with(tw_lm, y) # 目的変数の取得 (tw_lm$y でも可)
(Sy <- sum((Y-mean(Y))^2)) # 目的変数のばらつき
(S <- sum(resid(tw_lm)^2)) # 残差のばらつき
(Sr <- sum((fitted(tw_lm)-mean(Y))^2)) # 回帰のばらつき
S+Sr # Sy と同じになっている
```temp ~ solar + press
[1] 201.508
[1] 94.82269
[1] 106.6853
[1] 201.508決定係数によるモデルの比較
以下では tidyverse の package::broom を利用して， 分析結果の情報を整理する方法を示す．
分析結果の情報の取得
関数 broom::tidy() は関数 base::summary() で得られる情報の一部 (coef(summary(lmの出力)) と同様の内容)を tibble クラスのデータフレームとして出力する．
tidy(x, conf.int = FALSE, conf.level = 0.95, exponentiate = FALSE, ...)
#' x: stats::lm() の返値
#' conf.int: 信頼区間を含むか否か(既定値は含めない)
#' conf.level: 信頼係数(既定値は0.95)
#' exponentiate: 係数を指数変換するか否か(既定値はしない)
#' 詳細は ?broom::tidy.lm関数 broom::glance() も同様に関数 base::summary() で得られる情報の一部 (決定係数やその他の統計量)を tibble クラスのデータフレームとして出力する．
glance(x, ...)
#' x: stats::lm() の返値
#' 詳細は ?broom::glance.lm関数 broom::augument() では 関数 stats::resid() や stats::fitted() で得られる情報を tibble クラスのデータフレームとして取得できる．
augment(
x,
data = model.frame(x),
newdata = NULL,
se_fit = FALSE,
interval = c("none", "confidence", "prediction"),
conf.level = 0.95,
...
)
#' x: stats::lm() の返値
#' data: 分析対象のデータフレーム(既定値ではmodel.frame(x)を利用)
#' newdata: 予測対象のデータフレーム
#' se_fit: 標準誤差を含むか否か(既定値は含めない)
#' interval: 信頼区間，予測区間を指定(既定値は区間推定を行わない)
#' conf.level: 信頼係数(既定値は0.95)
#' 詳細は ?broom::augment.lm実行例
決定係数を取得する．
```{r}
library(broom)
glance(bw_lm2) # 全ての指標を tibble クラスで取得
glance(bw_lm2)[["r.squared"]] # 決定係数
glance(bw_lm2)[["adj.r.squared"]] # 自由度調整済み決定係数
bw_lm2 |> # パイプ演算子を用いた方法
glance() |>
pull(adj.r.squared) # 列をベクトルとして取得
bind_rows( # 2つのモデルの評価指標を結合
glance(bw_lm1),
glance(bw_lm2)) |>
select(r.squared,adj.r.squared) |> # 必要な2列を選択
mutate(name = c("モデル1","モデル2"), # モデルの名称
.before = 1) # 先頭に付加
#' 列番号がわかっている場合は select(1:2) でも良い
```# 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.828 0.796 0.287 26.4 0.0000000406 4 -1.80 15.6 23.4
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
[1] 0.8275278
[1] 0.7961692
[1] 0.7961692
# A tibble: 2 × 3
name r.squared adj.r.squared
<chr> <dbl> <dbl>
1 モデル1 0.257 0.195
2 モデル2 0.828 0.796問題
決定係数を用いてモデルの比較を行いなさい．
広告費と売上データ
sales ~ TV sales ~ radio sales ~ TV + radio東京の9月の気候データ
temp ~ solar temp ~ solar + press temp ~ solar + press + cloud
解答欄
```{r}
#' ここにコードを書く
```解答例
広告費と売上データでのモデルの比較は以下のとおりである．
Code
```{r}
#| code-fold: true
adv_formula1 <- sales ~ TV # モデル式を定義
adv_formula2 <- sales ~ radio
adv_formula3 <- sales ~ TV + radio
adv_lm1 <- lm(adv_formula1, data = adv_data) # 各モデルで分析
adv_lm2 <- lm(adv_formula2, data = adv_data)
adv_lm3 <- lm(adv_formula3, data = adv_data)
glance(adv_lm1)[["adj.r.squared"]] # 自由度調整済み決定係数
glance(adv_lm2)[["adj.r.squared"]] # model1より減少
glance(adv_lm3)[["adj.r.squared"]] # model1より上昇
```[1] 0.6099148
[1] 0.3286589
[1] 0.8961505予測値と実測値の比較は以下のようにすれば良い．
Code
```{r}
#| code-fold: true
adv_data |>
mutate(model1 = fitted(adv_lm1), # モデルごとに予測値を追加
model2 = fitted(adv_lm2),
model3 = fitted(adv_lm3)) |>
pivot_longer(starts_with("model"), # モデルをラベルとして予測値をまとめる
names_to = "model", # モデルのラベルを model 列
values_to = "fitted") |> # あてはめ値を fitted 列
ggplot(aes(x = sales, y = fitted)) + # 実測値をx軸，予測値をy軸で表示
geom_abline(slope = 1, intercept = 0, colour = "red") + # 基準線
geom_point(aes(colour = model, shape = model)) + # 予測値をモデル別に表示
labs(y = "fitted values")
```
東京の気候データも同様で，モデルの比較は以下のとおりである．
Code
```{r}
#| code-fold: true
tw_subset <- tw_data |> filter(month == 9) # 9月のデータの抽出
tw_formula1 <- temp ~ solar # モデル式を定義
tw_formula2 <- temp ~ solar + press
tw_formula3 <- temp ~ solar + press + cloud
tw_lm1 <- lm(tw_formula1, data = tw_subset) # 回帰分析
tw_lm2 <- lm(tw_formula2, data = tw_subset)
tw_lm3 <- lm(tw_formula3, data = tw_subset)
bind_rows( # 3つのモデルの評価指標を結合
glance(tw_lm1),
glance(tw_lm2),
glance(tw_lm3)) |>
select(r.squared,adj.r.squared) |> # 2種類の決定係数の列を抽出
mutate(name = paste0("モデル",1:3), # モデル名の列
.before = 1) # 1列目の前に付加
```# A tibble: 3 × 3
name r.squared adj.r.squared
<chr> <dbl> <dbl>
1 モデル1 0.0867 0.0541
2 モデル2 0.529 0.495
3 モデル3 0.556 0.505 比較表を作成するには，推定されたモデルの情報を取得して tibble 形式のデータフレームを作成する broom::glance() と package::gt を組み合わせて，例えば次のようにすれば良い．
Code
```{r}
#| code-fold: true
library(gt) # 表作成のためのパッケージ
bind_rows(glance(tw_lm1), # 3つのモデルの評価指標を結合
glance(tw_lm2),
glance(tw_lm3)) |>
select(1:2) |> # 決定係数と自由度調整済み決定係数を表示
mutate(model = map_vec(c(tw_formula1,
tw_formula2,
tw_formula3),
deparse),
.before = 1) |> # モデル式を文字列に変換して先頭に model 列として置く
set_names(c("モデル式","決定係数","自由度調整済み決定係数")) |>
gt() |> # 表に整形
fmt_number(decimals = 3) # 表示を3桁に
```| モデル式 | 決定係数 | 自由度調整済み決定係数 |
|---|---|---|
| temp ~ solar | 0.087 | 0.054 |
| temp ~ solar + press | 0.529 | 0.495 |
| temp ~ solar + press + cloud | 0.556 | 0.505 |
予測値と実測値の比較は以下のようになる．5
5 プロットのための作成する tibble クラスの列名を日本語にする場合の例となっている．オブジェクトは日本語(unicode)でも操作できるので，コードの可搬性・可読性を鑑みて適宜利用すると良い．
Code
```{r}
#| code-fold: true
tw_subset |>
mutate(model1 = fitted(tw_lm1), # モデルごとにあてはめ値をデータフレームに追加
model2 = fitted(tw_lm2),
model3 = fitted(tw_lm3)) |>
pivot_longer(starts_with("model"), # モデルをラベルとしてあてはめ値をまとめる
names_to = "モデル",
values_to = "あてはめ値") |>
ggplot(aes(x = temp, y = あてはめ値)) + # 気温の実測値をx軸，あてはめ値をy軸で表示
geom_abline(slope = 1,
intercept = 0,
colour = "red") + # 基準線
geom_point(aes(colour = モデル, shape = モデル)) + # あてはめ値をモデル別に表示
labs(x = "実際の気温")
```