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:
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)
surv <- survfit(Surv(days, status)~1)
#survival by group
kmsurv <- survfit(Surv(days,status) ~ strata(group), df)
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)
multinom <- multinom(group ~ demo1 + demo2 + demo3, df)
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)
#Cox model using multinomial propensity probabilities.
cox2 <- coxph(Surv(days,status) ~ group2x.ind + group3x.ind + X1 + X2, df)
Note: Comments removed and placed in Answer.