Is Chorussell really working correctly?

[a-torgovitsky]

This does not pass the sanity test.

library("lpinfer")
source("~/R/x86_64-pc-linux-gnu-library/3.6/lpinfer/example/dgp_mixedlogit.R")

set.seed(5)
dgp <- mixedlogit_dgp()
df <- mixedlogit_draw(dgp, n = 4000)

lpm <- lpmodel(A.obs = mixedlogit_Aobs(dgp),
               beta.obs = function(d) mixedlogit_betaobs(d, dgp),
               A.shp = rep(1, nrow(dgp$vdist)),
               beta.shp = 1,
               A.tgt = mixedlogit_Atgt_dfelast(dgp, w2eval = 1, eeval = -1))

set.seed(5)
r <- chorussell(data = df, lpmodel = lpm, ci = TRUE, beta.tgt = .2)
print(r)

lb <- sort(r$lb.bs[[1]])
ub <- sort(r$ub.bs[[1]])

print(min(lb))
print(max(ub))

print(mean(lb < .637))

produces

95%-confidence interval: [0.63708, Inf]                                                        
[1] 0.2177688
[1] 0.7402449
[1] 0.51

The largest bootstrap estimate of the upper bound is .74, yet the upper bound of the confidence interval is +Inf (or 1 if we use the logical confidence interval)? More than half of the bootstrap estimates are smaller than the lower bound of the confidence interval?

This can't be correct. Can you please walk through the calculation step-by-step here? Something is clearly wrong.

[conroylau]

I went through the implementation notes and the paper again. I am sorry that I think I have confused myself with the sign for the candidates of the upper bound earlier. The code has also been updated to reflect the change.

The updated confidence interval that I get from the above example [0.63392, 0.88549]. However, the same problem described above still appears and I am still trying to find the reason behind it. Please find a report that explains my calculations for your review here. Thank you!

[a-torgovitsky]

It's probably something wrong either with the computational approach I described in my note, or with the way it was implemented.

Maybe the best thing to do is try a simpler, brute force computational procedure first? The objective is to choose c{lb} and c{ub} to minimize the length of the resulting confidence interval (c{lb} + c{ub}). Why not just build a 2 dimensional grid. Iterate over all values in the grid. If (c{lb}, c{ub}) meets the constraints, then record their sum. Otherwise don't record it. Then just look at the smallest recorded sum.

Compare that to the result you get with the current, more complicated procedure.

If the results approximately match, then there's either something more fundamentally wrong about how we are interpreting Cho-Russell, or else there is something fundamentally wrong with their procedure. If the results do not match, then there's something wrong with our more complicated computational procedure.

[conroylau]

Thanks for the suggestions. I have just implemented this and added the results to the document here. Please refer to the last section "Obtaining the confidence interval without step 4" on page 4 to 5.

The confidence interval obtained with the brute force approach is the same as the ones obtained with the refinement step. Hence, the problems are still here. I think this rules out the case where the error is coming from the refinement step.

I will go back to read the paper again to see if I can get any hints from what might be causing the confusing results. Thanks!

[a-torgovitsky]

Ok, thanks, let me know what you find. It is possible that we both interpreted the paper incorrectly. It is also possible that their method is just weird.

By the way, an excellent unit test that you could add is to check whether you get the same answer with remove.const = FALSE as with remove.const = TRUE.

[conroylau]

I have not yet find out the reason behind the confusing results, but here are some updates on my attempts:

• I find that the Cho-Russell paper has imposed a threshold rule on the delta term (length of the identified set) in their paper. I tried to do the same as the delta term of our example is quite close to 0 as well. However, we still have the same problems.
• I tried other specifications like trying the 1-norm or setting another value of kappa, but this also do not help.

Please see here for the updated document with more details on my attempts (I added the new contents on page 6 and 7). I will continue try to see what's the reason behind the issue.

In addition, I have included an extra unit test to make sure the same answers are obtained in remove.const = FALSE as with remove.const = TRUE.

Thakns!

[a-torgovitsky]

What about considering a simpler example with 10 bootstraps and alpha = .10? That would help us understand exactly how c{lb} and c{ub} arise since we would only have 10 pairs of lower/upper bounds to look at.

[conroylau]

Thanks for the suggestions. I tried to run a simpler example with R = 10 and alpha = .1 and the problems still appear. Please find the details of the simulation and a brief analysis here. I find that the chorussell procedure is picking the smallest feasible lower and upper bounds of the confidence interval in this simulation exercise, although the smallest possible lower bound of the CI is still higher than half of bootstrap lower bounds (as seen from the figure).

[a-torgovitsky]

Ok, I think I am satisfied that this makes sense. I think this is something about their procedure in cases where lb and ub are approximately the same. But that's not our problem! So I'll close this issue now. Thanks!

[conroylau]

I see. Sure, thank you!