Chapter 2: SAS Code

Chapter 2: SAS Code

Sample dataset codebook:

treat = Binary indicator of treatment versus control group

x1-x5 = continuous confounders associated with Treat

cont_out = Continuous outcome of interest

bin_out = Binary outcome of interest

Estimating the propensity score in SAS with logistic regression

SAS> proc logistic data=mydata descending;

model treat=x1 x2 x3 x4 x5;

output out=mydata_ps prob=pscore;

run;

MATCHING

K:1 matching without replacement using GMATCH macro

/* Include %GMATCH

SAS> %include ‘gmatch.sas’;

/* Have previously estimated propensity score “pscore”

/* 1:1 nearest neighbor matching, matching on propensity score

SAS> %gmatch(data=mydata_ps, group=treat, id=id, mvars=pscore, seedca=111222, seedco=111333, out=matched);

/* 2:1 matching, matching on PS, caliper=0.10 units

SAS> %gmatch(data=mydata_ps, group=treat, id=id; mvars=pscore, dmaxk=0.10, ncontls=2, seedca=111222, seedco=111333, out=matched);

Caliper matching using GMATCH macro

/* Want to specify caliper in terms of SD of PS logit

/* Calculate standard deviation of logit of propensity score

SAS> proc means std data=mydata_ps;

var logit_ps;

output out=std_ps (keep=std) std=std;

run;

/* Create new variable that is 0.20 SD of logit of PS

SAS> data std_ps; set std_ps;

std=0.2*std;

run;

/*Create macro variable that specifies caliper size for matching

SAS> data _null_;

set std_ps;

call symput(‘stdcal’,std);

run;

/* 2:1 caliper matching, caliper=0.20 SD of logit of PS

SAS> %gmatch(data=mydata_ps, group=treat, id=id, mvars=logit_ps, dmaxk=&stdcal, ncontls=2, seedca=111222, seedco=111333, out=matched);

K:1 matching, with and without replacement using PSMatching macro

/* Include %PSMatching macro

SAS> %include ‘PSMatching.sas’;

/* Have previously estimated propensity score “pscore”

/* Separate datasets for treated and control, only contain ID and propensity score (or alternatively logit of PS)

SAS> data mydataC mydataT;

set mydata_ps;

if treat=0 then do;

idC=id; pscoreC=pscore; output mydataC;

end;

if treat=1 then do;

idT=id; pscoreT=pscore; output mydataT;

end;

/* 1:1 Nearest neighbor matching without replacement

SAS> %PSMatching(datatreatment=mydataT, datacontrol= mydataC, method=nn, numberofcontrols=1, replacement=no, out=matches);

Caliper matching using PSMatching macro

/* 2:1 caliper matching with replacement, 0.10 caliper

SAS> %PSMatching(datatreatment= mydataT, datacontrol= mydataC, method=caliper, caliper=0.10, numberofcontrols=2, replacement=yes, out=matches);

Radius matching using PSMatching macro

SAS> %PSMatching(datatreatment= mydataT, datacontrol= mydataC, method=radius, caliper=0.20, replacement=no, out=matches);

Optimal matching using DIST and VMATCH macros

/* Include %DIST and %VMATCH macro

SAS> %include ‘dist.sas’;

SAS> %include ‘vmatch.sas’;

/* If have already calculated distance matrix, can run vmatch

SAS> %vmatch(dist=mydata_dist, idca=id, a=2, b=2,lilm=num_of_control, n=num_of_treat, firstco=C_first,lastco=C_last);

/* Can run %VMATCH through %DIST macro

/* Does 1:1 optimal matching

SAS> %dist(data=mydata, group=treat, id=id, mvars=pscore,

wts=1, vmatch=Y, a=1, b=1, lilm= num_of_control, dmax=0.1,

outm=mp1_b, summatch=n, mergeout=mpropen);

/* Optimal matching; each treated matched with 1-3 controls

SAS> %dist(data=mydata, group=treat, id=id, mvars=pscore,

wts=1, vmatch=Y, a=1, b=3, lilm= num_of_control, dmax=0.1,

outm=mp1_b, summatch=n, mergeout=mpropen);

Balance diagnostics

/* Sort data by treatment groups

SAS> proc sort data = mydata_ps;

by treat;

/* Boxplot of PS by treatment group

SAS> proc boxplot data= mydata_ps;

symbol width = 2;

plot pscore*treat;

run;

/* Stacked histogram of PS by treatment group

SAS> ODS graphics on;

proc univariate noprint data= mydata_ps;

var pscore;

class treat;

histogram pscore / nrows=2 ;

run;

ODS graphics off;

PROPENSITY SCORE WEIGHTING, PARAMETRIC PS ESTIMATION

/* Estimate the propensity score with logistic regression

SAS> proc logistic data=mydata descending;

model treat = x1 x2 x3 x4 x5;

output out=mydata_ps prob=pscore;

run;

/* Calculate ATE and ATT propensity score weights

SAS> data mydata_ps.w;

set mydata_ps;

w_ate= treat/pscore + (1-treat)/(1-pscore);

w_att= treat + (1-treat)*(pscore/(1-pscore));

run;

/* Use ATE weights as probability weights in final analysis

SAS> proc genmod data=mydata_ps.w,

class id;

model cont_out = treat x1 x2 x3 x4 x5 / error=B;

weight w_ate;

run;

/* Use ATT weights as probability weights in final analysis

SAS> proc genmod data=mydata_ps.w,

class id;

model cont_out = treat x1 x2 x3 x4 x5 / error=B;

weight w_att;

run;

SUBCLASSIFICATION

Creating 5 propensity score subclasses

/* After generating propensity scores, create quintiles

SAS> proc rank data = mydata_ps groups=5 out= ps_strataranks;

var pscore;

ranks ps_quint;

run;

/* Sort data by quintiles

SAS> proc sort data = ps_strataranks;

by ps_quint;

Estimating subclass-specific and overall effect estimates

/* Binary outcome: Mantel-Haenszel stratified analysis

SAS> proc freq data= ps_strataranks;

table ps_quint*treat*bin_out / nocol cmh;

run;

/* Continuous outcome: can combine subclass-specific t-tests

SAS> proc ttest data=ps_strataranks;

by ps_quint;

class treat;

var cont_out;

ods output statistics = strata_out;

/* Calculate weights (weight is inverse of square of SE of group difference); weight group differences

SAS> data weights;

set strata_out;

if class = 'Diff (1-2)';

wt_i = 1/(StdErr**2);

wt_diff = wt_i*Mean;

/* Sum weighted means

SAS> proc means noprint data = weights;

var wt_i wt_diff;

output out = total sum = sum_wt sum_diff;

/* Calculate and output overall treatment difference, SE

SAS> data overall;

set total;

diff_mean = sum_diff/sum_wt;

diff_SE = SQRT(1/sum_wt);

SAS> proc print data = overall;

run;