diff --git a/R/test_attrition.R b/R/test_attrition.R index a68e43c..e87e2ae 100644 --- a/R/test_attrition.R +++ b/R/test_attrition.R @@ -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") @@ -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)) @@ -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( @@ -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") @@ -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") @@ -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") @@ -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) @@ -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. @@ -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\".") diff --git a/Stata/robust_attrition.ado b/Stata/robust_attrition.ado index 27f9204..0ac2c1c 100644 --- a/Stata/robust_attrition.ado +++ b/Stata/robust_attrition.ado @@ -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" @@ -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 { @@ -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 @@ -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)) /// @@ -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