UCI Machine Learning Repository の Obesity データの分析2 - 探索的データ分析: EDA (Exploratory Data Analysis) の実践その1

www.crosshyou.info

の続きです。今回は 探索的データ分析: EDA (Exploratory Data Analysis) という作業をしていきます。このデータの目的は、obesity の7つのカテゴリーを予測する、ということですから、obesity と他の変数の関係性を調べようと思います。

obesity と Gender: 性別の関係を見ます。

obesity, Gender ともにカテゴリカル変数なので、棒グラフで頻度をグラフにします。

Obesity_Type_II は、ほとんど男性、Obesity_TYpe_III は、ほとんど女性だとわかります。

また、obesity のタイプ別の数は、だいたい同じだということもわかります。

次は、obesity と Age: 年齢 です。Age は数値データなので、箱ひげ図にしてみます。

Age の平均値が高い obesity のタイプが一番上にくるように reorder() 関数を使っています。Insufficient_Weight が一番若いです。若い人は痩せている人が多いですね。

次は obesity と Height: 身長です。これも Age と同じく箱ひげ図を描きます。

Obesity_Type_II の身長が突出しています。これは、Obesity_Type_II がほとんど男性だからだと思います。

次は、obesity と Weight]: 体重です。

体重は明らかに obesity のタイプと関係していますね。というか体重と身長を基にして obesity のタイプを分類しているのかもしれません。実際に分類モデルを作るときは、Weight と Height は除いて分類モデルを作ることにします。

今回はここまでにして、残りの変数は次回以降にやろうと思います。

次回は

www.crosshyou.info

です。

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# EDA (Exploratory Data Analysis)]
#
# obesity と Gender
df |> 
  ggplot(aes(x = obesity, fill = Gender)) +
  geom_bar() +
  theme_minimal()
#
# obesity と Age
df |> 
  mutate(obesity = reorder(obesity, Age)) |> 
  ggplot(aes(x = Age, y = obesity, group = obesity)) +
  geom_boxplot(aes(fill = obesity)) +
  theme_minimal() +
  theme(legend.position = "none")
#
# obesity と Height
df |> 
  mutate(obesity = reorder(obesity, Height)) |> 
  ggplot(aes(x = Height, y = obesity, group = obesity)) +
  geom_boxplot(aes(fill = obesity)) +
  theme_minimal() +
  theme(legend.position = "none")
#
# obesity と Weight
df |> 
  mutate(obesity = reorder(obesity, Weight)) |> 
  ggplot(aes(x = Weight, y = obesity, group = obesity)) +
  geom_boxplot(aes(fill = obesity)) +
  theme_minimal() +
  theme(legend.position = "none")
#

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは Landscape of natural wildness grass field, allover covering yellow Garbera flowers, photo です。)

UCI Machine Learning Repository の Obesity データの分析1 - CSVファイルのデータをRに読み込ませる。

今回からしばらくは、UCI Machine Learning Repository の Obesity(肥満)のデータを使ってみたいと思います。

Estimation of Obesity Levels Based On Eating Habits and Physical Condition  [Dataset]. (2019). UCI Machine Learning Repository. https://doi.org/10.24432/C5H31Z.

ダウンロードした CSV ファイルは下の画像のようなものでした。

一番右の列の NObeyesdad が肥満タイプで、このカテゴリカル変数を分類する、というものです。

R で分析しますので、tidyverse と tidymodels パッケージを読み込みしておきます。

read_csv() 関数で CSVファイルを読み込みます。

relocate() 関数で、目的変数の NObeyesdad を一番左にもってきました。

glimpse() 関数でデータの様子を確認します。

問題なくデータは読み込まれました。NObeyesdad は 名前をわかりやすく、obesity に、family_history_with_overweight は FHWO と短縮します。これは rename() 関数で実行します。

NA の有無を確認します。sapply(), function(), sum(), is.na() 関数を使いました。

全部の変数で、0ですから、NA はありません。
次に文字列をファクター型にします。mutate(), across(), where(), as.factor() を使います。

summary() 関数でそれぞれの変数のサマリー統計値を確認します。

目的変数の obesity は7つのタイプがだいたい均等に分かれていることがわかります。

今回は以上です。

次回は

www.crosshyou.info

です。

今回のコードは以下になります。

#
# tidyverse, tidymodels の読み込み
library(tidyverse)
library(tidymodels)
#
# CSVファイルの読み込み
df_raw <- read_csv("obesity.csv") |> 
  relocate(NObeyesdad)
#
# データの確認
glimpse(df_raw)
#
# 変数名の変更
df <- df_raw |> 
  rename(obesity = NObeyesdad,
         FHWO = family_history_with_overweight)
