Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Vancouver_application_8-8-2018/submitted JASA data/Bryan Graham Stata software/bandwidth_selection_300/iptATE.ado
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
297 lines (237 sloc)
8.92 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
version 11.0 | |
capture mata mata drop m_ipt_ate() | |
capture mata mata drop m_ipt_logit() | |
capture mata mata drop m_ipt_logit_crit() | |
capture mata mata drop m_ipt_logit_phi() | |
mata: | |
void m_ipt_ate(string scalar varlist, string scalar touse, string scalar clusterid, string scalar weight, string scalar optroutine) | |
{ | |
real matrix YDh_X, h_X, t_X, D1_ipt_logit_results, D0_ipt_logit_results | |
real matrix m_delta1, m_delta0, m_hat, M11_hat, M22_hat, M31_hat, M32_hat, M33_hat, M_hat, OMEGA_hat, V_hat | |
real colvector Y, D, wgts, delta1_ipt, delta0_ipt, theta_ipt, v1_ipt, v0_ipt, pi1_ipt, pi0_ipt, m_gamma, selvec | |
real rowvector m_g | |
real scalar N, M, gamma_ipt, NG | |
string colvector t_X_names | |
string matrix pscore_D1_names, pscore_D0_names, coef_names | |
YDh_X = Y = D = h_X = . | |
st_view(YDh_X, ., tokens(varlist), touse) | |
st_subview(Y, YDh_X, ., 1) | |
st_subview(D, YDh_X, ., 2) | |
st_subview(h_X, YDh_X, ., (3\.)) | |
N = rows(h_X) | |
M = cols(h_X) | |
if (weight != "") { | |
st_view(wgts, ., tokens(weight), touse) | |
wgts = wgts :/ (colsum(wgts) :/ N) | |
} | |
else { | |
wgts = J(N,1,1) | |
} | |
t_X = (h_X,J(N,1,1)) | |
display("") | |
display("") | |
display("Computing treated subsample tilt") | |
D1_ipt_logit_results = m_ipt_logit(D, h_X, wgts, optroutine) | |
delta1_ipt = D1_ipt_logit_results[1,.]' | |
v1_ipt = t_X*delta1_ipt | |
pi1_ipt = (N^-1) :* (D :/ invlogit(v1_ipt)) | |
m_delta1 = D1_ipt_logit_results[1+2*(1+M)+1..1+2*(1+M)+N,1..1+M] | |
display("") | |
display("") | |
display("Computing control subsample tilt") | |
D0_ipt_logit_results = m_ipt_logit((1 :- D), h_X, wgts, optroutine) | |
delta0_ipt = -D0_ipt_logit_results[1,.]' | |
v0_ipt = t_X*delta0_ipt | |
pi0_ipt = (N^-1) :* ((1 :- D) :/ (1 :- invlogit(v0_ipt))) | |
m_delta0 = D0_ipt_logit_results[1+2*(1+M)+1..1+2*(1+M)+N,1..1+M] | |
gamma_ipt = sum(pi1_ipt :* Y - pi0_ipt :* Y) | |
m_gamma = wgts :* ((N :* pi1_ipt :* Y) - (N :* pi0_ipt :* (Y :+ gamma_ipt))) | |
m_hat = (m_delta1, m_delta0, m_gamma) | |
M11_hat = D1_ipt_logit_results[(1+1+M+1)..(1+2*(1+M)),.] | |
M22_hat = D0_ipt_logit_results[(1+1+M+1)..(1+2*(1+M)),.] | |
M31_hat = colsum(wgts :* ((D :* -exp(-v1_ipt)) :* Y ) :* t_X) | |
M32_hat = colsum(wgts :* (((1 :- D) :* exp(v0_ipt)) :* (Y :+ gamma_ipt)) :* t_X) | |
M33_hat = -N | |
M_hat = (M11_hat, J(1+M,1+M,0), J(1+M,1,0) \ J(1+M,1+M,0), M22_hat, J(1+M,1,0) \ M31_hat,M32_hat,M33_hat) | |
OMEGA_hat = J(2*(1+M)+1,2*(1+M)+1,0) | |
if (clusterid != "") { | |
st_view(G, ., tokens(clusterid), touse) | |
NG = max(G) | |
st_numscalar("nclust",NG) | |
for (g=1; g<=NG; g++) { | |
selvec = (G:==g) | |
m_g = colsum(m_hat :* selvec) | |
OMEGA_hat = OMEGA_hat + m_g'*m_g | |
} | |
V_hat = luinv(M_hat)*OMEGA_hat*luinv(M_hat)' | |
} | |
else { | |
OMEGA_hat = m_hat'*m_hat | |
V_hat = luinv(M_hat)*OMEGA_hat*luinv(M_hat)' | |
} | |
t_X_names = tokens(varlist) | |
t_X_names = t_X_names[.,3..1+1+M] | |
t_X_names = (t_X_names, "_cons")' | |
pscore_D1_names = (("delta1" :* J(1+M,1,1)), t_X_names) | |
pscore_D0_names = (("delta0" :* J(1+M,1,1)), t_X_names) | |
coef_names = (pscore_D1_names \ pscore_D0_names \ ("ate", "gamma")) | |
st_matrix("V_hat",V_hat) | |
st_matrixrowstripe("V_hat", coef_names) | |
st_matrixcolstripe("V_hat", coef_names) | |
theta_hat = (delta1_ipt' ,delta0_ipt', gamma_ipt) | |
st_matrix("theta_hat",theta_hat) | |
st_matrixcolstripe("theta_hat", coef_names) | |
st_numscalar("nobs",N) | |
} | |
real matrix m_ipt_logit(real colvector D, real matrix h_X, real colvector wgts, string scalar optroutine) | |
{ | |
real matrix vcov_ipt, hess_ipt, m2_ipt, A | |
real rowvector delta_ipt | |
t_Xt_X = cross(h_X,1,h_X,1) | |
if (rank(t_Xt_X) < cols(h_X) + 1) { | |
errprintf("singular or near singular matrix") | |
exit(499) | |
} | |
t_XD = cross(h_X,1,D,0) | |
delta_sv = cholsolve(t_Xt_X,t_XD)' | |
M = moptimize_init() | |
moptimize_init_evaluator(M, &m_ipt_logit_crit()) | |
moptimize_init_which(M, "max") | |
if (optroutine=="e2") { | |
moptimize_init_evaluatortype(M, "e2") | |
} | |
else { | |
moptimize_init_evaluatortype(M, "e1") | |
moptimize_init_technique(M, "bfgs") | |
} | |
moptimize_init_depvar(M, 1, D) | |
moptimize_init_eq_indepvars(M, 1, h_X) | |
moptimize_init_eq_coefs(M, 1, delta_sv) | |
moptimize_init_vcetype(M,"robust") | |
moptimize_init_userinfo(M,1,wgts) | |
moptimize(M) | |
delta_ipt = moptimize_result_eq_coefs(M) | |
vcov_ipt = moptimize_result_V(M,1) | |
hess_ipt = moptimize_result_Hessian(M) | |
m2_ipt = (moptimize_result_scores(M) :* h_X, moptimize_result_scores(M)) | |
A = (delta_ipt \ vcov_ipt \ hess_ipt \ m2_ipt) | |
return(A) | |
} | |
function m_ipt_logit_crit(transmorphic M, real scalar todo, | |
real rowvector delta, crit, foc, soc) | |
{ | |
real scalar N | |
real colvector t_Xdelta | |
real colvector D, phi, phi1, phi2, soc_i | |
real matrix phi012 | |
t_Xdelta = moptimize_util_xb(M, delta, 1) | |
D = moptimize_util_depvar(M, 1) | |
N = length(D) | |
phi012 = m_ipt_logit_phi(t_Xdelta, N) | |
phi = phi012[.,1] | |
wgts = moptimize_util_userinfo(M,1) | |
crit = moptimize_util_sum(M, wgts :* (D :* phi :- t_Xdelta)) | |
if (todo>=1) { | |
phi1 = phi012[.,2] | |
foc = wgts :* (D :* phi1 :- 1) | |
if (todo==2) { | |
phi2 = phi012[.,3] | |
soc_i = wgts :* (D :* phi2) | |
soc = moptimize_util_matsum(M, 1,1, soc_i, crit) | |
} | |
} | |
} | |
real matrix m_ipt_logit_phi(real colvector v, real scalar N) | |
{ | |
real scalar a, b, c, v_star | |
real colvector phi, phi1, phi2 | |
real matrix A | |
a = -(N - 1)*(1 + log(1/(N - 1)) + 0.5*(log(1/(N - 1)))^2) | |
b = N + (N - 1)*log(1/(N - 1)) | |
c = -(N-1) | |
v_star = log(1/(N - 1)) | |
phi = (v:>v_star) :* (v - exp(-v)) + (v:<=v_star) :* (a :+ (b:*v) + (0.5*c):*(v:^2)) | |
phi1 = (v:>v_star) :* (1 :+ exp(-v)) + (v:<=v_star) :* (b :+ (c:*v)) | |
phi2 = (v:>v_star) :* (-exp(-v)) + (v:<=v_star) :* c | |
A = phi,phi1,phi2 | |
return(A) | |
} | |
end | |
program define iptATE, eclass | |
version 11 | |
syntax varlist(numeric) [if] [in] [pweight] [, cluster(string) optroutine(string) BALance] | |
marksample touse | |
tokenize `varlist' | |
local outcomevar `1' | |
local treatdum `2' | |
macro shift | |
macro shift | |
local controlvars `*' | |
if "`balance'" != "" { | |
tempvar controldum | |
quietly gen double `controldum' = 1 - `treatdum' if `touse' | |
display "" | |
display "" | |
display "Mean covariate values by treatment status" | |
display "----------------------------------------------------------------------------------------" | |
display "- Variable T = 0 T = 1 Diff p-value -" | |
display "----------------------------------------------------------------------------------------" | |
foreach controlvar in `controlvars' { | |
quietly { | |
reg `controlvar' `controldum' `treatdum' [`pweight'] if `touse', nocons cluster(`cluster') robust | |
lincom `controldum' | |
local cm = r(estimate) | |
local cm_se = r(se) | |
lincom `treatdum' | |
local tm = r(estimate) | |
local tm_se = r(se) | |
lincom `treatdum' - `controldum' | |
local dm = r(estimate) | |
local dm_se = r(se) | |
test `treatdum' = `controldum' | |
local dm_pv = r(p) | |
} | |
display as result %25s abbrev("`controlvar'",25) " " %9.3f `cm' " " %9.3f `tm' " " %9.3f `dm' | |
display as result %25s "" " (" %9.3f `cm_se' ") (" %9.3f `tm_se' ") (" %9.3f `dm_se' ") " %9.3f `dm_pv' | |
display "" | |
} | |
display "----------------------------------------------------------------------------------------" | |
display "" | |
} | |
if "`weight'" != "" { | |
tempvar wgts | |
quietly gen double `wgts' `exp' if `touse' | |
} | |
else { | |
tempvar wgts | |
quietly gen double `wgts' = 1 if `touse' | |
} | |
if "`cluster'" != "" { | |
tempvar clusterid | |
egen `clusterid' = group(`cluster') if `touse' | |
mata: m_ipt_ate("`varlist'", "`touse'", "`clusterid'","`wgts'","`optroutine'") | |
} | |
else { | |
mata: m_ipt_ate("`varlist'", "`touse'", "`cluster'","`wgts'","`optroutine'") | |
} | |
display "" | |
display "" | |
display "Inverse probability tilting propensity score & ATE estimates" | |
display "Outcome variable : `outcomevar'" | |
display "Treatment indicator : `treatdum'" | |
display "Control variables : `controlvars'" | |
ereturn post theta_hat V_hat [`weight'], esample(`touse') | |
ereturn scalar N = nobs | |
ereturn display | |
if "`weight'" != "" { | |
display "Probability weights option specified: `weight' `exp'" | |
} | |
if "`cluster'" != "" { | |
display "Total number of primary sampling units: " nclust | |
} | |
display "Total number of observations: " nobs | |
display "" | |
display "Test of equality of two tilting coefficient vectors" | |
quietly: test [delta1]_cons = [delta0]_cons | |
test [delta1=delta0], common accum | |
end |