ENH: Return full dataset for subsequent analysis

[nickeubank]

As per discussion with @nehargupta via email, it'd be helpful if dame returned an analysis dataset. results is a start, but it drops treatment assignment and outcome, and the insertion of "*" values means it can't actually be used for analysis.

Basically, most people I know who do matching in the social sciences expect the matching package to return a data set with all the matched pairs. In 1 to 1 matching, that data set is just a subset of the original data set for which matches were found. In many to one matching, observations that are matched more than once end up being repeated in the data set but with weights that reflect that fact. you can see this in the matchit docs here: https://kosukeimai.github.io/MatchIt/articles/MatchIt.html#estimating-the-treatment-effect (you may have to scroll up a little for context).

Here's a super-crude implementation of what I'm looking for, where result_of_fit is the output of the match (the "result"), and model is my fitted model.

def get_dataframe(model, result_of_fit):

    # Get original data
    better = model.input_data.loc[result_of_fit.index]

   # Get match groups for clustering
   matched_data["match_group"] = np.nan
   for idx, group in enumerate(model.units_per_group):
       matched_data.loc[group, "match_group"] = idx

    # Weights
    better["weights"] = 1 / model.groups_per_unit

    if not model.repeats:
        assert (better["weights"] == 1).all()

    # Make sure right N!
    assert len(result_of_fit) == len(better)

    return better

EDIT: Oops. @nehargupta points out the result that the index in result_of_fit IS preserved, so edited with much much simpler solution.

[nickeubank]

@nehargupta you were right about index. Fixed solution.

[nickeubank]

Realizing the "weights" solution there is wrong -- repeats=False does not imply 1:1 matching. Gotta remember how weights are calculated when not doing 1:1 matching... It's the inverst for 1:1 matching with replacement, but here we can get odd sized groups of matches (e.g. 3 observations all matching).

[nickeubank]

This might implement weights correctly?

def get_dataframe(model, result_of_fit):

    # Get original data
    better = model.input_data.loc[result_of_fit.index]

    # Get match groups for clustering
    better["match_group"] = np.nan
    better["match_group_size"] = np.nan
    for idx, group in enumerate(model.units_per_group):
        better.loc[group, "match_group"] = idx
        better.loc[group, "match_group_size"] = len(group)

    # Get weights. I THINK this is right?! At least for with repeat=False?
    t = model.treatment_column_name
    better["t_in_group"] = better.groupby("match_group")[t].transform(np.sum)
    better["share_t_in_group"] = better["t_in_group"] / better["match_group_size"]
    better.loc[better[t] == 0, "prelim_weight"] = 1
    better.loc[better[t] == 1, "prelim_weight"] = 1 / better["share_t_in_group"]
    better["group_prelim_weight"] = better.groupby("match_group")[
        ["prelim_weight"]
    ].transform(sum)

    better["weights"] = better["prelim_weight"] * (
        better["match_group_size"] / better["group_prelim_weight"]
    )

    better = better.drop(
        ["prelim_weight", "group_prelim_weight", "t_in_group", "share_t_in_group"],
        axis="columns",
    )

    # Make sure right N!
    assert len(result_of_fit) == len(better)

    return better

But I'll admit my grasp of the stats of matching is limited -- I've never used it in research, so I've forgotten much of what I learned in class long ago. :)

[nehargupta]

Hi @nickeubank, I’m not sure I am familiar with the goal...is it the definition of weights here? https://www.rdocumentation.org/packages/MatchIt/versions/3.0.2/topics/matchit I initially thought it would be just the groups_per_unit parameter as it is, but that parameter is just the number of times a unit appears in a group (which is 1 for each unit if repeats=False, but that doesn't mean it's only matched to one unit, since groups can be any size 2->n where n is the number of units).

Sorry for the vocabulary confusion, I'm not that familiar with the other packages, but we welcome any suggestions that would help us to have more typical naming convention.

Glad you got the reindex issue figured out (we force a reindex if non-unique indexes because that would break it, we might want to add a user-warning in that case I guess, or possibly also non-integer if I remember correctly). Let me know if there's anything else!

[nickeubank]