#
# NAの有無
sapply(df, function(x) sum(is.na(x)))
#
# 文字列型はファクター型に変換する
df <- df |> 
  mutate(
    across(where(is.character), as.factor)
  )
df 
#
# データのサマリー
summary(df)
#

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは This is a beautiful landscape photograph. Under a blue sky, a range of snow-capped mountains stretches into the distance. In the foreground is a meadow, with a close-up of a red sweet pea in bloom. です。)

都道府県別の定期健康診断結果報告のデータの分析7 - 勾配ブースティングモデルでの回帰分析

www.crosshyou.info

の続きです。今回は勾配ブースティングモデルで所見率を回帰分析してみます。

xgboost パッケージを読み込みます。

説明変数(per_jushin, log_place, log_jushin)を行列に変換します。

shokenritsu を被説明変数として取り出します。

XGBoost用のDMatrixを作成します。

ハイパーパラメータの設定をします。

モデルの学習をします。

予測をします。

RMSEと相関係数を計算します。

さすが勾配ブースティングモデルですね。RMSEは三つのモデルの中で一番小さく、相関係数は一番大きいです。

変数の重要度を見てみます。

log_jushinのGainの値が 0.387で一番大きいので、log_jushinが一番重要ですね。

実際の値と予測値の散布図を描きます。

緑色のBが勾配ブースティングモデルでの結果です。散布図でもBが一番黒の直線に近いことがわかりますね。

今回は以上です。

肇から読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# 勾配ブースティング
# パッケージの読み込み
library(xgboost)
#
# 説明変数を行列に変換
X <- as.matrix(df[, c("per_jushin", "log_place", "log_jushin")])
#
# 被説明変数を取り出す
y <- df$shokenritsu
#
# XGBoost用のDMatrixを作成
dtrain <- xgb.DMatrix(data = X, label = y)
#
# ハイパーパラメータの設定
params <- list(
  objective = "reg:squarederror",
  max_depth = 1,      # 木の深さが1
  eta = 0.1,
  subsample = 0.8,
  colsample_bytree = 0.8
)
#
# モデルの学習
set.seed(1234)
xgb_mod <- xgb.train(
  params = params,
  data = dtrain,
  nrounds = 400 # 木の数は400本
)
#
# 予測
pred_xgb <- predict(xgb_mod, dtrain)
#
# RMSEとCorrelation
df_rmse <- df_rmse |> 
  rbind(
    tibble(
      Model = "Boosting",
      RMSE = sqrt(mean*1,
      Corr = cor(pred_xgb, df$shokenritsu) 
    )
  )
df_rmse
#
# 変数重要度
xgb.importance(model = xgb_mod)
#
# 実際の値と予測図のサンプル
x_range <- range(c(lm_pred, tree_pred, pred_xgb))
y_range <- range(df$shokenritsu)
plot(x_range, y_range, type = "n",
     xlab = "estimate", ylab = "actual",
     main = "Linear Model vs. Tree Model vs. Boosting")
points(predict(lm_mod4v5), df$shokenritsu, pch = "L", col = "red")
points(tree_pred, df$shokenritsu, pch = "T", col = "blue")
points(pred_xgb, df$shokenritsu, pch = "B", col = "green")
abline(0, 1, col = "black")
legend("topleft", c("Linear Model", "Tree Model","Boosting"),
       pch = c("L", "T", "B"), col = c("red", "blue", "green"), bty = "n")
#

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Natural landscape, small mountain is fully covered by red "Tsutsuji Flowers" under the blue sky, Photo です。)

*1:pred_xgb - df$shokenritsu)^2

読書記録 - 「オベリスクの門 <破壊された地球> 三部作」 N. K. ジェミシン 著 (創元SF文庫)

オベリスクの門 〈破壊された地球〉三部作 (創元SF文庫)

<破壊された地球>三部作の第2作目です。

前作は、エッスン、アマダ、サイアナイトの三人の話が交互に語られるという話でしたが、今作はエッスンとその娘のナッスンの物語です。

エッスンはナッスンのことをとても愛してるのに、ナッスンはそうでないことが悲しかったです。

オロジェン、石喰い、守護者の謎や月が何故なくなったかやオベリスクは誰が作ったかどはまだ解明されていないくて、最終作に持ち込みです。エッスンにとってハッピーエンドになって欲しいですが、難しいかもとも思います。

第5の季節を迎えた地球はますます人類が生き延びるには過酷な環境になっていきます。現実の地球も温暖化がどんどんと進み、暮らしていくには大変な環境になったらどうなるのかな、と思いました。

都道府県別の定期健康診断結果報告のデータの分析6 - 決定木モデルでの回帰分析

