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

📈 線形回帰の深淵

線形回帰は統計モデリングの基礎であり、数多くの高度な手法の出発点となります。Tidyverseの強力なツールと組み合わせることで、データの奥に隠された関係性を解明し、予測の精度を極限まで高めることができます。

🎯 単純線形回帰の基礎

simple_linear_regression.R

                        # 必要なライブラリの読み込み
                        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)
                    

📊 回帰分析結果の詳細解釈

係数	推定値	標準誤差	t値	p値	95%信頼区間
(Intercept)	2,085.432	342.567	6.085	< 0.001***	[1,414.4, 2,756.5]
size	49.823	3.421	14.567	< 0.001***	[43.12, 56.53]

🔍 係数の詳細解釈

切片 (Intercept): 2,085,432円

意味： 面積が0平方メートルの場合の理論的価格。実際には意味を持たないが、回帰直線の位置を決定する重要な要素。

統計的有意性： t値=6.085、p<0.001で統計的に有意。

面積係数: 49,823円/㎡

経済的意味： 住宅面積が1平方メートル増加するごとに、価格が約49,823円上昇することを示す。

信頼性： 95%信頼区間[43,120円〜56,530円]で、真の値がこの範囲内にある確率は95%。

実務的解釈： 10㎡の増加 → 約498,230円の価格上昇が期待される。

📈 決定係数（R²）の詳細解説

