jkcshea / ivmte

An R package for implementing the method in Mogstad, Santos, and Torgovitsky (2018, Econometrica).
GNU General Public License v3.0
18 stars 2 forks source link

Using custom weights to reproduce a LATE #220

Closed johnnybonney closed 2 years ago

johnnybonney commented 2 years ago

To preface - I'm not completely sure if this is a package error, or if I'm just specifying custom weights incorrectly (in theory or in practice).

I am trying to use custom weight/knot functions to estimate LATEs. As a starting point, I wanted to reproduce LATEs as estimated by ivmte---but I've been getting very different results. Consider the following example.

# fix a few arguments throughout
args <- list(data = AE,
             outcome = 'worked',
             propensity = morekids ~ samesex,
             m0 = ~u + I(u^2),
             m1 = ~u + I(u^2),
             solver = 'gurobi')

# estimate using package functionality
late_ivmte <- do.call(ivmte, c(args, list(target = 'late',
                                          late.from = list(samesex = 0),
                                          late.to = list(samesex = 1))))

yields bounds [-0.2995747, 0.1298702]. I then try to estimate the same LATE with custom weights, with the understanding that the weights should be $1[u \in (p(0), p(1)]] \times [p(1) - p(0)]^{-1}$. So, the weights I want are $[p(1) - p(0)]^{-1}$ for $u \in (p(0), p(1)]$ and 0 everywhere else:

# estimate manually
pmod <- glm(morekids ~ samesex, family = "binomial", AE)
p.from <- predict(pmod, newdata = data.table(samesex = 0), type = "response")
p.to <- predict(pmod, newdata = data.table(samesex = 1), type = "response")

# custom weights, knots
late.wt <- c(0, 1 / (p.to - p.from), 0)
late.knots <- c(p.from, p.to)

late_manual <- do.call(ivmte, c(args, list(target.weight0 = late.wt,
                                           target.weight1 = late.wt,
                                           target.knots0 = late.knots,
                                           target.knots1 = late.knots)))

but then the bounds I get are [0.7972684, 1.2264]---wildly different... Am I goofing something up here, or is ivmte reacting differently to two equivalent ways of articulating the same problem?

jkcshea commented 2 years ago

Hi @johnnybonney,

Sorry for the slow reply, especially because the solution is very simple. You just forgot to use the negative weight for the untreated group. Once you set target.weight0 = -late.wt, you should replicate the original LATE estimate.

johnnybonney commented 2 years ago

Ah, of course. I'm getting rusty... Thanks, Josh!

One follow-up question - if I specify target = "late" when covariates are included, the target is $E[LATE(X)]$, is that right? And not $E[Y{1} - Y{0} | U \in (p(0,X), p(1,X)]]$? I think you've answered this for me before, but it's been so long.

a-torgovitsky commented 2 years ago

It should be the unconditional LATE So $LATE = E[Y(1) - Y(0) \vert \text{cp}] = E[E[Y(1) - Y(0) \vert \text{cp}, X] | \text{cp}] \equiv E[LATE(X) \vert \text{cp}]$ i.e. it's going to adjust for compliance rates by $X$ See Section 4.3.1 in the paper: https://a-torgovitsky.github.io/shea-torgovitsky.pdf

[It would be useful/safe if you can check it in a simple example to make sure it works as expected.]

johnnybonney commented 2 years ago

I believe that it is returning $E[LATE(X)]$ and not $E[Y(1) - Y(0)|\text{cp}]$ (which is what I was referring to with $E[Y(1) - Y(0)|U \in (p(0,X), p(1,X)]]$).

I could be wrong(!!), but in your linked paper, the RHS of $LATE(z{0} \rightarrow z{1})$ (the second-to-last expression on page 15) actually appears to be equal to $E[LATE(X)]$. The term within that expectation is $LATE(X)$, I think... I believe the unconditional LATE would reweight by $[p(X,z{1}) - p(X, z{0})] / E[p(X,z{1}) - p(X, z{0})]$. But let me know if I'm mistaken.

I think the package is doing exactly what the paper says, and returning $E[LATE(X)]$. (But if my understanding of the theory is wrong above, then all bets are off.)

Take my initial example and modify by fully interacting the black covariate.

