第11章: 機械学習とtidymodels

# tidymodelsエコシステムの読み込み library(tidymodels) library(tidyverse) # コアパッケージの確認 tidymodels_packages() # サンプルデータの準備 library(palmerpenguins) data(penguins) # データの概要 penguins_clean <- penguins %>% drop_na() %>% select(-year) glimpse(penguins_clean) # データ分布の確認 penguins_clean %>% count(species, sort = TRUE)

# 初期データ分割（訓練:テスト = 80:20） set.seed(123) penguins_split <- initial_split( penguins_clean, prop = 0.8, strata = species # 層化サンプリング ) # 訓練データとテストデータの作成 penguins_train <- training(penguins_split) penguins_test <- testing(penguins_split) # データサイズの確認 cat("訓練データ:", nrow(penguins_train), "行\\n") cat("テストデータ:", nrow(penguins_test), "行\\n") # 層化サンプリングの確認 cat("\\n訓練データの種分布:\\n") penguins_train %>% count(species) %>% mutate(prop = round(n / sum(n), 3)) %>% print() cat("\\nテストデータの種分布:\\n") penguins_test %>% count(species) %>% mutate(prop = round(n / sum(n), 3)) %>% print()

# 10-fold クロスバリデーション penguins_folds <- vfold_cv( penguins_train, v = 10, strata = species ) print(penguins_folds) # Bootstrap リサンプリング penguins_boots <- bootstraps( penguins_train, times = 25, strata = species ) print(penguins_boots) # Time series用の分割（例：月次データ） ts_splits <- rolling_origin( penguins_train, initial = 150, assess = 50, skip = 25, cumulative = FALSE ) print(ts_splits)

# 基本的なレシピの作成 penguins_recipe <- recipe( species ~ ., # 目的変数と説明変数の指定 data = penguins_train ) %>% # 数値変数の正規化 step_normalize(all_numeric_predictors()) %>% # カテゴリ変数のダミー化 step_dummy(all_nominal_predictors()) %>% # ゼロ分散変数の除去 step_zv(all_predictors()) # レシピの確認 print(penguins_recipe) # レシピの訓練（prep） penguins_prep <- prep(penguins_recipe) # 変換後データの確認 penguins_baked <- bake(penguins_prep, new_data = NULL) glimpse(penguins_baked) # 正規化の確認 penguins_baked %>% select(where(is.numeric)) %>% summarise(across(everything(), list(mean = mean, sd = sd))) %>% pivot_longer(everything()) %>% print()

# 高度な前処理レシピ advanced_recipe <- recipe( species ~ ., data = penguins_train ) %>% # 外れ値の処理 step_filter( bill_length_mm < quantile(bill_length_mm, 0.99, na.rm = TRUE), bill_depth_mm < quantile(bill_depth_mm, 0.99, na.rm = TRUE) ) %>% # 新しい特徴量の作成 step_mutate( bill_ratio = bill_length_mm / bill_depth_mm, body_mass_log = log(body_mass_g), flipper_body_ratio = flipper_length_mm / body_mass_g ) %>% # 交互作用項の追加 step_interact(~ bill_length_mm:bill_depth_mm) %>% # 多項式特徴量の作成 step_poly(flipper_length_mm, degree = 2) %>% # 主成分分析 step_pca(all_numeric_predictors(), num_comp = 3) %>% # カテゴリ変数のエンコーディング step_dummy(all_nominal_predictors()) %>% # 低分散変数の除去 step_zv(all_predictors()) %>% # 高相関変数の除去 step_corr(all_numeric_predictors(), threshold = 0.9) # レシピの訓練と適用 advanced_prep <- prep(advanced_recipe) advanced_data <- bake(advanced_prep, new_data = NULL) cat("処理後の特徴量数:", ncol(advanced_data) - 1) # 目的変数を除く glimpse(advanced_data)

# recipesパッケージによる基本的な前処理パイプライン library(tidymodels) library(tidyverse) # 日本の住宅価格データセット例 set.seed(123) housing_data <- tibble( 面積 = rnorm(1000, 70, 20), 築年数 = sample(0:40, 1000, replace = TRUE), 駅距離 = rexp(1000, 0.3), 階数 = sample(1:15, 1000, replace = TRUE), エリア = sample(c("都心部", "住宅地", "郊外", "駅前"), 1000, replace = TRUE), 建物種別 = sample(c("マンション", "一戸建て", "アパート"), 1000, replace = TRUE), リフォーム済み = sample(c("済み", "未済"), 1000, replace = TRUE) ) %>% mutate( 価格 = 3000 + 50 * 面積 - 30 * 築年数 - 200 * 駅距離 + 10 * 階数 + rnorm(1000, 0, 300) ) # 【1】基本的なrecipe：標準化と正規化 basic_recipe <- recipe(価格 ~ ., data = housing_data) %>% # 数値変数の標準化（平均0、分散1） step_normalize(all_numeric_predictors()) %>% # カテゴリ変数のダミー化 step_dummy(all_nominal_predictors()) %>% # 分散がゼロに近い変数の除去 step_zv(all_predictors()) # recipeの準備と適用 prepped_recipe <- prep(basic_recipe) processed_data <- bake(prepped_recipe, new_data = NULL) cat("元データの列数：", ncol(housing_data), "\n") cat("処理後の列数：", ncol(processed_data), "\n") cat("追加された変数名：\n") print(setdiff(names(processed_data), names(housing_data)))

