Yes, with a Poisson GLM (log linear model) you can fit multinomial models. Hence multinomial logistic or log linear Poisson models are equivalent.
You need to see random counts $y_{ij}$ as Poisson random variables with means $μ_{ij}$ and specify the following the following log-linear model
$\log(μ_{ij}) = o + p_i + c_j + x_iβ_j$
To get a multinomial logit model the parameters are:
A parameter $p_i$ for each multinomial observation, for example individuals or group. This assures exact reproduction of the multinomial denominators and actually establishes the equivalence of Poisson and multinomial model. They are fixed in the multinomial likelihood, but random in the Poisson likelihood.
A parameter $c_j$ for each response category. This way the counts can be different for each response category and the margins can be non-uniform.
What you are really interested in are the interaction terms $x_iβ_j$ that represent the effects of $x_i$ on the log-odds of response $j$.
The log-odds can be simply calculated by $\log(μ_{ij}/μ_{ik}) = (c_j-c_k) +x_i(β_j-β_k)$. It is the log odds that observation i will fall in response category j relative to the response category $k$.
Then, the parameters in the multinomial logit model (denoted in latin letters) can be obtained as differences between the parameters in the corresponding log-linear model, i.e. $a_j = α_j-α_k$ and $b_j = β_j-β_k$.
The emmeans package now provides estimated marginal means for GEE models:
library(geepack)
library(emmeans)
warp.gee <- geeglm(breaks ~ tension, id=wool, family=gaussian, data=warpbreaks)
emmeans(warp.gee, ~tension)
tension lsmean SE df asymp.LCL asymp.UCL
L 36.38889 5.774705 Inf 25.07067 47.70710
M 26.38889 1.689200 Inf 23.07812 29.69966
H 21.66667 2.042753 Inf 17.66294 25.67039
plot(emmeans(warp.gee, ~tension), horizontal=FALSE, ylab="Estimated mean")
Best Answer
Try using the R package glmnet. It has very good documentation, it is the state of the art and certainly does multinomial regression. I have noticed some of these packages expect you to figure out cross validation yourself. In that sense, glmnet is nice. The cv function does it for you.
Besides I don't know what algorithm GEE uses (I haven't looked into it, they probably do the newer methods as is). Some older packages resort to Interior Point Methods which was state of the art at the time but has since become outdated.