args <- list(data = AE,
             outcome = 'worked',
             propensity = morekids ~ samesex + black + black:samesex,
             m0 = ~black + u + I(u^2) + black:u + black:I(u^2),
             m1 = ~black + u + I(u^2) + black:u + black:I(u^2),
             solver = 'gurobi')

# estimate using package functionality
late_ivmte <- do.call(ivmte, c(args, list(target = 'late',
                                          late.from = list(samesex = 0),
                                          late.to = list(samesex = 1))))

Package estimation yields bounds [-0.2565717, 0.08778715].

Knots should be covariate-specific and at the propensity score values. Per my understanding of the theory, the unconditional LATE weight for the applicable region should just be $1 / E[p(X,1) - p(X,0)]$ for the $d = 1$ MTR, and the same negative for $d = 0$.

# estimate manually
pmod <- glm(morekids ~ samesex + black + black:samesex, family = "binomial", AE)
p.from <- predict(pmod,
                  newdata = data.frame(black = 0:1, samesex = 0),
                  type = "response")
p.to <- predict(pmod,
                newdata = data.frame(black = 0:1, samesex = 1),
                type = "response")
pr.x <- c(1 - mean(AE$black), mean(AE$black))

# custom weights, knots to reflect unconditional LATE
late.wt <- c(0, 1 / sum(pr.x * (p.to - p.from)), 0)
knot1 <- function(black) p.from[(black + 1)]
knot2 <- function(black) p.to[(black + 1)]

late_manual <- do.call(ivmte, c(args, list(target.weight0 = -late.wt,
                                           target.weight1 = late.wt,
                                           target.knots0 = c(knot1, knot2),
                                           target.knots1 = c(knot1, knot2))))

This yields bounds [-0.2502527, 0.08424642]---close to the "canned" LATE bounds, but not the same. However, I can reproduce the package's LATE bounds by using the E[LATE(X)] weights:

# custom weights, knots to reflect E[LATE(X)]
exlate.wt0 <- function(black) -1 / (knot2(black) - knot1(black))
exlate.wt1 <- function(black) 1 / (knot2(black) - knot1(black))

exlate_manual <- do.call(ivmte, c(args, list(target.weight0 = c(0, exlate.wt0, 0),
                                             target.weight1 = c(0, exlate.wt1, 0),
                                             target.knots0 = c(knot1, knot2),
                                             target.knots1 = c(knot1, knot2))))

which yields bounds [-0.2565717, 0.08778715], just like ivmte originally produced.

a-torgovitsky commented 2 years ago

Whoops! Yes you're right about the theory. We need to fix that in the paper. @jkcshea did we code up both target parameters or just this "population-averaged" LATE? @johnnybonney I'm now remembering that we had some discussion about this way back...whether it was more interesting to look at population-averaged LATEs or the unconditional LATE. Do you remember this?

johnnybonney commented 2 years ago

Yeah, I vaguely remember that discussion, but I don't remember the conclusion... I tried to find it in the resolved GitHub issues, but no luck (so maybe this was an off-GitHub discussion).

It seems to me that the unconditional LATE would be (almost?) always more interesting to consider than the population-averaged LATEs - it's hard for me to come up with a case when the population-averaged LATEs would be more useful/interesting than the unconditional LATE.

The unconditional LATE also strikes me as closer to what a researcher might be expecting if they specify target = "late".

a-torgovitsky commented 2 years ago

I definitely agree with this:

The unconditional LATE also strikes me as closer to what a researcher might be expecting if they specify target = "late".

I think this is a mistake we need to fix @jkcshea .

My first instinct is also that the population-averaged LATE was more interesting...but I feel like in one of the applications we were working on were using both for some reason. I believe we called one of them the "natural" weighting and one of them the "population" weighting. Does that ring any bells? It's not a Github thing...it's somewhere in our files/results I think.

johnnybonney commented 2 years ago

OK, yes, that rings some bells. If I'm remembering correctly, the question was, what should we compare to the bounds? The TSLS estimate, a "population weighted average" of covariate-specific TSLS estimates, or the unconditional LATE? And I believe we decided when the instrument was multivalued, the population-weighted average of subgroup TSLS estimates was a better comparison point than the unconditional LATE.

But then when we wrote the first draft of the local LATEs paper, we just compared everything to the TSLS estimate anyway, to keep things uniform across applications.

Also

My first instinct is also that the population-averaged LATE was more interesting

You meant to write the unconditional LATE here, right?


