Rの関数:horvitz_thompson {estimatr}

# パッケージの読み込み
library(ggplot2)
library(dplyr)

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

# サンプルサイズ
N <- 1000

# 顧客データの生成

# 顧客タイプ: 30%が優良顧客、70%が一般顧客
type_preferred <- 0.3
type_general <- 0.7

df_sim <- data.frame(
  id = 1:N,
  type = sample(c("優良顧客", "一般顧客"), N, replace = TRUE, prob = c(type_preferred, type_general))
)

# 割り付け確率の設定
# 優良顧客は 80% の確率でクーポン配布 (Z=1)
# 一般顧客は 20% の確率でクーポン配布 (Z=1)
coupon_preferred <- 0.8
coupon_general <- 0.2

df_sim$prob_assign <- ifelse(df_sim$type == "優良顧客", coupon_preferred, coupon_general)

# 処置（クーポン配布）の実行
df_sim$treatment <- rbinom(N, size = 1, prob = df_sim$prob_assign)

# 潜在的な購入額（Potential Outcomes）の生成
# ベースライン購入額: 優良顧客は平均5000円、一般顧客は平均2000円
baseline_preferred <- 5000
baseline_general <- 2000
base_amount <- ifelse(df_sim$type == "優良顧客", baseline_preferred, baseline_general) + rnorm(N, 0, 500)

# 真の処置効果（クーポンの効果）: 一律 +1000円 とする
true_effect <- 1000

# 観測される購入額 Y
# 処置群 (Z=1) なら +1000円 される
df_sim$outcome <- base_amount + (true_effect * df_sim$treatment)

# 処置変数のラベル化
df_sim$treatment_label <- ifelse(df_sim$treatment == 1, "クーポンあり", "クーポンなし")

cat("--- データ概要 ---\n")
cat(sprintf("総観測数: %d 人\n", N))
cat(sprintf("割り付け確率: 優良顧客 = %d%%, 一般顧客 = %d%% \n", coupon_preferred * 100, coupon_general * 100))
cat(sprintf("ベースライン購入額: 優良顧客 = %d円, 一般顧客 = %d円\n", baseline_preferred, baseline_general))
cat(sprintf("真のクーポン効果: +%d 円\n\n", true_effect))

# 確率の違いを無視して差を計算
naive_model <- lm(outcome ~ treatment, data = df_sim)
naive_est <- coef(naive_model)["treatment"]

print(summary(naive_model))
cat("------------------------------------------------------------\n")
cat(sprintf("単純比較による推定効果: %.1f 円\n", naive_est))

# 割り付け確率 (condition_prs) を明示して補正を行う
ht_model <- horvitz_thompson(
  formula = outcome ~ treatment,
  data = df_sim,
  condition_prs = prob_assign, # 各個体の割り付け確率を指定
  se_type = "youngs"
)

# 結果の表示
print(ht_model)

cat("------------------------------------------------------------\n")

# 結果の抽出
ht_est <- ht_model$coefficients
ht_se <- ht_model$std.error
ht_ci_low <- ht_model$conf.low
ht_ci_high <- ht_model$conf.high

cat(sprintf("Horvitz-Thompson 推定値: %.1f 円\n", ht_est))
cat(sprintf("標準誤差: %.1f\n", ht_se))
cat(sprintf("95%%信頼区間: [%.1f, %.1f]\n", ht_ci_low, ht_ci_high))