Calculating IPW and SMR in SAS and Stata

Author: Bailey DeBarmore

Learning about a method in class, like inverse probability weighting, is different than implementing it in practice. 

This post will remind you why we might be interested in propensity scores to control for confounding - specifically inverse probability of treatment weights and SMR - and then show how to do so in SAS and Stata.

If you have corresponding code in R that you'd like to add to this post, please contact me.

A note about weighting versus multivariable regression:
Effect estimate interpretations when you use weighting are marginal effect in the target population. When you adjust for covariates in a regression model, you are interpreting a conditional effect, that is, the effect of the exposure holding (conditional on) the covariates being constant. 

Conditional estimates are troublesome with time-varying covariates because we run into collider bias and conditioning on mediators, thus weights are preferable. In simpler situations, using weights over multivariable regression can help with convergence issues .

Files to Download: .txt file with SAS and Stata code, as well as a PDF version of this post with code (perfect for students) available to download at the end of the post or at my github

Jump to:

• SMR in SASwww.baileydebarmore.com/epicode/calculating-ipw-and-smr-in-sas#smrsas
• SMR in Stata
• IPTW in SAS
• IPTW in Stata
• Files to download

propensity scores

A propensity score is a predicted probability that may be used to predict exposure (or treatment) status, but can also be used for censoring or missingness. 

How do we use propensity scores for confounding?
We can use propensity scores to generate WEIGHTS which, when applied to the final model, make the exposure independent from confounders (Figure 1B), usually by modeling the association between the exposure and confounders (instead of the main analysis where we model the outcome and exposure). 

Propensity scores can also control for confounding via covariate adjustment (I discourage you from this option), stratification, and matching, in addition to weighting.

The 2 types of weights I'll be discussing in this post are

• SMR (standardized mortality/morbidity ratio) weights
• IPTW (inverse probability of treatment weights)

Figure 1. Panel A shows the observed population in our data set, where the relationship between exposure and outcome (orange) is confounded by well, confounders. In B, we have removed the arrow from confounders to exposure. We can remove the arrow in several ways, including using propensity scores (of various types) to create a pseudopopulation where the exposure and confounder are no longer associated.

Figure 2. Panel A shows the usual multivariable model we run in our analyses - to estimate the association of the exposure with the outcome, controlled for confounders. When we want to use propensity scores, we first create the weights that we will later use in our final model by modeling the association of the confounders with the exposure - so we can remove that arrow.

SMR

When you use SMR weights, you're estimating the average treatment effect in the treated (ATT). In other words, you estimate the effect had the exposed group been exposed, versus the exposed group been unexposed. The pseudopopulation that you create has a covariate (confounder) distribution equal to that observed in the exposed group.

You can also generate your SMR weights where the unexposed group is the target of interest, and model the effect had the unexposed group been unexposed versus the unexposed group been exposed.

Figure 3. In Panel A, the target group is the exposed group, so we use SMR to model had the exposed group been exposed versus unexposed, with the covariate distribution of the exposed group (yellow). In Panel B, we use the covariate distribution of the unexposed group (green), and model had unexposed group been unexposed vs exposed.

Figure 4. Equations to calculate the SMR weights for given target population corresponding to Figure 3 examples. A indicates exposure status (=1 for exposed, =0 for unexposed). The numerator is the probabilitiy of exposure given the covariate distribution of the target group, and the denominator is the probability of being assigned to what was observed. You can see that in panel B, it switches when the target group is unexposed.

coding SMR in SAS

To use SMR, you'll be generating the probabilities for the numerator and denominator, and then creating the weight (numerator/denominator) and applying it in the appropriate statement.

