install.packages("https://github.com/grf-labs/sufrep/blob/master/sufrep_0.1.0.tar.gz?raw=true", repos = NULL, type = "source")
require(grf)
require(maq)
library(ggplot2)
require(sufrep)
library(forcats)
set.seed(20260409)
data <- read.csv("nhefs.csv")
# 結果変数の欠測データを取り除いた。
data <- data[!is.na(data$wt82_71),]
# 死亡有無を結果変数、禁煙を処置変数、それ以外の8変数(pretreatment)を特徴量とする
Y = data[,c("wt82_71")]
#Y = sign(Y)*log(1+abs(Y))
W = data[,c("qsmk")]
X = data[,c("sex", "race", "age", "education", "smokeintensity",
"smokeyrs", "exercise", "active", "wt71")]
#X$age2 = X$age**2
#X$smokeintensity2 = X$smokeintensity**2
#X$smokeyrs2 = X$smokeyrs**2
#X$wt71_2 = X$wt71**2
cont_val <- c("sex", "race", "age", "smokeintensity","smokeyrs", "wt71")
#"age2", "smokeintensity2","smokeyrs2", "wt71_2"
# factor 変数に対してmean encodingを行う
fact_val <- c("education","exercise", "active")
X$education <- as.factor(X$education)
X$exercise <- as.factor(X$exercise)
X$active <- as.factor(X$active)
X_newfactor <- fct_cross(X$education, X$exercise, X$active)
#means Encoding
encoder <- make_encoder(X[cont_val], X_newfactor, method="means")
X <- encoder(X[cont_val], X_newfactor)
# サンプル数
n = dim(data)[[1]]
p = dim(X)[[2]]
# 単純な2群の比較
summary(lm(Y ~ W))
# ランダムフォレストによる因果効果（喫煙が体重に与える影響）の推定結果
ps.full = probability_forest(X, as.factor(W), num.trees = 2000, min.node.size = 10)
llf_Y = ll_regression_forest(X, Y, num.trees = 2000, min.node.size = 5, enable.ll.split = TRUE, honesty = TRUE)
Y.hat = predict(llf_Y, X, linear.correction.variables = 1:p)
cf.full = causal_forest(X, Y, W, W.hat = ps.full$predictions[,2], Y.hat = Y.hat$predictions)
average_treatment_effect(cf.full)
# 推定された傾向スコアのヒストグラム
# オーバーラップが存在するために、0および1から離れていることが必要
hist(cf.full$W.hat, xlab = "Estimated propensity scores", main = "")
# 結果変数モデルの当てはまり
plot(cf.full$Y.hat + (W - cf.full$W.hat) * cf.full$predictions,Y)
# バイアスの診断プロット
par(mfcol=c(1,2))
hist((Y- cf.full$Y.hat - (W - cf.full$W.hat)*cf.full$predictions)[W==1]/(W-cf.full$W.hat)[W==1],
main = "Treated (W=1) Group",
xlab = "Inverse probability weighted Y-residual")
hist((Y- cf.full$Y.hat - (W - cf.full$W.hat)*cf.full$predictions)[W==0]/(W-cf.full$W.hat)[W==0],
main = "Untreated (W=0) Group",
xlab = "Inverse probability weighted Y-residual")
par(mfcol=c(1,1))
# 半分のサンプル（学習データ）で Causal Forest を学習させる
# 傾向スコアとアウトカム回帰モデルはランダムフォレスト
# Split data into a train and test sample.
train = sample(nrow(X), 0.6 * nrow(X))
test = -train
# Causal Forest の当てはめ
ps.train = probability_forest(X[train, ], as.factor(W[train]), num.trees = 2000, min.node.size = 10)
llf_Y.train = ll_regression_forest(X[train,], Y[train], num.trees = 2000, min.node.size = 5, enable.ll.split = TRUE, honesty = TRUE)
Y.hat.train = predict(llf_Y.train, X[train,], linear.correction.variables = 1:p)
cate.forest = causal_forest(X[train, ], Y[train], W[train],
W.hat = ps.train$predictions[,2], Y.hat = Y.hat.train$predictions)
# 半分のサンプル（評価用データ）で予測を計算する
X.test = X[test, ]
tau.hat.test = predict(cate.forest, X.test)$predictions
# CATEの推定量のヒストグラム
hist(tau.hat.test, xlab = "Estimated CATEs", main = "")
# CATE の異質性をうまく捉えられているかを確認する
num.groups = 3 # 4 for quartiles, 5 for quintiles, etc.
quartile = cut(tau.hat.test,
quantile(tau.hat.test, seq(0, 1, by = 1 / num.groups)),
labels = 1:num.groups,
include.lowest = TRUE)
# CATEの分位点ごとに評価用データを分割する
samples.by.quartile = split(seq_along(quartile), quartile)
# 各分位点ごとのCATE を計算するために、評価用データにCausal Forest を当てはめる
ps.eval = probability_forest(X.test, as.factor(W[test]), num.trees = 2000, min.node.size = 10)
llf_Y.eval = ll_regression_forest(X.test, Y[test], num.trees = 2000, min.node.size = 5, enable.ll.split = TRUE, honesty = TRUE)
Y.hat.eval = predict(llf_Y.train, X.test, linear.correction.variables = 1:p)
eval.forest = causal_forest(X.test, Y[test], W[test],
W.hat = ps.eval$predictions[,2],
Y.hat = Y.hat.eval$predictions)
# 各グループごとにATEを計算する
ate.by.quartile = lapply(samples.by.quartile, function(samples) {
average_treatment_effect(eval.forest, subset = samples)
})
# 各グループごとのATEと95%信頼区間を計算
df.plot.ate = data.frame(
matrix(unlist(ate.by.quartile), num.groups, byrow = TRUE, dimnames = list(NULL, c("estimate","std.err"))),
group = 1:num.groups
)
ggplot(df.plot.ate, aes(x = group, y = estimate)) +
geom_point() +
geom_errorbar(aes(ymin = estimate - 1.96 * std.err, ymax = estimate + 1.96 * std.err, width = 0.2)) +
xlab("Estimated CATE quantile") +
ylab("Average treatment effect")
# *** TOC/AUTOC を用いた異質性評価 ***
# TOC/AUTOC は分位点をなめらかに繋いだもの
# eval.forest を用いて、TOC/AUTOCに対する2重頑健推定量を計算する
rate.cate = rank_average_treatment_effect(
eval.forest,
tau.hat.test,
q = seq(0.05, 1, length.out = 100)
)
# Plot the TOC.
plot(rate.cate)
# An estimate and standard error of AUTOC.
print(rate.cate)
# Get a 2-sided p-value Pr(>|t|) for RATE = 0 using a t-value.
2 * pnorm(-abs(rate.cate$estimate / rate.cate$std.err))
# *** Qini曲線による介入の評価 ****
# Form a doubly robust estimate of the Qini curve, using R=200 bootstrap replications for std.errors.
cost = 1
qini = maq(tau.hat.test, cost, get_scores(eval.forest), R = 200)
# Form a baseline Qini curve that assigns treatment uniformly.
qini.baseline = maq(tau.hat.test, cost, get_scores(eval.forest), R = 200,
target.with.covariates = FALSE)
# 分位点ごとの介入効果の利得（平均的な増加量）
plot(qini, xlim = c(0, 1), xlab = "quantile")
plot(qini.baseline, add = TRUE, lty = 2, ci.args = NULL) # leave out CIs for legibility.
# Plot the Qini curve with 95% confidence intervals.
# This shows the expected gain from withholding units from deployments in order of predicted benefit.
# 禁煙が 2000 人の場合にどのようになるか？
# `maq_scale()` を用いて Qini 曲線を現実のスケールに拡大する
max.deployment = 2000
# Qini曲線をプロットする（スケール版）
scale_maq(qini, max.deployment) |>
plot(ylab = "Amount of Weights Inflation",
xlab = "Units held back from deployment")
# ベースラインをプロット（スケール版）
scale_maq(qini.baseline, max.deployment) |>
plot(add = TRUE, ci.args = NULL)
# 500人が禁煙した場合の体重増加量の合計を、曲線から推定する
average_gain(scale_maq(qini, max.deployment), 500)
# ランダムに禁煙した場合ど、禁煙による体重増加が多い人から禁煙した場合の差
difference_gain(scale_maq(qini, max.deployment),
scale_maq(qini.baseline, max.deployment), 500)
# *** CATEの推定関数と説明変数の関係性 ****
# Our `cate.forest` has given us some estimated function \hat tau().
# Let's have a closer look at how this function stratifies our sample in terms of "covariate" profiles.
# One way to do so is to look at histograms of our covariates by for example low / high CATE predictions.
# First, we'll use a simple heuristic to narrow down the number of predictors to look closer at.
# Here we use the variable importance metric of the fit CATE function to select 4 predictors to look closer at.
varimp.cate = variable_importance(cate.forest)
ranked.variables = order(varimp.cate, decreasing = TRUE)
top.varnames = colnames(X)[ranked.variables[1:4]]
print(top.varnames)
# Select the test set samples predicted to have low/high CATEs.
# [We could also have used the full sample for this exercise.]
low = samples.by.quartile[[1]]
high = samples.by.quartile[[num.groups]]
# Make some long format data frames for ggplot.
df.lo = data.frame(
covariate.value = unlist(as.vector(X.test[low, top.varnames])),
covariate.name = rep(top.varnames, each = length(low)),
cate.estimates = "Low"
)
df.hi = data.frame(
covariate.value = unlist(as.vector(X.test[high, top.varnames])),
covariate.name = rep(top.varnames, each = length(high)),
cate.estimates = "High"
)
df.plot.hist = rbind(df.lo, df.hi)
# Plot overlaid histograms of the selected covariates by low/high classification.
ggplot(df.plot.hist, aes(x = covariate.value, fill = cate.estimates)) +
geom_histogram(alpha = 0.7, position = "identity") +
facet_wrap(~ covariate.name, scales = "free", ncol = 2)
# *** Best linear projections (BLP) ****
# Select some potential effect modifier(s) we are interested in.
blp.vars = c("age", "smokeintensity", "smokeyrs","wt71")
# Estimate the best linear projection on our variables.
for(i in c(1:length(blp.vars))){
print(best_linear_projection(cf.full, X[, blp.vars[i]]))
}
print(best_linear_projection(cf.full, X[, blp.vars]))
cor(X[, blp.vars])