--- title: "統計データ解析" subtitle: "第8講 実習" date: "`r Sys.Date()`" format: html: toc: true html-math-method: katex self-contained: true grid: margin-width: 350px execute: echo: true warning: false reference-location: margin citation-location: margin tbl-cap-location: margin fig-cap-location: margin editor: visual editor_options: chunk_output_type: console --- ## 準備 以下で利用する共通パッケージを読み込む． ```{r} library(conflicted) # 関数名の衝突を警告 conflicts_prefer( # 優先的に使う関数を指定 dplyr::filter(), dplyr::select(), dplyr::lag(), ) library(tidyverse) library(MASS) library(broom) # 解析結果を tibble 形式に集約 library(gt) # 表の作成 library(gtsummary) # 分析結果の表の作成 library(ggfortify) # biplot表示のため ``` ## 線形判別 ### 問題 東京の気候データを用いて以下の分析を行いなさい． - 9月と10月の気温と湿度のデータを抽出する． ``` r tw_data <- read_csv("data/tokyo_weather.csv") tw_subset <- tw_data |> filter(month %in% c(9,10)) |> select(temp, humid, month) |> mutate(month = as_factor(month)) # 月を因子化 ``` - 半分のデータを用いて線形判別関数を構成し，残りのデータを用いて判別を行う． ``` r idx <- seq(2, 60, by = 2) tw_train <- tw_subset[ idx,] # 訓練データ tw_test <- tw_subset[-idx,] # 試験データ tw_lda <- lda(month ~ temp + humid, data = tw_train) # 線形判別関数の構成 tw_lda_fitted <- predict(tw_lda) # 判別関数によるクラス分類結果の取得 tw_lda_predict <- predict(tw_lda, newdata = tw_test) # 新しいデータの予測 ``` ```{r} ``` ### 解答例 必要な気象データを整理する． ```{r} tw_data <- read_csv("data/tokyo_weather.csv") tw_subset <- tw_data |> filter(month %in% c(9,10)) |> # 9,10月のデータ select(temp, humid, month) |> # 気温・湿度・月を選択 mutate(month = as_factor(month)) # 月を因子化 idx <- seq(1, nrow(tw_subset), by = 2) # データの分割(1つおき) tw_train <- tw_subset[ idx,] # 訓練データ(推定に用いる) tw_test <- tw_subset[-idx,] # 試験データ(評価に用いる) ``` ::: callout-tip ランダムに分割するには例えば以下のようにすれば良い． ``` r n <- nrow(tw_subset); idx <- sample(1:n, n/2) ``` ::: ::: callout-tip 関数 `lubridate::month()` を用いれば月を文字列(因子)とすることもできる． ``` r mutate(month = month(month, label = TRUE)) ``` 以下のようにして最初に月を因子化することもできるが，関数によっていは処理で使われない因子(例えば1月)も表示してしまうので注意が必要である． ```{r} read_csv("data/tokyo_weather.csv") |> mutate(month = as_factor(month)) ``` また複数の項目を因子化するには例えば以下のようにすればよい． ``` r mutate(across(c(year, month, day), as_factor)) ``` ::: 気温と湿度の散布図を作成して，データを視覚化する． ```{r} tw_subset |> # ggplot(aes(x = temp, y = humid)) + geom_point(aes(colour = month)) + # 月ごとに点の色を変える labs(x = "Temperature", y = "Humidity", title = "September & October") ``` 訓練データを用いて線形判別関数(等分散性を仮定)を作成する． ```{r} tw_model <- month ~ temp + humid tw_lda <- lda(formula = tw_model, data = tw_train) ``` 判別関数によるクラス分類結果を確認する． ```{r} tw_lda_fitted <- predict(tw_lda) # table(true = tw_train[["month"]], # 真値 predict = tw_lda_fitted[["class"]]) # 予測値の比較(混同行列) tibble(true = tw_train[["month"]]) |> mutate(predict = tw_lda_fitted[["class"]]) |> View() # 真値と予測値の対応表を作成して直接比較 as_tibble(tw_lda_fitted[["x"]]) |> # リストから判別関数値(配列)を取得 mutate(month = tw_train[["month"]]) |> # 月(真値)を追加 ggplot(aes(x = LD1, fill = month)) + geom_histogram(aes(y = after_stat(density)), bins = 30) + facet_grid(month ~ .) # 月ごとに表示(y方向に並べる) ``` 試験データによる評価を行う ```{r} tw_lda_predict <- predict(tw_lda, newdata = tw_test) table(true = tw_test[["month"]], # 混同行列 predict = tw_lda_predict[["class"]]) tibble(true = tw_test[["month"]]) |> # 比較表 mutate(predict = tw_lda_predict[["class"]]) |> View() as_tibble(tw_lda_predict[["x"]]) |> # 判別関数値の視覚化 mutate(month = tw_test[["month"]]) |> ggplot(aes(x = LD1, fill = month)) + geom_histogram(aes(y = after_stat(density)), bins = 30) + facet_grid(month ~ .) ``` 推定された線形判別関数の識別境界を図示する． ```{r} range_x <- range(tw_subset[["temp"]]) # 気温の値域 range_y <- range(tw_subset[["humid"]]) # 湿度の値域 grid_x <- pretty(range_x, 100) # 気温の値域の格子点を作成 grid_y <- pretty(range_y, 100) # 湿度の値域の格子点を作成 grid_xy <- expand.grid(temp = grid_x, humid = grid_y) # 2次元の格子点を作成 tw_lda_grid <- predict(tw_lda, # 格子点上の判別関数値を計算 newdata = grid_xy) #' 判別関数により予測されるラベルの図示 p <- as_tibble(grid_xy) |> mutate(predict = tw_lda_grid[["class"]]) |> ggplot(aes(x = temp, y = humid)) + geom_tile(aes(fill = predict), alpha = 0.3) + labs(x = "Temperature", y = "Humidity", title = "Linear Discriminant Analysis") print(p) # 判別境界の表示 p + # データ点の重ね描き geom_point(data = tw_subset, aes(x = temp, y = humid, colour = month)) ``` 判別関数値の図示 (参考) ```{r} p2 <- as_tibble(grid_xy) |> mutate(LD = tw_lda_grid[["x"]][,"LD1"]) |> # LD1列のみ利用 ggplot(aes(x = temp, y = humid)) + geom_raster(aes(fill = LD), alpha = 0.5) + scale_fill_gradientn(colours=c("red","white","blue")) + labs(x = "Temperature", y = "Humidity", title = "Scores of Discriminant Variables") print(p2) # 判別関数値の表示 p2 + # データ点の重ね描き geom_point(data = tw_subset, aes(x = temp, y = humid, colour = month)) ``` ## 2次判別 ### 問題 - 東京の気候データを用いて以下の分析を行いなさい - 前問と同様な設定で2次判別を行いなさい #+begin_src R :eval no :tangle no tw_qda <- qda(month ~ temp + humid, data = tw_train) # 2次判別関数の構成 tw_qda_fitted <- predict(tw_qda) # 判別関数によるクラス分類結果の取得 tw_qda_predict <- predict(tw_qda, newdata = tw_test) # 新しいデータの予測 #+end_src - 別の月や変数を用いて判別分析を行いなさい ```{r} ``` ### 解答例 #+begin_src R :exports none #' --------------------------------------------------------------------------- #' @practice 2次判別 #+end_src #+begin_src R :eval no :exports none #' 東京の気象データによる判別分析 #' 前の練習問題で作成した 'tw_train', 'tw_test' を利用 #' 訓練データで2次判別関数(等分散性を仮定しない)を作成 #' 線形判別関数と同じモデル式を用いる tw_qda <- qda(formula = tw_model, data = tw_train) tw_qda_fitted <- predict(tw_qda) # 判別関数によるクラス分類結果の取得 table(true = tw_train[["month"]], # 真値 predict = tw_qda_fitted[["class"]]) # 予測値の比較(混同行列) tibble(true = tw_train[["month"]]) |> mutate(predict = tw_qda_fitted[["class"]]) |> View() # 真値と予測値の対応表を作成して直接比較 #' 試験データによる評価 tw_qda_predict <- predict(tw_qda, newdata = tw_test) table(true = tw_test[["month"]], # 混同行列 predict = tw_qda_predict[["class"]]) tibble(true = tw_test[["month"]]) |> # 比較表 mutate(predict = tw_qda_predict[["class"]]) |> View() #' 推定された2次判別関数の図示 tw_qda_grid <- predict(tw_qda, # 格子点上の判別関数値を計算 newdata = grid_xy) #' 判別関数により予測されるラベルの図示 p <- as_tibble(grid_xy) |> mutate(predict = tw_qda_grid[["class"]]) |> ggplot(aes(x = temp, y = humid)) + geom_tile(aes(fill = predict), alpha = 0.3) + labs(x = "Temperature", y = "Humidity", title = "Quadratic Discriminant Analysis") print(p) # 判別境界の表示 p + # データ点の重ね描き geom_point(data = tw_subset, aes(x = temp, y = humid, colour = month)) #+end_src #+begin_src R :exports none #' --------------------------------------------------------------------------- #+end_src ## 多値判別 ### 問題 - 東京の気候データを用いて以下の分析を行いなさい - 9月，10月，11月の気温と湿度のデータを用いて判別関数を作成しなさい． #+begin_src R :eval no :tangle no tw_subset <- tw_data |> filter(month %in% c(9,10,11)) |> select(temp, humid, month) |> mutate(month = as_factor(month)) #+end_src - 別の月や変数を用いて判別分析を行いなさい #+begin_src R :eval no :tangle no #' 雨の有無を識別する例 tw_subset2 <- tw_data |> mutate(rain = factor(rain > 0), month = as_factor(month)) |> # 雨の有無でラベル化する select(rain, temp, solar, wind, month) tw_lda2 <- lda(rain ~ ., data = tw_subset2) # 'rain' をそれ以外で判別 #+end_src ### 解答例 #+begin_src R :exports none #' --------------------------------------------------------------------------- #' @practice 多値判別 #+end_src #+begin_src R :eval no :exports none #' 東京の気象データによる判別分析 #' 3値判別のためのデータの整理 tw_subset3 <- tw_data |> filter(month %in% c(9,10,11)) |> select(temp, humid, month) |> mutate(month = as_factor(month)) #' 線形判別関数(3値)を作成 tw_lda3 <- lda(formula = tw_model, data = tw_subset3) tw_lda3_fitted <- predict(tw_lda3) # 判別関数によるクラス分類結果の取得 table(true = tw_subset3[["month"]], # 真値 predict = tw_lda3_fitted[["class"]]) # 予測値の比較(混同行列) #' 推定された線形判別関数(3値)により予測されるラベルの図示 #' 格子点は再計算する必要がある range_x <- range(tw_subset3[["temp"]]) # 気温の値域 range_y <- range(tw_subset3[["humid"]]) # 湿度の値域 grid_x <- pretty(range_x, 100) # 気温の値域の格子点を作成 grid_y <- pretty(range_y, 100) # 湿度の値域の格子点を作成 grid_xy <- expand.grid(temp = grid_x, humid = grid_y) # 2次元の格子点を作成 tw_lda3_grid <- predict(tw_lda3, # 格子点上の判別関数値を計算 newdata = grid_xy) as_tibble(grid_xy) |> mutate(predict = tw_lda3_grid[["class"]]) |> ggplot(aes(x = temp, y = humid)) + geom_tile(aes(fill = predict), alpha = 0.3) + geom_point(data = tw_subset3, aes(x = temp, y = humid, colour = month)) + labs(x = "Temperature", y = "Humidity", title = "Multi-label Discriminant Analysis") #' 3値判別の場合には2つの判別関数を構成するので #' これを用いてデータ点の散布図を作成することができる p <- bind_cols(tw_lda3_fitted[["x"]], # 判別関数値(LD1,LD2) tw_subset3["month"]) |> # データフレームとして抽出 ggplot(aes(x = LD1, y = LD2)) + geom_point(aes(colour = month)) # 月ごとに色を変更 print(p) #' 関数geom_tile()が座標軸に沿った格子点を想定しているため #' 判別関数値の散布図上で判別境界とデータ点を表示するには工夫が必要 #' ここでは関数geom_point()で代用した簡便な例を示す p + geom_point(data = # 判別関数値と予測ラベルのデータフレームを作成 bind_cols(tw_lda3_grid[["x"]], predict = tw_lda3_grid[["class"]]), aes(colour = predict), alpha = 0.1) #' 12ヶ月分のデータを用いる #' 数が多いのでサンプリングする idx <- sample(nrow(tw_data), 100) tw_subset12 <- slice(tw_data, idx) |> mutate(month = as_factor(month)) tw_lda12 <- lda(month ~ temp + solar + wind + humid, data = tw_subset12) tw_lda12_fitted <- predict(tw_lda12) table(true = tw_subset12[["month"]], # 混同行列 predict = tw_lda12_fitted[["class"]]) bind_cols(tw_lda12_fitted[["x"]], # 判別関数値の散布図を作成 tw_subset12["month"]) |> GGally::ggpairs(aes(colour = month)) #' 判別関数は説明変数の数までしか作成できないので，精度はあまり高くないことがわかる #' 雨の有無を識別する例 tw_rain <- tw_data |> mutate(rain = factor(rain > 0), # 雨の有無でラベル化する month = as_factor(month)) # 月ごとの気候の違いの補正のため tw_rain_lda <- lda(rain ~ temp + solar + wind + month, data = tw_rain, subset = idx) # 一部のデータで推定，12ヶ月分の例とは別の指定の仕方 tw_rain_lda_fitted <- predict(tw_rain_lda) as_tibble(tw_rain_lda_fitted[["x"]]) |> # 判別関数値の視覚化 mutate(rain = tw_rain[["rain"]][idx]) |> ggplot(aes(x = LD1, fill = rain)) + geom_histogram(aes(y = after_stat(density)), bins = 30) + facet_grid(rain ~ .) #' 全データを予測 tw_rain_lda_predict <- predict(tw_rain_lda, newdata = tw_rain) table(true = tw_rain[["rain"]][idx], # 推定に用いたデータの混同行列 predict = tw_rain_lda_predict[["class"]][idx]) table(true = tw_rain[["rain"]][-idx], # 未知データに対する予測の混同行列 predict = tw_rain_lda_predict[["class"]][-idx]) #+end_src #+begin_src R :exports none #' --------------------------------------------------------------------------- #+end_src