# 【2】多項式特徴量とカテゴリエンコーディング polynomial_recipe <- recipe(価格 ~ ., data = housing_data) %>% # 面積の2次・3次項を追加（非線形関係をキャプチャ） step_poly(面積, degree = 3) %>% # 築年数の対数変換（減衰効果をモデル化） step_log(築年数, offset = 1) %>% # 駅距離の平方根変換（距離効果の非線形性） step_sqrt(駅距離) %>% # 頻度エンコーディング（カテゴリの出現頻度を数値化） step_count(エリア, 建物種別) %>% # ターゲットエンコーディング（カテゴリごとの目的変数平均） step_lencode_glm(エリア, outcome = vars(価格)) %>% # 交互作用項の作成 step_interact(~ 面積:階数) %>% step_normalize(all_numeric_predictors()) # 【3】外れ値処理とデータクリーニング outlier_recipe <- recipe(価格 ~ ., data = housing_data) %>% # 四分位範囲による外れ値除去 step_rm_outliers(面積, 価格, method = "IQR", multiplier = 1.5) %>% # 欠損値の処理：数値変数は中央値で補完 step_impute_median(all_numeric_predictors()) %>% # カテゴリ変数は最頻値で補完 step_impute_mode(all_nominal_predictors()) %>% # Winsorization：極値を95%・5%点でクリップ step_range(all_numeric_predictors(), min = 0.05, max = 0.95) %>% step_normalize(all_numeric_predictors()) # 【4】時系列特徴量の作成 # 日付データを含むサンプル作成 timeseries_data <- housing_data %>% mutate( 売却日 = sample( seq(as.Date("2020-01-01"), as.Date("2023-12-31"), by = "day"), 1000, replace = TRUE ) ) datetime_recipe <- recipe(価格 ~ ., data = timeseries_data) %>% # 日付から年・月・曜日を抽出 step_date(売却日, features = c("year", "month", "dow")) %>% # 四半期情報の抽出 step_mutate(四半期 = lubridate::quarter(売却日)) %>% # 季節性の表現（三角関数変換） step_harmonic(売却日_month, frequency = 12, cycle_size = 12) %>% # 元の日付列は削除 step_rm(売却日) %>% step_dummy(all_nominal_predictors()) %>% step_normalize(all_numeric_predictors())

# 【5】テキストデータを含む不動産データ text_housing_data <- housing_data %>% mutate( 物件説明 = sample(c( "駅近で便利な立地です", "閑静な住宅街にあります", "リフォーム済みで綺麗です", "日当たり良好で明るいです", "ペット可物件です", "南向きで眺望良好", "新築同様の美品です", "駐車場完備しています", "商業施設が近くにあります" ), 1000, replace = TRUE) ) # textrecipesパッケージが必要 # install.packages("textrecipes") library(textrecipes) text_recipe <- recipe(価格 ~ ., data = text_housing_data) %>% # テキストのトークン化 step_tokenize(物件説明) %>% # 低頻度単語の除去（5回未満の単語を削除） step_tokenfilter(物件説明, min_times = 5) %>% # TF-IDF変換 step_tfidf(物件説明) %>% # 数値変数の標準化 step_normalize(all_numeric_predictors()) %>% # カテゴリ変数のダミー化 step_dummy(all_nominal_predictors()) # より高度なテキスト前処理 advanced_text_recipe <- recipe(価格 ~ ., data = text_housing_data) %>% # 文字数・単語数の特徴量を作成 step_mutate( 説明文字数 = nchar(物件説明), 説明単語数 = stringr::str_count(物件説明, "\\S+") ) %>% # キーワードの有無を特徴量化 step_mutate( 駅近フラグ = ifelse(stringr::str_detect(物件説明, "駅近"), 1, 0), リフォームフラグ = ifelse(stringr::str_detect(物件説明, "リフォーム"), 1, 0), 日当たりフラグ = ifelse(stringr::str_detect(物件説明, "日当たり"), 1, 0) ) %>% step_rm(物件説明) %>% step_normalize(all_numeric_predictors()) %>% step_dummy(all_nominal_predictors())

