\n", "Implements several types of propensity-score matching, balance diagnostics for group characteristics, average treatment effect on the treated (ATT) estimates, and bootstraping to estimate standard errors of the estimated ATT.\n", "\n", "**NB: one should proceed with caution using bootstrap standard errors... there's no guarantee it will work. See [\"On the Failure of the Bootstrap For Matching Estimators\" (Abadie and Imbens, 2008)](http://scholar.harvard.edu/imbens/files/on_the_failure_of_the_bootstrap_for_matching_estimators.pdf).**\n", "
\n",
"Propensity scoring methods:
\n",
"One-to-one: Matches individuals in the treatment group to a single individual in the control group. Variations include whether or not controls are matched with replacement, whether or not a caliper should be used, and the scale of the caliper.\n",
"
\n",
"One-to-many: Matches individuals in the treatment group to as many individuals in the control group as possible, subject to some criteria. Variations include whether or not controls are matched with replacement, caliper methods as in one-to-one matching, and k nearest neighbor matching."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import math\n",
"import numpy as np\n",
"import scipy\n",
"from scipy.stats import binom, hypergeom, gaussian_kde\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.linear_model import LogisticRegression\n",
"from ModelMatch import binByQuantiles\n",
"import statsmodels.api as sm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Propensity score estimation\n",
"\n",
"We use a logistic regression of treatment on covariates to estimate the propensity score. This is not the only way to estimate the propensity score, but it is the most common."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def computePropensityScore(formula, data, verbosity=1):\n",
" '''\n",
" Compute propensity scores \n",
" \n",
" Inputs:\n",
" formula = string of the form 'Treatment~ covariate1 + covariate2 + ...', where these are column names in data\n",
" data = matrix-like object with columns corresponding to terms in the formula\n",
" verbosity = whether or not to print glm summary\n",
" \n",
" dependencies: LogisticRegression from sklearn.linear_model\n",
" statsmodels as sm\n",
" '''\n",
" \n",
" ####### Using LogisticRegression from sklearn.linear_model \n",
" #propensity = LogisticRegression()\n",
" #propensity.fit(predictors, groups)\n",
" #return propensity.predict_proba(predictors)[:,1]\n",
" \n",
" ####### Using sm.GLM\n",
" #predictors = sm.add_constant(predictors, prepend=False)\n",
" #glm_binom = sm.GLM(groups, predictors, family=sm.families.Binomial())\n",
"\n",
" ####### Using sm.formula.glm with formula call\n",
" glm_binom = sm.formula.glm(formula = formula, data = data, family = sm.families.Binomial())\n",
" res = glm_binom.fit()\n",
" if verbosity:\n",
" print res.summary()\n",
" return res.fittedvalues"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Matching\n",
"\n",
"\n",
"We implement several variants of matching: one-to-one matching, one-to-many matching, with or without a caliper, and without or without replacement. Variants of the methods are examined in the following paper. \n",
"\n",
"Austin, P. C. (2014), A comparison of 12 algorithms for matching on the propensity score. Statist. Med., 33: 1057–1069. doi: 10.1002/sim.6004"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def Match(groups, propensity, caliper = 0.05, caliper_method = \"propensity\", replace = False):\n",
" ''' \n",
" Implements greedy one-to-one matching on propensity scores.\n",
" \n",
" Inputs:\n",
" groups = Array-like object of treatment assignments. Must be 2 groups\n",
" propensity = Array-like object containing propensity scores for each observation. Propensity and groups should be in the same order (matching indices)\n",
" caliper = a numeric value, specifies maximum distance (difference in propensity scores or SD of logit propensity) \n",
" caliper_method = a string: \"propensity\" (default) if caliper is a maximum difference in propensity scores,\n",
" \"logit\" if caliper is a maximum SD of logit propensity, or \"none\" for no caliper\n",
" replace = Logical for whether individuals from the larger group should be allowed to match multiple individuals in the smaller group.\n",
" (default is False)\n",
" \n",
" Output:\n",
" A series containing the individuals in the control group matched to the treatment group.\n",
" Note that with caliper matching, not every treated individual may have a match.\n",
" '''\n",
"\n",
" # Check inputs\n",
" if any(propensity <=0) or any(propensity >=1):\n",
" raise ValueError('Propensity scores must be between 0 and 1')\n",
" elif not(0<=caliper<1):\n",
" if caliper_method == \"propensity\" and caliper>1:\n",
" raise ValueError('Caliper for \"propensity\" method must be between 0 and 1')\n",
" elif caliper<0:\n",
" raise ValueError('Caliper cannot be negative')\n",
" elif len(groups)!= len(propensity):\n",
" raise ValueError('groups and propensity scores must be same dimension')\n",
" elif len(groups.unique()) != 2:\n",
" raise ValueError('wrong number of groups: expected 2')\n",
" \n",
" \n",
" # Transform the propensity scores and caliper when caliper_method is \"logit\" or \"none\"\n",
" if caliper_method == \"logit\":\n",
" propensity = log(propensity/(1-propensity))\n",
" caliper = caliper*np.std(propensity)\n",
" elif caliper_method == \"none\":\n",
" caliper = 0\n",
" \n",
" # Code groups as 0 and 1\n",
" groups = groups == groups.unique()[0]\n",
" N = len(groups)\n",
" N1 = groups[groups == 1].index; N2 = groups[groups == 0].index\n",
" g1, g2 = propensity[groups == 1], propensity[groups == 0]\n",
" # Check if treatment groups got flipped - the smaller should correspond to N1/g1\n",
" if len(N1) > len(N2):\n",
" N1, N2, g1, g2 = N2, N1, g2, g1\n",
" \n",
" \n",
" # Randomly permute the smaller group to get order for matching\n",
" morder = np.random.permutation(N1)\n",
" matches = {}\n",
"\n",
" \n",
" for m in morder:\n",
" dist = abs(g1[m] - g2)\n",
" if (dist.min() <= caliper) or not caliper:\n",
" matches[m] = dist.argmin() # Potential problem: check for ties\n",
" if not replace:\n",
" g2 = g2.drop(matches[m])\n",
" return (matches)\n",
"\n",
"\n",
"\n",
"def MatchMany(groups, propensity, method = \"caliper\", k = 1, caliper = 0.05, caliper_method = \"propensity\", replace = True):\n",
" ''' \n",
" Implements greedy one-to-many matching on propensity scores.\n",
" \n",
" Inputs:\n",
" groups = Array-like object of treatment assignments. Must be 2 groups\n",
" propensity = Array-like object containing propensity scores for each observation. Propensity and groups should be in the same order (matching indices)\n",
" method = a string: \"caliper\" (default) to select all matches within a given range, \"knn\" for k nearest neighbors,\n",
" k = an integer (default is 1). If method is \"knn\", this specifies the k in k nearest neighbors\n",
" caliper = a numeric value, specifies maximum distance (difference in propensity scores or SD of logit propensity) \n",
" caliper_method = a string: \"propensity\" (default) if caliper is a maximum difference in propensity scores,\n",
" \"logit\" if caliper is a maximum SD of logit propensity, or \"none\" for no caliper\n",
" replace = Logical for whether individuals from the larger group should be allowed to match multiple individuals in the smaller group.\n",
" (default is True)\n",
" \n",
" Output:\n",
" A series containing the individuals in the control group matched to the treatment group.\n",
" Note that with caliper matching, not every treated individual may have a match within calipers.\n",
" In that case we match it to its single nearest neighbor. The alternative is to throw out individuals with no matches, but then we'd no longer be estimating the ATT.\n",
" '''\n",
"\n",
" # Check inputs\n",
" if any(propensity <=0) or any(propensity >=1):\n",
" raise ValueError('Propensity scores must be between 0 and 1')\n",
" elif not(0<=caliper<1):\n",
" if caliper_method == \"propensity\" and caliper>1:\n",
" raise ValueError('Caliper for \"propensity\" method must be between 0 and 1')\n",
" elif caliper<0:\n",
" raise ValueError('Caliper cannot be negative')\n",
" elif len(groups)!= len(propensity):\n",
" raise ValueError('groups and propensity scores must be same dimension')\n",
" elif len(groups.unique()) != 2:\n",
" raise ValueError('wrong number of groups: expected 2')\n",
" \n",
" \n",
" # Transform the propensity scores and caliper when caliper_method is \"logit\" or \"none\"\n",
" if method == \"caliper\":\n",
" if caliper_method == \"logit\":\n",
" propensity = log(propensity/(1-propensity))\n",
" caliper = caliper*np.std(propensity)\n",
" elif caliper_method == \"none\":\n",
" caliper = 0\n",
" \n",
" # Code groups as 0 and 1\n",
" groups = groups == groups.unique()[0]\n",
" N = len(groups)\n",
" N1 = groups[groups == 1].index; N2 = groups[groups == 0].index\n",
" g1, g2 = propensity[groups == 1], propensity[groups == 0]\n",
" # Check if treatment groups got flipped - the smaller should correspond to N1/g1\n",
" if len(N1) > len(N2):\n",
" N1, N2, g1, g2 = N2, N1, g2, g1\n",
" \n",
" \n",
" # Randomly permute the smaller group to get order for matching\n",
" morder = np.random.permutation(N1)\n",
" matches = {}\n",
" \n",
" for m in morder:\n",
" dist = abs(g1[m] - g2)\n",
" dist.sort()\n",
" if method == \"knn\":\n",
" caliper = dist.iloc[k-1]\n",
" # PROBLEM: when there are ties in the knn. \n",
" # Need to randomly select among the observations tied for the farthest eacceptable distance\n",
" keep = np.array(dist[dist<=caliper].index)\n",
" if len(keep):\n",
" matches[m] = keep\n",
" else:\n",
" matches[m] = [dist.argmin()]\n",
" if not replace:\n",
" g2 = g2.drop(matches[m])\n",
" return (matches)\n",
"\n",
"\n",
" \n",
"def whichMatched(matches, data, many = False, unique = False):\n",
" ''' \n",
" Simple function to convert output of Matches to DataFrame of all matched observations\n",
" Inputs:\n",
" matches = output of Match\n",
" data = DataFrame of covariates\n",
" many = Boolean indicating if matching method is one-to-one or one-to-many\n",
" unique = Boolean indicating if duplicated individuals (ie controls matched to more than one case) should be removed\n",
" '''\n",
"\n",
" tr = matches.keys()\n",
" if many:\n",
" ctrl = [m for matchset in matches.values() for m in matchset]\n",
" else:\n",
" ctrl = matches.values()\n",
" # need to remove duplicate rows, which may occur in matching with replacement\n",
" temp = pd.concat([data.ix[tr], data.ix[ctrl]])\n",
" if unique == True:\n",
" return temp.groupby(temp.index).first()\n",
" else:\n",
" return temp"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def getWeights(matches, groups):\n",
" ''' computes weights for mean & regression according to how many times a control was matched in one-many matching'''\n",
" \n",
" ctrl = [m for matchset in matches.values() for m in matchset]\n",
" weights = groups.copy()\n",
" for c in ctrl:\n",
" weights[c] += 1\n",
" return weights\n",
" \n",
" \n",
"def whichMatched(matches, data, many = False, unique = False):\n",
" ''' \n",
" Simple function to convert output of Matches to DataFrame of all matched observations\n",
" Inputs:\n",
" matches = output of Match\n",
" data = DataFrame of covariates\n",
" many = Boolean indicating if matching method is one-to-one or one-to-many\n",
" unique = Boolean indicating if duplicated individuals (ie controls matched to more than one case) should be removed\n",
" '''\n",
"\n",
" tr = matches.keys()\n",
" if many:\n",
" ctrl = [m for matchset in matches.values() for m in matchset]\n",
" else:\n",
" ctrl = matches.values()\n",
" # need to remove duplicate rows, which may occur in matching with replacement\n",
" temp = pd.concat([data.ix[tr], data.ix[ctrl]])\n",
" if unique == True:\n",
" return temp.groupby(temp.index).first()\n",
" else:\n",
" return temp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Estimate average treatment effect on the treated\n",
"\n",
"Using the matches we found above, we can estimate ATT. We use two estimation strategies below: taking the difference in mean outcomes between treated and control groups and OLS regression."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def averageTreatmentEffect(groups, response, matches):\n",
" '''\n",
" Computes ATT using difference in means.\n",
" The data passed in should already have unmatched individuals and duplicates removed.\n",
"\n",
" Inputs:\n",
" groups = Series containing treatment assignment. Must be 2 groups\n",
" response = Series containing response measurements. Indices should match those of groups.\n",
" matches = output of Match or MatchMany\n",
" '''\n",
" if len(groups.unique()) != 2:\n",
" raise ValueError('wrong number of groups: expected 2')\n",
" \n",
" groups = (groups == groups.unique()[0])\n",
" response = response.groupby(response.index).first()\n",
" response1 = []; response0 = []\n",
" for k in matches.keys():\n",
" response1.append(response[k])\n",
" response0.append( (response[matches[k]]).mean() ) # Take mean response of controls matched to treated individual k\n",
" return np.array(response1).mean() - np.array(response0).mean()\n",
"\n",
"\n",
"def regressAverageTreatmentEffect(groups, response, covariates, matches=None, verbosity = 0):\n",
" '''\n",
" Computes ATT by regression.\n",
" This works for one-to-one matching. The data passed in should already have unmatched individuals removed.\n",
" Weights argument will be added later for one-many matching\n",
" \n",
" Inputs:\n",
" groups = Series containing treatment assignment. Must be 2 groups\n",
" response = Series containing response measurements. Indices should match those of groups.\n",
" covariates = DataFrame containing the covariates to include in the linear regression\n",
" matches = optional: if using one-many matching, should be the output of MatchMany.\n",
" Use None for one-one matching.\n",
" \n",
" Dependencies: statsmodels.api as sm, pandas as pd\n",
" '''\n",
" if len(groups.unique()) != 2:\n",
" raise ValueError('wrong number of groups: expected 2')\n",
" \n",
" weights = pd.Series(data = np.ones(len(groups)), index = groups.index)\n",
" if matches:\n",
" ctrl = [m for matchset in matches.values() for m in matchset] \n",
" matchcounts = pd.Series(ctrl).value_counts()\n",
" for i in matchcounts.index:\n",
" weights[i] = matchcounts[i]\n",
" if verbosity:\n",
" print weights.value_counts(), weights.shape\n",
" X = pd.concat([groups, covariates], axis=1)\n",
" X = sm.add_constant(X, prepend=False)\n",
" linmodel = sm.WLS(response, X, weights = weights).fit()\n",
" return linmodel.params[0], linmodel.bse[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def sampleWithinGroups(groups, data):\n",
" '''\n",
" To use in bootstrapping functions. \n",
" Sample with replacement from each group, return bootstrap sample dataframe\n",
" \n",
" Inputs:\n",
" groups = Series containing treatment assignment. Must be 2 groups\n",
" data = Dataframe containing observations from which to create bootstrap sample.\n",
" '''\n",
" bootdata = pd.DataFrame()\n",
" for g in groups.unique():\n",
" sample = np.random.choice(data.index[data.groups==g], sum(groups == g), replace = True)\n",
" newdata =(data[data.groups==g]).ix[sample]\n",
" bootdata = bootdata.append(newdata)\n",
" bootdata.index = range(len(groups))\n",
" return bootdata\n",
"\n",
"\n",
"def bootstrapATT(groups, response, propensity, many=True, B = 500, method = \"caliper\", k = 1, caliper = 0.05, caliper_method = \"propensity\", replace = False):\n",
" '''\n",
" Computes bootstrap standard error of the average treatment effect on the treated.\n",
" Sample observations with replacement, within each treatment group. Then match them and compute ATT.\n",
" Repeat B times and take standard deviation.\n",
" \n",
" Inputs:\n",
" groups = Series containing treatment assignment. Must be 2 groups\n",
" response = Series containing response measurements\n",
" propensity = Series containing propensity scores\n",
" many = Boolean: are we using one-many matching?\n",
" B = number of bootstrap replicates. Default is 500\n",
" caliper, caliper_method, replace = arguments to pass to Match or MatchMany\n",
" method, k = arguments to pass to MatchMany\n",
" '''\n",
" if len(groups.unique()) != 2:\n",
" raise ValueError('wrong number of groups: expected 2')\n",
"\n",
" data = pd.DataFrame({'groups':groups, 'response':response, 'propensity':propensity})\n",
" boot_ate = np.empty(B)\n",
" for i in range(B):\n",
" bootdata = sampleWithinGroups(groups, data)\n",
" if many:\n",
" pairs = MatchMany(bootdata.groups, bootdata.propensity, method = method, k = k, caliper = caliper, caliper_method = caliper_method, replace = replace)\n",
" else:\n",
" pairs = Match(bootdata.groups, bootdata.propensity, caliper = caliper, caliper_method = caliper_method, replace = replace)\n",
" boot_ate[i] = averageTreatmentEffect(bootdata.groups, bootdata.response, matches = pairs)\n",
" return boot_ate.std()\n",
"\n",
"\n",
"def bootstrapRegression(groups, response, propensity, covariates, many = True, B = 500, method = \"caliper\", k = 1, caliper = 0.05, caliper_method = \"propensity\", replace = False):\n",
" '''\n",
" Computes bootstrap standard error of the ATT.\n",
" Sample observations with replacement, within each treatment group. Then match them and compute ATT.\n",
" Repeat B times and take standard deviation.\n",
" \n",
" Inputs:\n",
" groups = Series containing treatment assignment. Must be 2 groups\n",
" response = Series containing response measurements\n",
" propensity = Series containing propensity scores\n",
" covariates = DataFrame containing covariates to use in regression\n",
" many = Boolean: are we using one-many matching?\n",
" B = number of bootstrap replicates. Default is 500\n",
" caliper, caliper_method, replace = arguments to pass to Match or MatchMany\n",
" method, k = arguments to pass to MatchMany\n",
" '''\n",
" if len(groups.unique()) != 2:\n",
" raise ValueError('wrong number of groups: expected 2')\n",
"\n",
" data = pd.DataFrame({'groups':groups, 'response':response, 'propensity':propensity})\n",
" data = pd.concat([data, covariates], axis=1)\n",
" boot_ate = np.empty(B)\n",
" for i in range(B):\n",
" bootdata = sampleWithinGroups(groups, data)\n",
" if many:\n",
" pairs = MatchMany(bootdata.groups, bootdata.propensity, method = method, k = k, caliper = caliper, caliper_method = caliper_method, replace = replace)\n",
" matched = whichMatched(pairs, bootdata, many = True, unique = True)\n",
" boot_ate[i] = regressAverageTreatmentEffect(matched.groups, matched.response, matched.ix[:,3:], matches=pairs, verbosity = 0)[0]\n",
" else:\n",
" pairs = Match(bootdata.groups, bootdata.propensity, caliper = caliper, caliper_method = caliper_method, replace = replace)\n",
" matched = whichMatched(pairs, bootdata, many = False)\n",
" boot_ate[i] = regressAverageTreatmentEffect(matched.groups, matched.response, matched.ix[:,3:], matches=None, verbosity = 0)[0]\n",
" return boot_ate.std()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Balance diagnostics\n",
"\n",
"After matching, we hope that the treatment group and matched controls are *close* in distribution with respect to each of their covariates."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def Balance(groups, covariates):\n",
" '''\n",
" Computes absolute difference of means and standard error for covariates by group\n",
" '''\n",
" means = covariates.groupby(groups).mean()\n",
" dist = abs(means.diff()).ix[1]\n",
" std = covariates.groupby(groups).std()\n",
" n = groups.value_counts()\n",
" se = std.apply(lambda(s): np.sqrt(s[0]**2/n[0] + s[1]**2/n[1]))\n",
" return dist, se\n",
"\n",
"def plotScores(groups, propensity, matches, many=True):\n",
" '''\n",
" Plot density of propensity scores for each group before and after matching\n",
" \n",
" Inputs: groups = treatment assignment, pre-matching\n",
" propensity = propensity scores, pre-matching\n",
" matches = output of Match or MatchMany\n",
" many = indicator - True if one-many matching was done (default is True), otherwise False\n",
" '''\n",
" pre = pd.DataFrame({'groups':groups, 'propensity':propensity}) \n",
" post = whichMatched(matches, pre, many = many, unique = False)\n",
" \n",
" plt.figure(1)\n",
" plt.subplot(121)\n",
" density0 = scipy.stats.gaussian_kde(pre.propensity[pre.groups==0])\n",
" density1 = scipy.stats.gaussian_kde(pre.propensity[pre.groups==1])\n",
" xs = np.linspace(0,1,1000)\n",
" #density0.covariance_factor = lambda : 0.5\n",
" #density0._compute_covariance()\n",
" #density1.covariance_factor = lambda : 0.5\n",
" #density1._compute_covariance()\n",
" plt.plot(xs,density0(xs),color='black')\n",
" plt.fill_between(xs,density1(xs),color='gray')\n",
" plt.title('Before Matching')\n",
" plt.xlabel('Propensity Score')\n",
" plt.ylabel('Density')\n",
" \n",
" plt.subplot(122)\n",
" density0_post = scipy.stats.gaussian_kde(post.propensity[post.groups==0])\n",
" density1_post = scipy.stats.gaussian_kde(post.propensity[post.groups==1])\n",
" xs = np.linspace(0,1,1000)\n",
" #density0.covariance_factor = lambda : 0.5\n",
" #density0._compute_covariance()\n",
" #density1.covariance_factor = lambda : 0.5\n",
" #density1._compute_covariance()\n",
" plt.plot(xs,density0_post(xs),color='black')\n",
" plt.fill_between(xs,density1_post(xs),color='gray')\n",
" plt.title('After Matching')\n",
" plt.xlabel('Propensity Score')\n",
" plt.ylabel('Density')\n",
" plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,
"nbformat_minor": 0
}