Rで線形回帰における各種区間の算出

library(ggplot2)
library(dplyr)

# === データ生成 (シナリオに基づいたシミュレーション) ===
seed <- 20250623
set.seed(seed)

# シナリオ設定
n <- 50 # データ収集日数
true_intercept <- 40000 # 経験0年のバリスタのみの場合の基礎売上 (4万円)
true_slope <- 5000 # 平均経験年数が1年増えるごとの売上増加額 (5千円)
error_sd <- 8000 # 日々のランダムな売上変動の標準偏差 (8千円)

# 架空のデータを生成
experience <- runif(n, 0.5, 8) # バリスタの平均経験年数を0.5年から8年の間でランダムに生成
true_sales <- true_intercept + true_slope * experience
sales <- true_sales + rnorm(n, mean = 0, sd = error_sd) # ランダムな誤差を加える

# データフレームにまとめる
cafe_data <- data.frame(experience, sales)


# === 線形回帰モデルの構築 ===
model <- lm(sales ~ experience, data = cafe_data)
model_summary <- summary(model)


# === 4つの区間の算出 ===

# 1. 母回帰係数（傾きと切片）の信頼区間
coef_ci <- confint(model, level = 0.95)


# 予測用のx軸データを生成 (グラフを滑らかにするため)
pred_grid <- data.frame(experience = seq(min(experience), max(experience), length.out = 100))

# 2. 応答の平均値の信頼区間 (t分布)
mean_conf_interval <- predict(model, newdata = pred_grid, interval = "confidence", level = 0.95) %>%
  as.data.frame() %>%
  rename(fit_conf = fit, lwr_conf = lwr, upr_conf = upr)

# 3. 目的変数の予測区間
pred_interval <- predict(model, newdata = pred_grid, interval = "prediction", level = 0.95) %>%
  as.data.frame() %>%
  rename(fit_pred = fit, lwr_pred = lwr, upr_pred = upr)

# 4. ワーキング-ホテリング信頼帯 (F分布)
f_value <- qf(0.95, 2, model$df.residual)
wh_multiplier <- sqrt(2 * f_value)
se_fit <- predict(model, newdata = pred_grid, se.fit = TRUE)$se.fit
wh_lwr <- mean_conf_interval$fit_conf - wh_multiplier * se_fit
wh_upr <- mean_conf_interval$fit_conf + wh_multiplier * se_fit

# グラフ描画用にすべてのデータを結合
plot_data <- bind_cols(pred_grid, mean_conf_interval, pred_interval) %>%
  mutate(wh_lwr = wh_lwr, wh_upr = wh_upr)


# === ggplotによる可視化 ===
ggplot() +
  # 予測区間 (一番外側・最も広い)
  geom_ribbon(data = plot_data, aes(x = experience, ymin = lwr_pred, ymax = upr_pred, fill = "Prediction"), alpha = 0.2) +
  # ワーキング-ホテリング信頼帯 (中間)
  geom_ribbon(data = plot_data, aes(x = experience, ymin = wh_lwr, ymax = wh_upr, fill = "Working-Hotelling"), alpha = 0.3) +
  # 応答の平均値の信頼区間 (一番内側・最も狭い)
  geom_ribbon(data = plot_data, aes(x = experience, ymin = lwr_conf, ymax = upr_conf, fill = "Confidence"), alpha = 0.5) +

  # 回帰直線
  geom_line(data = plot_data, aes(x = experience, y = fit_conf), color = "red", linewidth = 1) +

  # 元データの散布図
  geom_point(data = cafe_data, aes(x = experience, y = sales), color = "darkred", size = 2) +

  # ラベルとデザイン
  scale_fill_manual(
    name = "95% Interval / Band",
    values = c("Prediction" = "lightcoral", "Working-Hotelling" = "green", "Confidence" = "navy"),
    labels = c("応答の平均値の信頼区間", "予測区間", "ワーキング-ホテリング信頼帯")
    #    labels = c("Confidence (Mean)", "Prediction (New Day)", "W-H Band (Whole Line)")
  ) +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title = "カフェの売上とバリスタの経験年数の関係",
    subtitle = "3種類の区間の比較",
    x = "バリスタの平均経験年数 (年)",
    y = "一日の総売上 (円)"
  ) +
  theme_bw() +
  theme(legend.position = "top")