3-to-1 propensity score matching and sensitivity analysis

M

msarin

New Member
#1
#1
I’m doing a 3-to-1 propensity score matching to compare general (control) and non-general (treatment) anesthesia groups in lumbar spine surgery. I've matched on a number of factors, and I've calculated SMD values for continuous variables and standardized differences for categorical variables. T tests could not be used due to the violation of the assumption of independence between the compared groups. The c statistic is .679, which is within the acceptable range for PS matching. I was given a couple of suggestions to improve the impact since the database is not very granular and PS matching can't adjust unmeasured bias. Also, the post-match results demonstrated that non-GA patients were fully matched -- none was unmatched. That rarely occurs even with PS matching. This is due to highly disproportionate numbers. In this case the generalizability could be an issue because the matched GA patients might have been very special / rare cases but remained in the final analysis simply because there were a lot of GA patients. Therefore, I have been asked to add the results of a sensitivity analysis (I was referred to https://academic.oup.com/aje/article/174/3/345/247053) and to also add the results of a multivariate analysis using PS as a confounding factor (such as Outcome ~ GA/non-GA + PS) so that we can show the effect of the non-GA/GA variable as a coefficient or odds ratio. After my initial matching I put together the following:

Code: 
* Perform sensitivity analysis;

ods graphics off;

proc logistic data=data09_19 (keep = anest_grpn gender outpatient ASA_1-ASA_5 ASA_NA CPT22102 CPT62380 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0

c63057_1 c63057_ge2 c0275T age1 bmi readm complication tothlos where=(gender NE . and bmi NE . and tothlos NE .));

class anest_grpn gender outpatient ASA_1-ASA_5 CPT22102 CPT62380 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1 c63057_ge2 c0275T;

model anest_grpn (ref='1') = age1 gender bmi outpatient ASA_1-ASA_5 CPT22102 CPT62380 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1

c63057_ge2 c0275T;

output out=scores_s p=pscore;

proc psmatch data=scores_s region=allobs (psmax=1);

class anest_grpn;

psdata treatvar=anest_grpn (Treated='2') ps=pscore;

match method=greedy (k=3) distance=lps caliper=0.25;

assess lps var = (age1 gender bmi outpatient ASA_1-ASA_5 CPT22102 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 CPT62380 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1 c63057_ge2 c0275T)

/ stddev=pooled (allobs=no) stdbinvar=no varinfo;

output out (obs=match) = outgss lps=_lps matchid=_matchID;

ods output varinfo=var_info_s;

ods output stddiff=smd_s;

ods graphics on;

* Perform conditional logistic regression for specified discrete outcome variable;

%macro clr(v);

proc logistic data=outgss;

strata _matchID;

model &v (event='1') = anest_grpn;

ods output oddsratios=or_&v;

proc print data=or_&v;

proc export data=or_&v

outfile="E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\or_&v..xls"

dbms=xls replace;

run;

%mend clr;

%clr(readm);

%clr(complication);

run;

* Perform generalized random block design with fixed block effects for continuous outcome variable (LOS);

ods graphics on / discretemax=2300;

proc glm data=outgss;

class _matchID anest_grpn;

model tothlos = _matchID anest_grpn _matchID*anest_grpn;

means anest_grpn _matchID*anest_grpn;

ods output modelanova=manova_los;

run;

ods graphics off;

proc print data=manova_los;

proc export data=manova_los

outfile='E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\manova_los.xls'

dbms=xls replace;

run;

* Perform multivariate analysis of specified outcome variable using propensity score as a confounding factor;

%macro glmps(v);

proc glm data=outgss outstat=glmout_&v;

class anest_grpn;

model &v = anest_grpn pscore / ss3;

lsmeans anest_grpn;

ods output overallanova=anovao_&v modelanova=anovam_&v;

proc print data=anovao_&v;

proc export data=anovao_&v

outfile="E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\anovao_&v..xls"

dbms=xls replace;

