I am utilizing gamlss currently because of its flexibility in specifying 'rarer' data families but run into an error when trying to compare 2-by-2 differences in mu effects for my fixed model factors.

This is the code:

m3=gamlss(y ~ x1 + x2 + x3 + x1 : x2 + x1 : x3 + x2 : x3,



emmeans(m3, x1)
emmip(m3, x1 ~ x2)

However, both of these commands related to the emmeans package result into the same error, namely:

> emmeans(m3, x1)
Error in V[idx, idx, drop = FALSE] : subscript out of bounds
> emmip(m3, x1 ~ x2) #look at 2-by-2 means with tukey HSD
Error in V[idx, idx, drop = FALSE] : subscript out of bounds

After looking for solutions online, I believe it is likely related to something in the gamlss baseline code, but cannot find where the issue lies exactly. Can anyone help with understanding what I am doing wrong? Solutions proposed in Post hoc analysis for gamlss model in R and https://stackoverflow.com/questions/56253376/emmeans-error-error-in-match-argtype-arg-should-be-one-of-link-respon are not fixing the issue.

(ggemmeans from 'ggeffects' package also returns the same error)

> ggemmeans(m3,"x1")
Can't compute estimated marginal means, 'emmeans::emmeans()' returned an error.

Reason: subscript out of bounds
You may try 'ggpredict()' or 'ggeffect()'.

  • $\begingroup$ Support for gamlss is pretty sketchy. The error you see is usually due to some mismatch in dimensions, e.g. returning some extra unexpected coefficients. It might be possible to cobble together something via qdrg(): by obtaining the correct coef, vcov, etc., examining each carefully that they match up and match the formula. See ? qdrg in emmeans $\endgroup$
    – Russ Lenth
    Feb 10, 2022 at 3:47

1 Answer 1


If anyone would still be looking for an answer, I found the solution (given people in my lab ran into similar issues).

As per usual, it is actually really stupid. Emmeans seems to not be able to read outputs from GAMLSS if your initial dataframe has ordered factors in it or things that were manipulated with dplyr on forehand. However, when you remove all the fancy things from your dataframe ,it runs smoothly and gives you proper estimates.

The fastest way to do that for me was saving it as .csv after manipulation and reading it back in with header = T.

write.csv(dfbasstress, "C:/Users/SReynaert/Desktop/R_analysis/FULLYJOINED.csv", row.names=T) TESTCSV=read.csv("C:/Users/SReynaert/Desktop/R_analysis/FULLYJOINED.csv",header=T)

Hope this is helps!


Not the answer you're looking for? Browse other questions tagged or ask your own question.