Post hoc analysis for gamlss model in R

Question

I used a gamlss model for my data and would like to do some post hoc analysis afterwards. I tried the packages emmeans and ggemmeans but both of them give me an error: "Error in match.arg(type) : 'arg' should be one of “link”, “response”, “terms”"

Here is the code:

library(gamlss)
model1 <- gamlss(y ~ x1 + x2 + x3, data=na.omit(Dataset), family=ZAGA)

library(ggemmeans)
ggemmeans(model1, terms = "x1")

results in

Can't compute marginal effects, 'emmeans::emmeans()' returned an error.
Reason: 'arg' should be one of “link”, “response”, “terms”
You may try 'ggpredict()' or 'ggeffect()'.

and

library(emmeans)
emmeans(model1, "x1")

results in

[15] ERROR:
'arg' should be one of “link”, “response”, “terms”

From this question I get that this may have nothing to do with the packages but with the way my data is organized. However the solution to the question in the linked post (using ggemmeans instead of emmeans) does not work for me.

Can anyone help me to understand how I should structure my data in order to do a post hoc analysis for a gamlss model in R? Any other suggestion on how to do a posthoc test on this model is very welcome.

Here is my data:

x1,x2,x3,y
a,K1,6696,0
a,K6,274,0
a,K3,10605,0
a,K5,90,0
a,K2,4986,0
a,K4,4999,0
a,K6,878,0
a,K3,3426,0
a,K3,718,0
a,K4,6445,0
a,K2,3778,0
a,K5,7601,0
a,K2,3041,0
a,K4,6808,0
a,K4,3601,0
a,K2,2646,0
a,K4,3938,0
a,K2,5569,0
a,K4,2741,0
a,K5,5886,0
a,K2,6983,0
a,K2,3200,0
a,K3,2643,0
a,K3,737,0
a,K3,11406,0
a,K3,9340,0
a,K1,8322,0
a,K2,714,0
a,K5,7168,0
a,K1,4469,0
a,K4,4381,0
a,K1,9943,0
a,K5,8066,0
a,K3,9696,0
a,K2,3058,0
a,K6,1038,0
a,K6,1014,0
a,K1,1866,0
a,K6,7696,0
a,K5,2056,0
a,K2,7339,0
a,K2,7339,0
a,K6,6860,0
a,K4,4720,0
a,K2,4726,0
a,K5,2700,0
a,K5,3587,0
a,K5,4880,0
a,K3,4190,0
a,K5,9836,0
a,K5,7728,0
a,K3,4880,0
a,K6,4080,0
a,K1,774,0
a,K6,671,0
a,K6,4200,0
a,K2,3221,0
a,K6,3836,0
a,K4,5919,0
a,K4,5006,0
a,K3,4254,0
a,K6,4400,0
a,K4,9829,0
a,K3,3162,0
a,K2,2410,0
a,K5,44946,0
a,K5,3662,0
a,K1,5124,0
a,K2,5348,0
a,K6,15103,0
a,K6,9783,0.04
a,K4,2439,0.05
a,K5,2068,0.05
a,K3,3015,0.05
a,K5,2885,0.05
a,K6,19781,0.05
a,K3,2164,0.054298642533937
a,K6,10865,0.054545454545455
a,K3,11297,0.066666666666667
a,K6,4080,0.066666666666667
a,K2,1567,0.083333333333333
a,K5,2523,0.1
a,K4,5683,0.1
a,K6,545,0.1
a,K4,5666,0.1
a,K4,1198,0.1
a,K4,45903,0.1
a,K5,20303,0.1
a,K6,6645,0.1
a,K6,200,0.1
a,K6,299,0.1
a,K4,6255,0.1
a,K5,6165,0.1
a,K4,5000,0.1
a,K6,11440,0.1
a,K3,112,0.117647058823529
a,K3,11586,0.142857142857143
a,K6,4367,0.15
a,K6,2349,0.15
a,K5,9818,0.193548387096774
a,K2,1620,0.2
a,K5,15600,0.2
a,K3,5928,0.2
a,K6,547,0.2
a,K3,3584,0.2
a,K4,4545,0.2
a,K3,13904,0.214285714285714
a,K6,5681,0.25
a,K4,23824,0.25
a,K4,2560,0.25
a,K4,9,0.25
a,K2,6970,0.25
a,K3,6607,0.266666666666667
a,K3,1797,0.3
a,K4,831,0.3
a,K3,7532,0.3
a,K2,1695,0.3
a,K2,4482,0.3
a,K2,2953,0.4
a,K3,10053,0.444444444444444
a,K6,22121,0.45
a,K3,7062,0.5
a,K6,8406,0.5
a,K6,18044,0.5
a,K2,6650,0.5
a,K6,7675,0.5
a,K2,3215,0.529100529100529
a,K3,7134,0.533333333333333
a,K3,23513,0.588235294117647
a,K3,4212,0.615384615384615
a,K2,10883,0.666666666666667
a,K1,6412,0.666666666666667
a,K1,1949,0.666666666666667
a,K2,6575,0.666666666666667
a,K2,1068,0.8
a,K4,9921,1
a,K1,10499,1.33333333333333
a,K1,1990,1.33333333333333
a,K2,3253,2
a,K2,22,2.14285714285714
a,K1,13763,3
b,K2,3784,0
b,K1,501,0
b,K3,6624,0
b,K4,474,0
b,K2,837,0
b,K2,1876,0
b,K5,2762,0
b,K4,39269,0
b,K1,1205,0
b,K2,4995,0
b,K2,5299,0
b,K2,304,0
b,K3,293,0
b,K6,7.65075785824725,0
b,K3,3822,0
b,K2,4794,0
b,K2,21065,0
b,K2,5958,0
b,K4,5157,0
b,K6,7544,0
b,K6,4492,0
b,K2,1614,0
b,K4,1062,0
b,K1,431,0
b,K3,575,0
b,K2,1223,0
b,K3,3664,0
b,K6,234,0
b,K2,6437,0
b,K2,6059,0
b,K3,2311,0
b,K3,5279,0
b,K1,4258,0
b,K3,4004,0
b,K6,3939,0
b,K4,4478,0
b,K1,4311,0
b,K6,9054,0
b,K6,1302,0
b,K5,3708,0
b,K3,6435,0
b,K2,1485,0
b,K4,2314,0
b,K6,6026,0
b,K3,3291,0
b,K6,623,0
b,K1,691,0
b,K3,22614,0
b,K1,6922,0
b,K4,4623,0
b,K2,12253,0
b,K4,304,0
b,K3,9245,0
b,K2,35,0
b,K6,160,0
b,K2,6163,0
b,K2,6040,0
b,K2,279,0
b,K3,5425,0
b,K1,7036,0
b,K1,10872,0
b,K3,34,0.025
b,K6,4018,0.045454545454546
b,K3,6601,0.049180327868853
b,K5,831,0.05
b,K5,2175,0.05
b,K5,10854,0.05
b,K4,5016,0.05
b,K4,1911,0.05
b,K3,7444,0.0625
b,K4,1995,0.085714285714286
b,K6,240,0.092307692307692
b,K4,5127,0.1
b,K6,4489,0.1
b,K6,2615,0.1
b,K3,6263,0.111111111111111
b,K3,17412,0.111111111111111
b,K6,5573,0.12
b,K3,2198,0.133333333333333
b,K2,8877,0.142857142857143
b,K5,3878,0.15
b,K3,8698,0.15
b,K6,2213,0.15
b,K3,3852,0.157894736842105
b,K5,3917,0.16
b,K6,0,0.162162162162162
b,K3,9006,0.166666666666667
b,K3,0,0.181818181818182
b,K6,1244,0.2
b,K5,7898,0.2
b,K5,2645,0.2
b,K4,27566,0.2
b,K6,11435,0.2
b,K4,34,0.2
b,K2,5668,0.25
b,K2,900,0.285714285714286
b,K3,1586,0.3
b,K6,620,0.3
b,K2,11576,0.333333333333333
b,K6,2315,0.35
b,K2,3076,0.4
b,K5,916,0.4
b,K6,13595,0.4
b,K6,0,0.4
b,K6,7675,0.4
b,K3,2311,0.4
b,K4,9288,0.4
b,K4,2664,0.428571428571429
b,K3,9413,0.470588235294118
b,K3,1637,0.476190476190476
b,K6,14400,0.5
b,K1,1025,0.5
b,K4,3,0.6
b,K6,2467,0.6
b,K5,1359,0.6
b,K5,916,0.6
b,K2,4.3369083067469,0.666666666666667
b,K1,4947,0.666666666666667
b,K1,10735,0.666666666666667
b,K3,2534,0.666666666666667
b,K2,7912,0.8
b,K2,6040,0.857142857142857
b,K6,5681,0.9
b,K1,1751,1
b,K1,0,1
b,K1,5937,1
b,K3,1797,1.16666666666667
b,K2,2661,1.2
b,K4,11826,1.3
b,K2,10229,1.4
b,K1,3124,2
b,K1,5265,2
b,K2,7720,2.22222222222222
b,K2,20245,2.34375
b,K1,3438,3.11111111111111
b,K5,34318,3.3
b,K1,11290,3.33333333333333
b,K1,1227,3.5
b,K1,5335,3.6
b,K1,1819,8
b,K2,19431,8.25

