第10章: 統計モデリング基礎

# 必要なライブラリの読み込み library(tidyverse) library(broom) library(modelr) # サンプルデータの生成 set.seed(42) housing_data <- tibble( size = rnorm(200, 100, 30), age = sample(1:50, 200, replace = TRUE), location_score = runif(200, 1, 10), price = 50 * size + 20 * location_score - 10 * age + rnorm(200, 0, 500) + 2000 ) # データの概要確認 glimpse(housing_data) # 単純線形回帰モデル simple_model <- lm(price ~ size, data = housing_data) # モデルの概要をtidyに model_summary <- tidy(simple_model) print(model_summary) # モデル適合度 model_stats <- glance(simple_model) print(model_stats)

# 重回帰モデルの構築 multiple_model <- lm( price ~ size + age + location_score, data = housing_data ) # モデル比較 model_comparison <- bind_rows( glance(simple_model) %>% mutate(model = "Simple"), glance(multiple_model) %>% mutate(model = "Multiple") ) %>% select(model, r.squared, adj.r.squared, AIC, BIC) print(model_comparison) # 予測精度の評価 predictions <- housing_data %>% add_predictions(multiple_model) %>% add_residuals(multiple_model) %>% mutate( abs_error = abs(resid), squared_error = resid^2 ) # 予測精度メトリクス performance_metrics <- predictions %>% summarise( MAE = mean(abs_error), RMSE = sqrt(mean(squared_error)), MAPE = mean(abs(resid / price)) * 100 ) print(performance_metrics)

# 包括的な統計分析と結果解釈 library(tidyverse) library(broom) library(performance) library(ggplot2) # サンプルデータ：不動産価格予測 set.seed(42) real_estate <- tibble( area = rnorm(200, 100, 30), age = sample(1:50, 200, replace = TRUE), distance_to_station = rexp(200, 0.5), floor = sample(1:15, 200, replace = TRUE) ) %>% mutate( price = 50 * area - 20 * age - 10 * distance_to_station + 5 * floor + rnorm(200, 0, 500) ) # 重回帰モデルの構築 full_model <- lm( price ~ area + age + distance_to_station + floor, data = real_estate ) # 【1】回帰係数とその統計的有意性 coefficients_table <- tidy(full_model) %>% mutate( significance = case_when( p.value < 0.001 ~ "***（極めて有意）", p.value < 0.01 ~ "**（とても有意）", p.value < 0.05 ~ "*（有意）", p.value < 0.1 ~ ".（限界的有意）", TRUE ~ "（非有意）" ), effect_interpretation = case_when( term == "area" ~ "面積1㎡増加で価格増", term == "age" ~ "築年数1年増で価格減", term == "distance_to_station" ~ "駅距離1km増で価格減", term == "floor" ~ "階数1階増で価格増", TRUE ~ "切片（基準値）" ) ) %>% select(term, estimate, std.error, statistic, p.value, significance, effect_interpretation) print("【回帰係数の解釈表】") print(coefficients_table)

# A/Bテストデータの生成 set.seed(123) ab_test_data <- tibble( group = rep(c("A", "B"), each = 100), conversion_rate = c( rbinom(100, 1, 0.15), # Group A rbinom(100, 1, 0.22) # Group B ), revenue = ifelse(conversion_rate == 1, rnorm(200, 50, 15), 0) ) # グループ別統計量 group_stats <- ab_test_data %>% group_by(group) %>% summarise( n = n(), conversion_rate = mean(conversion_rate), avg_revenue = mean(revenue), sd_revenue = sd(revenue), .groups = 'drop' ) print(group_stats) # t検定の実行 revenue_test <- t.test( revenue ~ group, data = ab_test_data, var.equal = FALSE ) # 結果をtidyフォーマットに test_results <- tidy(revenue_test) print(test_results) # 効果量の計算（Cohen's d） cohens_d <- ab_test_data %>% group_by(group) %>% summarise(mean_rev = mean(revenue), .groups = 'drop') %>% summarise( effect_size = abs(diff(mean_rev)) / sd(ab_test_data$revenue) ) print(cohens_d)