***********************************************
* Calculating SMR weights where exposure = 0, 1
**************************************(********;

&let data=<data>;
&let y=<outcome>;
&let x=<exposure>;
&let id=<id>;

*Estimate the predicted probability given covariates;

​proc logistic data=&data desc;
model &x=<covariates>;
output out=pred p=p1;
run;

*Generate the weights by exposure status, for exposed group = target their weight will be 1;

data <newdata>;
set pred;
p0 = 1-p1;
odds = p1/p0;
if &x=1 then wt=1;
else wt=odds;
run;

*Final weighted analysis;
proc logistic data=<newdata> desc;
weight wt;
model &y = &x;
run;


***********************************************
* Calculating SMR weights where exposure = categorical
**************************************(********;

&let data=<data>;
&let y=<outcome>;
&let x=<exposure>;
&let id=<id>;

*Estimate the predicted probability given covariates;

​proc logistic data=&data desc;
model &x=<covariates> /LINK= glogit;
output out=pred p=p1;
run;

*Generate the weights by exposure status, for exposed group = target their weight will be 1;

data <newdata>;
set pred;
p0 = 1-p1;
odds = p1/p0;
if &x=1 then wt=1;
else wt=odds;
run;

*Final weighted analysis;
proc logistic data=<newdata> desc;
weight wt;
model &y = &x;
run;

What is &let?
So that you can easily adapt my code, I use &let statements at the beginning of my code blocks. After the equals sign, you would replace <data> with your dataset name, <exposure> with your exposure variable, and <outcome> with your outcome variable. The code as written will then run with those chosen variables. Note that you do need to replace <covariates> in the model statement with your confounders.

If you don't want to use &let statements, simply go through the code and anywhere you see &<text>, replace both the & and the text with your regular code.

coding SMR in Stata

The approach to calculating SMR is different in Stata. You'll be using the built in teffects command (reference manual [TE] here) and using options to specify SMR versus IPW. 

Since using SMR gives us the average treatment effect in the treated, we'll use the option atet (average treatment effect on the treated) instead of ate (average treatment effect) which we'll use for IPTW.

******************************************
* Calculating SMR weights
*****************************************;
​* Syntax for teffects statement

*teffects ipw (<outcome>) (<exposure> <covariates>), atet
*where <outcome> is your outcome variable, <exposure> is your exposure variable, and <covariates> is a list of your covariates to generate your weights. 

*Example: Binary
*Outcome = lowbirthwt
*Exposure = maternalsmoke
*Covariates = maternalage nonwhite

*Use the teffects statement to generate your weights and then apply them in a logistic (default) model all in 1 step

teffects ipw (lowbirthwt) (maternalsmoke maternalage nonwhite), atet

*If your outcome is continuous, you can specify a probit model

*Example: Continuous
*Outcome = birthwt
*Exposure = maternalsmoke
*Covariates = maternalage nonwhite

teffects ipw (birthwt) (maternalsmoke maternalage nonwhite, probit), atet
​

inverse probability of treatment weights (IPTW)

In contrast to SMR weights, when you use IPTW weights you are estimating the average treatment effect (ATE), that is the treatment effect in a population with covariate distribution equal to the entire observed study population, not just the exposed or unexposed.

​Thus, you're modeling the complete counterfactual. In other words, you're estimating the effect had the entire population been exposed versus the entire population been unexposed.

You can use unstabilized or stabilized IPTW in your final model. Choosing one over the other doesn't change your ultimate interpretation but affects the width of your confidence interval (wider when you use unstabilized). 

Figure 5. Panel A shows what happens in the pseudopopulation when we use unstabilized weights. We duplicate our N, creating 2 new populations with the covariate distribution of the entire observed population, but had everyone been exposed versus everyone unexposed. In B we apply the covariate distribution within strata of exposure, maintaining the same N. (E with a bar (macron) indicates unexposed).

Figure 6. In Panel A we have the equation for unstabilized IPTW for exposed and unexposed, and in Panel B, the stabilized IPTW for exposed and unexposed. You'll see additional notes that show you that you can also calculate the weight components for unexposed by subtracting the exposed probabilities from 1. You'll see this used in our code below.

Unstabilized

Create a pseudopopulation 2x the size of our observed - one where everyone is exposed and one where everyone is unexposed.

Stabilized

Create a pseudopopulation maintaining the original population size, but we adjust the covariate distribution within each strata of exposure group by upweighting and downweighting people to match the overall covariate distribution. 

You can see in Figure 5 versus Figure 3 that we are using different covariate distributions. With SMR where we match the covariate distribution of the exposed or unexposed group (whichever is target). This is why IPTW generates an average treatment effect while SMR generates an average treatment effect in the treated.

coding IPTW in SAS

******************************************
* Calculating IPTW
*****************************************;
​
&let data=<data>;
&let y=<outcome>;
&let x=<exposure>;
&let id=<id>;

*Estimate denominator - output a dataset with results of regression called denom, with the resulting probabilities stored in variable d;
​
proc logistic data=&data desc;
model &x = <covariates>;
output out=denom p=d;
run;

*Generate numerator for stabilized weights - output a dataset with results of regression called num, with the resulting probabilities stored in variable n - note that there is nothing on the right side of the equation because the numerator will simply be P(A=a), where a = observed exposure status;

proc logistic data=&data desc;
model &x=;
output out=num p=n;
run;

*Generate stabilized and unstabilized weights by merging the datasets with regression output (merge on the unique identifier in your dataset, &id);

data <newdata>;
merge &data denom num;
by &id;
if &x=1 then do;
uw = 1/d; 
sw = n/d;
end;

*Remember we can use 1 - P(exposed) for the unexposed weight components;
​
else if &x=0 then do;
uw=1/(1-d);
sw=(1-n)/(1-d);
end;
​run;

*Check the distribution of your IPTW - the mean should be 1. Is the sum for uw twice the sum of sw? why? is the range of uw greater than sw? why?;

proc means data=<newdata> mean sum min max;
var uw sw;
run;

*You can check to see if your exposure and covariates are associated in your new pseudopopulation (<newdata>);

proc logistic data=<newdata> desc;
weight sw;
model &x=<covariates>;
run;

*Now you can run your main analyses and apply the weights using the weight statement - use sw variable for stabilized weights, and use uw for unstabilized weights - you can use proc genmod, glm, logistic, etc. I'll show you below with logistic you can see now we're using &y and &x - and we don't need the covariates because the confounder -> x arrow is encompassed in the sw weight statement;

proc logistic data=<newdata> desc;
weight sw;
model &y = &x;
run;

coding IPTW in Stata

The approach to calculating IPTW is different in Stata than SAS. You'll be using the built in teffects command (reference manual [TE] here) and using options to specify IPTW.

Since using IPTW results in the average treatment effect overall, we'll use option ate (average treatment effect) instead of  option atet (average treatment effect on the treated) that we use for SMR. 

******************************************
* Calculating IPTW
*****************************************;
​
​* Syntax for teffects statement

*teffects ipw (<outcome>) (<exposure> <covariates>), ate
*where <outcome> is your outcome variable, <exposure> is your exposure variable, and <covariates> is a list of your covariates to generate your weights. 

*Example: Binary
*Outcome = lowbirthwt
*Exposure = maternalsmoke
*Covariates = maternalage nonwhite

*Use the teffects statement to generate your weights and then apply them in a logistic (default) model all in 1 step

teffects ipw (lowbirthwt) (maternalsmoke maternalage nonwhite), ate

*If your outcome is continuous, you can specify a probit model

*Example: Continuous
*Outcome = birthwt
*Exposure = maternalsmoke
*Covariates = maternalage nonwhite

teffects ipw (birthwt) (maternalsmoke maternalage nonwhite, probit), ate

* Syntax to manually create IPTW for binary exposure (treatment)
logistic treatment vars
predict p
gen iptw = 1/p if treatment==1
replace iptw=1/(1-p) if treatment==0

*To calculate stabilized IPTW, run tab treatment and use the proportions for X below
tab treatment
gen siptw = X/p if treatment==1
replace siptw = (1-X)/(1-p) if treatment==0

---

UPDATE: Feb 22 2021
If your treatment (exposure) is categorical you'll need to change things up a bit. I'm writing this blurb in response to a comment below, and in response to this topic coming up several times in just the past week! 

Figure 7. Calculating IPTW for a categorical A.

​In your software of choice, you need to calculate the probability of being treated given your covariates, save those probabilities, and then assign the probability given each observations observed status.

In SAS, you can use code like the sample below (thanks, Paul!) which uses a generalized logit model for a treatment variable <treatment> with multiple categories. 

The output data set will have probabilities for every value of treatment. For example, a person with treatment=1 will have values for treatment=0,1,2,3 so you need to do some more coding to tell the computer you only want the predicted value for treatment if treatment=treatment.

​PROC LOGISTIC DATA=...;
    MODEL treatment = vars / LINK=glogit;
    OUTPUT OUT=denom_ipw P=d;
RUN;

In Stata, you can run mlogit to estimate the probability of <treatment> adjusted for variables <vars> and then use the postestimation commands, predict, to create a denominator variable for each value of treatment (0, 1, 2, 3). You use the option "outcome" to tell Stata this. 

mlogit treatment vars
predict p0, outcome(0)
predict p1, outcome(1)
predict p2, outcome(2)
predict p3, outcome(3)

Then create your weight variable, iptw below.

gen iptw=.
replace iptw=1/p0 if treatment==0
replace iptw=1/p1 if treatment==1
replace iptw=1/p2 if treatment==2
replace iptw=1/p3 if treatment==3

If you want to create stabilized weights, you can run a tab to get the proportion in each category, and then calculate your weights. For example, let's say my groups are distributed as P(A=0)=0.6, P(A=1)=0.14, P(A=2)=0.2, and P(A=3)=0.06, I would make a stabilized iptw variable, siptw, below.

gen siptw=.
replace siptw=0.6/p0 if treatment==0
replace siptw=0.14/p1 if treatment==1
replace siptw=0.2/p2 if treatment==2
replace siptw=0.06/p3 if treatment==3

---

Download resources here, or get the most up-to-date versions on my GitHub

Calculating SMR and IPW - SAS - 2021 - EPICODE
File Size:	8 kb
File Type:	txt

Download File

---

Calculating SMR and IPW in Stata - 2022
File Size:	6 kb
File Type:	txt

Download File

---

Coding IPW and SMR in SAS and Stata - PDF for teachers
File Size:	551 kb
File Type:	pdf

Download File

---

---

Suggested citation

DeBarmore BM. “Coding IPW and SMR in SAS and Stata”. 2019. Updated 2021. Retrieved from http://www.baileydebarmore.com/epicode/calculating-ipw-and-smr-in-sas

About the author

Bailey DeBarmore is a doctoral student at the University of North Carolina at Chapel Hill studying epidemiology. Find her on Twitter @BaileyDeBarmore and blogging for the American Heart Association on the Early Career Voice blog. ​​​