Skip to content

Add interaction plot to robust_attrition() in R and Stata #129

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jrosen48 merged 6 commits into from
Aug 6, 2026
Merged

Add interaction plot to robust_attrition() in R and Stata #129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
a2fd196
Add interaction plot to robust_attrition() in R and Stata
JihoonChoi26 Jul 24, 2026
7b09b50
Update test_attrition.R
JihoonChoi26 Jul 27, 2026
359fa9d
Update test_attrition.R
JihoonChoi26 Jul 27, 2026
0a6a62f
Update test_attrition.R
JihoonChoi26 Jul 27, 2026
1229293
Update test_attrition.R
JihoonChoi26 Jul 29, 2026
b78ef19
robust_attrition() update
JihoonChoi26 Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
159 changes: 129 additions & 30 deletions R/test_attrition.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -376,26 +376,19 @@ robust_attrition <- function(
corr_tvalidatep <- (a_root - j_root * corr_rlimp) /
sqrt(z_root - w_root * corr_rlimp^2)
corr_tvalidaten <- (a_root - j_root * corr_rlimn) /
sqrt(z_root - w_root * corr_rlimn^2)

sqrt(z_root - w_root * corr_rlimn^2)
# ===================================================================
# Output
# Printed report (shared by to_return = "print" and "plot")
# ===================================================================
if (to_return == "print") {

# -----------------------------------------------------------
# Header
# -----------------------------------------------------------
print_report <- function() {
cat(crayon::bold("[BETA] Robustness to Differential Attrition\n\n"))

# -----------------------------------------------------------
# Verbose opening (context)
# -----------------------------------------------------------
if (isTRUE(verbose)) {
cat("This analysis quantifies what would have to be true in the\n")
cat("attritted (missing) data to nullify an inference based on the\n")
cat("observed data. Two complementary blocks are reported:\n")
cat(" (1) Nonparametric: the treatment effect in the missing data\n")
cat(" (1) Nonparametric: the treatment estimate in the missing data\n")
cat(" required to bring the combined-sample effect to the\n")
cat(" significance threshold.\n")
cat(" (2) Correlation-based: the correlation (and implied\n")
Expand All @@ -404,9 +397,6 @@ robust_attrition <- function(
cat(" from the attrition x treatment interaction.\n\n")
}

# -----------------------------------------------------------
# Observed inputs
# -----------------------------------------------------------
cat(sprintf("Observed treatment mean: %.4f\n", yobt))
cat(sprintf("Observed control mean: %.4f\n", yobc))
cat(sprintf("Estimated effect (observed): %.4f\n", esteffect))
Expand All @@ -418,9 +408,6 @@ robust_attrition <- function(
ntreatmi, 100 * alphat,
ncontrolmi, 100 * alphac))

# -----------------------------------------------------------
# Verbose notes on internally-computed defaults
# -----------------------------------------------------------
if (isTRUE(verbose)) {
if (is.null(R2xz)) {
cat(sprintf(
Expand All @@ -438,9 +425,6 @@ robust_attrition <- function(

cat("\n")

# -----------------------------------------------------------
# Nonparametric block
# -----------------------------------------------------------
cat("Nonparametric:\n")
if (isTRUE(verbose)) {
cat(" (effect in missing data required to pull the\n")
Expand All @@ -452,9 +436,6 @@ robust_attrition <- function(
nonpar_deltacomb))
cat("\n")

# -----------------------------------------------------------
# Correlation-based block
# -----------------------------------------------------------
cat("Correlation-based:\n")
if (isTRUE(verbose)) {
cat(" (correlation in missing data required to nullify,\n")
Expand All @@ -474,9 +455,6 @@ robust_attrition <- function(
))
cat("\n")

# -----------------------------------------------------------
# Footer
# -----------------------------------------------------------
if (isTRUE(verbose)) {
cat("For all computed quantities, run with to_return = \"raw_output\".\n")
cat("See ?robust_attrition for help.\n")
Expand All @@ -485,6 +463,14 @@ robust_attrition <- function(
cat("For all computed quantities, run with to_return = \"raw_output\".\n")
}
cat("\n")
}

# ===================================================================
# Output
# ===================================================================
if (to_return == "print") {

print_report()

# -----------------------------------------------------------
# Invisible return (same structure as raw_output)
Expand Down Expand Up @@ -617,6 +603,10 @@ robust_attrition <- function(

} else if (to_return == "plot") {

# Print the same report shown by to_return = "print", honoring verbose,
# then build and return the two plots.
print_report()

# x-window: centered on the observed pooled mean (yob), the mean
# assumed for the missing data in the original calculation. Half-width
# defaults to pad_frac (15%) of |yob|; raise pad_frac for a longer line.
Expand Down Expand Up @@ -691,11 +681,120 @@ robust_attrition <- function(
ggplot2::theme_bw() +
ggplot2::theme(
plot.title.position = "plot",
plot.title = ggplot2::element_text(size = 14),
plot.subtitle = ggplot2::element_text(size = 12)
plot.title = ggplot2::element_text(size = 16),
plot.subtitle = ggplot2::element_text(size = 14),
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 13),
strip.text = ggplot2::element_text(size = 13),
legend.text = ggplot2::element_text(size = 13),
legend.position = "bottom"
)

# ---- Second plot: cell means by group and attrition status ----
# Four cell means from Y = b0 + b1*M + b2*T + b3*M*T. The observed
# (M = 0) cells are shared; the attritted (M = 1) cells differ between
# the mean-difference and correlation-based (parametric) frameworks.
means_df <- data.frame(
framework = factor(
rep(c("Mean-difference", "Correlation-based"), each = 4),
levels = c("Mean-difference", "Correlation-based")
),
group = factor(
rep(c("Treatment", "Treatment", "Control", "Control"), times = 2),
levels = c("Treatment", "Control")
),
M = rep(c(0, 1, 0, 1), times = 2),
Y = c(
yobt, nonpar_ytreatmi, yobc, nonpar_ycontrolmi,
yobt, corr_ytreatmi, yobc, corr_Ycontrolmi
)
)
# push each value label vertically away from the lines: the upper
# point in a column goes up, the lower point goes down. Within each
# (framework, M) column, compare the two group means.
means_df$vj_lab <- ave(
means_df$Y, means_df$framework, means_df$M,
FUN = function(v) ifelse(v == max(v), -0.8, 1.8)
)

# delta indicators: observed effect (shared) and the framework-specific
# attritted effect, drawn just outside the M = 0 and M = 1 columns.
gap_df <- data.frame(
framework = factor(
rep(c("Mean-difference", "Correlation-based"), each = 2),
levels = c("Mean-difference", "Correlation-based")
),
x = rep(c(-0.08, 1.08), times = 2),
ylo = c(pmin(yobt, yobc), pmin(nonpar_ytreatmi, nonpar_ycontrolmi),
pmin(yobt, yobc), pmin(corr_ytreatmi, corr_Ycontrolmi)),
yhi = c(pmax(yobt, yobc), pmax(nonpar_ytreatmi, nonpar_ycontrolmi),
pmax(yobt, yobc), pmax(corr_ytreatmi, corr_Ycontrolmi)),
val = c(esteffect, nonpar_deltami, esteffect, corr_rbetami),
sym = c("ob", "mi", "ob", "mi")
)
gap_df$hj <- ifelse(gap_df$x < 0, 1.15, -0.15)
gap_df$txt <- sprintf("hat(delta)^{%s} == %.3f", gap_df$sym, gap_df$val)

p_means <- ggplot2::ggplot(
means_df,
ggplot2::aes(x = M, y = Y, color = group, group = group)
) +
# treatment-control gap (double-headed) behind the lines
ggplot2::geom_segment(
data = gap_df, inherit.aes = FALSE,
ggplot2::aes(x = x, xend = x, y = ylo, yend = yhi),
color = "grey45", linewidth = 0.5,
arrow = grid::arrow(ends = "both",
length = grid::unit(0.05, "in"))
) +
ggplot2::geom_text(
data = gap_df, inherit.aes = FALSE, parse = TRUE,
ggplot2::aes(x = x, y = (ylo + yhi) / 2, label = txt, hjust = hj),
color = "grey20", size = 4.2
) +
ggplot2::geom_line(linewidth = 1.1) +
ggplot2::geom_point(size = 2.6) +
ggplot2::geom_text(
ggplot2::aes(label = sprintf("%.3f", Y), vjust = vj_lab),
color = "grey20", size = 4.2, show.legend = FALSE
) +
ggplot2::facet_wrap(~ framework) +
ggplot2::scale_color_manual(
values = c("Treatment" = "#1F78B4", "Control" = "#33A02C")
) +
ggplot2::scale_x_continuous(
breaks = c(0, 1),
labels = c("Observed\n(M = 0)", "Attritted\n(M = 1)"),
limits = c(-0.55, 1.55)
) +
ggplot2::scale_y_continuous(
expand = ggplot2::expansion(mult = 0.09)
) +
ggplot2::labs(
title = expression(
paste("Treatment estimates (", hat(delta)^{mi},
") in missing (attritted) data necessary to nullify the inference")
),
subtitle = paste0(
"Observed estimate on the left side (M = 0) and missing right side (M = 1).\n",
"The left side is the observed treatment estimate; the right side is the treatment estimate \n",
"in the missing data necessary to nullify the inference."
),
x = "Missing status, M",
y = "Outcome mean, Y",
color = NULL
) +
ggplot2::theme_bw() +
ggplot2::theme(
plot.title.position = "plot",
plot.title = ggplot2::element_text(size = 16),
plot.subtitle = ggplot2::element_text(size = 14),
legend.position = "bottom"
)

return(p)
print(p)
print(p_means)
return(invisible(list(robustness = p, means = p_means)))

} else {
stop("to_return must be one of \"print\", \"raw_output\", or \"plot\".")
Expand Down
105 changes: 98 additions & 7 deletions Stata/robust_attrition.ado
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -354,7 +354,7 @@ program define robust_attrition, rclass
di as txt "This analysis quantifies what would have to be true in the"
di as txt "attritted (missing) data to nullify an inference based on"
di as txt "the observed data. Two complementary blocks are reported:"
di as txt " (1) Nonparametric: the treatment effect in the missing"
di as txt " (1) Nonparametric: the treatment estimate in the missing"
di as txt " data required to bring the combined-sample effect"
di as txt " to the significance threshold."
di as txt " (2) Correlation-based: the correlation (and implied"
Expand Down Expand Up @@ -433,9 +433,12 @@ program define robust_attrition, rclass
di ""

if (`graph') {
di as txt "See the figure: effect in the missing data required to nullify"
di as txt "the inference, across assumed missing-data means. The red point"
di as txt "marks the original calculation (Y# = observed pooled mean)."
di as txt "See the two figures. The first shows the effect in the missing"
di as txt "data required to nullify the inference, across assumed"
di as txt "missing-data means; the red point marks the original"
di as txt "calculation (Y# = observed pooled mean). The second shows the"
di as txt "outcome means by group and attrition status, under both the"
di as txt "nonparametric and correlation-based frameworks."
di ""
}
else {
Expand Down Expand Up @@ -505,6 +508,33 @@ program define robust_attrition, rclass
local s_slope = strtrim(string(`b_slp', "%12.3f"))
local lgd3 = "Y# = `s_center', delta_mi = `s_delta', slope = `s_slope'"

* ---- values for the second plot (cell means) ----
* Four cell means from Y = b0 + b1*M + b2*T + b3*M*T. The observed
* (M = 0) cells are shared; the attritted (M = 1) cells differ between
* the mean-difference and correlation-based (parametric) frameworks.
local m_nptre = `np_ytreatmi'
local m_npctl = `np_ycontrolmi'
local m_cotre = `corr_ytreatmi'
local m_coctl = `corr_Ycontrolmi'
local m_dob = `esteffect'
local m_dnp = `np_deltami'
local m_dco = `corr_rbetami'

* shared y-window across both panels, with a small margin
local m_ymin = min(`yobt', `yobc', `m_nptre', `m_npctl', `m_cotre', `m_coctl')
local m_ymax = max(`yobt', `yobc', `m_nptre', `m_npctl', `m_cotre', `m_coctl')
local m_pad = 0.09 * (`m_ymax' - `m_ymin')
if (`m_pad' >= . | `m_pad' <= 0) local m_pad = 1
local m_ylo = `m_ymin' - `m_pad'
local m_yhi = `m_ymax' + `m_pad'

local m_sob = strtrim(string(`m_dob', "%12.3f"))
local m_snp = strtrim(string(`m_dnp', "%12.3f"))
local m_sco = strtrim(string(`m_dco', "%12.3f"))

* delta with a combining circumflex (U+0302), i.e. delta-hat
local dhat = ustrunescape("\u03B4\u0302")

quietly {
preserve
clear
Expand All @@ -514,11 +544,11 @@ program define robust_attrition, rclass

twoway ///
(line yline_y xline_x, lcolor(navy) lwidth(medthick)) ///
(pci `dmi' `x_lo' `dmi' `center', lpattern(dash) lcolor(red) lwidth(thin)) ///
(pci `y_winlo' `center' `dmi' `center', lpattern(dash) lcolor(red) lwidth(thin)) ///
(pci `dmi' `x_lo' `dmi' `x_hi', lpattern(dash) lcolor(red) lwidth(thin)) ///
(scatteri `dmi' `center', mcolor(red) msymbol(O) msize(medium)) ///
, ///
yline(0, lcolor(gs11) lwidth(thin)) ///
xline(`center', lpattern(dash) lcolor(red) lwidth(thin)) ///
yscale(range(`y_winlo' `y_winhi') noextend) ///
xscale(range(`x_lo' `x_hi') noextend) ///
ylabel(, angle(horizontal)) ///
Expand All @@ -527,14 +557,75 @@ program define robust_attrition, rclass
subtitle("Effect the missing data would need to nullify the observed inference," "by assumed mean outcome in the missing data", size(medsmall)) ///
xtitle("Assumed mean outcome in the missing data (Y#)") ///
ytitle("Required effect in the missing data (delta_mi)") ///
text(`dmi' `center' "Original calculation", placement(e) color(red) size(small)) ///
text(`dmi' `center' "Original calculation", placement(se) color(red) size(small)) ///
text(`lgd_y' `lgd_x' "Red point: original calculation (Y# = pooled mean)" "Line slope: rise in delta_mi per unit of Y#" "`lgd3'", placement(se) box fcolor(white) lcolor(black) margin(small) size(small) just(left)) ///
legend(off) ///
graphregion(color(white)) ///
name(attrition_plot, replace)

* ---- Second plot: cell means by group and attrition status ----
clear
set obs 4
gen byte fw = cond(_n <= 2, 1, 2)
gen double mstat = cond(mod(_n, 2) == 1, 0, 1)
label define _fwlbl 1 "Mean-difference" 2 "Correlation-based", replace
label values fw _fwlbl

gen double y_t = `yobt' if mstat == 0
gen double y_c = `yobc' if mstat == 0
replace y_t = `m_nptre' if fw == 1 & mstat == 1
replace y_c = `m_npctl' if fw == 1 & mstat == 1
replace y_t = `m_cotre' if fw == 2 & mstat == 1
replace y_c = `m_coctl' if fw == 2 & mstat == 1

* treatment-control gap, drawn just outside each column
gen double gap_x = cond(mstat == 0, -0.08, 1.08)
gen double gap_lo = min(y_t, y_c)
gen double gap_hi = max(y_t, y_c)
gen double gap_mid = (gap_lo + gap_hi) / 2
gen byte gappos = cond(mstat == 0, 9, 3)

gen str50 gaplab = "`dhat'{superscript:ob} = `m_sob'" if mstat == 0
replace gaplab = "`dhat'{superscript:mi} = `m_snp'" if fw == 1 & mstat == 1
replace gaplab = "`dhat'{superscript:mi} = `m_sco'" if fw == 2 & mstat == 1

* nudge value labels off the crossing lines
gen byte pos_t = cond(mstat == 0, 12, 6)
gen byte pos_c = cond(mstat == 0, 6, 12)
gen str12 lab_t = string(y_t, "%9.3f")
gen str12 lab_c = string(y_c, "%9.3f")

twoway ///
(rcap gap_hi gap_lo gap_x, lcolor(gs7) lwidth(medium)) ///
(scatter gap_mid gap_x, msymbol(none) mlabel(gaplab) ///
mlabcolor(black) mlabsize(small) mlabvposition(gappos)) ///
(line y_t mstat, lcolor("31 120 180") lwidth(medthick)) ///
(line y_c mstat, lcolor("51 160 44") lwidth(medthick)) ///
(scatter y_t mstat, mcolor("31 120 180") msymbol(O) msize(medium) ///
mlabel(lab_t) mlabcolor(black) mlabsize(small) mlabvposition(pos_t)) ///
(scatter y_c mstat, mcolor("51 160 44") msymbol(O) msize(medium) ///
mlabel(lab_c) mlabcolor(black) mlabsize(small) mlabvposition(pos_c)) ///
, ///
by(fw, ///
title("Treatment estimates (`dhat'{superscript:mi}) in missing (attritted)" "data necessary to nullify the inference", size(medsmall)) ///
subtitle("Observed estimate on the left side (M = 0) and missing right side (M = 1)." "The left side is the observed treatment estimate; the right side is the treatment estimate" "in the missing data necessary to nullify the inference.", size(small)) ///
note("") ///
iscale(*1.1) ///
legend(position(6)) ///
) ///
legend(order(5 "Treatment" 6 "Control") position(6) rows(1)) ///
yscale(range(`m_ylo' `m_yhi')) ///
ylabel(, angle(horizontal) labsize(small)) ///
xscale(range(-0.45 1.45)) ///
xlabel(0 "Observed (M = 0)" 1 "Attritted (M = 1)", labsize(small)) ///
xtitle("Missing status, M", size(small)) ///
ytitle("Outcome mean, Y", size(small)) ///
graphregion(color(white)) ///
name(attrition_means, replace)
restore
}

graph display attrition_plot
graph display attrition_means
}
end
Loading