# マーケティングキャンペーンデータ campaign_data <- tibble( age_group = sample( c("18-30", "31-45", "46-60", "60+"), 500, replace = TRUE ), response = sample( c("Yes", "No"), 500, replace = TRUE, prob = c(0.3, 0.7) ) ) # クロス集計表の作成 contingency_table <- campaign_data %>% count(age_group, response) %>% pivot_wider(names_from = response, values_from = n, values_fill = 0) print(contingency_table) # カイ二乗検定の実行 chi_test <- chisq.test( table(campaign_data$age_group, campaign_data$response) ) chi_results <- tidy(chi_test) print(chi_results) # Cramer's Vの計算（効果量） cramers_v <- sqrt(chi_test$statistic / (sum(table(campaign_data$age_group, campaign_data$response)) * (min(length(unique(campaign_data$age_group)), length(unique(campaign_data$response))) - 1))) cat("Cramer's V (効果量):", round(cramers_v, 3))

# 顧客離反予測データ set.seed(456) churn_data <- tibble( customer_id = 1:1000, tenure_months = sample(1:72, 1000, replace = TRUE), monthly_charges = runif(1000, 20, 120), total_charges = tenure_months * monthly_charges + rnorm(1000, 0, 100), support_tickets = rpois(1000, 2), contract_type = sample( c("Monthly", "One year", "Two year"), 1000, replace = TRUE ) ) %>% mutate( # 離反確率の設定（複雑なロジック） churn_prob = plogis( -2 + -0.05 * tenure_months + 0.02 * monthly_charges + 0.3 * support_tickets + case_when( contract_type == "Monthly" ~ 1, contract_type == "One year" ~ 0, contract_type == "Two year" ~ -0.5 ) ), churn = rbinom(1000, 1, churn_prob) ) %>% select(-churn_prob) # ロジスティック回帰モデル logit_model <- glm( churn ~ tenure_months + monthly_charges + support_tickets + contract_type, data = churn_data, family = binomial() ) # モデル係数の解釈 logit_results <- tidy(logit_model) %>% mutate( odds_ratio = exp(estimate), conf_low = exp(estimate - 1.96 * std.error), conf_high = exp(estimate + 1.96 * std.error) ) print(logit_results) # 予測とモデル評価 predictions <- churn_data %>% mutate( pred_prob = predict(logit_model, type = "response"), pred_class = ifelse(pred_prob > 0.5, 1, 0) ) # 混同行列 confusion_matrix <- table(predictions$churn, predictions$pred_class) print(confusion_matrix) # 精度指標 accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix) precision <- confusion_matrix[2,2] / sum(confusion_matrix[,2]) recall <- confusion_matrix[2,2] / sum(confusion_matrix[2,]) f1_score <- 2 * (precision * recall) / (precision + recall) cat("Accuracy:", round(accuracy, 3), "\n") cat("Precision:", round(precision, 3), "\n") cat("Recall:", round(recall, 3), "\n") cat("F1-score:", round(f1_score, 3), "\n")

# クロスバリデーション用のライブラリ library(rsample) # k-fold クロスバリデーション cv_folds <- vfold_cv(housing_data, v = 5) # 各foldでのモデル評価関数 evaluate_fold <- function(split) { train_data <- training(split) test_data <- testing(split) # モデル訓練 model <- lm(price ~ size + age + location_score, data = train_data) # テストデータで予測 predictions <- predict(model, test_data) # 評価指標の計算 rmse <- sqrt(mean((test_data$price - predictions)^2)) mae <- mean(abs(test_data$price - predictions)) r2 <- cor(test_data$price, predictions)^2 tibble(RMSE = rmse, MAE = mae, R_squared = r2) } # 各foldで評価実行 cv_results <- cv_folds %>% mutate( metrics = map(splits, evaluate_fold) ) %>% unnest(metrics) # クロスバリデーション結果の要約 cv_summary <- cv_results %>% summarise( mean_RMSE = mean(RMSE), sd_RMSE = sd(RMSE), mean_MAE = mean(MAE), sd_MAE = sd(MAE), mean_R2 = mean(R_squared), sd_R2 = sd(R_squared) ) print(cv_summary) # Bootstrap信頼区間 bootstrap_results <- bootstraps(housing_data, times = 100) %>% mutate( models = map(splits, ~ lm(price ~ size, data = training(.x))), slope = map_dbl(models, ~ tidy(.x)$estimate[2]) ) # 95%信頼区間 slope_ci <- quantile(bootstrap_results$slope, c(0.025, 0.975)) cat("面積係数の95%信頼区間:", round(slope_ci, 2))