library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse) 第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 を用いた ボルドーワインの価格と気候の関係に関する回帰分析は以下で実行できる．
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 rowsbw_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 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 問題 (広告費と売上データ)
データセット https://www.statlearning.com/s/Advertising.csv は 広告費(TV,radio,newspapers)と売上の関係を調べたもの 2 である． このデータセットを用いて以下の回帰式を推定しなさい．
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解答欄
解答例
データの読み込みを行う．
adv_data <- read_csv("https://www.statlearning.com/s/Advertising.csv")インターネット上にあるデータをローカルに保存しておくには， 例えば以下のようにすればよい． 2回目以降はローカルファイルを参照することになる．
adv_url <- "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_url) # リモートファイルを読み込む
write_csv(adv_data, file = adv_file) # ローカルファイルに保存する
}3 標準設定では図は横長で描画される．この例では縦横が同じ大きさ(6inch)になるように指定している．
GGally パッケージを利用して散布図を描く．3
library(GGally)
ggpairs(adv_data, columns = 2:5)
“TV” の宣伝費で売上を説明する．
adv_lm1 <- lm(sales ~ TV, data = adv_data)
adv_lm1
Call:
lm(formula = sales ~ TV, data = adv_data)
Coefficients:
(Intercept) TV
7.03259 0.04754 #' ()は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)
“radio” の宣伝費で売上を説明する．
adv_lm2 <- lm(sales ~ radio, data = adv_data)
adv_lm2
Call:
lm(formula = sales ~ radio, data = adv_data)
Coefficients:
(Intercept) radio
9.3116 0.2025 adv_data |>
ggplot(aes(x = radio, y = sales)) +
geom_point(colour = "orange") +
geom_smooth(method = lm, se = FALSE)
両者の宣伝費で売上を説明する．
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次元の散布図を描くには package::scatterplot3d を利用することができる．
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解答欄
解答例
データの読み込みを行う．
tw_data <- read_csv("data/tokyo_weather.csv")9月の”気温”を目的変数，“日射量・気圧”を説明変数とするモデルの推定を行う．
tw_formula <- temp ~ solar + press # モデル式の定義
class(tw_formula) # formula クラスであることを確認[1] "formula"tw_lm <- lm(tw_formula, # 回帰係数の推定
data = tw_data,
subset = month == 9) # 9月のデータの抽出
tw_lm
Call:
lm(formula = tw_formula, data = tw_data, subset = month == 9)
Coefficients:
(Intercept) solar press
386.2285 0.3020 -0.3602 tw_df <- model.frame(tw_lm) |> # 推定に用いたデータフレームの抽出(詳しくは後述)
as_tibble() # tibble クラスに変換
tw_df# A tibble: 30 × 3
temp solar press
<dbl> <dbl> <dbl>
1 29.2 24.0 1012.
2 29.6 22.1 1010.
3 29.1 18.6 1011.
4 26.1 7.48 1008.
5 29.3 22.6 1005.
6 27.5 13.2 1004.
7 27 11.0 1008.
8 21.9 2.1 1008.
9 24.8 8.81 1007.
10 27.8 17.6 1009.
# ℹ 20 more rowsggpairs(tw_df) # 散布図
散布図と回帰式の定める平面の描画の3次元プロットを行う．
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 R をでは以下に示す関数で，個別に情報を取り出すことができる．
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 = (X^TY)^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}
- あてはめ値と残差が直交すること．
- 回帰式が標本平均を通ること．
解答欄
解答例
広告費と売上データの回帰分析において， 回帰係数と正規方程式の解の一致は以下の計算で確認できる．
beta <- coef(adv_lm3) # 推定された回帰係数
beta(Intercept) TV radio
2.92109991 0.04575482 0.18799423 X <- model.matrix(adv_lm3) # デザイン行列
Y <- model.frame(adv_lm3)[[1]] # 目的変数 (データフレームの1列目のベクトル)
solve(crossprod(X)) %*% crossprod(X, Y) # 正規方程式の解 [,1]
(Intercept) 2.92109991
TV 0.04575482
radio 0.18799423あてはめ値と残差の直交性はベクトルの内積を利用すればよい．
yhat <- fitted(adv_lm3) # あてはめ値
ehat <- resid(adv_lm3) # 残差
yhat %*% ehat # 直交すれば内積はO(数値誤差は考慮する必要がある) [,1]
[1,] -2.337233e-13回帰式が標本平均を通ることは説明変数の平均のあてはめ値を求めて比較すればよい．
colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値 [,1]
[1,] 14.0225mean(Y) # 目的変数の標本平均 [1] 14.0225東京の気候データにおいても同様に確認できる．
beta <- coef(tw_lm) # 推定された回帰係数
beta(Intercept) solar press
386.2284847 0.3019861 -0.3602392 X <- model.matrix(tw_lm) # デザイン行列
X |> head(n = 10) # 先頭10行だけ表示 (Intercept) solar press
244 1 24.01 1012.1
245 1 22.07 1010.3
246 1 18.64 1010.6
247 1 7.48 1007.5
248 1 22.58 1005.2
249 1 13.17 1003.6
250 1 11.01 1007.9
251 1 2.10 1007.8
252 1 8.81 1006.8
253 1 17.57 1009.1Y <- model.frame(tw_lm)[[1]] # 目的変数 (データフレームの1列目のベクトル)
Y |> head(n = 10) # 先頭10成分だけ表示 [1] 29.2 29.6 29.1 26.1 29.3 27.5 27.0 21.9 24.8 27.8solve(crossprod(X)) %*% crossprod(X, Y) # 正規方程式の解(回帰係数と一致する) [,1]
(Intercept) 386.2284847
solar 0.3019861
press -0.3602392yhat <- fitted(tw_lm) # あてはめ値
ehat <- resid(tw_lm) # 残差
yhat %*% ehat # 直交すれば内積はO(数値誤差は考慮する必要がある) [,1]
[1,] 1.123776e-14colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値 [,1]
[1,] 26.69333mean(Y) # 目的変数の標本平均 [1] 26.69333関数 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}
解答欄
解答例
広告費と売上データでは以下のように確認できる．
summary(adv_lm3) # 用いるモデル3の概要を表示
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-16Y <- model.frame(adv_lm3)[[1]] # 目的変数の取得
Sy <- sum((Y-mean(Y))^2) # 目的変数のばらつき
Sy[1] 5417.149S <- sum(resid(adv_lm3)^2) # 残差のばらつき
S[1] 556.914Sr <- sum((fitted(adv_lm3)-mean(Y))^2) # 回帰のばらつき
Sr[1] 4860.235S+Sr # Sy = S+Sr であることを確認[1] 5417.149東京の気候データも同様に確認できる． 以下の例では目的変数を推定結果に含める lm() のオプションを利用している．
tw_formula # モデルの確認temp ~ solar + presstw_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) # 目的変数のばらつき
Sy[1] 151.0787S <- sum(resid(tw_lm)^2) # 残差のばらつき
S[1] 55.66168Sr <- sum((fitted(tw_lm)-mean(Y))^2) # 回帰のばらつき
Sr[1] 95.41699S+Sr # Sy = S+Sr であることを確認[1] 151.0787決定係数によるモデルの比較
以下では tidyverse の package::broom を利用して， 分析結果の情報を整理する方法を示す
分析結果の情報の取得 (
tidyverse)
関数 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::augumet() では 関数 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
実行例
決定係数を取得する．
library(broom)
glance(bw_lm2) # 全ての指標を tibble クラで取得# 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>glance(bw_lm2)[["r.squared"]] # 決定係数[1] 0.8275278glance(bw_lm2)[["adj.r.squared"]] # 自由度調整済み決定係数[1] 0.7961692bw_lm2 |> # パイプ演算子を用いた方法
glance() |>
pull(adj.r.squared) # 列をベクトルとして取得[1] 0.7961692bind_rows( # 3つのモデルの評価指標を結合
glance(bw_lm1),
glance(bw_lm2)) |>
select(r.squared,adj.r.squared) |> # 必要な3列を選択
mutate(name = c("モデル1","モデル2"), # モデルの名称
.before = 1) # 先頭に付加# 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#' 列番号がわかっている場合は select(1:2) でも良い問題
それぞれのデータに対して，決定係数を用いて以下のモデルの比較を行いなさい．
広告費と売上データ
sales ~ TV sales ~ radio sales ~ TV + radio東京の9月の気候データ
temp ~ solar temp ~ solar + press temp ~ solar + press + cloud
解答欄
解答例
広告費と売上データでのモデルの比較は以下のとおりである．
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"]] # 自由度調整済み決定係数[1] 0.6099148glance(adv_lm2)[["adj.r.squared"]] # model1より減少[1] 0.3286589glance(adv_lm3)[["adj.r.squared"]] # model1より上昇[1] 0.8961505予測値と実測値の比較は以下のようにすれば良い．
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")
東京の気候データも同様で，モデルの比較は以下のとおりである．
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.414 0.393
2 モデル2 0.632 0.604
3 モデル3 0.633 0.591比較表を作成するには，推定されたモデルの情報を取得して tibble クラスのデータフレームを作成する broom::glance() と package::gt を組み合わせて，例えば次のようにすれば良い．
library(gt) # 表作成のためのパッケージ
bind_rows( # 3つのモデルの評価指標を結合
glance(tw_lm1),
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.414 | 0.393 |
| temp ~ solar + press | 0.632 | 0.604 |
| temp ~ solar + press + cloud | 0.633 | 0.591 |
予測値と実測値の比較は以下のようになる．
tw_subset |>
mutate(model1 = fitted(tw_lm1), # モデルごとに予測値をデータフレームに追加
model2 = fitted(tw_lm2),
model3 = fitted(tw_lm3)) |>
pivot_longer(starts_with("model"), # モデルをラベルとして予測値をまとめる
names_to = "model",
values_to = "fitted") |>
ggplot(aes(x = temp, y = fitted)) + # 気温の実測値をx軸，予測値をy軸で表示
geom_abline(slope = 1, intercept = 0, colour = "red") + # 基準線
geom_point(aes(colour = model, shape = model)) + # 予測値をモデル別に表示
labs(x = "temperature", y = "fitted values")