Upload case study on causal IV analyses

[joonho112]

Hello, I've uploaded a case study document titled "Instrumental Variables Analysis of Randomized Experiments with One-Sided Noncompliance." This study is a collaborative work with Profs. Avi Feller and Sophia Rabe-Hesketh from UC Berkeley. The document is a result of our work on the IES project, which you can review here: IES Project Link. I would greatly appreciate it if you could approve this merge request. Thank you.

[bob-carpenter]

Hi, @joonho112 and thanks for submitting. I think all we're waiting for is to clean up the checked-in html.

[WardBrian]

I don’t know if we have anyone who is familiar with the subject matter to do a more thorough review.

One small thing, but it would be nice to use the new array syntax in the Stan code. RStan just updated this morning to allow it, so it should be reasonable to update the case study and keep using RStan for it, while making the code usable by cmdstanr/py users

[bob-carpenter]

I just gave this a read. I'm not an expert on causal inference, but this looks fine to me. Everything below is mostly about some Stan details and formatting. I'm OK letting the authors fix as much or as little of this as they want then posting it.

I would recommend is a bit more inter-line spacing in sequences of equations (i.e., putting in \\[4pt] or something similar for extra vertical space.

For Stan code, we recommend always including spaces after commas, like in lower=0, upper=1 and always including space before and after parens, like in if (Z[n].

You can make the Stan code more efficient by avoiding recomputation of log(pi_c) and log(1 - pi_c) (which should be coded as log1m(pi_c) in Stan for arithmetic stability). Just define local variables outside of the loop and define

real log_pi_c = log(pi_c);
real log1m_pi_ci = log1m(pi_c);

It can be sped up a lot by partitioning the data so that the treatment/control and compiler/never-taker sets are split so that the target increments can be vectorized and the Bernoullis are unfolded with prior computation of the 1 and 0 outcome probabilities. But if it's already fast enough, I wouldn't worry about any of this other than making sure to use log1m(pi_c) for arithmetic stability. I covered all of this in the efficiency chapter of the User's Guide.

Personally, I'm often inclined to transform these models from the natural parameterization in terms of probabilities to log odds to make hierarchical modeling easier, but that'd take you way outside of just replicating what Rubin and Imbens did.

You could've used sampling notation instead of target +=:

target += log(pi_c) + bernoulli_lpmf(Y[n] | eta_c1);

can be replaced with the equivalent code

1 ~ bernoulli(pi_c);
Y[n] ~ bernoulli(eta_c1);

We have a built-in for mixtures, that can simplify

target += log_sum_exp(
        log(1 - pi_c) + bernoulli_lpmf(Y[n] | eta_n),  // Never-taker
        log(pi_c) + bernoulli_lpmf(Y[n] | eta_c0)      // Complier
      ); 

to

target += log_mix(pi_c, bernoulli_lpmf(Y[n] | eta_c0), bernoulli_lpmf(Y[n] | eta_n));

Note that I reversed the order to match pi_c.

Rather than this

  real CACE; 
  CACE = (eta_c1 - eta_c0)*10^3;

they can be combined into one line,

  real CACE = (eta_c1 - eta_c0) * 10^3;

Same for the later example that includes NACE.

We recommend spaces around all operators other than ^.

The model fit very nicely with NUTS drawing roughly independent samples.

Aki's current recommendation on R-hat is something like 1.05 or 1.01. It's a moving target, but 1.1 is definitely too high (though we discount the ESS when R-hat is high).

The replicated figure is very nice.

While a 90% interval of (1.28, 5.17) does exclude zero, what do the tails look like? What's the 99% interval (that would require many more MCMC runs to estimate accurately, though)? Can you say anything about the scale of 1 to 5 in terms of the effect. It's positive, but does it matter? Do we know what inter-person variation is?

It'd be nice to have the Stan program names like Model_02_CACE_without_Exclusion_Restriction.stan point to the GitHub repo or to an appendix that spells them out in full. I'd strongly recommend trying to keep file names shorter. I'd have gone with something like cace_wo_exclusion.stan and I definitely wouldn't have included the Model_02, which appears to be an attempt to get the files to sort in some order under ls (I'd suggest just putting all the models in a subdirectory called stan).

I would have run the stan_fit_noER until the R-hat on lp__ was below 1.01. You can't really trust the R-hat on the other parameters if the lp__ variable hasn't converged (it's a non-linear function of all the other variables).

I find the red on grey very hard to read.

I would've preferred disks with alpha << 1 for the scatterplot rather than crosses.

[joonho112]

Hello @bob-carpenter and @WardBrian,

Thank you so much for your invaluable comments and feedback. I have carefully reviewed and incorporated Bob's suggestions into the updated case study. For your convenience, I am attaching the revised case study with this comment. Additionally, I've taken steps to remove the extraneous case-studies.html file as you pointed out.

Please do let me know if there are any other aspects that require my attention before you proceed with the publication of this case study. I truly appreciate your guidance in refining this work. cace_one-sided.html.zip

[joonho112]

I've updated the Stan models using the new array syntax.

[joonho112]

I have made the adjustments to the Stan code, moving the codes from the "transformed parameters" section to the "generated quantities" section. Additionally, I have removed all .DS_Store files from the repository.

Could you proceed to merge this Pull Request and publish this case study? I am looking forward to the final sign-off from @WardBrian or @bob-carpenter.

[bob-carpenter]

I just merged @joonho112 (sorry if this gives you a double notice---I wasn't sure you'd see the merge).