# 準備 以下で利用する共通パッケージを読み込む．[^1] [^1]: : 必要なパッケージは適宜追加して良い． ```{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 } ``` # 回帰係数の推定 ## 線形回帰のための関数 回帰係数は関数 `stats::lm()` を用いて推定することができる． ``` r 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 # データの表示(一部) ``` ```{r} bw_lm1 <- lm(formula = LPRICE2 ~ WRAIN + HRAIN, # 冬と収穫期の雨で価格を予測 data = bw_data) bw_lm1 # 推定結果の簡単な表示 ``` ```{r} bw_lm2 <- lm(formula = LPRICE2 ~ . - VINT, # VINTを除く全てで価格を予測 data = bw_data) bw_lm2 # 推定結果の簡単な表示 ``` ::: column-margin パイプ (`|>`) を利用して， "データフレーム \|\> 操作" という形で統一的に書きたいのであれば， `_` (プレースフォルダ; placeholder) を利用して以下のように書くことができる． ``` r bw_fit <- bw_data |> lm(formula = LPRICE2 ~ . - VINT, data = _) ``` ::: ## 問題 (広告費と売上データ) データセット [^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. [参考](https://www.statlearning.com) ``` r formula = sales ~ TV formula = sales ~ radio formula = sales ~ TV + radio ``` ## 解答欄 ```{r} #' ここにコードを書く ``` ## 解答例 データの読み込みを行う． ``` r adv_data <- read_csv("https://www.statlearning.com/s/Advertising.csv") ``` ::: column-margin インターネット上にあるデータをローカルに保存しておくには， 例えば以下のようにすればよい． 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) # ローカルファイルに保存する } ``` ::: GGally パッケージを利用して散布図を描く．[^3] [^3]: : 標準設定では図は横長で描画される．この例では縦横が同じ大きさ(7inch)になるように指定している． ```{r} #| code-fold: true #| fig-width: 6 #| fig-height: 6 library(GGally) ggpairs(adv_data, columns = 2:5) ``` "TV" の宣伝費で売上を説明する． ```{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) ``` "radio" の宣伝費で売上を説明する． ```{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) ``` 両者の宣伝費で売上を説明する． ```{r} #| code-fold: true adv_lm3 <- lm(sales ~ TV + radio, data = adv_data) adv_lm3 ``` ::: column-margin 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]: : [気象庁のサイト](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php) ``` r formula = temp ~ solar + press ``` ## 解答欄 ```{r} #' ここにコードを書く ``` ## 解答例 データの読み込みを行う． ```{r} #| code-fold: true tw_data <- read_csv("data/tokyo_weather.csv") ``` 9月の"気温"を目的変数，"日射量・気圧"を説明変数とするモデルの推定を行う． ```{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 ``` 推定に用いた変数の関係を図示する． ```{r} #| code-fold: true #| fig-width: 6 #| fig-height: 6 tw_df <- model.frame(tw_lm) |> # 推定に用いたデータフレームの抽出(詳しくは後述) as_tibble() # tibble クラスに変換 tw_df ggpairs(tw_df) # 散布図 ``` ::: column-margin 散布図と回帰式の定める平面の描画の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()` の出力には様々な情報が含まれる． ``` 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` クラスで表示することができる． ``` r library(broom) # 解析結果を tibble クラスに集約 broom::tidy(lmの出力) # 推定結果のまとめ．coef(summary(lmの出力)) と同様 broom::glance(lmの出力) # 評価指標(統計量)のまとめ．決定係数やF値などを整理 broom::augment(lmの出力) # 入出力データのまとめ．あてはめ値・残差などを整理 ``` ## ベクトルと行列の計算 データフレーム以外の重要なデータ構造として以下の2つがある． - ベクトル (`vector`) : 1次元の配列 - 行列 (`matrix`) : 2次元の同じデータ型の配列 データフレーム(の一部)は 必要であればベクトルと行列に明示的に変換できる． ``` r as.vector(データフレーム[列名]) # base::as.vector() as_vector(データフレーム[列名]) # purrr::as_vector() as.matrix(データフレーム) # base::as.matrix() #' データフレームを行列に変換する場合は全て同じデータ型でなくてはならない ``` 代表的な行列とベクトルの計算は以下のようにする． - 行列 $A,B$ の積 $AB$ および ベクトル $a,b$ の内積 $a\cdot b$ ``` r A %*% B # 行列の大きさは適切である必要がある a %*% b # ベクトルは同じ長さである必要がある ``` - 正方行列 $A$ の逆行列 $A^{-1}$ ``` r solve(A) # 他にもいくつか関数はある ``` - $X^{\mathsf{T}}Y$ および $X^{\mathsf{T}}X$ の計算 ``` r crossprod(X, Y) # cross product の略 #' X: 行列 (またはベクトル) #' Y: 行列 (またはベクトル) crossprod(X) # 同じものを掛ける場合は引数は1つで良い ``` 行列の転置を計算する関数 `t()` と積を組み合わせてもよい． ``` r 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} #' ここにコードを書く ``` ## 解答例 広告費と売上データの回帰分析において， 回帰係数と正規方程式の解の一致は以下の計算で確認できる． ```{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) # 正規方程式の解 ``` あてはめ値と残差の直交性はベクトルの内積を利用すればよい． ```{r} #| code-fold: true yhat <- fitted(adv_lm3) # あてはめ値 ehat <- resid(adv_lm3) # 残差 yhat %*% ehat # 直交すれば内積はO(数値誤差は考慮する必要がある) ``` 回帰式が標本平均を通ることは説明変数の平均のあてはめ値を求めて比較すればよい． ```{r} #| code-fold: true colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値 mean(Y) # 目的変数の標本平均 ``` 東京の気候データにおいても同様に確認できる． ```{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) # 正規方程式の解(回帰係数と一致する) ``` ```{r} #| code-fold: true yhat <- fitted(tw_lm) # あてはめ値 ehat <- resid(tw_lm) # 残差 yhat %*% ehat # 直交すれば内積はO(数値誤差は考慮する必要がある) ``` ```{r} #| code-fold: true colMeans(X) %*% beta # 説明変数の標本平均のあてはめ値 mean(Y) # 目的変数の標本平均 ``` ::: column-margin 関数 `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` は全て数値なので行列に変換して計算できる． ``` r 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} #' ここにコードを書く ``` ## 解答例 広告費と売上データでは以下のように確認できる． ```{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 と同じになっている ``` 東京の気候データも同様に確認できる． 以下の例では目的変数を推定結果に含める `lm()` のオプションを利用している． ```{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 と同じになっている ``` # 決定係数によるモデルの比較 以下では `tidyverse` の `package::broom` を利用して， 分析結果の情報を整理する方法を示す． ## 分析結果の情報の取得 関数 `broom::tidy()` は関数 `base::summary()` で得られる情報の一部 (`coef(summary(lmの出力))` と同様の内容)を `tibble` クラスのデータフレームとして出力する． ``` r 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` クラスのデータフレームとして出力する． ``` r glance(x, ...) #' x: stats::lm() の返値 #' 詳細は ?broom::glance.lm ``` 関数 `broom::augument()` では 関数 `stats::resid()` や `stats::fitted()` で得られる情報を `tibble` クラスのデータフレームとして取得できる． ``` r 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) でも良い ``` ## 問題 決定係数を用いてモデルの比較を行いなさい． - 広告費と売上データ ``` r sales ~ TV sales ~ radio sales ~ TV + radio ``` - 東京の9月の気候データ ``` r temp ~ solar temp ~ solar + press temp ~ solar + press + cloud ``` ## 解答欄 ```{r} #' ここにコードを書く ``` ## 解答例 広告費と売上データでのモデルの比較は以下のとおりである． ```{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より上昇 ``` 予測値と実測値の比較は以下のようにすれば良い． ```{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") ``` 東京の気候データも同様で，モデルの比較は以下のとおりである． ```{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列目の前に付加 ``` 比較表を作成するには，推定されたモデルの情報を取得して `tibble` 形式のデータフレームを作成する `broom::glance()` と `package::gt` を組み合わせて，例えば次のようにすれば良い． ```{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桁に ``` 予測値と実測値の比較は以下のようになる．[^5] [^5]: プロットのための作成する `tibble` クラスの列名を日本語にする場合の例となっている．オブジェクトは日本語(unicode)でも操作できるので，コードの可搬性・可読性を鑑みて適宜利用すると良い． ```{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 = "実際の気温") ```