Can multiple outcomes be regressed on the same linear model?

[nicklausmillican]

From your book, you mention that a multinomial can be decomposed into several Poisson regression, where each output category gets its own linear model. Say that, instead of 2-3 category outcomes, there are many more (e.g., K=100). Is there a way to regress K-1 outcome categories on the same linear model, i.e., without having to write K-1 linear models?

I'm imagining something like:

Y <- rep(c("A","B","C"), times=c(10,8,6))
Y_ind <- 1:3
x <- c(rnorm(n=10, mean=0, sd=1),
rnorm(n=8, mean=2, sd=1),
rnorm(n=6, mean=5, sd=1))
one_hot <- model.matrix(~ Y - 1)
Y_counts <- colSums(one_hot)

m <- ulam(
alist(
Y_counts[Y_ind] <- dpois(lambda[Y_ind]),
log(lambda[Y_ind]) <- a[Y_ind] + b[Y_ind]*x,
a[Y_ind] ~ dnorm(0,1),
b[Y_ind] ~ dnorm(0,1)
), data=list(Y_counts=Y_counts, Y_ind=Y_ind)
)