get_coef.gamlss bug

[vincentarelbundock]

Reported here first:

strengejacke/ggeffects#297 (comment)

[vincentarelbundock]

remotes::install_github("vincentarelbundock/marginaleffects")

Restart your R session. Then:

library(gamlss)
library(marginaleffects)
dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv", na.strings = "")
dat <- dat |>
    transform(prop = rBE(nrow(dat), mu = 0.5, sigma = 0.2)) |> 
    na.omit()
mod <- gamlss::gamlss(
        prop ~ sex * body_mass_g + year + re(random = list(~ 1 | species, ~ 1 | island)),
        family = BE(),
        data = dat,
        trace = FALSE)
avg_comparisons(mod, what = "mu")
# Warning in vcov.gamlss(model, what = dots$what): Additive terms exists in the  mu formula. 
#   Standard errors for the linear terms maybe are not appropriate
# 
#         Term      Contrast  Estimate Std. Error       z Pr(>|z|)   S     2.5 %    97.5 %
#  body_mass_g +1             6.34e-06   5.39e-06    1.18    0.240 2.1 -4.23e-06  1.69e-05
#  sex         male - female -1.92e-02   3.56e-02   -0.54    0.589 0.8 -8.90e-02  5.06e-02
#  year        +1            -1.29e-02   3.69e-05 -348.60   <0.001 Inf -1.29e-02 -1.28e-02
# 
# Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
# Type:  response

[vincentarelbundock]

Here strengejacke/ggeffects#297 (comment) @DHLocke asks:

What about a what = 'all' argument that combines the mu and sigma?

marginaleffects only supports values that are supported by the underlying predict() function. If that function does not support all, you will have to call marginaleffects::avg_comparisons() several times. The good news is that avg_comparisons() returns simple data frames, and they are easy to rbind() and manipulate.

[DHLocke]

Amazing, thanks so much. sorry for a novice question, what is the appropriate way to combined those then?

plot_predictions(mod, what = 'mu', condition = c('body_mass_g', 'sex')) # excellent!
plot_predictions(mod, what = 'mu', condition = c('sex', 'body_mass_g')) # fabulous

great, but

Warning message:

In vcov.gamlss(model, what = dots$what) :

Additive terms exists in the mu formula.

Standard errors for the linear terms maybe are not appropriate

try to add predictions together

preds_mu <- predictions(mod, what = 'mu')
preds_sigma <- predictions(mod, what = 'sigma')

library(tidyverse)
preds_mu |>
bind_cols(
preds_sigma |>
select(estimate_sigma = estimate, conf.low.sigma = conf.low, conf.high.sigma = conf.high)
) |>
tibble() |>
mutate(pred = estimate + estimate_sigma
, ll = conf.low - conf.low.sigma
, ul = conf.high + conf.high.sigma
) |>
ggplot(aes(pred, body_mass_g, color = sex)) +
geom_point()

??

[vincentarelbundock]

I'm not a specialist of these models, so I can't tell you if it makes sense to just add them up like this. But if you are just asking how to do it mechanically, then what you are doing seems fine. You probably want to check the rowid and term columns to make sure they match (a proper join is always safer than binding columns).