統計データ解析

第8講 サンプルコード

準備

以下で利用する共通パッケージを読み込む．

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月の気温と湿度のデータを抽出する．

  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)) # 月を因子化

• 半分のデータを用いて線形判別関数を構成し，残りのデータを用いて判別を行う．

  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) # 新しいデータの予測

解答例

必要なデータを整理する．

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,] # 試験データ(評価に用いる)

ランダムに分割するには例えば以下のようにすれば良い．

n <- nrow(tw_subset); idx <- sample(1:n, n/2)

関数 lubridate::month() を用いれば月を文字列(因子)とすることもできる．

mutate(month = month(month, label = TRUE))

以下のようにして最初に月を因子化することもできるが，関数によっていは処理で使われない因子(例えば1月)も表示してしまうので注意が必要である．

read_csv("data/tokyo_weather.csv") |>
    mutate(month = as_factor(month))

また複数の項目をそれぞれ因子化するには例えば以下のようにすればよい．

mutate(across(c(year, month, day), as_factor))

気温と湿度の散布図を作成して，データを視覚化する．

tw_subset |> # 
    ggplot(aes(x = temp, y = humid)) + 
    geom_point(aes(colour = month)) + # 月ごとに点の色を変える
    labs(x = "Temperature", y = "Humidity",
         title = "September & October")

訓練データを用いて線形判別関数(等分散性を仮定)を作成する．

tw_model <- month ~ temp + humid
tw_lda <- lda(formula = tw_model, data = tw_train)

判別関数によるクラス分類結果を確認する．

tw_lda_fitted <- predict(tw_lda) # 
tibble(true = tw_train[["month"]]) |>             # 真値と予測値の対応表を作成
    mutate(predict = tw_lda_fitted[["class"]]) |>
    table()                                       # 混同行列

    predict
true  9 10
  9  14  1
  10  1 15

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方向に並べる)

この例ではデータフレームを作成しているが，関数 table() のみで混同行列を作ることはできる．

table(true = tw_train[["month"]],
      predict = tw_lda_fitted[["class"]])

試験データによる評価を行う．

tw_lda_predict <- predict(tw_lda, newdata = tw_test) 
tibble(true = tw_test[["month"]]) |>               # 真値と予測値の対応表を作成
    mutate(predict = tw_lda_predict[["class"]]) |>
    table()                                        # 混同行列

    predict
true  9 10
  9  14  1
  10  0 15

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 ~ .) 

推定された線形判別関数の識別境界を図示する．

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)) 

判別関数値を色づけして図示するには以下のようにすればよい．

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次判別を行いなさい．

  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) # 新しいデータの予測

• 別の月や変数を用いて判別分析を行いなさい．

解答例

前の練習問題で作成した tw_train および tw_test を利用する．

まず，訓練データで2次判別関数(等分散性を仮定しない)を作成する． 線形判別関数と同じモデル式を用いる．

tw_qda <- qda(formula = tw_model, data = tw_train)
tw_qda_fitted <- predict(tw_qda)                  # 判別関数によるクラス分類結果の取得
tibble(true = tw_train[["month"]]) |>             # 真値と予測値の対応表を作成
    mutate(predict = tw_qda_fitted[["class"]]) |>
    table()                                       # 予測値の比較(混同行列) 

    predict
true  9 10
  9  14  1
  10  1 15

試験データによる評価を行う．

tw_qda_predict <- predict(tw_qda, newdata = tw_test) 
table(true = tw_test[["month"]], # 混同行列
      predict = tw_qda_predict[["class"]])

    predict
true  9 10
  9  14  1
  10  0 15

tibble(true = tw_test[["month"]]) |> # 比較表
    mutate(predict = tw_qda_predict[["class"]])

# A tibble: 30 × 2
   true  predict
   <fct> <fct>  
 1 9     9      
 2 9     9      
 3 9     9      
 4 9     9      
 5 9     9      
 6 9     9      
 7 9     9      
 8 9     9      
 9 9     9      
10 9     9      
# ℹ 20 more rows

推定された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))

多値判別

問題

東京の気候データを用いて以下の分析を行いなさい．

• 9月，10月，11月の気温と湿度のデータを用いて判別関数を作成しなさい．

  tw_subset  <- tw_data |>
      filter(month %in% c(9,10,11)) |>
      select(temp, humid, month) |>
      mutate(month = as_factor(month))

• 別の月や変数を用いて判別分析を行いなさい．

  #' 雨の有無を識別する例
  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' をそれ以外で判別

解答例

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"]]) # 予測値の比較(混同行列)

    predict
true  9 10 11
  9  28  2  0
  10  1 26  4
  11  0 10 20

推定された線形判別関数(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"]])

    predict
true  1  2  3  4  5  6  7  8  9 10 11 12
  1  12  0  0  0  0  0  0  0  0  0  0  0
  2   2  0  0  0  0  0  0  0  0  0  0  2
  3   0  0  5  1  0  0  0  0  0  0  1  1
  4   0  0  1  6  2  0  0  0  0  1  2  0
  5   0  0  0  3  3  0  0  0  0  0  0  0
  6   0  0  0  0  0  3  0  0  0  1  0  0
  7   0  0  0  0  0  0  5  4  1  0  0  0
  8   0  0  0  0  0  0  0 13  0  0  0  0
  9   0  0  0  0  0  0  1  1  1  0  0  0
  10  0  0  0  0  0  0  1  0  0 10  0  0
  11  1  0  1  0  0  0  0  0  0  3  1  4
  12  2  0  0  0  0  0  0  0  0  0  3  2

bind_cols(tw_lda12_fitted[["x"]],           # 判別関数値の散布図を作成
          tw_subset12["month"]) |>
    GGally::ggpairs(aes(colour = month),
                    upper = list(continuous = GGally::wrap("cor", size = 1.5)))

判別関数は説明変数の数までしか作成できないので，精度はあまり高くないことがわかる．

雨の有無を識別する例は以下のようになる．

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)                       # 一部のデータで推定する例
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])

       predict
true    FALSE TRUE
  FALSE    69    6
  TRUE     11   14

table(true = tw_rain[["rain"]][-idx],                 # 未知データに対する予測の混同行列
      predict = tw_rain_lda_predict[["class"]][-idx])

       predict
true    FALSE TRUE
  FALSE   178   12
  TRUE     38   37