While we are talking about LATEs, @cblandhol and I have one more question (rather than open a new issue...) In the simple cases (such as my original example), LATEs are nonparametrically point-identified, so we don't need to extrapolate. They used to be point-identified in earlier versions of ivmte as well, but not any more. Is that to be expected, given the regression-based approach?

a-torgovitsky commented 2 years ago

You meant to write the unconditional LATE here, right?

Yes, sorry.


While we are talking about LATEs, @cblandhol and I have one more question (rather than open a new issue...) In the simple cases (such as my original example), LATEs are nonparametrically point-identified, so we don't need to extrapolate. They used to be point-identified in earlier versions of ivmte as well, but not any more. Is that to be expected, given the regression-based approach?

I think it depends. What moments were you matching before? (The moment-matching approach where you specify ivlike should not have changed since previous versions, so presumably you will still get the same thing.) Were you able to "exactly match" the moments before (minimum criterion of 0)?

johnnybonney commented 2 years ago

I am just talking about simple cases, no covariates, where you use the sharp information set, E[Y|D,Z]. Consider my example from the initial post, where ivmte gave bounds of [-0.2995747, 0.1298702] via the regression-based approach. If I add args$ivlike <- list(worked ~ morekids + samesex + morekids:samesex) and re-run (via the moment-matching approach), I get [-0.08484221, -0.08484221] with a minimum criterion of 0.

a-torgovitsky commented 2 years ago

Can you collate the entire example for me. I'm not able to reproduce your initial numbers:

library("ivmte")

# fix a few arguments throughout
args <- list(data = AE,
             outcome = 'worked',
             propensity = morekids ~ samesex,
             m0 = ~u + I(u^2),
             m1 = ~u + I(u^2),
             solver = 'gurobi')

# estimate using package functionality
late_ivmte <- do.call(ivmte, c(args, list(target = 'late',
                                          late.from = list(samesex = 0),
                                          late.to = list(samesex = 1))))
print(late_ivmte)

I'm getting

Bounds on the target parameter: [-0.2758405, 0.1061561]
Audit terminated successfully after 1 round 
jkcshea commented 2 years ago

@jkcshea did we code up both target parameters or just this "population-averaged" LATE?

No, unfortunately not. This is my fault---sorry! It is very easy to include both the unconditional LATE and the population-average LATE. Is that something we want to do?

a-torgovitsky commented 2 years ago

I think we want to do this

jkcshea commented 2 years ago

Okay, this has now been fixed. But I'll push it after we resolve the discussion in #221. Thanks so much for catching this, @johnnybonney.

Something to mention: The bounds I get for the unconditional LATE are [-0.250187, 0.08424317]. This differs slightly from the late_manual bounds I get when running @johnnybonney's code here, which are [-0.2503054, 0.08424609]. This discrepancy arises because the propensity scores estimated by ivmte differ from those manually estimated by < 1e-16. I'm surprised such a small difference in the objective leads to a perceptible difference in the bounds.

Note also my late_manual bounds differ from those of @johnnybonney's. Hopefully that will be resolved as part of #221.

a-torgovitsky commented 2 years ago

I would say that's a barely perceptible difference. It's probably because the propensity score is in the denominator, so small differences are going to get blown up. And small differences in the propensity score are because the logit is probably being solved slightly differently on different machines.

johnnybonney commented 2 years ago

Edit: this probably has nothing to do with this issue, given that the behavior disappears when presolve is turned off, but I'm leaving this comment here just in case.

Here's another example that we were discussing in #225.

dt <- AE[AE$yob < 55,]

args <- list(data = dt,
             outcome = "worked",
             propensity = morekids ~ factor(samesex) * factor(yob),
             m1 = ~factor(yob)*u + factor(yob):I(u^2),
             m0 = ~factor(yob)*u + factor(yob):I(u^2),
             solver = "gurobi",
             solver.options = list(presolve = 0),
             direct = "lp",
             initgrid.nu = 50,
             audit.nu = 200)  

# package estimation, LATE(yob = 53)
est_ivmte <- do.call(
  ivmte,
  c(args, list(target = "late",
               late.to = list(samesex = 1),
               late.from = list(samesex = 0),
               late.X = list(yob = 53)))
)
print(est_ivmte)
Bounds on the target parameter: [-0.2190124, -0.2614573]
Audit terminated successfully after 2 rounds 

