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

library(tidyverse) library(broom) library(modelr) # サンプルデータの作成 set.seed(123) house_data <- tibble( area = runif(200, 50, 200), # 面積（平方メートル） age = runif(200, 0, 30), # 築年数 station_distance = runif(200, 1, 20), # 駅距離（分） price = 30 + area * 0.8 - age * 2 - station_distance * 1.5 + rnorm(200, 0, 15) ) %>% mutate(price = pmax(price, 20)) # 最低価格を20万円に設定 print(head(house_data)) # 面積と価格の単純線形回帰 simple_model <- lm(price ~ area, data = house_data) # broomパッケージで結果を整理 model_summary <- tidy(simple_model) model_metrics <- glance(simple_model) model_predictions <- augment(simple_model) print("回帰係数:") print(model_summary) print("モデル評価指標:") print(model_metrics) # 予測値と残差の分析 prediction_analysis <- model_predictions %>% select(area, price, .fitted, .resid) %>% mutate( abs_error = abs(.resid), percent_error = abs(.resid) / price * 100 ) print("予測分析（上位10件）:") print(head(prediction_analysis, 10)) # モデルの解釈 area_coefficient <- model_summary$estimate[2] r_squared <- model_metrics$r.squared print(paste("面積1平方メートル増加あたりの価格上昇:", round(area_coefficient, 2), "万円")) print(paste("決定係数 R²:", round(r_squared, 3))) print(paste("面積で説明できる価格の変動:", round(r_squared * 100, 1), "%"))

# 重回帰モデルの構築 multiple_model <- lm(price ~ area + age + station_distance, data = house_data) # 結果の比較 multiple_summary <- tidy(multiple_model) multiple_metrics <- glance(multiple_model) print("重回帰係数:") print(multiple_summary) print("重回帰評価指標:") print(multiple_metrics) # モデル比較 model_comparison <- tibble( model = c("単純回帰", "重回帰"), r_squared = c(model_metrics$r.squared, multiple_metrics$r.squared), adj_r_squared = c(model_metrics$adj.r.squared, multiple_metrics$adj.r.squared), rmse = c(model_metrics$sigma, multiple_metrics$sigma), aic = c(model_metrics$AIC, multiple_metrics$AIC) ) print("モデル比較:") print(model_comparison) # 残差分析 multiple_predictions <- augment(multiple_model) residual_analysis <- multiple_predictions %>% summarise( mean_residual = mean(.resid), sd_residual = sd(.resid), rmse = sqrt(mean(.resid^2)), mae = mean(abs(.resid)), outliers_count = sum(abs(.resid) > 2 * sd(.resid)) ) print("残差統計:") print(residual_analysis) # 予測精度の比較 accuracy_comparison <- house_data %>% mutate( simple_pred = predict(simple_model, .), multiple_pred = predict(multiple_model, .), simple_error = abs(price - simple_pred), multiple_error = abs(price - multiple_pred) ) %>% summarise( simple_mae = mean(simple_error), multiple_mae = mean(multiple_error), improvement = (simple_mae - multiple_mae) / simple_mae * 100 ) print("予測精度の改善:") print(accuracy_comparison)