proc print data=anovam_&v;

proc export data=anovam_&v

outfile="E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\anovam_&v..xls"

dbms=xls replace;

run;

%mend glmps;

%glmps(readm);

%glmps(complication);

%glmps(tothlos);

run;
I have a few questions:

1. Is each of these steps correct?

2. Is the generalized random block design with fixed block effects the correct choice for the analysis of the continuous outcome variable?

3. In each of the 4 steps above, what are the appropriate statistical results to report?
 
Last edited by a moderator:
hlsmith

hlsmith

Less is more. Stay pure. Stay poor.
#2
#2
Can you provide the code with color coding which would have been in the syntax and provide some basic formatting (indentations, capitalizations of functions, etc.), so it doesn't look like a big hot mess of text!

I personally am not a fan of PS matching. Using inverse PS as weights helps ensure the use of more observations - typically. Given this, did they want your sensitivity analyses to be the subset of matched patients with PS as a coefficient? I would ignore that and use the PSs as a weight not a coefficient.

Love plots and examination of overlap in PSs via histograms is the standard for reporting balance. Some may say Love plots are moot and the PSs overlap is more important.
 
hlsmith

hlsmith

Less is more. Stay pure. Stay poor.
#3
#3
Just skimmed the referenced AJE paper. What is you dependent variable and how is it formatted? Does totLOS represent length of stay? The AJE approach appears to require some contextual knowledge and assumption on your part to quantify the magnitude of an outstanding confounder. If your outcome was binary, which I am guessing it is not - you could get away with calculating an Evalue (Vander Weele). Otherwise the AJE approach would require a little tweaking perhaps since the presented example was for a binary dependent variable. Also, the approach was applied to PSs not used in matching, right? I didn't look at their supplemental documents.
 
M

msarin

New Member
#4
#4
* Perform sensitivity analysis;

ods graphics off;

PROC LOGISTICS data=data09_19
(keep = anest_grpn gender outpatient ASA_1-ASA_5 ASA_NA CPT22102 CPT62380 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1 c63057_ge2 c0275T age1 bmi readm complication tothlos where=(gender NE . and bmi NE . and tothlos NE .));

CLASS anest_grpn gender outpatient ASA_1-ASA_5 CPT22102 CPT62380 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1 c63057_ge2 c0275T;

MODEL anest_grpn (ref='1') = age1 gender bmi outpatient ASA_1-ASA_5 CPT22102 CPT62380 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1 c63057_ge2 c0275T;

OUTPUT out=scores_s p=pscore;

PROC PSMATCH data=scores_s region=allobs (psmax=1);

CLASS anest_grpn;

PSDATA treatvar=anest_grpn (Treated='2') ps=pscore;

MATCH method=greedy (k=3) distance=lps caliper=0.25;

ASSESS lps var = (age1 gender bmi outpatient ASA_1-ASA_5 CPT22102 CPT63005 CPT63011 CPT63012 CPT63017 CPT63030 CPT63047 CPT63056 CPT62380 c63035_0 c63035_1 c63035_ge2 c63048_0 c63048_1 c63048_ge2 c63057_0 c63057_1 c63057_ge2 c0275T) / stddev=pooled (allobs=no) stdbinvar=no varinfo;

OUTPUT out (obs=match) = outgss lps=_lps matchid=_matchID;

ods output varinfo=var_info_s;

ods output stddiff=smd_s;

ods graphics on;


* Perform conditional logistic regression for specified discrete outcome variable;

%macro clr(v);

PROC LOGISTIC data=outgss;

STRATA _matchID;

MODEL &v (event='1') = anest_grpn;

ods output oddsratios=or_&v;

PROC PRINT data=or_&v;

PROC EXPORT data=or_&v

outfile="E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\or_&v..xls"

dbms=xls replace;

run;

%mend clr;

%clr(readm);

%clr(complication);

run;

* Perform generalized random block design with fixed block effects for continuous outcome variable (LOS);