Answer 1

The way the emmeans support currently works is that it extracts the right stuff for the given model, then passes that to mgcv::predict.gam(..., type = "lpmatrix"). That is done in hopes of supporting the possibility that the model involves smoothers.

Theoretically, this should work because:

> class(model1)
[1] "gamlss" "gam"    "glm"    "lm"

This suggests that model1 inherits from mgcv::gam. However, it apparently doesn't inherit from gam. Actually, it doesn't have all the right slots for lm or glm either, but it is easier to doctor-up. I modified the code so it passes off to the glm support instead, and it now works for this example:

> emmeans(model1, "x1")
 x1 emmean     SE  df asymp.LCL asymp.UCL
 a   -1.09 0.1031 Inf     -1.30    -0.892
 b   -0.58 0.0917 Inf     -0.76    -0.400

Results are averaged over the levels of: x2 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

As it stands now, it will not work if there are smoothers in the model, but it will work more reliably for other models. And it works as well for alternative modes besides the dafault "mu" mode. This model is not a great demo for that because we only have an intercept for components other than mu; but we can test that at least:

> emmeans(model1, "1", what="sigma")
 1       emmean     SE  df asymp.LCL asymp.UCL
 overall -0.207 0.0524 Inf     -0.31    -0.104

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

> emmeans(model1, "1", what = "sigma", type = "response")
 1       response     SE  df asymp.LCL asymp.UCL
 overall    0.813 0.0426 Inf     0.734     0.901

Confidence level used: 0.95 
Intervals are back-transformed from the log scale

I have uploaded the revised package to the emmeans GitHub repository