4
$\begingroup$

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

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Can you confirm that the propensity score will be used to balance treatment assignment and that you have more than two treatments? $\endgroup$
    – Todd D
    Aug 19, 2016 at 2:49
  • $\begingroup$ Yes That is right! $\endgroup$
    – Vincent
    Aug 19, 2016 at 22:36
  • $\begingroup$ 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! $\endgroup$
    – Vincent
    Aug 20, 2016 at 7:51
  • $\begingroup$ 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? $\endgroup$
    – Vincent
    Aug 31, 2016 at 10:24
  • $\begingroup$ 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. $\endgroup$
    – Todd D
    Sep 1, 2016 at 0:19

2 Answers 2

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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:

Lechner, Michael, Identification and Estimation of Causal Effects of Multiple Treatments Under the Conditional Independence Assumption (September 1999). Available at SSRN: https://ssrn.com/abstract=177089 or http://dx.doi.org/10.2139/ssrn.177089

Improve this answer
$\endgroup$
2
  • $\begingroup$ which paper? Link is broken $\endgroup$
    – Papayapap
    Oct 19, 2023 at 0:32
  • 1
    $\begingroup$ @Papayapap Fixed the link and added the full citation. $\endgroup$
    – dimitriy
    Oct 19, 2023 at 4:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.