2
$\begingroup$

https://i.stack.imgur.com/fDm4x.png

I want to analyse the association between the outcome "Other CTR-CVD" and the independent variables would be "anthracyclines", "Her2", "VEGF", "TKI, "Prot Inh", RAF/MEK Inh, "ICI" and Fluoropirimidines"

(this are the 8 first rows of the table I have with the date of 296 patients)

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ Is your data structured so that it is one row per person or one row per condition, i.e. of the 5 rows of tachyarrythmia, one of those cases may also belong to syncope? Also do you have recurrence of specific conditions, i.e. does one row represent at least one occurrence of "X", or if recurrence occurs, are there two or three or more rows as needed? Also are there any negative controls without any heart failure conditions? $\endgroup$
    – AdamO
    May 12 at 16:05
  • $\begingroup$ Do you mean that the outcome is one of three (or more) categories rather than one of two categories? $\endgroup$
    – Dave
    May 12 at 16:21
  • $\begingroup$ Thanks AdamO, the data is structured so that it is one row per person. Each patient had one of the "Other CTR-CVD" events and was exposed to either one or more of the drugs listed there (each drug is one column). $\endgroup$
    – Maria Sol
    May 13 at 10:57
  • $\begingroup$ "Do you mean that the outcome is one of three (or more) categories rather than one of two categories?" – Yes, it is one of 9 actually $\endgroup$
    – Maria Sol
    May 13 at 10:58
  • $\begingroup$ if you are interested in pure categorization and don't need things like odds ratios, consider discriminant analysis, although that has assumptions as well $\endgroup$
    – Ralph Winters
    May 13 at 13:55

2 Answers 2

Reset to default
2
$\begingroup$

When you want to do something like logistic regression but with $3+$ outcome categories instead of two, the $y$ is multinomial instead of binomial. Consequently, the analogous model is multinomial logistic regression, sometimes called polytomous logistic regression or softmax regression.

In R, this can be performed by using multinom in the nnet package.

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ I do wonder, however, whether the assumption of "independence of irrelevant alternatives" would hold for this type of outcome data. I haven't thought that through, but the OP might need to address it. (+1) $\endgroup$
    – EdM
    May 13 at 12:30
  • $\begingroup$ @EdM Given that the snippet of data suggests the data might be mostly zeros, I'm wondering whether the idea of fitting a regression model is reasonable to begin with. $\endgroup$
    – dipetkov
    May 13 at 17:12
1
$\begingroup$

If each patient was exposed to either one or more of the drugs listed, categories are not mutually exclusive. This can be coded as a multinomial outcome, as suggested in an earlier answer. Depending on the number of drug combination patterns observed, however, this might yield a large number of categories ($2^9$ at worst). A multinomial regression will estimate a coefficient for each predictor and category (plus an intercept), yielding at worst $2^9 \times (p+1)$ (where $p$ = number of predictor variables) coefficients being estimated.

One can also model a multivariate binary outcome; this limits the number of regression coefficients to $9$ (number of drug types) $\times p$. Can either be analysed with data in wide form using a structural equation model (e.g., in R using package lavaan). Or with data in long form, using a multilevel model (e.g., in R using package lme4).

An example (but different number of predictors and binary outcomes):

## generate 3 predictor variables
set.seed(42)
x1 <- rnorm(1000) 
x2 <- rnorm(1000)
x3 <- rnorm(1000)

## generate 5 binary outcomes
y1 <- rbinom(1000, 1, prob = 1 / (1+exp(-(x1))) )
y2 <- rbinom(1000, 1, prob = 1 / (1+exp(-(x2))) )
y3 <- rbinom(1000, 1, prob = 1 / (1+exp(-(x3))) )
y4 <- rbinom(1000, 1, prob = 1 / (1+exp(-(.5*(x1+x2)))) )
y5 <- rbinom(1000, 1, prob = 1 / (1+exp(-(.5*(x2+x3)))) )

table(y1, y2) ## not mutually exclusive

## Create multinomial response
y <- factor(paste0(y1, y2, y3, y4, y5))
length(unique(y)) ## 2^5 = 32 categories
data <- data.frame(x1, x2, x3, y1, y2, y3, y4, y5, y)

## Fit multinomial (penalized) regression
library("glmnet")
pglm <- glmnet(x = data[,1:3], y = data$y, lambda = .01, family = "multinomial")
t(do.call("cbind", coef(pglm))) ## 32x4 = 128 estimated coefficients

## Fit wide format multivariate probit regression
sem_mod <- '
  y1 ~ x1 + x2 + x3
  y2 ~ x1 + x2 + x3
  y3 ~ x1 + x2 + x3
  y4 ~ x1 + x2 + x3
  y5 ~ x1 + x2 + x3
' 
sem_fit <- lavaan(sem_mod, data = data, ordered = c(paste0("y", 1:5)))
coef(sem_fit) ## 3x5 = 15 estimated coefficients
Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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