Bootstrap weights

[alanlujan91]

Add weights kwarg to handle survey weights

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Files with missing lines	Coverage Δ
src/estimagic/bootstrap.py	99.06% <100.00%> (ø)
src/estimagic/bootstrap_helpers.py	100.00% <100.00%> (ø)
src/estimagic/bootstrap_outcomes.py	100.00% <100.00%> (ø)
src/estimagic/bootstrap_samples.py	100.00% <100.00%> (ø)
src/estimagic/msm_weighting.py	98.27% <ø> (ø)

[timmens]

Hey Alan! Thanks for opening another pull request.

Since this time, the PR is not related to a minor fix but an extension, could you extend the PR description? In particular, could you state the problem (i.e. the missing feature), how this is called in the bootstrap literature, and how you plan to implement it. If you know other libraries that implement this feature, this would also be great to know.

[timmens]

Nice work so far! Thank you!

As mentioned above, I think we need to go through a few things first. If implementing these changes will take too much of your time, you can also re-open this pull request as an issue (feature request). Then we will think about how to implement the feature.

[timmens]

You only check that the get_bootstrap_indices function works when called with the new weights arguments. But I think you can do more:

1. I would rename "wts" to "weights", just to be explicit
2. If the weights are uniform (as in your example), the bootstrap indices with and without the weights argument should be the same. Could you test that? You will need to create two random number generators with the same state.
3. Additionally, it would be very helpful to see if the weights work. An extreme case to check this would be a weighting vector where one entry is 1, and the other is 0. Then, you could (probabilistically) check that the index corresponding to the 0 weight is never chosen.

[timmens]

It would be more readable if you compute the vector of probabilities outside of these if conditions. If weights is None, the vector of probabilities could be None. In this case, you could also write the same code for both cases:

probs = _calculate_bootstrap_indices_probs(data, weight_by, cluster_by) 

bootstrap_indices = list(
    rng.choice(
        n_obs,
        size=(n_draws, n_obs),
        replace=True,
        p=probs,
    )
)

This will be slower than np.integers(n_obs, ...) but I think this is a trade-off we would want to make here.

[timmens]

This has to be updated. Especially the line clusters[clusters_by] is very confusing. I would think that clusters is equal to something like data[cluster_by].

Similarly as above, it would be nice to pre-compute the probabilities.

[alanlujan91]

Hey Alan! Thanks for opening another pull request.

Since this time, the PR is not related to a minor fix but an extension, could you extend the PR description? In particular, could you state the problem (i.e. the missing feature), how this is called in the bootstrap literature, and how you plan to implement it. If you know other libraries that implement this feature, this would also be great to know.

Will try to address all of this in the coming days! Thanks!

[timmens]

Hey @alanlujan91, thanks again for working on this feature. We are planning a major refactoring and want to close as many PRs as possible beforehand.

Are you planning to work on this PR in the near future?

[alanlujan91]

hi @timmens

I know this feature as proportional sampling, which is something we have done in Econ-ARK but I don't have a good link to the literature.

https://en.wikipedia.org/wiki/Probability-proportional-to-size_sampling

I can work on this more if you think this is a desired feature, I've just been pretty busy so this hasn't been a priority.

[janosg]

Hi Alan, thanks for the PR. I think this would be nice to have in estimagic. As Tim said, we are planning major refactorings and would like to have few PRs open when we start. Do you think you can finish within the next three weeks? Otherwise it might be a good idea to convert this to an issue and postpone the actual implementation until you have more free time.

[alanlujan91]

@janosg I can definitely finish this within 3 weeks

[alanlujan91]

Modified code according to above comments from @timmens, except for

If the weights are uniform (as in your example), the bootstrap indices with and without the weights argument should be the same. Could you test that? You will need to create two random number generators with the same state.

This one is complicated because of the following note from numpy

"Setting user-specified probabilities through p uses a more general but less efficient sampler than the default. The general sampler produces a different sample than the optimized sampler even if each element of p is 1 / len(a)."

It can be done, but requires generating dummy p's as p = np.ones(n_obs)/n_obs when weights_by is None. Let me know if this is desirable

[timmens]

The changes look really good, thank you! I only have a few minor remarks and two small request, and then we can merge 🎉

Request 1

The function _get_probs_for_bootstrap_indices is not tested yet and the docstring is not that informative. Can you

1. Extend the docstring by explaining what to expect in the different cases (i.e., weight_by is not None, cluster_by is not None, etc.)
2. Write unit tests for these cases. Looking at the function, it is not directly clear to me that the case where both weight_by and cluster_by are not None is correctly handled. With an example in a unit test, this would be much easier to see.

Request 2

It can be done, but requires generating dummy p's as p = np.ones(n_obs)/n_obs when weights_by is None. Let me know if this is desirable

Ah, I see; very good point! In this case, can you add a structural test like this:

def test_get_bootstrap_indices_heterogeneous_weights():
    data = pd.DataFrame({"id": [0, 1], "w_homogenous": [0.5, 0.5], "w_heterogenous": [0.1, 0.9]})
    
    res_homogenous = get_bootstrap_indices(data, weight_by="w_homogenous", n_draws=1_000, rng=get_rng(seed=0))
    res_heterogenous = get_bootstrap_indices(data, weight_by="w_heterogenous", n_draws=1_000, rng=get_rng(seed=0))
    
    # Given the weights, the first sample mean should be close to 0.5, while the second one should be close to 0.9
    assert np.mean(res_homogenous) < 0.75 < np.mean(res_heterogenous)

[janosg]

@alanlujan91 do you think it makes sense to keep this PR open or should we revisit the topic some other time?

[alanlujan91]

Sorry I took so long to get back to this. I finally addressed all of the comments, but now it seems there's some unrelated issue to this PR making the checks fail

[janosg]

Thanks @alanlujan91! Yes, this is unrelated and comes from a new release of DFO-LS. I'll fix it soon.