Rで機械学習:ランダムフォレスト

# 再現性のための乱数シード設定
seed <- 20250812
set.seed(seed)

# データ生成
# スマートファームで収穫されたトマトのセンサーデータを生成します。
# 特徴量は「糖度」「酸度」「硬度」「色相」「重量」の5つです。
n <- 250 # 各クラスのデータ数
# 高品質トマトのデータ
premium_tomato <- data.frame(
  sugar_content = rnorm(n, mean = 8.5, sd = 0.8),
  acidity = rnorm(n, mean = 0.4, sd = 0.05),
  firmness = rnorm(n, mean = 12, sd = 1.0),
  color_hue = rnorm(n, mean = 25, sd = 3),
  weight = rnorm(n, mean = 150, sd = 15),
  quality = "高品質"
)
# 標準品質トマトのデータ
standard_tomato <- data.frame(
  sugar_content = rnorm(n, mean = 6.5, sd = 1.0),
  acidity = rnorm(n, mean = 0.5, sd = 0.08),
  firmness = rnorm(n, mean = 10, sd = 1.5),
  color_hue = rnorm(n, mean = 35, sd = 5),
  weight = rnorm(n, mean = 130, sd = 20),
  quality = "標準品質"
)

tomato_data <- rbind(premium_tomato, standard_tomato)
tomato_data$quality <- as.factor(tomato_data$quality)

# 各特徴量の分布をバイオリンプロットで比較
tomato_data_jp <- tomato_data %>%
  rename(
    糖度 = sugar_content,
    酸度 = acidity,
    硬度 = firmness,
    色相 = color_hue,
    重量 = weight,
    品質 = quality
  )

# データをプロットしやすい「ロング形式」に変換
tomato_long <- tomato_data_jp %>%
  pivot_longer(
    cols = c(糖度, 酸度, 硬度, 色相, 重量),
    names_to = "特徴量",
    values_to = "値"
  )

# バイオリンプロットの作成
p_violin <- ggplot(tomato_long, aes(x = 品質, y = 値, fill = 品質)) +
  geom_violin(trim = FALSE, alpha = 0.6) + # バイオリンプロット
  geom_boxplot(width = 0.15, fill = "white", outlier.shape = NA) + # 箱ひげ図を重ねる
  facet_wrap(~特徴量, scales = "free_y", ncol = 5) + # 特徴量ごとにプロットを分割
  scale_fill_manual(values = c("高品質" = "#F8766D", "標準品質" = "#00BFC4")) +
  labs(
    title = "トマトの品質別 センサーデータ分布",
    subtitle = "各センサー値について「高品質」と「標準品質」の分布を比較",
    x = "トマトの品質",
    y = "センサー測定値"
  ) +
  theme_bw() +
  theme(
    legend.position = "none", # 凡例は不要
    axis.text.x = element_text(angle = 0, vjust = 0),
    strip.text = element_text(size = 12)
  ) # 各プロットのタイトルサイズ

print(p_violin)

# GGallyパッケージのggpairs関数でペアプロットを作成
library(GGally)
p_pairs <- ggpairs(
  tomato_data_jp,
  columns = 1:5, # 1から5列目（特徴量）を使用
  mapping = aes(color = 品質, shape = 品質, alpha = 0.7),
  title = "トマトのセンサーデータ相関分析（品質別）",
  upper = list(continuous = wrap("cor", size = 3)), # 上半分は相関係数を表示
  lower = list(continuous = wrap("points", size = 1.5)), # 下半分は散布図を表示
  diag = list(continuous = wrap("densityDiag", alpha = 0.5)) # 対角線は密度プロットを表示
) +
  scale_fill_manual(values = c("高品質" = "#F8766D", "標準品質" = "#00BFC4")) +
  scale_shape_manual(values = c("高品質" = 1, "標準品質" = 2)) +
  theme_bw() +
  theme(strip.text = element_text(size = 10))

print(p_pairs)

# ランダムフォレストモデルを500本の決定木で学習させます
# importance=TRUE にして特徴量の重要度を計算させる
rf_model <- randomForest(quality ~ ., data = tomato_data, ntree = 500, importance = TRUE)

# モデルの概要を出力
print(rf_model)

# どの特徴量がトマトの品質判定に重要だったかを可視化します

# 特徴量の重要度データを取得
importance_data <- as.data.frame(importance(rf_model))
importance_data$feature <- rownames(importance_data)
# 日本語の特徴量名を追加
importance_data$feature_jp <- factor(c("糖度", "酸度", "硬度", "色相", "重量"))

# MeanDecreaseGini（ジニ不純度の減少量平均）を使ってプロット
p_importance <- ggplot(importance_data, aes(x = reorder(feature_jp, MeanDecreaseGini), y = MeanDecreaseGini)) +
  geom_col(fill = "skyblue", color = "black") +
  coord_flip() + # 横向きの棒グラフにする
  labs(
    title = "特徴量の重要度（トマトの品質判定）",
    subtitle = "値が大きいほど、品質の分類に大きく貢献した特徴量であることを示す",
    x = "特徴量",
    y = "重要度 (Mean Decrease Gini)"
  ) +
  theme_bw()

print(p_importance)

# 決定木の数を増やしていくとOOB誤差率がどのように変化するか確認します

# OOB誤差率のデータを抽出してデータフレームに変換
oob_error_data <- as.data.frame(rf_model$err.rate)
oob_error_data$num_trees <- 1:nrow(oob_error_data)

# データをプロットしやすい形式（ロング形式）に変換
oob_error_long <- oob_error_data %>%
  pivot_longer(cols = -num_trees, names_to = "error_type", values_to = "error_rate") %>%
  mutate(error_type_jp = case_when(
    error_type == "OOB" ~ "全体 (OOB)",
    error_type == "高品質" ~ "高品質クラス",
    error_type == "標準品質" ~ "標準品質クラス"
  ))

# OOB誤差率の推移をプロット
p_oob <- ggplot(oob_error_long, aes(x = num_trees, y = error_rate, color = error_type_jp)) +
  geom_line(linewidth = 1) +
  labs(
    title = "決定木の数とOOB誤差率の関係",
    subtitle = "木の数が増えるにつれて誤差率が収束し、モデルが安定していく様子",
    x = "決定木の数",
    y = "OOB誤差率",
    color = "誤差の種類"
  ) +
  scale_y_continuous(labels = scales::percent_format()) + # y軸をパーセント表示に
  theme_bw() +
  theme(legend.position = "bottom")

print(p_oob)