how to get multinomial regression coefficients from gamlss

Question

I am trying to fit a mixed effects multinomial model for my 4 level categorical outcome variable (y). I am using the gamlss function from gamlss package. The call looks like:

library(mixcat)
data("schizo")
head(schizo)

gamlss::gamlss(
formula = y ~ wk + gamlss::random(factor(trt)),
sigma.formula = ~1,
family = gamlss.dist::MN4,
data = schizo
)

The output I get is like this.

Family:  c("MN4", "Multinomial 4 levels") 
Fitting method: RS() 

Call:  gamlss::gamlss(formula = y ~ wk + gamlss::random(factor(trt)),
   sigma.formula = ~1, family = gamlss.dist::MN4,  
    data = schizo) 

Mu Coefficients:
                (Intercept)                           wk  gamlss::random(factor(trt))  
                    -2.5082                       0.4665                           NA  
Sigma Coefficients:
(Intercept)  
    -0.1063  
Nu Coefficients:
(Intercept)  
    -0.2464  

 Degrees of Freedom for the fit: 4 Residual Deg. of Freedom   1599 
Global Deviance:     4077.31 
            AIC:     4085.31 
            SBC:     4106.82 

I am confused about this output from the gamlssfunction. I was expecting 3 sets of regression coefficients one for each level of y except the baseline, k-1 sets of regression coefficients. Instead I get one set of coefficients , shown above.

Does the gamlss function generate k-1 sets of regression coefficients? If so, how do I extract these coefficients?

Answer 1

MN4 has 4 levels and hence 3 distribution parameters, mu, sigma and nu. (See Rigby et al. (2019) page 528-529.)

You would need to put the model

wk + gamlss::random(factor(trt))

into all 3 parameters.

Some thought is needed to interpret this model.