How do I do a logistic regression model in R for an outcome with multiple values? 
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)
 A: 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.
A: 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