# 【6】高次元データの作成（多くの特徴量を持つデータ） high_dim_data <- housing_data %>% mutate( # 周辺環境指標を追加 学校距離 = rexp(1000, 0.5), 病院距離 = rexp(1000, 0.3), 商店街距離 = rexp(1000, 0.4), 公園距離 = rexp(1000, 0.6), コンビニ距離 = rexp(1000, 0.8), # 建物詳細指標 バルコニー面積 = rnorm(1000, 8, 3), 収納面積 = rnorm(1000, 5, 2), 天井高 = rnorm(1000, 2.4, 0.2), 部屋数 = sample(1:5, 1000, replace = TRUE) ) # PCAによる次元削減recipe pca_recipe <- recipe(価格 ~ ., data = high_dim_data) %>% # カテゴリ変数のダミー化 step_dummy(all_nominal_predictors()) %>% # 分散ゼロの変数を除去 step_zv(all_predictors()) %>% # 標準化（PCAの前に必須） step_normalize(all_numeric_predictors()) %>% # PCA：分散の90%を保持する成分数を自動選択 step_pca(all_numeric_predictors(), threshold = 0.9) %>% # 高相関変数の除去 step_corr(all_numeric_predictors(), threshold = 0.8) # ICAによる独立成分分析 ica_recipe <- recipe(価格 ~ ., data = high_dim_data) %>% step_dummy(all_nominal_predictors()) %>% step_zv(all_predictors()) %>% step_normalize(all_numeric_predictors()) %>% # 独立成分分析：5成分を抽出 step_ica(all_numeric_predictors(), num_comp = 5) # UMAP/t-SNEによる非線形次元削減 # embed パッケージが必要 umap_recipe <- recipe(価格 ~ ., data = high_dim_data) %>% step_dummy(all_nominal_predictors()) %>% step_normalize(all_numeric_predictors()) %>% # UMAP：2次元に削減 step_umap(all_numeric_predictors(), num_comp = 2)

# 【7】不均衡データの処理 # themisパッケージが必要（install.packages("themis")） library(themis) # 分類問題用の不均衡データを作成 imbalanced_data <- housing_data %>% mutate( 高額物件 = ifelse(価格 > quantile(価格, 0.8), "Yes", "No") %>% as.factor() ) %>% select(-価格) resampling_recipe <- recipe(高額物件 ~ ., data = imbalanced_data) %>% # SMOTE：合成マイノリティオーバーサンプリング step_smote(高額物件) %>% # または、単純なアップサンプリング # step_upsample(高額物件) %>% # または、ダウンサンプリング # step_downsample(高額物件) %>% step_normalize(all_numeric_predictors()) %>% step_dummy(all_nominal_predictors()) # 【8】カスタム変換関数の作成 # 独自の前処理関数を定義 価格帯分類 <- function(価格) { case_when( 価格 < 2000 ~ "低価格帯", 価格 < 4000 ~ "中価格帯", 価格 < 6000 ~ "高価格帯", TRUE ~ "超高級" ) } custom_recipe <- recipe(価格 ~ ., data = housing_data) %>% # カスタム関数を使った特徴量作成 step_mutate( # 面積と築年数の比率 面積築年数比 = 面積 / (pmax(築年数, 1)), # 駅距離の逆数（近いほど大きい値） 駅近度 = 1 / (1 + 駅距離), # 条件分岐による特徴量 新築フラグ = ifelse(築年数 <= 3, 1, 0), 大型物件フラグ = ifelse(面積 >= 100, 1, 0) ) %>% # 外れ値の処理：Yeo-Johnson変換 step_YeoJohnson(all_numeric_predictors()) %>% # Box-Cox変換（正の値のみ） # step_BoxCox(all_numeric_predictors()) %>% step_normalize(all_numeric_predictors()) %>% step_dummy(all_nominal_predictors()) # 【9】役割分担による複雑なrecipe comprehensive_recipe <- recipe(価格 ~ ., data = housing_data) %>% # 【段階1】データクリーニング step_filter(面積 > 0, 築年数 >= 0) %>% step_impute_knn(all_numeric_predictors(), neighbors = 5) %>% # 【段階2】特徴量エンジニアリング step_mutate(log_面積 = log(面積 + 1)) %>% step_interact(~ 面積:エリア) %>% # 【段階3】エンコーディング step_novel(all_nominal_predictors()) %>% step_other(all_nominal_predictors(), threshold = 0.01) %>% step_dummy(all_nominal_predictors()) %>% # 【段階4】正規化とスケーリング step_zv(all_predictors()) %>% step_normalize(all_numeric_predictors()) %>% # 【段階5】最終調整 step_corr(all_numeric_predictors(), threshold = 0.9) # 全てのrecipeの確認 cat("=== Recipe処理結果の比較 ===\n") recipes_list <- list( 基本処理 = basic_recipe, 多項式処理 = polynomial_recipe, 外れ値処理 = outlier_recipe, PCA処理 = pca_recipe, 包括処理 = comprehensive_recipe ) purrr::imap(recipes_list, function(recipe, name) { prepped <- prep(recipe) processed <- bake(prepped, new_data = NULL) cat(sprintf("%s: %d列 → %d列\n", name, ncol(housing_data), ncol(processed))) })