# 顧客満足度データの作成 set.seed(456) customer_data <- tibble( customer_id = 1:300, service_type = sample(c("ベーシック", "プレミアム", "エンタープライズ"), 300, replace = TRUE, prob = c(0.5, 0.3, 0.2)), region = sample(c("関東", "関西", "その他"), 300, replace = TRUE), satisfaction = case_when( service_type == "ベーシック" ~ rnorm(n(), 6.5, 1.2), service_type == "プレミアム" ~ rnorm(n(), 7.5, 1.0), service_type == "エンタープライズ" ~ rnorm(n(), 8.2, 0.8) ), monthly_usage = case_when( service_type == "ベーシック" ~ rpois(n(), 15), service_type == "プレミアム" ~ rpois(n(), 25), service_type == "エンタープライズ" ~ rpois(n(), 40) ) ) %>% mutate(satisfaction = pmax(1, pmin(10, satisfaction))) # 1-10の範囲に制限 # 記述統計 descriptive_stats <- customer_data %>% group_by(service_type) %>% summarise( count = n(), mean_satisfaction = mean(satisfaction), sd_satisfaction = sd(satisfaction), mean_usage = mean(monthly_usage), sd_usage = sd(monthly_usage), .groups = 'drop' ) print("サービス別記述統計:") print(descriptive_stats) # t検定（2群比較） basic_premium_test <- customer_data %>% filter(service_type %in% c("ベーシック", "プレミアム")) %>% t.test(satisfaction ~ service_type, data = .) print("ベーシック vs プレミアムのt検定:") print(tidy(basic_premium_test)) # 分散分析（ANOVA） anova_model <- aov(satisfaction ~ service_type, data = customer_data) anova_results <- tidy(anova_model) print("サービス別満足度の分散分析:") print(anova_results) # 多重比較（Tukey HSD） tukey_results <- TukeyHSD(anova_model) tukey_tidy <- tidy(tukey_results) print("多重比較結果:") print(tukey_tidy) # カイ二乗検定（地域とサービス種別の関連性） contingency_table <- customer_data %>% count(region, service_type) %>% pivot_wider(names_from = service_type, values_from = n, values_fill = 0) print("地域別サービス利用状況:") print(contingency_table) # カイ二乗検定の実行 chi_square_test <- customer_data %>% select(region, service_type) %>% table() %>% chisq.test() print("カイ二乗検定結果:") print(tidy(chi_square_test)) # 効果量の計算（Cohen's d） cohens_d_calculation <- customer_data %>% filter(service_type %in% c("ベーシック", "プレミアム")) %>% group_by(service_type) %>% summarise( mean_sat = mean(satisfaction), sd_sat = sd(satisfaction), n = n(), .groups = 'drop' ) # Cohen's dの手動計算 mean_diff <- cohens_d_calculation$mean_sat[2] - cohens_d_calculation$mean_sat[1] pooled_sd <- sqrt(((cohens_d_calculation$n[1] - 1) * cohens_d_calculation$sd_sat[1]^2 + (cohens_d_calculation$n[2] - 1) * cohens_d_calculation$sd_sat[2]^2) / (cohens_d_calculation$n[1] + cohens_d_calculation$n[2] - 2)) cohens_d <- mean_diff / pooled_sd print(paste("Cohen's d (ベーシック vs プレミアム):", round(cohens_d, 3)))

# 顧客離脱データの作成 set.seed(789) churn_data <- customer_data %>% mutate( contract_length = sample(1:36, n(), replace = TRUE), support_calls = rpois(n(), 2), # 離脱確率を複数要因で決定 churn_prob = plogis( -2 + ifelse(service_type == "ベーシック", 1, 0) + (10 - satisfaction) * 0.3 + support_calls * 0.2 - contract_length * 0.05 ), churned = rbinom(n(), 1, churn_prob) ) %>% select(-churn_prob) # 離脱率の確認 churn_summary <- churn_data %>% group_by(service_type) %>% summarise( total_customers = n(), churned_customers = sum(churned), churn_rate = mean(churned), .groups = 'drop' ) print("サービス別離脱率:") print(churn_summary) # ロジスティック回帰モデル logistic_model <- glm( churned ~ satisfaction + monthly_usage + support_calls + contract_length + service_type, data = churn_data, family = binomial() ) # モデル結果の整理 logistic_summary <- tidy(logistic_model, exponentiate = TRUE) logistic_metrics <- glance(logistic_model) print("ロジスティック回帰結果（オッズ比）:") print(logistic_summary) # 予測と評価 predictions <- churn_data %>% mutate( predicted_prob = predict(logistic_model, type = "response"), predicted_class = ifelse(predicted_prob > 0.5, 1, 0) ) # 混同行列 confusion_matrix <- predictions %>% count(churned, predicted_class) %>% pivot_wider(names_from = predicted_class, values_from = n, values_fill = 0) %>% rename(predicted_0 = "0", predicted_1 = "1") print("混同行列:") print(confusion_matrix) # モデル評価指標 TP <- confusion_matrix %>% filter(churned == 1) %>% pull(predicted_1) TN <- confusion_matrix %>% filter(churned == 0) %>% pull(predicted_0) FP <- confusion_matrix %>% filter(churned == 0) %>% pull(predicted_1) FN <- confusion_matrix %>% filter(churned == 1) %>% pull(predicted_0) model_evaluation <- tibble( metric = c("Accuracy", "Precision", "Recall", "F1-Score", "Specificity"), value = c( (TP + TN) / (TP + TN + FP + FN), # Accuracy TP / (TP + FP), # Precision TP / (TP + FN), # Recall 2 * (TP / (TP + FP)) * (TP / (TP + FN)) / ((TP / (TP + FP)) + (TP / (TP + FN))), # F1 TN / (TN + FP) # Specificity ) ) print("モデル評価指標:") print(model_evaluation) # 特徴量重要度（オッズ比の解釈） feature_importance <- logistic_summary %>% filter(term != "(Intercept)") %>% mutate( importance = abs(log(estimate)), direction = ifelse(estimate > 1, "離脱促進", "離脱抑制") ) %>% arrange(desc(importance)) print("特徴量重要度:") print(feature_importance)