www.crosshyou.info

の続きです。前回は lm() 関数を使って、線形モデルで重回帰分析をしました。R-squared は 0.22 ということで残念ながら線形モデルでは、所見率は上手く説明できないようでした。そこで今回は決定木モデルを使ってみます。

はじめに rpart パッケージ、rpart.plot パッケージの読み込みをします。

モデルを作成します。

cp の確認をします。

cp = 0.01 のときが一番、xerror が小さくなります。

最適な cp を保存して、この cp で木の選定をします。

作成した決定木をグラフにします。

一番左の枝に注目すると、log_jushin が 12 以上で、さらに per_jushin が 118 以上だと、shokenritsu は 52 となる、ということです。

このモデルの予想値と実際の値を散布図にしてみます。

前回の線形モデルと今回の決定木モデルを比較してみます。RMSEと相関係数を比較します。

RMSEは線形モデルが 3.18 なのに対して、決定木モデルは 2.51 と 21% も低下しました。相関係数は、線形モデルが 0.473 なのに対して、決定木モデルは 0.720 と 52% も改善しています。shokenritsu は per_jushin, log_place, log_jushin に対しては線形の関係ではないことがわかりますね。

最後に線形モデルと決定木モデルの予測値と実際の値の散布図を描いてみます。

青い T が決定木モデルの予測値、赤い L が線形モデルの予測値です。黒い直線が傾き 1, 切片 0 の直線です。

今回は以上です。

次回は

www.crosshyou.info

