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
 A: 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

#survival adjusted to group effect
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)

Note: Comments removed and placed in Answer.
A: 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.
