Step 3: Applying PS covariate adjustment
Next, use the produced propensity score to apply propensity score covariate adjustment method. First try doing this on your own using the simulated data sets.
Reveal example code below if you need a hint or want to check your work.
Example code (click to reveal)
* Obtaining adjusted odds ratio using the output in cell 15E, can also use the trimmed dataset from exercise 2;
proc logistic data=PS_out;
class Exposure (ref='0') / param=ref;
model i_diab (event='1') = Exposure ps;
run;
Step 4: Applying PS matching
Next, apply propensity score matching method once with greedy matching with caliper and once with optimal matching without caliper. You can estimate the propensity score in a PROC LOGISTIC and then match on propensity score. Instead, we can use PROC PSMATCH and produce some diagnostics automatically. First try doing this on your own.
Reveal example code below if you need a hint or want to check your work.
Example code (click to reveal)
Example code: Greedy matching with calipers
****************** Matching method: GREEDY, WITH CALIPERS******************;
ods graphics on;
proc psmatch data=cohort region=treated;
class exposure 4_drug 4_diag htn cholest pvd chf loop;
psmodel exposure (Treated='1')= year 4_drug 4_diag htn cholest pvd chf loop;
match distance=lps method=greedy(k=1 order=descending) CALIPER=0.20;
assess lps var=(year class 4_drug 4_diag htn cholest pvd chf loop)
/ plots(nodetails)=all;
output out(obs=match)=PS_out matchid=_MatchID_ ;
run;
**********************************************************************************************
*********************** Adjusted odds ratio estimate, ATT*******************************
**********************************************************************************************;
proc logistic data=PS_out;
class exposure(ref='0');
model i_diab(event='1') = exposure;
ods select OddsRatios;
run;
**********************************************************************************************
*****A better approach is to take matching into account in the standard errors*****
***********e.g. this is a GEE model that treats matched pairs as clusters************
***********************************************************************************************;
ods trace on / listing;
proc genmod data=PS_out;
class exposure(ref='0') _MatchID_;
model i_diab (event='1') = exposure / link=logit dist=bin cl;
REPEATED subject = _MatchID_ / type=exch corrw;
LSMEANS TRT / diff ilink exp cl;
ods select ResponseProfile GEEEmpPEst GEEExchCorr LSMeans Diffs;
run;
***************** Matching method: OPTIMAL, WITHOUT CALIPERS*************;
proc psmatch data=cohort region=treated;
class exposure 4_drug 4_diag htn cholest pvd chf loop;
psmodel exposure(Treated='1')= year 4_drug 4_diag htn cholest pvd chf loop;
match distance=lps method=OPTIMAL(k=1) caliper=.;
assess lps var=(year class 4_drug 4_diag htn cholest pvd chf loop);
output out(obs=match)=PS_out_optimal matchid=_MatchID_ ;
run;
***************************Adjusted odds ratio estimate, ATT************************;
proc logistic data=PS_out_optimal;
class exposure(ref='0');
model i_diab (event='1') = exposure;
ods select OddsRatios;
run;
[END]
Step 5: Applying PS stratification
Apply propensity score stratification method using quintiles with 5 strata. Produce both ATT and ATE when pooling the stratum-specific estimates. Instead of estimating propensity scores in a PROC LOGISTIC and then stratifying, we can use PROC PSMATCH.
First try doing this on your own. Reveal example code below if you need a hint or want to check your work.
Example code (click to reveal)
**************************Estimate ATE "stratumwgt=total"****************************;
PROC PSMATCH data=Cohort REGION=ALLOBS ;
class exposure 4_drug 4_diag htn cholest pvd chf loop;
psmodel exposure (Treated='1')= year 4_drug 4_diag htn cholest pvd chf loop;
STRATA nstrata=5 key=total stratumwgt=total;
assess ALLCOV PS / plots=(barchart STDDIFFPLOT) ;
ods exclude BarChart;
output out(obs=all)=PS_out ;
run;
********Obtaining the strata-specific OR estimates, then pooling************;
proc freq data=PS_out;
tables _STRATA_*exposure*i_diab / nopercent nocol relrisk cmh;
run;
*******************Estimate ATT "stratumwgt=treated"**********************;
PROC PSMATCH data=Cohort REGION=ALLOBS ;
class exposure 4_drug 4_diag htn cholest pvd chf loop;
psmodel exposure (Treated='1')= year 4_drug 4_diag htn cholest pvd chf loop;
STRATA nstrata=5 key=total stratumwgt=treated;
assess ALLCOV PS / plots=(barchart STDDIFFPLOT) ;
ods exclude BarChart;
output out(obs=all)=PS_out ;
run;
********Obtaining the strata-specific OR estimates, then pooling************;
proc freq data=PS_out;
tables _STRATA_*exposure*i_diab / nopercent nocol relrisk cmh;
run;
Step 6: Applying IPTW method
You can simply apply a logistic model to produce propensity scores and then calculate IPTW and use a weighted outcome model.
First try doing this on your own. Reveal example code below if you need a hint or want to check your work.
Example code (click to reveal)
**********************Step 1: estimate PS, using the PS output from 15E****
**********************Step 2: calculate the weights******************************;
data ps_weight;
set PS_out;
if exposure = ""1"" then ps_weight=1/ps;
if exposure = ""0"" then ps_weight=1/(1 - ps);
run;
************************Step 3: weighted outcome model*****************************;
proc logistic data=ps_weight;
class exposure(ref='0') / param=ref;
model i_diab(event='1') = exposure;
weight ps_weight / normalize;
ods select MoldeInfo NObs OddsRatios;
run;
**********Instead, you can use proc PSmatch to replace steps 1 and 2**********;
PROC PSMATCH data=cohort REGION=ALLOBS ;
class exposure 4_drug 4_diag htn cholest pvd chf loop;
psmodel exposure (Treated='1')= year 4_drug 4_diag htn cholest pvd chf loop;
psweight weight=_ATEWgt_ ;
assess lps ALLCOV / plots=all ;
output out(obs=REGION)=PS_out atewgt=_ATEWgt_ ;
run;
********************************Weighted outcome model****************************
** Obtaining adjusted ATE odds ratio;
proc logistic data=PS_out;
class exposure(ref='0') / param=ref;
model i_diab(event='1') = exposure;
weight _ATEWgt_ / normalize;
ods select MoldeInfo NObs OddsRatios;
run;
*******************************************************************************************
Note: Using normalized weights (with the option normalize) means we don't need the step of stabilizing weights, but here is an example code to stabilize weights;
proc sql noprint;
select mean(_ATEWgt_)
into :ATEWgt_mu
from PS_out ;
quit;
data PS_out; set PS_out;
stdzd_weight = _ATEWgt_ / &ATEWgt_mu;
run;"