です。

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# rpartパッケージで tree model で regression
#
library(rpart)
library(rpart.plot)
#
# モデルを作成
set.seed(277)
rpart_mod <- rpart(
  shokenritsu ~ per_jushin + log_place + log_jushin,
  data = df,
  method = "anova"
)
#
# cp の確認
printcp(rpart_mod)
plotcp(rpart_mod)
#
# 最適 cp で剪定
best_cp <- rpart_mod$cptable[which.min(rpart_mod$cptable[,"xerror"]), "CP"]
best_cp
pruned_rpart <- prune(rpart_mod, cp = best_cp)
#
# 木の可視化
rpart.plot(pruned_rpart)
#
# 予測
tree_pred <- predict(pruned_rpart, newdata = df)
#
# 予測値と実際の値の散布図
tibble(
  estimate = tree_pred,
  actual = df$shokenritsu
) |> 
  ggplot(aes(x = estimate, y = actual)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red") +
  theme_minimal()
#
# Linear Modelの予測とTree ModelのRMSE, 相関の比較
lm_pred <- predict(lm_mod4v5, newdata = df)
df_rmse <- tibble(
  Model = c("Linear", "Tree"),
  RMSE = c(sqrt(mean*1,
           sqrt(mean*2
)
df_rmse
#
# 散布図での比較
x_range <- range(c(lm_pred, tree_pred))
y_range <- range(df$shokenritsu)
plot(x_range, y_range, type = "n",
     xlab = "estimate", ylab = "actual",
     main = "Linear Model vs. Tree Model")
points(predict(lm_mod4v5), df$shokenritsu, pch = "L", col = "red")
points(tree_pred, df$shokenritsu, pch = "T", col = "blue")
abline(0, 1, col = "black")
legend("topleft", c("Linear Model", "Tree Model"),
       pch = c("L", "T"), col = c("red", "blue"), bty = "n")
#

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Natural landscape of white Jasmin flower fields, under the blue sky, there are so much flowers, Photo です。)

*1:lm_pred - df$shokenritsu)^2

*2:tree_pred - df$shokenritsu)^2))),
  Corr = c(cor(lm_pred, df$shokenritsu),
           cor(tree_pred, df$shokenritsu

都道府県別の定期健康診断結果報告のデータの分析5 - 所見率を被説明変数にして重回帰分析をする。

www.crosshyou.info

の続きです。今回は、shokenritsu: 所見率を被説明変数にして重回帰分析をしてみたいと思います。

まずはじめに、per_jushin: 1事業所当たりの受診者数の数、log_place: 事業所の数の対数変換値、log_jushin: 受信者数の対数変換値を説明変数にして重回帰分析をしてみました。

p-value: 1.268e-05 とあります。これは、0.00001268 ですので、統計学的には有意なモデルです。ただ、各変数の p値は0.05よりも大きいです。R-squared が0.1978ですから、約20%くらいしか shokenritsu を説明できません。

このモデルに2乗項と交差項を加えてみます。

log_place の2乗項と log_jushin の2乗項の p値が0.05よりも小さくなりました。R-squareは0.2557と少し上昇しました。

step() 関数で不要な説明変数を削除してみます。

per_jushin: log_place: log_jushon の p値が 0.57626 ですから、これも削除できそうですね。他にも0.05以上のものがあるので、update() 関数で手動で削除していきます。

まだまだ削除できます。

I(per_jushin^2) を削除します。

これで、モデルの中の全ての変数の p値が0.05よりも小さな値になりました。R-squaredは0.2242です。このモデルで shokenritsu を 22% くらいは説明できるということですね。これだけ頑張っても 22% ですから、shokenritsu は 他の要因の影響が大きいということですね。

この最終的なモデル、lm_mod4v5 の予測値と実際の shokenritsu の散布図を描いてみます。

横軸がモデルの予測値、縦軸が実際の shokenritsu です。

今度は、このモデルに year と pref を加えます。

R-squared が 0.9654 に跳ね上がりました！都道府県と年と説明変数に加えると一気にモデルの説明力が増しました。

予測と実際の値の散布図を描いてみます。

さきほどの散布図と比べると、予測精度が改善されていることがわかります。

year と pref だけが説明変数のモデルも試してみます。

R-squared は 0.9642 ですので、ほとんど変わらないです。散布図を描きます。

散布図もほとんど変わらないですね。

今回は以上です。

次回は、

www.crosshyou.info

です。

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# shokenritsuをper_jushin, log_place, log_jushinで回帰分析
lm_mod3 <- lm(shokenritsu ~ per_jushin + log_place + log_jushin, data = df)
summary(lm_mod3)
#
# modelに2乗項と交差項を加える
lm_mod4 <- lm(shokenritsu ~ per_jushin * log_place * log_jushin +
                I(per_jushin^2) + I(log_place^2) + I(log_jushin^2),
              data = df)
summary(lm_mod4)
#
# lm_mod4 から不要な変数を削除
lm_mod4v2 <- step(lm_mod4)
summary(lm_mod4v2)
#
# 手動で不要な変数を削除 1回目
lm_mod4v3 <- update(lm_mod4v2, . ~ . - per_jushin:log_place:log_jushin)
summary(lm_mod4v3)
#
# 手動で不要な変数を削除 2回目
lm_mod4v4 <- update(lm_mod4v3, . ~ . - per_jushin:log_jushin)
summary(lm_mod4v4)
#
# 手動で不要な変数を削除 3回目
lm_mod4v5 <- update(lm_mod4v4, . ~ . - I(per_jushin^2))
summary(lm_mod4v5)
#
# 実際のshokenritsuとlm_mod4v5の予測値の散布図
tibble(
  actual = df$shokenritsu,
  estimate = predict(lm_mod4v5)
) |> 
  ggplot(aes(x = estimate, y = actual)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal()
#
# lm_mod4v5にyearとprefを追加する
lm_mod5 <- update(lm_mod4v5, ~ . + year + pref, data = df)
summary(lm_mod5)$r.squared
#
# 実際のshokenritsuとlm_mod5の予測値の散布図
tibble(
  actual = df$shokenritsu,
  estimate = predict(lm_mod5)
) |> 
  ggplot(aes(x = estimate, y = actual)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal()
#
# shokenritsuをyearとprefだけで回帰分析
lm_mod6 <- lm(shokenritsu ~ year + pref, data = df)
summary(lm_mod6)$r.squared
#
# 実際のshokenritsuとlm_mod6の予測値の散布図
tibble(
  actual = df$shokenritsu,
  estimate = predict(lm_mod6)
) |> 
  ggplot(aes(x = estimate, y = actual)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal()
#

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Natural landscape, there are white, red and pink Rhododendron transiens all over the field, under the blue sky, photo です。)

読書記録 - 「第五の季節 <破壊された地球> 三部作」N. K. ジェミシン 著 (創元SF文庫)

第五の季節 〈破壊された地球〉三部作 (創元ＳＦ文庫)

<破壊された地球>三部作の第1作目です。

この三部作で3年連続してヒューゴー賞を受賞したということで、期待をもって読み始めました。読み終わった感想は、期待どおりの面白さで、早く続きが読みたい、という気持ちともう一度最初から読み直したい、という気持ちが二つあります。

石喰いって何？オベリスクは誰が作った？などの謎がどのように解明されるのか楽しみです。一部の人間には「オロジェニー」と呼ばれる大地を操る能力が生まれながらに備わっていて、そういう人たちは、「オロジェン」と呼ばれています。その能力で回りの人を殺してしまうこともあるが故に、差別や迫害を受けているという世界です。

優れた能力を持つ少数の人間がその能力ゆえに、普通の人間社会からは迫害される、という設定はこの本だけでなくいろいろな本で見られる設定です。自分たちと異質なものを排除しようとするのは動物としての本能なのかな、と思いました。