# Generating a propensity score for multiple treatment using multinomial logistic regression

I am analyzing data from a representative cohort (>10,000 persons, 10 years follow-up) and I would like to perform a retrospective cohort study comparing the effect of a treatment on the outcomes. Although the population is representative, applying regression models solely might introduce a selection bias because these are historical data and people were assigned to different treatments for different reasons. Certainly controlling my regression models for a full set of covariates should work just fine but still I would like to perform a sensitivity analysis adjusting for a calculated propensity score (regression adjustment rather than matching). However, I would normally calculate the PS with a logit model when the treatment variable is a binary variable, while in this case I have three categories. Would it work just fine to calculate the PS using a multinomial regression model? If not, could you suggest any other way (I use Stata as statistical software)?

Thanks a lot

• Can you confirm that the propensity score will be used to balance treatment assignment and that you have more than two treatments? – Todd D Aug 19 '16 at 2:49
• Yes That is right! – Vincent Aug 19 '16 at 22:36
• Hi Todd, thanks for this! Indeed I came across those two topics as well but I might have overlooked them but none really advises about the best way to calculate the score with multiple treatments! – Vincent Aug 20 '16 at 7:51
• Hi Todd. Thanks again for this! I read with attention the circulation paper, and I found it brilliant. My only question is regarding the use of the three propensity scores in the Cox model. If each score would exclude one of the three groups per time (1 vs 2; 2 vs 3; 1 vs 3), that means that each time one of the groups would have missing value. Wouldn't that be a problem when including the three variables in the Cox model due to missing values? – Vincent Aug 31 '16 at 10:24
• I think you will be OK as each of the 3 logistic models will predict some non-zero probability of 1 vs 2, 2 vs 3, or 3 vs 1. Each cell should possess some tiny chance of membership in one or the other group, which avoids the problem you describe. Essentially, the logistic model should not fully "exclude" a group assignment. If there is not some tiny chance, then your group assignment is too highly correlated with a measured variable to lend itself to propensity analysis. – Todd D Sep 1 '16 at 0:19

I referenced the article in this answer: Propensity Score Matching for more than 2 groups to compile the present answer.

The Circulation article states: "These [propensity] scores were developed from binary logistic regression models and were represented in the final Cox regression models by three variables consisting of the linear score or logit from each of the three logistic models (ie, CABG versus PTCA, CABG versus medical therapy, PTCA versus medical therapy)." They proceeded with a Cox model stratified for treatment.

You could proceed with any number of models by including the propensity as either: 1) 2 probabilities: p(group = 1) and p(group=2) fitted from a multinomial model of treatment assignment, or 2) 3 logit probabilities from each of 3 binary models for each of 1 vs 2, 2 vs 3 and 1 vs 3.

While I cannot provide a quantitative proof, the intuition of this approach seems sound. I compiled the following R' code which, reassuringly, provides the same treatment effects and log likelihood functions for either the binary logistic or multinomial approach:

library("survival")
require("survival")
library("nnet")
require("nnet")

set.seed(42)
days <- rpois(100, 3)
group <- sample(c(1,2,3), 100, replace=TRUE)
status <- rbinom(100,1,0.65)
demo1 <- rnorm(100,100,25)
demo2 <- rpois(100,10)
demo3 <- rbinom(100,1,0.67)

df <- data.frame(days, status, group, demo1, demo2,demo3)

#overall survival
surv <- survfit(Surv(days, status)~1)
summary(surv)
plot(surv)

#survival by group
kmsurv <- survfit(Surv(days,status) ~ strata(group), df)
plot(kmsurv)

#propensity
df12 <- subset(df ,group != 3)
df12$group <- df12$group == 2
model12 <- glm(group ~demo1 + demo2 + demo3, df12, family= "binomial")
df$pred12 <- predict(model12,df ,type="response") df23 <- subset(df ,group != 1) df23$group <- df23$group == 3 model23 <- glm(group ~demo1 + demo2 + demo3, df23, family= "binomial") df$pred23 <- predict(model23,df ,type="response")

df13 <- subset(df ,group != 2)
df13$group <- df13$group == 1
model13 <- glm(group ~demo1 + demo2 + demo3, df13, family= "binomial")
df$pred13 <- predict(model13,df ,type="response") #dummy variable for group df$group2x.ind <- df$group == 2 df$group3x.ind <- df\$group == 3

cox <- coxph(Surv(days,status) ~ group2x.ind + group3x.ind +pred12 +pred23 +pred13, df)
summary(cox)

#multinomial propensity
multinom <- multinom(group ~ demo1 + demo2 + demo3, df)
summary(multinom)
df <- data.frame(fitted(multinom), df)

#Cox model using logit propensity probabilities.
options(digits = 5)
cox1 <- coxph(Surv(days,status) ~ group2x.ind + group3x.ind + pred12 + pred23 + pred13, df)
format(cox1, scientific=F)
summary(cox1)

#Cox model using multinomial propensity probabilities.
cox2 <- coxph(Surv(days,status) ~ group2x.ind + group3x.ind + X1 + X2, df)
summary(cox2)


In Stata, you can try teffects multinomial.
Alternatively, you can also feed propensity scores estimated with multinomial commands to user-written commands like psmatch2` as in this paper.