ods graphics on / discretemax=2300;

PROC GLM data=outgss;

CLASS _matchID anest_grpn;

MODEL tothlos = _matchID anest_grpn _matchID*anest_grpn;

MEANS anest_grpn _matchID*anest_grpn;

ods output modelanova=manova_los;

run;

ods graphics off;

PROC PRINT data=manova_los;

PROC EXPORT data=manova_los

outfile='E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\manova_los.xls'

dbms=xls replace;

run;

* Perform multivariate analysis of specified outcome variable using propensity score as a confounding factor;

%macro glmps(v);

PROC GLM data=outgss outstat=glmout_&v;

CLASS anest_grpn;

MODEL &v = anest_grpn pscore / ss3;

LSM anest_grpn;

ods output overallanova=anovao_&v modelanova=anovam_&v;

PROC PRINT data=anovao_&v;

PROC EXPORT data=anovao_&v

outfile="E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\anovao_&v..xls"

dbms=xls replace;

PROC PRINT data=anovam_&v;

PROC EXPORT data=anovam_&v

outfile="E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\anovam_&v..xls"

dbms=xls replace;

run;

%mend glmps;

%glmps(readm);

%glmps(complication);

%glmps(tothlos);

run;


In answer to your question, I have provided the exact description of what I've been asked to do. I've also produced these plots:

1630602968842.png

1630603003469.png 1630602968842.png 1630603003469.png
 
M

msarin

New Member
#5
#5
hlsmith said:
Just skimmed the referenced AJE paper. What is you dependent variable and how is it formatted? Does totLOS represent length of stay? The AJE approach appears to require some contextual knowledge and assumption on your part to quantify the magnitude of an outstanding confounder. If your outcome was binary, which I am guessing it is not - you could get away with calculating an Evalue (Vander Weele). Otherwise the AJE approach would require a little tweaking perhaps since the presented example was for a binary dependent variable. Also, the approach was applied to PSs not used in matching, right? I didn't look at their supplemental documents.
Click to expand...

* Perform generalized random block design with fixed block effects for continuous outcome variable (LOS);
ods graphics on / discretemax=2300;
proc glm data=outgss;
class _matchID anest_grpn;
model tothlos = _matchID anest_grpn _matchID*anest_grpn;
means anest_grpn _matchID*anest_grpn;
ods output modelanova=manova_los;
run;
ods graphics off;
proc print data=manova_los;
proc export data=manova_los
outfile='E:\LogonData\UserFolders\sarinm\anesthesia_lumbar\manova_los.xls'
dbms=xls replace;
run;

I think you're referring to this section. The dependent variable here is tothlos, length of stay, which is continuous. I got this approach from https://people.math.carleton.ca/~sm... Some Advanced Experimental Designs.htm#GRBDf, where I used the section GRBD with fixed block effects.
 
hlsmith

hlsmith

Less is more. Stay pure. Stay poor.
#6
#6
What is on the y-axis of your waterfall plot?

What is the purpose of the analyses in Post #5. Seems to be a cross-level interaction - that would reveal if the effect differed by match group. What would this add to the analysis? I am not familiar with the approach.

Thanks.
 
M

msarin

New Member
#7
#7
The y axis contains propensity scores; the x axis represents frequencies. As for post 5, I had been searching for the type of analysis that should be used for sensitivity analysis involving a continuous outcome variable. I had found this in https://people.math.carleton.ca/~sm...de for Some Advanced Experimental Designs.htm:
  • RCBD with fixed block effects, with subsampling
  • proc glm data=yourdata;
    class block tx;
    model y = block tx block*tx;
    test h=tx e=block*tx;
    run;
 
M

msarin

New Member
#8
#8
Actually, this is the one that I used.

  • GRBD with fixed block effects
  • proc glm data=yourdata;
    class block tx;
    model y = block tx block*tx;
    means tx block*tx;
    run;
 
You must log in or register to reply here.