while manually specifying weights/knots gives

# manual estimation, LATE(yob = 53)
p1_x53 <- mean(dt[dt$samesex == 1 & dt$yob == 53,]$morekids)
p0_x53 <- mean(dt[dt$samesex == 0 & dt$yob == 53,]$morekids)
Pr_x53 <- mean(dt$yob == 53)

late_wt_m1 <- function(yob) as.integer(yob == 53) / ((p1_x53 - p0_x53) * Pr_x53)
late_wt_m0 <- function(yob) -1 * late_wt_m1(yob)

est_manual <- do.call(
  ivmte,
  c(args, list(target.weight0 = c(0, late_wt_m0, 0),
               target.weight1 = c(0, late_wt_m1, 0),
               target.knots0 = c(p0_x53, p1_x53),
               target.knots1 = c(p0_x53, p1_x53)))
)
print(est_manual)
Bounds on the target parameter: [-0.2190124, -0.3006759]
Audit terminated successfully after 2 rounds 

However, if I specify solver.options = list(presolve = 0), then I get the same bounds for both ways of specifying the problem: [-0.2520733, -0.258758].

a-torgovitsky commented 2 years ago

Thanks let's leave the comment. @jkcshea just as a reminder to check whether covariates have been incorporated correctly when adjusting the LATE target parameters. The late would be $E[Y(1) - Y(0) \vert \text{cp}, V = v]$ While the avglate would be $E[E[Y(1) - Y(0) \vert \text{cp}, X, V = v] | V = v]$ where $V = v$ is the late.X thing (I'm using the notation in the paper here)

jkcshea commented 2 years ago

Thanks for the example conditional on late.X. Using the updated code and having presolve = 0, my bounds match those of @johnnybonney's. (But as in @johnnybonney's example, the bounds match whether we choose the late or avglate target in ivmte since we are not integrating over covariates after conditioning on yob == 53. Let me know if I'm misunderstanding, though.)

I'll look into what's going on with presolve. We've incorporated a simple switch for it in the past (i.e., solver.presolve = FALSE), but never dived into what exactly it does. I know one of Gurobi's developers has written a paper on presolve for MILPs. Presolve.pdf

a-torgovitsky commented 2 years ago

Ah right it's not a good example for that reason. When you adjusted the code though did you check that it would work with late.X in the correct way?

jkcshea commented 2 years ago

I did not... :( And sadly it's not working with late.X. I'll keep working on this.

a-torgovitsky commented 2 years ago

Ok, it should just be a matter of getting the conditioning correct in the weights. Let me know if you want me to sketch it out.

jkcshea commented 2 years ago

Yes, I think I spotted where my mistake was. I didn't correctly condition on late.X at a certain point in the code.

To make sure I'm doing things correctly, here are the weights I'll be implementing: conditional-late.pdf Let me know if I got something wrong!

a-torgovitsky commented 2 years ago

I think you've switched the role of X and Z from what we have in the paper? https://a-torgovitsky.github.io/shea-torgovitsky.pdf Section 4.3.1 Can you switch it back? Otherwise I'm surely going to get confused checking the algebra

jkcshea commented 2 years ago

Ah yes, I'm very sorry. Here's the corrected version conditional-late.pdf .

a-torgovitsky commented 2 years ago

Thanks @jkcshea I read through your note and it all looks right to me.

I guess we will also want $E[LATE(X) \vert V = v]$ for the avglate parameter, but the weights for this should be more straightforward (no Bayes rule trickery).

jkcshea commented 2 years ago

Done!

Here are some examples to check the implementation of $\mathbb{E}[\mathit{LATE}(X) \mid V = v]$ for the late and avglate parameters. conditional-late.zip

I haven't merged this code to the master branch, though. More than once, while playing with the other issues, I was confused why I wasn't able to replicate the examples in the thread. The reason was because I was on the branch used to debug this issue. And on that branch, target = 'late' no longer had the same meaning as in the code from the examples. So I'll withhold merging the code until this batch of stability issues are resolved. But if someone wants to use this code, it's available on the branch bugfix/issue220.

I'm leaving this branch open as a reminder for myself to merge that branch into the master branch.

jkcshea commented 2 years ago

Since the other issues seem to be related to #227/stability, I've merged the changes from this issue into the master branch and pushed them.