Hahaha -- don't worry -- this is not a domain in which I'm expert, so any errors are probably the same as mine.

When the weights are right, then running a weighted least squares should give us our treatment effect estimate. e.g.:

smf.wls(
    "annual_earnings ~ has_college", matched_data, weights=matched_data["weights"]
).fit().summary()

But now you get statistical significance, can add coefficients for non-exact match corrections, cluster standard errors, etc.

OK, the more current docs for matchit have a longer definition (https://www.rdocumentation.org/packages/MatchIt/versions/4.1.0/topics/matchit), but also one I think is incomplete...

How Matching Weights Are Computed

Matching weights are computed in one of two ways depending on whether matching was done with replacement or not.

For matching without replacement, each unit is assigned to a subclass, which represents the pair they are a part of (in the case of k:1 matching) or the stratum they belong to (in the case of exact matching, coarsened exact matching, full matching, or subclassification). The formula for computing the weights depends on the argument supplied to estimand. A new stratum "propensity score" (p) is computed as the proportion of units in each stratum that are in the treated group, and all units in that stratum are assigned that propensity score. Weights are then computed using the standard formulas for inverse probability weights: for the ATT, weights are 1 for the treated units and p/(1-p) for the control units; for the ATC, weights are (1-p)/p for the treated units and 1 for the control units; for the ATE, weights are 1/p for the treated units and 1/(1-p) for the control units.

For matching with replacement, units are not assigned to unique strata. For the ATT, each treated unit gets a weight of 1. Each control unit is weighted as the sum of the inverse of the number of control units matched to the same treated unit across its matches. For example, if a control unit was matched to a treated unit that had two other control units matched to it, and that same control was matched to a treated unit that had one other control unit matched to it, the control unit in question would get a weight of 1/3 + 1/2 = 5/6. For the ATC, the same is true with the treated and control labels switched. The weights are computed using the match.matrix component of the matchit output object.

In each treatment group, weights are divided by the mean of the nonzero weights in that treatment group to make the weights sum to the number of units in that treatment group. If sampling weights are included through the s.weights argument, they will be included in the matchit output object but not incorporated into the matching weights. match.data, which extracts the matched set from a matchit object, combines the matching weights and sampling weights.

[nickeubank]

If weights are just 1 for treated units, and then 1/(num controls matched to a treatment) for treatments in a match group, then this would do it I think?

def get_dataframe(model, result_of_fit):

    # Get original data
    better = model.input_data.loc[result_of_fit.index]
    if not better.index.is_unique:
        raise ValueError("Need index values in input data to be unique")

    # Get match groups for clustering
    better["match_group"] = np.nan
    better["match_group_size"] = np.nan
    for idx, group in enumerate(model.units_per_group):
        better.loc[group, "match_group"] = idx
        better.loc[group, "match_group_size"] = len(group)

    # Get weights. I THINK this is right?! At least for with repeat=False?
    t = model.treatment_column_name
    better["t_in_group"] = better.groupby("match_group")[t].transform(np.sum)

    # Make weights
    better["weights"] = np.nan
    better.loc[better[t] == 1, "weights"] = 1  # treaments are 1

    # Controls start as proportional to num of treatments
    # each observation is matched to.
    better.loc[better[t] == 0, "weights"] = better["t_in_group"] / (
        better["match_group_size"] - better["t_in_group"]
    )

    # Then re-normalize for num unique control observations.
    control_weights = better[better[t] == 0]["weights"].sum()

    num_control_obs = len(better[better[t] == 0].index.drop_duplicates())
    renormalization = num_control_obs / control_weights
    better.loc[better[t] == 0, "weights"] = (
        better.loc[better[t] == 0, "weights"] * renormalization
    )
    assert better.weights.notnull().all()

    better = better.drop(["t_in_group"], axis="columns")

    # Make sure right length and values!
    assert len(result_of_fit) == len(better)
    assert better.loc[better[t] == 0, "weights"].sum() == num_control_obs

    return better

[nehargupta]

Oh! I didn't realize you wanted the weights to put the dataset in a WLS regression, now I see! Ok, thank you. I will need to think about this more before proceeding...

[nickeubank]

(nudge)

[nickeubank]

Dropping another nudge here.