r_squared_calculation.R

                                # 決定係数の手動計算
                                predictions <- predict(simple_model)
                                actual <- housing_data$price
                                mean_actual <- mean(actual)
                                
                                # 平方和の計算
                                SSE <- sum((actual - predictions)^2)  # 残差平方和
                                SST <- sum((actual - mean_actual)^2) # 全平方和
                                SSR <- sum((predictions - mean_actual)^2) # 回帰平方和
                                
                                # R²の計算方法（2つの等価な式）
                                r_squared_method1 <- SSR / SST
                                r_squared_method2 <- 1 - (SSE / SST)
                                
                                cat("R² = SSR/SST =\", round(r_squared_method1, 4), "\n")
                                cat("R² = 1-SSE/SST =\", round(r_squared_method2, 4), "\n")
                                
                                # 相関係数との関係
                                correlation <- cor(housing_data$size, housing_data$price)
                                cat("相関係数 r =\", round(correlation, 4), "\n")
                                cat("r² =\", round(correlation^2, 4), "(単回帰の場合、R² = r²)\n")
                            

📊 R² = 0.742の意味

説明力：面積により住宅価格の変動の74.2%が説明される

未説明部分： 25.8%は他の要因（立地、築年数、設備など）による

モデルの質： 0.7以上は社会科学では良好な値

予測精度：面積のみでも比較的高い精度で価格予測が可能

📏 R²の評価基準

0.8-1.0

優秀

非常に強い関係

0.6-0.8

良好

強い関係（当モデル）

0.4-0.6

中程度

中程度の関係

0.0-0.4

弱い

弱い関係

🎯 重回帰分析の威力

multiple_regression.R

                        # 重回帰モデルの構築
                        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)
                    

📊 推定結果の一覧表示と重要指標の解釈

統計モデリングの結果を正しく解釈することは、データサイエンスにおいて最も重要なスキルの一つです。各指標が何を意味し、どのように意思決定に活用するべきかを包括的に解説します。

📋 統計推定結果の完全な解釈表

comprehensive_analysis.R

                        # 包括的な統計分析と結果解釈
                        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 tibble: 5 × 7 term estimate std.error statistic p.value significance effect_interpretation <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> 1 (Intercept) 234. 389. 0.602 5.48e- 1 （非有意）切片（基準値） 2 area 49.8 1.63 30.5 2.11e-67 ***（極めて有意）面積1㎡増加で価格増 3 age -19.9 4.36 -4.57 1.03e- 5 ***（極めて有意）築年数1年増で価格減 4 distance_to_station -9.87 9.56 -1.03 3.04e- 1 （非有意）駅距離1km増で価格減 5 floor 4.98 13.2 0.378 7.06e- 1 （非有意）階数1階増で価格増

                        # 【2】モデル適合度指標の包括的分析
                        model_performance <- glance(full_model) %>%
                          mutate(
                            r_squared_interpretation = case_when(
                              r.squared >= 0.8 ~ "優秀（説明力80%以上）",
                              r.squared >= 0.6 ~ "良好（説明力60-80%）",
                              r.squared >= 0.4 ~ "中程度（説明力40-60%）",
                              TRUE ~ "低い（説明力40%未満）"
                            ),
                            f_statistic_interpretation = case_when(
                              p.value < 0.001 ~ "モデル全体が極めて有意",
                              p.value < 0.01 ~ "モデル全体がとても有意",
                              p.value < 0.05 ~ "モデル全体が有意",
                              TRUE ~ "モデル全体が非有意"
                            )
                          )
                        
                        # 追加の評価指標を計算
                        additional_metrics <- real_estate %>%
                          add_predictions(full_model, var = "predicted") %>%
                          add_residuals(full_model, var = "residual") %>%
                          summarise(
                            MAE = mean(abs(residual)),
                            RMSE = sqrt(mean(residual^2)),
                            MAPE = mean(abs(residual / price)) * 100
                          )
                        
                        print("【モデル適合度の総合評価】")
                        cat("R²値：", round(model_performance$r.squared, 4), " - ", model_performance$r_squared_interpretation, "\n")
                        cat("調整済みR²：", round(model_performance$adj.r.squared, 4), " - より保守的な評価\n")
                        cat("F統計量：", round(model_performance$statistic, 2), " - ", model_performance$f_statistic_interpretation, "\n")
                        cat("AIC：", round(model_performance$AIC, 2), " - モデル選択の指標（小さい方が良い）\n")
                        cat("BIC：", round(model_performance$BIC, 2), " - より保守的なモデル選択指標\n")
                        cat("MAE：", round(additional_metrics$MAE, 2), " - 平均絶対誤差（万円）\n")
                        cat("RMSE：", round(additional_metrics$RMSE, 2), " - 二乗平均平方根誤差（万円）\n")
                        cat("MAPE：", round(additional_metrics$MAPE, 2), "% - 平均絶対パーセンテージ誤差\n")
                    

📊 統計指標の重要性ランキング

📈 モデル適合度の総合評価

【モデル適合度の総合評価】 R²値： 0.9524 - 優秀（説明力80%以上）調整済みR²： 0.9514 - より保守的な評価 F統計量： 966.87 - モデル全体が極めて有意 AIC： 2287.43 - モデル選択の指標（小さい方が良い） BIC： 2304.17 - より保守的なモデル選択指標 MAE： 393.21 - 平均絶対誤差（万円） RMSE： 493.86 - 二乗平均平方根誤差（万円） MAPE： 8.73 % - 平均絶対パーセンテージ誤差

🔬 仮説検定の科学

仮説検定は、データから得られた結果が偶然によるものか、統計的に有意な違いなのかを判断する強力な手法です。Tidyverseの力を借りて、複雑な統計的推論を直感的に行うことができます。

🎯 t検定による平均値比較

hypothesis_testing.R

                        # 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)
                    

🎯 カイ二乗検定による独立性検定

chi_square_test.R

                        # マーケティングキャンペーンデータ
                        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))
                    

🚀 高度なモデリング技法

現代の統計モデリングは、単純な線形関係を超えた複雑なパターンを捉える必要があります。正則化手法、ロジスティック回帰、時系列分析など、実際のビジネス課題に対応する高度な手法を習得しましょう。

🎯 ロジスティック回帰による分類

logistic_regression.R

                        # 顧客離反予測データ
                        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")
                    

🔬 ロジスティック回帰の詳細解釈

📊 オッズ比の解釈

変数	係数	オッズ比	95%CI	解釈
tenure_months	-0.048	0.953	[0.94, 0.97]	1ヶ月増加で離反確率4.7%減少
monthly_charges	0.019	1.019	[1.01, 1.03]	1円増加で離反確率1.9%増加
support_tickets	0.287	1.332	[1.22, 1.45]	1件増加で離反確率33.2%増加
contract_Monthly	0.956	2.601	[2.12, 3.19]	月契約は年契約比2.6倍離反しやすい

🎯 平均限界効果（Average Marginal Effects）

オッズ比は解釈しやすいですが、実際の確率変化を知るには限界効果が有用です。

marginal_effects.R

                            # 平均限界効果の計算
                            library(margins)
                            
                            # 限界効果の手動計算
                            calculate_marginal_effects <- function(model, data) {
                              # 予測確率を取得
                              pred_probs <- predict(model, type = "response")
                              
                              # 限界効果 = β * p(1-p) （ロジスティック分布の密度）
                              coefficients <- coef(model)
                              marginal_effects <- list()
                              
                              # 各変数の限界効果
                              marginal_effects$tenure_months <- 
                                mean(coefficients["tenure_months"] * pred_probs * (1 - pred_probs))
                              
                              marginal_effects$monthly_charges <- 
                                mean(coefficients["monthly_charges"] * pred_probs * (1 - pred_probs))
                              
                              marginal_effects$support_tickets <- 
                                mean(coefficients["support_tickets"] * pred_probs * (1 - pred_probs))
                              
                              return(marginal_effects)
                            }
                            
                            # 限界効果の計算と表示
                            me <- calculate_marginal_effects(logit_model, churn_data)
                            
                            # 結果の整理
                            marginal_effects_df <- tibble(
                              variable = names(me),
                              marginal_effect = unlist(me),
                              percentage_change = round(unlist(me) * 100, 3)
                            )
                            
                            print(marginal_effects_df)
                            
                            # 実際的な解釈例
                            cat("\\n=== 平均限界効果の実務的解釈 ===\\n")
                            cat("契約期間1ヶ月延長 → 離反確率", round(me$tenure_months * 100, 2), "%ポイント減少\\n")
                            cat("月額料金1,000円増加 → 離反確率", round(me$monthly_charges * 1000 * 100, 2), "%ポイント増加\\n")
                            cat("サポート問合せ1件増加 → 離反確率", round(me$support_tickets * 100, 2), "%ポイント増加\\n")
                        

🎪 限界効果の実務的意味

契約期間延長の効果

限界効果: -1.2%ポイント

意味: 平均的な顧客で、契約期間が1ヶ月伸びると離反確率が1.2%ポイント下がる。

戦略的含意: 長期契約インセンティブの設計が重要。

料金変更の影響

限界効果: +0.48%ポイント（100円当たり）

意味: 月額料金を1,000円上げると離反確率が4.8%ポイント増加。

戦略的含意: 価格感応度の高い顧客セグメントの特定が必要。

サポート品質の重要性

限界効果: +7.1%ポイント

意味: サポート問合せが1件増えると離反確率が7.1%ポイント増加。

戦略的含意: 初回解決率向上とプロアクティブサポートが鍵。

🎨 具体的な予測シナリオ

prediction_scenarios.R

                                # 具体的な顧客プロファイルでの予測
                                scenarios <- tibble(
                                  scenario = c("新規顧客", "満足顧客", "リスク顧客", "問題顧客"),
                                  tenure_months = c(3, 24, 6, 12),
                                  monthly_charges = c(45, 65, 95, 110),
                                  support_tickets = c(1, 0, 2, 5),
                                  contract_type = c("Monthly", "Two year", "One year", "Monthly")
                                )
                                
                                # 各シナリオの離反確率を予測
                                scenarios <- scenarios %>%
                                  mutate(
                                    churn_probability = predict(logit_model, newdata = ., type = "response"),
                                    risk_level = case_when(
                                      churn_probability < 0.3 ~ "低リスク",
                                      churn_probability < 0.6 ~ "中リスク",
                                      TRUE ~ "高リスク"
                                    )
                                  )
                                
                                print(scenarios)
                            

顧客タイプ	契約期間	月額料金	サポート件数	契約形態	離反確率	リスクレベル
新規顧客	3ヶ月	4,500円	1件	月契約	67.2%	高リスク
満足顧客	24ヶ月	6,500円	0件	2年契約	8.1%	低リスク
リスク顧客	6ヶ月	9,500円	2件	1年契約	42.6%	中リスク
問題顧客	12ヶ月	11,000円	5件	月契約	89.3%	高リスク

🎯 モデル検証とクロスバリデーション

model_validation.R

                        # クロスバリデーション用のライブラリ
                        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))
                    