Rの関数:lm_robust {estimatr}

# パッケージの読み込み
library(ggplot2)
library(dplyr)
library(texreg) # 結果比較用

# 乱数シードの固定
seed <- 20251218
set.seed(seed)

# 設定: 小売チェーン店の広告効果分析
# N: 総観測数
# n_stores: 店舗数（クラスター数）
n_stores <- 50
n_obs_per_store <- 20
N <- n_stores * n_obs_per_store

# データフレーム作成
df_sim <- data.frame(
  store_id = rep(1:n_stores, each = n_obs_per_store)
)

# 1. 説明変数 X (広告費): 一様分布
df_sim$ad_spend <- runif(N, min = 10, max = 100)

# 2. クラスター効果 (店舗ごとの固有要因)
# 店舗ごとにベースの売上が異なる
store_effect <- rnorm(n_stores, mean = 0, sd = 10)
df_sim$store_intercept <- store_effect[df_sim$store_id]

# 3. 誤差項 (不均一分散 + クラスター相関)
# 不均一分散: 広告費が多いほど、売上のブレ（分散）が大きくなる
# クラスター相関: 同じ店舗内では誤差の方向が似る（ここでは簡易的にstore_interceptで表現済みだが、さらにノイズを加える）
# heteroskedastic_error: Xが大きいほどSDが大きい
hetero_error <- rnorm(N, mean = 0, sd = df_sim$ad_spend * 0.5)

# 4. 被説明変数 Y (売上)
# True Model: Sales = 500 + 5 * Ad_Spend + Store_Effect + Error
true_beta <- 5
df_sim$sales <- 500 + (true_beta * df_sim$ad_spend) + df_sim$store_intercept + hetero_error

cat("--- データ概要 ---\n")
cat(sprintf("観測数: %d (店舗数: %d)\n", N, n_stores))

# モデル1: 通常の lm (等分散・独立を仮定)
model_ols <- lm(sales ~ ad_spend, data = df_sim)

# モデル2: lm_robust (デフォルト = HC2: 不均一分散に対処)
# クラスターは指定せず、個体ごとの不均一分散のみ考慮
model_robust_hc2 <- lm_robust(sales ~ ad_spend, data = df_sim, se_type = "HC2")

# モデル3: lm_robust (クラスターロバスト = CR2: クラスター相関に対処)
# 店舗IDをクラスターとして指定
model_robust_cluster <- lm_robust(sales ~ ad_spend, data = df_sim, clusters = store_id, se_type = "CR2")

cat("--- 推定結果の比較 ---")
# texregを使って見やすく並べる
screenreg(
  list(model_ols, model_robust_hc2, model_robust_cluster),
  custom.model.names = c("通常のlm", "lm_robust(HC2)", "lm_robust(Cluster)"),
  include.ci = FALSE,
  digits = 3
)

cat("\n--- 標準誤差の比較 ---\n")

# 係数と標準誤差の抽出
se_ols <- summary(model_ols)$coefficients["ad_spend", "Std. Error"]
se_hc2 <- summary(model_robust_hc2)$coefficients["ad_spend", "Std. Error"]
se_cl <- summary(model_robust_cluster)$coefficients["ad_spend", "Std. Error"]

cat(sprintf("1. 通常のlmの標準誤差       : %.4f\n", se_ols))
cat(sprintf("2. 頑健(HC2)の標準誤差      : %.4f\n", se_hc2))
cat(sprintf("3. クラスター頑健の標準誤差 : %.4f\n\n", se_cl))

# fixed_effects 引数を利用
model_fe <- lm_robust(sales ~ ad_spend, data = df_sim, fixed_effects = ~store_id, clusters = store_id)

print(summary(model_fe))