Panel data and multinomial logistic regression in R

I have panel data of varying length for some 400+ companies. Each company is identified by a code.

date         code    var1   var2  var3    category
2016-01-01   AAA        1     2      3    2
2016-02-01   AAA        2     3      3    3
2016-01-01   BBB        1     2      3    1
2016-02-01   BBB        2     3      3    3


where the category is 1, 2 or 3

I want to do a regression to see which variables affect the Category.

So far I digged it down that I need to use pglm function from the pglm package which is the Panel Estimators for Generalized Linear Models.

Information on pglm is quite limited and is basically just examples.

Though I wanted to run the model with fixed effects (at the date and company level), it does not seem to allow this (Why pglm fails for within model?)

So I ran it with "pooling"

formula_lm1 <- category ~ var1 +var2 +var3
f_pglm <- pglm(
formula_lm1,
data = test,
family = ordinal('logit'),
model = "pooling",
index = c('code', 'date'),
print.level = 0,
method = 'nr'
)
summary(f_pglm)


Questions:

1. Am I using the correct approach - correct function for my task?

2. Does anyone know if using the index = c('code', 'date') allows to have fixed effects?

3. With my regression, how do I adjust standard errors?

4. Are there any alternative packages/functions to use for my task?