Auxiliary-to-Study Tilting for Estimation of the Average Treatment Effect on the Treated (ATT)

Bryan Graham - University of California - Berkeley

May 2016 (Python 3.6 updated July 2018)

Code citation:

Graham, Bryan S. (2016). "Auxiliary-to-Study Tilting for Estimation of the Average Treatment Effect on the Treated (ATT) Python Jupyter Notebook," (Version 1.0) [Computer program]. Available at http://bryangraham.github.io/econometrics/ (Accessed 04 May 2016)

Paper citation:

Graham, Bryan S., Cristine Pinto and Daniel Egel. (2016). “Efficient estimation of data combination models by the method of auxiliary-to-study tilting (AST),” Journal of Business and Economic Statistics 31 (2): 288 - 301. Available at http://www.tandfonline.com/doi/abs/10.1080/07350015.2015.1038544

This notebook illustrates how to use the Python 3.6 package ipt to estimate the Average Treatment Effect on the Treated (ATT) estimand under exogeneity (e.g., Imbens (2004, Review of Economics and Statistics). This package is registered on PyPi, with a GitHub repository at https://github.com/bryangraham/ipt.

The ipt package has the following dependencies: numpy, pandas, numpy.linalg, scipy, scipy.optimize and scipy.stats. These are standard libraries and are included in most scientific Python distributions. For example they are included in the highly recommended Anaconda distribution of Python. If you are using the Anaconda distribution of Python, then you can follow the (straightforward but tedious) instructions here to learn how install the ipt package from PyPi and make it available in Anaconda using the "conda" package manager. For users who anticipate only infrequent use, permanent installation of the ipt package may not be worth the trouble. One possibility is to just clone (ie., copy) the GitHub repository https://github.com/bryangraham/ipt, which contains the latest version of ipt. Then append the path pointing to the location of the ipt package on your local machine to your sys directory. This is what is done in the next snippet of code.


In [1]:
# Append location of ipt module base directory to system path
# NOTE: only required if permanent install not made (see comments above)
import sys
sys.path.append('/Users/bgraham/Dropbox/Sites/software/ipt/')

# Load ipt module
import ipt as ipt

Hide Runtime warning messages to clean up output.


In [2]:
import warnings

def fxn():
    warnings.warn("runtime", RuntimeWarning)

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fxn()

Now load additional libraries that will be used below.


In [3]:
# Direct Python to plot all figures inline (i.e., not in a separate window)
%matplotlib inline

# Load additional libraries
import numpy as np
import pandas as pd
import seaborn as sns 
import matplotlib.pyplot as plt
import statsmodels.api as sm

Introducing ipt.att()

The ipt package includes a Python 3.6 implementation of the ATT estimator introduced in Graham, Pinto and Egel (2016, Journal of Business and Economic Statistics). An implementation of the ATE estimator introduced in Graham, Pinto and Egel (2012, Review of Economic Studies) is planned for the near future.

The basic notation is as follows. $D \in \left \{ 0,1 \right \}$ is the binary treatment indicator. It equals one for treated units and zero for control units. The potential outcome under treatment is $Y$ and the potential outcome under control is $X$. This notation is admittedly non-standard but aligns with that used in the JBES paper. The paper considers a more general class of data combination problems of which ATT estimation is a special case.

The vector of always-observed pre-treatment control variables is $W$. Let $r(W)$ be a vector of known functions of $W$ such that the propensity score is a logit function with an index linear in $r(W)$. Let $\dim \left \{ r(W) \right \} = 1 + L$. The vector of known functions of $W$, $t(W)$, are those for which "exact balance" is imposed (let $\dim \left \{ t(W) \right \} = 1 + M$). Both $r(W)$ and $t(W)$ are assumed to contain a constant.

The treated and control units are re-weighted such that after re-weighting the mean of $t(W)$ in each of these subsamples equals its corresponding semiparametrically efficient estimate:

$$\bar{t}_{\textrm{eff}} = \frac{\sum_{i=1}^{N} G(r(W_{i})'\hat{\delta}) \times t(W_{i})} { \sum_{i=1}^{N} G(r(W_{i})'\hat{\delta})}$$

Here $\hat{\delta}$ is the maximum likelihood estimate of the propensity score parameter and $G(v)$ is the sigmoid/logit function.

Note that in the special case where all the elements in $t(W)$ are also contained in $r(W)$, the efficient estimate $\bar{t}_{\textrm{eff}}$ is numerically identical to the sub-sample mean of $t(W)$ across treated units alone (this is a peculiarity of the the structure of the logit score vector and would not be true if, for example, the propensity score took a probit form). In this special case, only the auxiliary/control units will need to be re-weighted to impose exact balance. However when $t(W)$ has at least some elements not in common with $r(W)$, then both the treatment and control subsamples will need to be re-weighted. In the typical empirical application the re-weighting of the treated units (relative to the empirical measure) tends to be quite modest, while that of the controls is more pronounced.

The various parameters needed by att() are listed in the function header. The function allows for survey sampling weights and clustered standard errors (although these two features have yet to be fully tested).


In [4]:
help(ipt.att)


Help on function att in module ipt.att:

att(D, Y, r_W, t_W, study_tilt=True, rlgrz=1, s_wgt=None, nocons=False, c_id=None, silent=False)
    AUTHOR: Bryan S. Graham, UC - Berkeley, bgraham@econ.berkeley.edu
    DATE: Python 2.7 code on 26 May 2016, updated for Python 3.6 on 15 July 2018        
    
    This function estimates the average treatment effect on the treated (ATT)
    using the "auxiliary-to-study tilting" (AST) method described by 
    Graham, Pinto and Egel (2016, Journal of Business and Economic Statistics). 
    The notation below mirrors that in the paper where possible. The Supplemental 
    Web Appendix of the paper describes the estimation algorithm implemented here 
    in detail. A copy of the paper and all supplemental appendices can be found 
    online at http://bryangraham.github.io/econometrics/
    
    INPUTS
    ------
    D         : N x 1 pandas.Series with ith element equal to 1 if ith unit in the merged
                sample is from the study population and zero if from the auxiliary
                population (i.e., D is the "treatment" indicator)
    Y         : N x 1  pandas.Series of observed outcomes                  
    r_W       : r(W), N x 1+L pandas.DataFrame of functions of always observed covariates
                -- these are the propensity score basis functions
    t_W       : t(W), N x 1+M pandas.DataFrame of functions of always observed covariates
                -- these are the balancing functions     
    study_tilt: If True compute the study sample tilt. This should be set to False 
                if all the elements in t(W) are also contained in h(W). In that 
                case the study_tilt coincides with its empirical measure.This 
                measure is returned in the pi_s vector when  study_tilt = False.
    rlgrz     : Regularization parameter. Should positive and less than or equal 
                to one. Smaller values correspond to less regularizations, but 
                may cause underflow problems when overlap is poor. The default 
                value will be adequate for most applications.
    s_wgt     : N x 1 pandas.Series of sampling weights variable (optional)
    nocons    : If True, then do NOT add constant to h_W and t_W matrix
                (only set to True if user passes in dataframes with constants included)
    c_id      : N X 1 pandas.Series of unique `cluster' id values (assumed to be integer valued) (optional)
                NOTE: Default is to assume independent observations and report heteroscedastic robust 
                      standard errors
                NOTE: Data are assumed to be pre-sorted by groups.
    silent    : if silent = True display less optimization information and use
                lower tolerance levels (optional)
    
    OUTPUTS
    -------
    gamma_ast         : AST estimate of gamma (the ATT)
    vcov_gamma_ast    : estimated large sample variance of gamma
    pscore_tests      : list of [study_test, auxiliary_test] where      
                        study_test     : ChiSq test statistic of H0 : lambda_s = 0; list with 
                                         [statistic, dof, p-val]
                                         NOTE: returns [None, None, None] if study_tilt = False
                        auxiliary_test : ChiSq test statistic of H0 : lambda_a = 0; list with 
                                         [statistic, dof, p-val]
    tilts             : numpy array with pi_eff, pi_s & pi_a as columns, sorted according
                        to the input data, and where                                     
                        pi_eff : Semiparametrically efficient estimate of F_s(W) 
                        pi_s   : Study sample tilt
                        pi_a   : Auxiliary sample tilt 
    exitflag          : 1 = success, 2 = can't compute MLE of p-score, 3 = can't compute study/treated tilt,
                        4 = can't compute auxiliary/control tilt
    
    FUNCTIONS CALLED  : logit()                             (...logit_logl(), logit_score(), logit_hess()...)
    ----------------    ast_crit(), ast_foc(), ast_soc()    (...ast_phi()...)

For problems where overlap is poor, the "rlgrz" regularization parameter can be an important user choice. Solving for the study/treatment and auxiliary/control "tilts" involves the solution to a convex problem over restricted domain. When $\bar{t}_{eff}$ is near the boundary of the convex hull of the scatter of $t(W)$ points across control units, solving for the auxiliary/control tilt can be numerically challenging and numerical underflow issues may arise. The default value of "rlgrz" (which is one) is generally sufficient for quick and reliable computation of the two tilts for most problems. For problems where overlap is poor sometimes the "rlgrz" parameter needs to be set to a smaller value to ensure that a valid solution is found. This corresponds to less regularizations, but also a harder optimization problem for the solver. Validity of the solution can be easily checked by ensuring that the balancing conditions are satisfied in practice (which are shown in the output). Note that in situations where overlap is very poor there will be no valid tilt of the control/auxiliary sample. In such cases the AST estimate is undefined. My general advice is to not proceed in such situations or to use parametric estimation procedures where researcher reliance on functional form assumptions is transparent. However, there are some interesting recent alternative proposals that may be worth carefully considering when overlap is poor. See for example the working paper by Athey, Imbens and Wagner (2016) on the Arxiv http://arxiv.org/abs/1604.07125.

Illustration with NSW evaluation dataset

Following Dehejia and Wahba (1999, Journal of American Statistical Association) I work with the subset of the NSW evaluation dataset where two pretreatment earnings levels are observed (earnings in both 1974 and 1975). Earnings in 1978, which is post-treatment, is the outcome of interest. I construct two data frames, the first is the subset of NSW participants with two years of pre-treatment earnings ($N = 445$). Of these units $185$ correspond to treated units, with the remaining $260$ being controls.

Following Dehejia and Wahba I also create an "observational" dataset consisting of the $185$ NSW treated units and a set of $2,490$ non-treated control units from the PSID. These two datasets are put into the nsw and nsw_obs pandas data frames. Some basic statistics on the "observational" dataset are reported below.

These datasets can be imported into a pandas dataframe directly from Rajeev Dehejia's webpage (http://users.nber.org/~rdehejia/index.html)


In [5]:
# Read in LaLonde NSW data as DataFrame directly from Rajeev Dehejia's webpage
nsw=pd.read_stata("http://www.nber.org/~rdehejia/data/nsw_dw.dta")

# Read in LaLonde PSID pseudo-controls as DataFrame directly from Rajeev Dehejia's webpage
psid = pd.read_stata("http://www.nber.org/~rdehejia/data/psid_controls.dta") 

# Create pseudo-observational dataset of NSW treated subjects and PSID "controls"
frames = [nsw[(nsw.treat != 0)], psid]
nsw_obs = pd.concat(frames)

# Make some adjustments to variable definitions in observational dataframe
nsw_obs['constant'] = 1                    # Add constant to observational dataframe
nsw_obs['age']      = nsw_obs['age']/10    # Rescale age to be in decades

nsw_obs['re74']        = nsw_obs['re74']/1000            # Recale earnings to be in thousands
nsw_obs['re75']        = nsw_obs['re75']/1000            # Recale earnings to be in thousands
nsw_obs['re74_X_re75'] = nsw_obs['re74']*nsw_obs['re75'] # Prior earnings interaction
nsw_obs['re74_sq']     = nsw_obs['re74']**2              # Squares of prior earnings
nsw_obs['re75_sq']     = nsw_obs['re75']**2              # Squares of prior earnings

nsw_obs['ue74']        = 1*(nsw_obs['re74']==0)          # Dummy for unemployed in 1974
nsw_obs['ue75']        = 1*(nsw_obs['re75']==0)          # Dummy for unemployed in 1975
nsw_obs['ue74_X_ue75'] = nsw_obs['ue74']*nsw_obs['ue75'] # Prior unemployment interaction


nsw_obs.describe()


Out[5]:
treat age education black hispanic married nodegree re74 re75 re78 constant re74_X_re75 re74_sq re75_sq ue74 ue75 ue74_X_ue75
count 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.0 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000 2675.000000
mean 0.069159 3.422579 11.994392 0.291589 0.034393 0.819439 0.333084 18.229973 17.850929 20502.371094 1.0 488.345612 520.561951 511.173492 0.129346 0.134579 0.101308
std 0.253759 1.049983 3.053559 0.454579 0.182268 0.384722 0.471402 13.722255 13.877773 15632.524414 0.0 825.004883 846.587952 890.807251 0.335645 0.341338 0.301793
min 0.000000 1.700000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
25% 0.000000 2.500000 10.000000 0.000000 0.000000 1.000000 0.000000 8.816700 7.605290 9243.400879 1.0 53.216442 77.734200 57.840645 0.000000 0.000000 0.000000
50% 0.000000 3.200000 12.000000 0.000000 0.000000 1.000000 0.000000 17.437475 17.008064 19432.103516 1.0 271.637817 304.065552 289.274261 0.000000 0.000000 0.000000
75% 0.000000 4.350000 14.000000 1.000000 0.000000 1.000000 1.000000 25.470469 25.583711 28815.667969 1.0 623.526764 648.744751 654.526581 0.000000 0.000000 0.000000
max 1.000000 5.500000 17.000000 1.000000 1.000000 1.000000 1.000000 137.148682 156.653229 121173.578125 1.0 21484.783203 18809.761719 24540.234375 1.000000 1.000000 1.000000

Experimental benchmark estimate of the ATT

A simple comparision of post-treatment earnings across NSW treated and control units provides an unbiased estimate of the ATT. Such an estimate is computed below. Treated individuals earned $1,794 more in 1978 than non-treated individuals.


In [6]:
Y = nsw['re78']
D = nsw['treat']

nsw_reg1=sm.OLS(Y,sm.add_constant(D)).fit(cov_type='HC0')
print('------------------------------------------------------------------------------')
print('- NSW Experimental Benchmark                                                 -')
print('------------------------------------------------------------------------------')
print('')
print(nsw_reg1.summary())


------------------------------------------------------------------------------
- NSW Experimental Benchmark                                                 -
------------------------------------------------------------------------------

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   re78   R-squared:                       0.018
Model:                            OLS   Adj. R-squared:                  0.016
Method:                 Least Squares   F-statistic:                     7.187
Date:                Tue, 30 Oct 2018   Prob (F-statistic):            0.00762
Time:                        11:17:07   Log-Likelihood:                -4542.7
No. Observations:                 445   AIC:                             9089.
Df Residuals:                     443   BIC:                             9098.
Df Model:                           1                                         
Covariance Type:                  HC0                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const       4554.8011    339.438     13.419      0.000    3889.514    5220.088
treat       1794.3424    669.315      2.681      0.007     482.508    3106.176
==============================================================================
Omnibus:                      282.071   Durbin-Watson:                   2.064
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             3430.892
Skew:                           2.547   Prob(JB):                         0.00
Kurtosis:                      15.613   Cond. No.                         2.47
==============================================================================

Warnings:
[1] Standard Errors are heteroscedasticity robust (HC0)

Observational estimate of the ATT using PSID controls

Next I attempt estimate the ATT using the PSID "observational" controls instead of the experimental ones. Even if we are willing to assume unconfoundedness holds in this setting, computing the ATT in this context is challenging. To show why this is the case, the next snippet of code plots the convex hull of joint empirical distribution of 1974 and 1975 earnings across the PSID controls. It also plots the mean earnings across these two years for both NSW treatment and PSID control units.

Note that mean earnings across treated units is very near the boundary of the convex hull of earnings values across control units. This means that any re-weighting of control units that equates their mean earnings with those observed across treated units will need to place a large amount of mass on a small number of units and close to zero mass on many units. The required tilt will be "far" from the empirical measure of the control sample. This tilt may be hard to compute.


In [7]:
# Directory where figures generated below are saved
workdir =  '/Users/bgraham/Dropbox/Sites/software/ipt/Notebooks/'

from scipy.spatial import ConvexHull

# Extract pre-treatment earnings for NSW treated units and PSID controls
earnings_treatment = np.asarray(nsw_obs[(nsw_obs.treat == 1)][['re74','re75']])
earnings_control   = np.asarray(nsw_obs[(nsw_obs.treat == 0)][['re74','re75']])

# Calculate convex hull of PSID control units' earnings realizations
hull = ConvexHull(earnings_control)

# Create a figure object to plot the convex hull
convexhull_fig = plt.figure(figsize=(6, 6))    

# Scatter plot of pre-treatment earnings in 1974 and 1975 for PSID controls
ax = convexhull_fig.add_subplot(1,1,1)
sns.regplot(x="re74", y="re75", data=nsw_obs[(nsw_obs.treat == 0)], \
            fit_reg=False, color='#FDB515')
plt.title('Convex Hull of Earnings for PSID Controls', fontsize=12)
plt.xlabel('Earnings in 1974')
plt.ylabel('Earnings in 1975')

# Set plot range and tick marks
plt.ylim([0,175])
plt.xlim([0,175])    
plt.xticks([0, 25, 50, 75, 100, 125, 150], fontsize=10)
plt.yticks([0, 25, 50, 75, 100, 125, 150], fontsize=10)

# Plot mean earnings for NSW treated units and PSID controls
plt.plot(np.mean(earnings_control[:,0]), np.mean(earnings_control[:,1]), \
         color='#00B0DA', marker='o', markersize=10)
plt.plot(np.mean(earnings_treatment[:,0]), np.mean(earnings_treatment[:,1]), \
         color='#EE1F60', marker='s', markersize=10)

# Plot convex hull
for simplex in hull.simplices:
    plt.plot(earnings_control[simplex, 0], earnings_control[simplex, 1], 'k-')

# Clean up the plot, add frames, remove gridlines etc.
ax = plt.gca()
ax.patch.set_facecolor('gray')               # Color of background
ax.patch.set_alpha(0.15)                     # Translucency of background
ax.grid(False)                               # Remove gridlines from plot

# Add frame around plot
for spine in ['left','right','top','bottom']:
    ax.spines[spine].set_visible(True)
    ax.spines[spine].set_color('k')
    ax.spines[spine].set_linewidth(2)

# Add legend to the plot
import matplotlib.lines as mlines

psid_patch      = mlines.Line2D([], [], color='#FDB515', marker='o', linestyle='None',\
                                markersize=5, label='PSID controls')
psid_mean_patch = mlines.Line2D([], [], color='#00B0DA', marker='o', linestyle='None',\
                                markersize=10, label='PSID mean')
nsw_mean_patch  = mlines.Line2D([], [], color='#EE1F60', marker='s', linestyle='None',\
                                markersize=10, label='NSW mean')

lgd = plt.legend(handles=[psid_patch, psid_mean_patch, nsw_mean_patch], \
                          loc='upper left', fontsize=12, ncol=2, numpoints = 1) 
        
# Render & save plot
plt.tight_layout()
plt.savefig(workdir+'Fig_LaLonde_Convex_Hull.png')


Computing the AST estimates of the ATT using the PSID controls

This next snippet of code computes the ATT estimate using the PSID controls. This shows how use the ipt.att() function. The vector $ t(W) $ is specified to include earnings in both 1974 and 1975, as well as their squares and cross product. This choice ensures exact balancing of the first two moments of the pre-treatment earnings distributions across treated and controls. I additionally include in $t(W)$ age, education and the black and hispanic dummies. I rescale some of these variables to make computation easier (which is good to do in practice). The constructed auxiliary tilt satisfies the exact balancing conditioning to (at least) four decimal places. I also set $r(W) = t(W)$


In [8]:
# Treatment indicator
D = nsw_obs['treat']

# Propensity score model
h_W = nsw_obs[['constant','black','hispanic','age','married','nodegree',\
               're74','re75','re74_X_re75','ue74','ue75','ue74_X_ue75']]

# Balancing moments
t_W = h_W

# Outcome
Y = nsw_obs['re78']

# Compute AST estimate of ATT
[gamma_ast, vcov_gamma_ast, pscore_tests, tilts, exitflag] = \
                                            ipt.att(D, Y, h_W, t_W, study_tilt = False, rlgrz = 1/2, \
                                                    s_wgt=None, nocons=True, c_id=None, silent=False)

# Compute log of auxiliary/control sample normalized weights to construct histogram below
pi_a = tilts[:,2]               # Auxiliary tilt is the last columnn in np.array tilts
b = np.where(1-D)[0]            # find indices of control units
Na = len(pi_a[b])
log_wgts = np.log(Na*pi_a[b])   # log of normalized weights to undo right skewness for plotting purposes


--------------------------------------------------------------
- Computing propensity score by MLE                          -
--------------------------------------------------------------
Fisher-Scoring Derivative check (2-norm): 3.07242496
Value of -logL = 1380.537270,  2-norm of score = 311857.625104
Value of -logL = 1011.727881,  2-norm of score = 136618.736764
Value of -logL = 781.262441,  2-norm of score = 58983.641036
Value of -logL = 641.725866,  2-norm of score = 25343.082588
Value of -logL = 559.194658,  2-norm of score = 10686.495400
Value of -logL = 513.715972,  2-norm of score = 4322.399768
Value of -logL = 491.765974,  2-norm of score = 1681.029336
Value of -logL = 453.692345,  2-norm of score = 1225.888561
Value of -logL = 400.103785,  2-norm of score = 597.192586
Value of -logL = 373.444453,  2-norm of score = 740.225515
Value of -logL = 372.932112,  2-norm of score = 172.461940
Value of -logL = 301.339647,  2-norm of score = 252.270987
Value of -logL = 278.364415,  2-norm of score = 17963.682831
Value of -logL = 276.591512,  2-norm of score = 423.681909
Value of -logL = 276.574884,  2-norm of score = 176.858042
Value of -logL = 265.972550,  2-norm of score = 211.811965
Value of -logL = 247.883268,  2-norm of score = 208.495592
Value of -logL = 247.879412,  2-norm of score = 65.265454
Value of -logL = 238.543679,  2-norm of score = 62.150336
Value of -logL = 224.826973,  2-norm of score = 131.938638
Value of -logL = 224.825497,  2-norm of score = 32.790992
Value of -logL = 221.198232,  2-norm of score = 35.091565
Value of -logL = 215.412924,  2-norm of score = 46.967811
Value of -logL = 215.362248,  2-norm of score = 21.094903
Value of -logL = 213.354464,  2-norm of score = 19.625374
Value of -logL = 209.895365,  2-norm of score = 19.581332
Value of -logL = 205.808250,  2-norm of score = 42.759714
Value of -logL = 205.808135,  2-norm of score = 8.018157
Value of -logL = 205.803938,  2-norm of score = 5.810941
Value of -logL = 205.799213,  2-norm of score = 2.785580
Value of -logL = 205.614151,  2-norm of score = 2.715314
Value of -logL = 205.309112,  2-norm of score = 3.572834
Value of -logL = 205.042375,  2-norm of score = 7.260134
Value of -logL = 205.042371,  2-norm of score = 1.209296
Value of -logL = 205.032485,  2-norm of score = 0.933703
Value of -logL = 205.014639,  2-norm of score = 0.779689
Value of -logL = 204.986267,  2-norm of score = 0.780225
Value of -logL = 204.981685,  2-norm of score = 0.393355
Value of -logL = 204.973319,  2-norm of score = 0.321765
Value of -logL = 204.960669,  2-norm of score = 0.176436
Value of -logL = 204.956000,  2-norm of score = 0.030700
Value of -logL = 204.955944,  2-norm of score = 0.004989
Value of -logL = 204.955940,  2-norm of score = 0.000455
Value of -logL = 204.955940,  2-norm of score = 0.000000
Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 204.955940
         Iterations: 44
         Function evaluations: 127
         Gradient evaluations: 159
         Hessian evaluations: 45

----------------------------------------------------------------------
- Tilt of study sample not requested by user (study_tilt = False).   -
- Validity of this requires all elements of t(W) to be elements of   -
- h(W) as well. User is advised to verify this condition.            -
----------------------------------------------------------------------


--------------------------------------------------------------
- Computing auxiliary/control sample tilt                    -
--------------------------------------------------------------
Auxiliary sample tilt derivative check (2-norm): 0.00054355
Value of ast_crit = 0.718308,  2-norm of ast_foc = 33.459844
Value of ast_crit = 0.717327,  2-norm of ast_foc = 12.279000
Value of ast_crit = 0.716970,  2-norm of ast_foc = 4.487388
Value of ast_crit = 0.716841,  2-norm of ast_foc = 1.626349
Value of ast_crit = 0.716796,  2-norm of ast_foc = 0.593654
Value of ast_crit = 0.710410,  2-norm of ast_foc = 0.639011
Value of ast_crit = 0.659198,  2-norm of ast_foc = 0.443667
Value of ast_crit = 0.589830,  2-norm of ast_foc = 1.004049
Value of ast_crit = 0.588063,  2-norm of ast_foc = 0.282596
Value of ast_crit = 0.585411,  2-norm of ast_foc = 0.163608
Value of ast_crit = 0.580366,  2-norm of ast_foc = 0.148367
Value of ast_crit = 0.571659,  2-norm of ast_foc = 0.123351
Value of ast_crit = 0.562478,  2-norm of ast_foc = 0.116339
Value of ast_crit = 0.562409,  2-norm of ast_foc = 0.030090
Value of ast_crit = 0.559840,  2-norm of ast_foc = 0.023407
Value of ast_crit = 0.556342,  2-norm of ast_foc = 0.016828
Value of ast_crit = 0.555690,  2-norm of ast_foc = 0.005597
Value of ast_crit = 0.555342,  2-norm of ast_foc = 0.001345
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000199
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000001
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.555308,  2-norm of ast_foc = 0.000000
Optimization terminated successfully.
         Current function value: 0.555308
         Iterations: 25
         Function evaluations: 32
         Gradient evaluations: 56
         Hessian evaluations: 25

-------------------------------------------------------------------------------------------
- Auxiliary-to-Study (AST) estimates of the ATT                                           -
-------------------------------------------------------------------------------------------
ATT:  2354.973392
     (746.472279)

-------------------------------------------------------------------------------------------
NOTE: Outcome variable   = re78
      Heteroscedastic-robust standard errors reported
      N1 = 185, N0 = 2490

-------------------------------------------------------------------------------------------
- Maximum likelihood estimates of the p-score                                             -
-------------------------------------------------------------------------------------------

Independent variable       Coef.    ( Std. Err.) 
-------------------------------------------------------------------------------------------
constant                   0.608334 (  0.615931)
black                      2.575545 (  0.334231)
hispanic                   2.155130 (  0.624783)
age                       -1.101892 (  0.152297)
married                   -1.541664 (  0.271927)
nodegree                   0.695969 (  0.268419)
re74                       0.007705 (  0.039750)
re75                      -0.311054 (  0.071125)
re74_X_re75                0.002299 (  0.000485)
ue74                       2.533851 (  0.674380)
ue75                      -2.634600 (  0.888030)
ue74_X_ue75                1.848389 (  1.008136)

-------------------------------------------------------------------------------------------
- Tilting parameter estimates                                                             -
-------------------------------------------------------------------------------------------

--------------------------------------------------------
- NOTE: Study tilt not computed (study_tilt = False).  -
-       Components of t(W) assumed to be also in h(W). -
--------------------------------------------------------


CONTROL (auxiliary) sample tilt
-------------------------------------------------------------------------------------------

Independent variable       Coef.    ( Std. Err.) 
-------------------------------------------------------------------------------------------
constant                   0.648563 (  0.615931)
black                     -0.872087 (  0.334231)
hispanic                  -0.675904 (  0.624783)
age                        0.085623 (  0.152297)
married                    1.176204 (  0.271927)
nodegree                   0.045673 (  0.268419)
re74                       0.008373 (  0.039750)
re75                       0.099005 (  0.071125)
re74_X_re75               -0.018313 (  0.000485)
ue74                      -1.587588 (  0.674380)
ue75                      -0.276814 (  0.888030)
ue74_X_ue75                1.583684 (  1.008136)

Specification test for p-score (H_0 : lambda_a = 0)
-------------------------------------------------------------------------------------------
chi-square(12) = 100968.821378   p-value:  0.000000

Summary statistics auxiliary/control re-weighting
-------------------------------------------------------------------------------------------
Kish's effective auxiliary/control sample size = 24

Percentiles of N_a * pi_a distribution
 1 percentile = 0.0000
 5 percentile = 0.0000
10 percentile = 0.0000
25 percentile = 0.0003
50 percentile = 0.0034
75 percentile = 0.0648
90 percentile = 0.6238
95 percentile = 2.3329
99 percentile = 13.0890

Means & standard deviations of t_W (pre-balance)                                           
-------------------------------------------------------------------------------------------
                            Treated (D = 1)        Control (D = 0)        Norm. Diff.      
-------------------------------------------------------------------------------------------
constant                   1.0000 (  0.0000)      1.0000 (  0.0000)         nan
black                      0.8432 (  0.3636)      0.2506 (  0.4334)      1.4816
hispanic                   0.0595 (  0.2365)      0.0325 (  0.1774)      0.1288
age                        2.5816 (  0.7136)      3.4851 (  1.0439)     -1.0104
married                    0.1892 (  0.3917)      0.8663 (  0.3404)     -1.8453
nodegree                   0.7081 (  0.4546)      0.3052 (  0.4605)      0.8805
re74                       2.0956 (  4.8734)     19.4287 ( 13.4042)     -1.7187
re75                       1.5321 (  3.2105)     19.0633 ( 13.5942)     -1.7750
re74_X_re75               13.1186 ( 50.6452)    523.6540 (844.2237)     -0.8537
ue74                       0.7081 (  0.4546)      0.0863 (  0.2809)      1.6454
ue75                       0.6000 (  0.4899)      0.1000 (  0.3000)      1.2309
ue74_X_ue75                0.5892 (  0.4920)      0.0651 (  0.2466)      1.3469

Means and standard deviations of t_W (post-balance)                                        
-------------------------------------------------------------------------------------------
                            Treated (D = 1)        Control (D = 0)        Efficient (D = 1)
-------------------------------------------------------------------------------------------
constant                   1.0000 (  0.0000)      1.0000 (  0.0000)      1.0000 (  0.0000)    
black                      0.8432 (  0.3636)      0.8432 (  0.3636)      0.8432 (  0.3636)    
hispanic                   0.0595 (  0.2365)      0.0595 (  0.2365)      0.0595 (  0.2365)    
age                        2.5816 (  0.7136)      2.5816 (  0.7605)      2.5816 (  0.7865)    
married                    0.1892 (  0.3917)      0.1892 (  0.3917)      0.1892 (  0.3917)    
nodegree                   0.7081 (  0.4546)      0.7081 (  0.4546)      0.7081 (  0.4546)    
re74                       2.0956 (  4.8734)      2.0956 (  4.7190)      2.0956 (  4.7439)    
re75                       1.5321 (  3.2105)      1.5321 (  3.2098)      1.5321 (  3.1509)    
re74_X_re75               13.1186 ( 50.6452)     13.1186 (199.4132)     13.1186 (199.5765)    
ue74                       0.7081 (  0.4546)      0.7081 (  0.4546)      0.7081 (  0.4546)    
ue75                       0.6000 (  0.4899)      0.6000 (  0.4899)      0.6000 (  0.4899)    
ue74_X_ue75                0.5892 (  0.4920)      0.5892 (  0.4920)      0.5892 (  0.4920)    
/Users/bgraham/Dropbox/Sites/software/ipt/ipt/att.py:658: RuntimeWarning: invalid value encountered in true_divide
  NormDif_t    = (mu_t_D1 - mu_t_D0)/np.sqrt((mu_t_D1_std**2 + mu_t_D0_std**2)/2)

Observe that the test of $H_{0} : \lambda_{a} = 0$ is decisively rejected. This should be unsurprising. The covariate distribution in the PSID is so different from that observed across the NSW treated units that it would be extraordinary if the p-score took a parsimonious form (recall from Graham, Pinto and Egel (2016, JBES), and work on discriminant analysis more generally, that the propensity score is related to the ratio of the covariate densities in the treated and control samples).

The ATT estimate is close to its experimental benchmark, with a somewhat larger standard error. Despite the PSID including 2,490 controls the standard error goes up relative to the experimental benchmark (which is based on a much smaller number of controls). This reflects poor overlap.

Experimental benchmark estimate of the ATT with covariate adjustment by AST

Next I compute the AST estimate of the ATT using the NSW experimental controls and the same specification for $t(W)$ given above. However, I set $r(W)$ to include a constant only (corresponding to the fact that we know a priori that the propensity score is constant due the RCT structure of the data). Note in this case that both the study and auxiliary tilts need to be computed.


In [9]:
# Make some adjustments to variable definitions in experimental dataframe
nsw['constant'] = 1                # Add constant to observational dataframe
nsw['age']      = nsw['age']/10    # Rescale age to be in decades

nsw['re74']        = nsw['re74']/1000          # Recale earnings to be in thousands
nsw['re75']        = nsw['re75']/1000          # Recale earnings to be in thousands
nsw['re74_X_re75'] = nsw['re74']*nsw['re75']   # Prior earnings interaction
nsw['re74_sq']     = nsw['re74']**2            # Squares of prior earnings
nsw['re75_sq']     = nsw['re75']**2            # Squares of prior earnings

nsw['ue74']        = 1*(nsw['re74']==0)        # Dummy for unemployed in 1974
nsw['ue75']        = 1*(nsw['re75']==0)        # Dummy for unemployed in 1975
nsw['ue74_X_ue75'] = nsw['ue74']*nsw['ue75']   # Prior unemployment interaction

# Treatment indicator
D = nsw['treat']

# Balancing moments
t_W = nsw[['constant','black','hispanic','age','married','nodegree',\
           're74','re75','re74_X_re75','ue74','ue75','ue74_X_ue75']]

# Propensity score variables
h_W = nsw[['constant']]

# Outcome
Y = nsw['re78']

# Compute AST estimate of ATT
[gamma_ast, vcov_gamma_ast, pscore_tests, tilts, exitflag] = \
                                            ipt.att(D, Y, h_W, t_W, study_tilt = True, rlgrz = 1/2, \
                                                    s_wgt=None, nocons=True, c_id=None, silent=False)

# Compute log of auxiliary/control sample normalized weights to construct histogram below
     
pi_a_nsw = tilts[:,2]                  # Auxiliary tilt is the last columnn in np.array tilts
b = np.where(1-D)[0]                   # find indices of control units
Na = len(pi_a_nsw[b])
log_wgts_nsw = np.log(Na*pi_a_nsw[b])  # log of normalized weights to undo right skewness for plotting purposes


--------------------------------------------------------------
- Computing propensity score by MLE                          -
--------------------------------------------------------------
Fisher-Scoring Derivative check (2-norm): 0.00000133
Value of -logL = 302.100574,  2-norm of score = 0.351080
Value of -logL = 302.100004,  2-norm of score = 0.000095
Value of -logL = 302.100004,  2-norm of score = 0.000000
Value of -logL = 302.100004,  2-norm of score = 0.000000
Optimization terminated successfully.
         Current function value: 302.100004
         Iterations: 4
         Function evaluations: 5
         Gradient evaluations: 8
         Hessian evaluations: 4

--------------------------------------------------------------
- Computing study/treated sample tilt                        -
--------------------------------------------------------------
Study sample tilt derivative check (2-norm): 0.00090181
Value of ast_crit = 0.725526,  2-norm of ast_foc = 0.284082
Value of ast_crit = 0.717273,  2-norm of ast_foc = 0.292519
Value of ast_crit = 0.706193,  2-norm of ast_foc = 0.623114
Value of ast_crit = 0.706108,  2-norm of ast_foc = 0.109790
Value of ast_crit = 0.693197,  2-norm of ast_foc = 0.090100
Value of ast_crit = 0.675873,  2-norm of ast_foc = 0.114986
Value of ast_crit = 0.675581,  2-norm of ast_foc = 0.033087
Value of ast_crit = 0.672829,  2-norm of ast_foc = 0.028932
Value of ast_crit = 0.668197,  2-norm of ast_foc = 0.023720
Value of ast_crit = 0.663593,  2-norm of ast_foc = 0.039369
Value of ast_crit = 0.663568,  2-norm of ast_foc = 0.003718
Value of ast_crit = 0.663522,  2-norm of ast_foc = 0.003292
Value of ast_crit = 0.663443,  2-norm of ast_foc = 0.002381
Value of ast_crit = 0.663357,  2-norm of ast_foc = 0.000592
Value of ast_crit = 0.663356,  2-norm of ast_foc = 0.000010
Value of ast_crit = 0.663356,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.663356,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.663356,  2-norm of ast_foc = 0.000000
Optimization terminated successfully.
         Current function value: 0.663356
         Iterations: 18
         Function evaluations: 25
         Gradient evaluations: 42
         Hessian evaluations: 18

--------------------------------------------------------------
- Computing auxiliary/control sample tilt                    -
--------------------------------------------------------------
Auxiliary sample tilt derivative check (2-norm): 0.00016040
Value of ast_crit = 0.216762,  2-norm of ast_foc = 0.165042
Value of ast_crit = 0.212013,  2-norm of ast_foc = 0.133727
Value of ast_crit = 0.203322,  2-norm of ast_foc = 0.113471
Value of ast_crit = 0.190292,  2-norm of ast_foc = 0.109981
Value of ast_crit = 0.185279,  2-norm of ast_foc = 0.085831
Value of ast_crit = 0.185276,  2-norm of ast_foc = 0.020239
Value of ast_crit = 0.183035,  2-norm of ast_foc = 0.016466
Value of ast_crit = 0.179977,  2-norm of ast_foc = 0.024188
Value of ast_crit = 0.179834,  2-norm of ast_foc = 0.003591
Value of ast_crit = 0.179574,  2-norm of ast_foc = 0.002963
Value of ast_crit = 0.179191,  2-norm of ast_foc = 0.001995
Value of ast_crit = 0.179072,  2-norm of ast_foc = 0.000783
Value of ast_crit = 0.179071,  2-norm of ast_foc = 0.000028
Value of ast_crit = 0.179071,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.179071,  2-norm of ast_foc = 0.000000
Value of ast_crit = 0.179071,  2-norm of ast_foc = 0.000000
Optimization terminated successfully.
         Current function value: 0.179071
         Iterations: 16
         Function evaluations: 23
         Gradient evaluations: 38
         Hessian evaluations: 16

-------------------------------------------------------------------------------------------
- Auxiliary-to-Study (AST) estimates of the ATT                                           -
-------------------------------------------------------------------------------------------
ATT:  1681.596850
     (689.815792)

-------------------------------------------------------------------------------------------
NOTE: Outcome variable   = re78
      Heteroscedastic-robust standard errors reported
      N1 = 185, N0 = 260

-------------------------------------------------------------------------------------------
- Maximum likelihood estimates of the p-score                                             -
-------------------------------------------------------------------------------------------

Independent variable       Coef.    ( Std. Err.) 
-------------------------------------------------------------------------------------------
constant                  -0.340326 (  0.099124)

-------------------------------------------------------------------------------------------
- Tilting parameter estimates                                                             -
-------------------------------------------------------------------------------------------

TREATED (study) sample tilt
-------------------------------------------------------------------------------------------

Independent variable       Coef.    ( Std. Err.) 
-------------------------------------------------------------------------------------------
constant                   0.740409 (  0.609634)
black                     -0.177025 (  0.409981)
hispanic                  -0.854413 (  0.558865)
age                        0.069895 (  0.151407)
married                    0.334960 (  0.302610)
nodegree                  -0.742990 (  0.259169)
re74                      -0.005430 (  0.040595)
re75                       0.093901 (  0.076997)
re74_X_re75               -0.006890 (  0.004309)
ue74                      -0.139205 (  0.464843)
ue75                      -1.019283 (  0.893028)
ue74_X_ue75                0.881073 (  0.931125)

Specification test for p-score (H_0 : lambda_s = 0)
-------------------------------------------------------------------------------------------
chi-square(12) =  28.701609   p-value:  0.004362

Summary statistics study/treated re-weighting
-------------------------------------------------------------------------------------------
Kish's effective study/treated sample size = 171

Percentiles of N_s * pi_s distribution
 1 percentile = 0.5850
 5 percentile = 0.6669
10 percentile = 0.6984
25 percentile = 0.7843
50 percentile = 0.9510
75 percentile = 1.1906
90 percentile = 1.2351
95 percentile = 1.5798
99 percentile = 1.9644

CONTROL (auxiliary) sample tilt
-------------------------------------------------------------------------------------------

Independent variable       Coef.    ( Std. Err.) 
-------------------------------------------------------------------------------------------
constant                   0.826404 (  0.727495)
black                     -0.326680 (  0.391582)
hispanic                  -0.923673 (  0.549713)
age                        0.091238 (  0.164059)
married                    0.138614 (  0.334748)
nodegree                  -0.619579 (  0.271710)
re74                      -0.017249 (  0.070147)
re75                       0.066937 (  0.093629)
re74_X_re75               -0.003591 (  0.008353)
ue74                      -0.315436 (  0.479376)
ue75                      -0.831875 (  0.916401)
ue74_X_ue75                0.806391 (  0.970829)

Specification test for p-score (H_0 : lambda_a = 0)
-------------------------------------------------------------------------------------------
chi-square(12) =  13.481717   p-value:  0.335020

Summary statistics auxiliary/control re-weighting
-------------------------------------------------------------------------------------------
Kish's effective auxiliary/control sample size = 251

Percentiles of N_a * pi_a distribution
 1 percentile = 0.7250
 5 percentile = 0.7575
10 percentile = 0.8068
25 percentile = 0.8961
50 percentile = 0.9425
75 percentile = 1.0641
90 percentile = 1.2457
95 percentile = 1.3853
99 percentile = 1.5891

Means & standard deviations of t_W (pre-balance)                                           
-------------------------------------------------------------------------------------------
                            Treated (D = 1)        Control (D = 0)        Norm. Diff.      
-------------------------------------------------------------------------------------------
constant                   1.0000 (  0.0000)      1.0000 (  0.0000)         nan
black                      0.8432 (  0.3636)      0.8269 (  0.3783)      0.0440
hispanic                   0.0595 (  0.2365)      0.1077 (  0.3100)     -0.1749
age                        2.5816 (  0.7136)      2.5054 (  0.7044)      0.1075
married                    0.1892 (  0.3917)      0.1538 (  0.3608)      0.0939
nodegree                   0.7081 (  0.4546)      0.8346 (  0.3715)     -0.3047
re74                       2.0956 (  4.8734)      2.1070 (  5.6770)     -0.0022
re75                       1.5321 (  3.2105)      1.2669 (  3.0970)      0.0841
re74_X_re75               13.1186 ( 50.6452)     14.5303 ( 60.8245)     -0.0252
ue74                       0.7081 (  0.4546)      0.7500 (  0.4330)     -0.0944
ue75                       0.6000 (  0.4899)      0.6846 (  0.4647)     -0.1772
ue74_X_ue75                0.5892 (  0.4920)      0.6577 (  0.4745)     -0.1417

Means and standard deviations of t_W (post-balance)                                        
-------------------------------------------------------------------------------------------
                            Treated (D = 1)        Control (D = 0)        Efficient (D = 1)
-------------------------------------------------------------------------------------------
constant                   1.0000 (  0.0000)      1.0000 (  0.0000)      1.0000 (  0.0000)    
black                      0.8337 (  0.3723)      0.8337 (  0.3723)      0.8337 (  0.3723)    
hispanic                   0.0876 (  0.2828)      0.0876 (  0.2828)      0.0876 (  0.2828)    
age                        2.5371 (  0.7000)      2.5371 (  0.7019)      2.5371 (  0.7092)    
married                    0.1685 (  0.3743)      0.1685 (  0.3743)      0.1685 (  0.3743)    
nodegree                   0.7820 (  0.4129)      0.7820 (  0.4129)      0.7820 (  0.4129)    
re74                       2.1023 (  5.3452)      2.1023 (  5.3569)      2.1023 (  5.3576)    
re75                       1.3771 (  3.1262)      1.3771 (  3.0947)      1.3771 (  3.1474)    
re74_X_re75               13.9434 ( 55.5826)     13.9434 ( 56.2502)     13.9434 ( 56.8189)    
ue74                       0.7326 (  0.4426)      0.7326 (  0.4426)      0.7326 (  0.4426)    
ue75                       0.6494 (  0.4771)      0.6494 (  0.4771)      0.6494 (  0.4771)    
ue74_X_ue75                0.6292 (  0.4830)      0.6292 (  0.4830)      0.6292 (  0.4830)    
/Users/bgraham/Dropbox/Sites/software/ipt/ipt/att.py:658: RuntimeWarning: invalid value encountered in true_divide
  NormDif_t    = (mu_t_D1 - mu_t_D0)/np.sqrt((mu_t_D1_std**2 + mu_t_D0_std**2)/2)

The next snippet of code plots histograms of $N_{a} \times \pi_{a}$ for both the observational and experimental estimates. Under perfect covariate balance these normalized weights should be symmetrically centered around 1 (and tightly so in larger samples). When the PSID controls are used the distribution of these weights is very right skewed: most weights are close to zero, reflecting the fact that these units do not have any similar treated units. A small number of weights are very large (~100). This reflects the fact that, in practice, only a small number of PSID control units meaningfully contribute to the value of the ATT estimate.

We can use Leslie Kish's formula from the book Survey Sampling to compute an estimate of the "effective" number of control units after re-weighting (ignoring any clustering effects). This gives

$$ \textrm{Effective number of controls} = \frac{(\sum_{i=1}^{N_a} \pi_{a,i})^{2}} { \sum_{i=1}^{N_a} \pi_{a,i}^{2}} $$

Note that the numerator sums to 1 by contruction (and hence does not need to be calculated). In the present case Kish's formula gives an effective sample size of $24$ for the PSID controls. Clearly the $2,490 $ units in the PSID are not a very good set of controls; our inferences would be just as accurate as if we instead had a random sample of $24$ units from the NSW treated population to use as controls. Kish's effective control sample size is reported in the ipt.att() estimation output above.


In [10]:
# Create a figure object to plot histograms of auxiliary weights
auxiliary_tilt_fig = plt.figure(figsize=(8, 3.5))    

# ------------------------------------- #
# - EXPERIMENTAL NSW CONTROLS SUBPLOT - #
# ------------------------------------- #

# Histogram of Na*pi_a (with no selection Na*pi_a should be close to one)
ax1 = auxiliary_tilt_fig.add_subplot(1,2,1)
sns.distplot(log_wgts_nsw, kde=False, color='#FDB515', norm_hist = True)
plt.plot((0, 0), (0, 7.55), 'k-') # plot vertical line at x = 1

# Set axis limits and tick marks (log scale with tick labels in levels)
plt.ylim([0,7.5])
plt.xlim(np.log([0.5,2]))

ax1.set_xticks(np.log([0.5, 0.75, 1.0, 1.5, 2.0]))
ax1.set_xticklabels([0.5, 0.75, 1.0, 1.5, 2.0])

# Add title and axis labels to histogram
plt.title('Experimental: NSW Controls', fontsize=12)
plt.xlabel("Normalized weights, " r"$N_{a} \times \pi_{a} $")
plt.ylabel('Density')

# Clean up the plot, add frames etc.
ax1 = plt.gca()
ax1.patch.set_facecolor('gray')               # Color of background
ax1.patch.set_alpha(0.15)                     # Translucency of background
ax1.grid(False)                               # Remove gridlines from plot

# Add frame around plot
for spine in ['left','right','top','bottom']:
    ax1.spines[spine].set_visible(True)
    ax1.spines[spine].set_color('k')
    ax1.spines[spine].set_linewidth(2)

# ---------------------------------------- #
# - OBSERVATIONAL PSID CONTROLS SUBPLOT  - #
# ---------------------------------------- #

# Histogram of Na*pi_a (with no selection Na*pi_a should be close to one)
ax2 = auxiliary_tilt_fig.add_subplot(1,2,2)
sns.distplot(log_wgts, kde=False, color='#FDB515', norm_hist = True)
plt.plot((0, 0), (0, 0.126), 'k-') # plot vertical line at x = 1

# Set axis limits and tick marks (log scale with tick labels in levels)
plt.ylim([0,0.125])
plt.xlim(np.log([1e-11,1000]))

ax2.set_xticks(np.log([1e-11,1e-9,1e-7,1e-5,1e-3,1e-1,1,10,100]))
ax2.set_xticklabels([1e-11,1e-9,1e-7,1e-5,1e-3,1e-1,1,10,100])

# Add title and axis labels to histogram
plt.title('Observational: PSID Controls', fontsize=12)
plt.xlabel("Normalized weights, " r"$N_{a} \times \pi_{a} $")
plt.ylabel('Density')

# Clean up the plot, add frames etc.
ax2 = plt.gca()
ax2.patch.set_facecolor('gray')               # Color of background
ax2.patch.set_alpha(0.15)                     # Translucency of background
ax2.grid(False)                               # Remove gridlines from plot

# Add frame around plot
for spine in ['left','right','top','bottom']:
    ax2.spines[spine].set_visible(True)
    ax2.spines[spine].set_color('k')
    ax2.spines[spine].set_linewidth(2)
    
# Render & save plot
plt.tight_layout()
plt.savefig(workdir+'Fig_LaLonde_Weights_Hist.png')


Small Monte Carlo Study

The next snippet of code executes a small Monte Carlo study. It is based upon the "dense" design used in Athey, Imbens and Wagner (2016), but with the number of balancing constraints small relative to the sample size (i.e., with $ \dim \left \{ t(W) \right \} = 1 + M $ much less than $N$). Note that in this design the propensity score is not logit with an index linear in the balancing covariates (although it is fairly well-approximated by such a model). Hence the specification test for the propensity score based on the auxiliary tilting parameters should (hopefully!) reject with frequency greater than size.

For a design with $M = 50$ and $N = 1000$ the AST procedure is computationally reliable, approximately mean and median unbiased, and produces a confidence interval with actual coverage close to target.


In [11]:
B      = 1000
N      = 1000
P      = 50

# normalize matrix to store Monte Carlo results
Simulation_Results = np.zeros((B,5))

# generate parameter vectors for covariate and outcome DGPs
delta_dense = []
beta_dense  = []
beta_harm   = []
X_names     = []

for p in range(0,P):
    delta_dense.append(4/np.sqrt(N))
    beta_dense.append(1/np.sqrt(p+1))
    X_names.append(str(p+1))

# rescale beta vector to have length 10
beta_dense =  10*(beta_dense/np.linalg.norm(beta_dense))
beta_dense = np.reshape(beta_dense,(-1,1)) # Turn into two-dimensional array, P x 1

# ------------------------------------ #
# PERFORM B MONTE CARLO REPLICATIONS - #
# ------------------------------------ #

for b in range(0,B):
    # generate treatment indicator
    W    = np.random.binomial(1, 0.5, (N,1))
    
    # generate covariate vector (with 2-cluster structure)   
    C0   = np.random.binomial(1, 0.8, (N,1))
    C1   = np.random.binomial(1, 0.2, (N,1))

    mu_X = (1 - W) * (C0 * np.zeros((1,P)) + (1 - C0) * delta_dense) \
           + W * (C1 * np.zeros((1,P)) + (1 - C1) * delta_dense) 
    X    = mu_X + np.random.multivariate_normal(np.zeros(P),np.identity(P),N)
    
    # generate outcome vector
    Y    = np.dot(X,beta_dense) + W*10 + np.random.normal(0, 1, (N,1))

    # add constant to X matrix for estimation purposes
    X    = np.concatenate((np.ones((N,1)), X), axis=1)                         
    
    # Convert W, Y and X to pandas objects
    W    = pd.Series(np.ravel(W), name="W")
    Y    = pd.Series(np.ravel(Y), name="Y")
    X    = pd.DataFrame(X)
    X.rename(columns = lambda x: str(x), inplace=True)   # Ensures column names are strings
    
    # compute AST estimate of the ATT
    try:
        [gamma_ast, vcov_gamma_ast, pscore_tests, tilts, exitflag] = \
                                                        ipt.att(W, Y, X, X, study_tilt=False, rlgrz=1, s_wgt=None, nocons=True, \
                                                                c_id=None, silent=True)
    
        Simulation_Results[b,0] = (exitflag==1)           # Record whether AST computation was successful
        Simulation_Results[b,1] = gamma_ast               # Record ATT estimate
        Simulation_Results[b,2] = np.sqrt(vcov_gamma_ast) # Record standard error estimates 
            
        # Check coverage of 95 percent Wald asymptotic confidence interval
        Simulation_Results[b,3] = (gamma_ast - 1.96*Simulation_Results[b,2] <= 10) & \
                                  (gamma_ast + 1.96*Simulation_Results[b,2] >= 10)
        # Check size/power properties of auxiliary tilt specification test
        Simulation_Results[b,4] = (pscore_tests[1][2] <= 0.05)
    except:        
        Simulation_Results[b,0] = False
                     
# ------------------------------------------------ #
# - SUMMARIZE RESULTS OF MONTE CARLO SIMULATIONS - #
# ------------------------------------------------ #
        
print("")
print("Number of times AST estimate of the ATT was successfully computed")
print(np.sum(Simulation_Results[:,0], axis = 0)) 
                      
# find simulation replicates where computation was successful
b = np.where(Simulation_Results[:,0])[0]
SimRes = Simulation_Results[b,:]
                      
# create Pandas dataframe with AST Monte Carlo results
SR=pd.DataFrame({'ATT'      :  SimRes[:,1], 'StdErr'         :  SimRes[:,2],\
                 'Coverage' :  SimRes[:,3], 'Size of p-test' :  SimRes[:,4],})                      

print("")
print("Monte Carlo summary statistics for AST")
print(SR.describe())

Q = SR[['ATT']].quantile(q=[0.05,0.95])
print("")
print('Robust estimate of standard deviation of ATT')
print((Q['ATT'][0.95]-Q['ATT'][0.05])/(2*1.645))


Number of times AST estimate of the ATT was successfully computed
1000.0

Monte Carlo summary statistics for AST
               ATT       StdErr     Coverage  Size of p-test
count  1000.000000  1000.000000  1000.000000     1000.000000
mean      9.999604     0.074090     0.948000        0.241000
std       0.074871     0.002884     0.222138        0.427904
min       9.709353     0.065687     0.000000        0.000000
25%       9.946656     0.072050     1.000000        0.000000
50%      10.000158     0.073936     1.000000        0.000000
75%      10.051836     0.075893     1.000000        0.000000
max      10.219297     0.084270     1.000000        1.000000

Robust estimate of standard deviation of ATT
0.0733624474606

In [12]:
# This imports an attractive notebook style from Github
from IPython.display import HTML
from urllib.request import urlopen
html = urlopen('http://bit.ly/1Bf5Hft')
HTML(html.read().decode('utf-8'))


Out[12]: