#StackBounty: #r #post-hoc #lsmeans #gamlss Post hoc analysis for gamlss model in R

Bounty: 50

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


Get this bounty!!!

#StackBounty: #mixed-model #lme4-nlme #gamlss Standard deviation of outcome in gamlss model with random intercepts in mean

Bounty: 50

I’m interested in a simple random-intercepts model:
$$Y_{ij} = alpha_0 + gamma_i + epsilon_{ij}$$
where $gamma_i sim N(0, sigma_{gamma}^2)$ independently of $epsilon_{ij} sim N(0, sigma_{epsilon}^2)$.

In nlme:

library(nlme)
data(ergoStool)
l1 <- lme(effort ~ 1, data = ergoStool, random = ~ 1|Subject, method="ML")

and gives estimates of (what I assume are) $sigma_gamma$ and $sigma_epsilon$:

summary(l1)

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:   0.9090556 2.020727

In gamlss, the same model can be fit, with the random intercepts included in a couple of ways:

library(gamlss)
t1 <- gamlss(effort ~ 1 + random(Subject), data = ergoStool)
t2 <- gamlss(effort ~ 1 + re(random = ~ 1|Subject), data = ergoStool)

The gamlss models give the same fit, and each has three variance ‘components’, which agree with one another:

c(sigb = getSmo(t1)$sigb, sige = getSmo(t1)$sige, exp(t1$sigma.coefficients))

       sigb        sige (Intercept) 
  0.9090593   1.0610962   1.9043752

summary(t2)

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:   0.9090594 1.061096

exp(t2$sigma.coefficients)

(Intercept) 
   1.904375

The first is obviously the same $hat{sigma}_gamma$, but I can’t work out how the other two (the residual standard deviation from the random intercept part, and the estimate of the standard deviation of the outcome) relate to $hat{sigma}_epsilon$. Is there a relationship between these estimates?


Get this bounty!!!

#StackBounty: #splines #gamlss Create Spline from Coefficients and Knots in GAMLSS

Bounty: 50

In the R package GAMLSS, it is possible to model a random variable $Y$ as a smoothed non-parametric function of some predictor $x$.

One option for such a function is the penalised spline using y~pb(x). This outputs a list of coefficients and knots which, combined with a set of basis splines, results in a smooth function of $x$.

How can one recreate the spline function, given the coefficients and the knots? (and without having to write my own b-spline generating function).

For example:

library(gamlss)
set.seed(9876)
xstart <- 0
xend <- 100
datan <- 20000

seq <- seq(xstart, xend)
mean <- sapply(seq, function(x){0.5+0.2*sin(x/10)})
xs <- ceiling(runif(datan, xstart, xend))
ys <- sapply(xs, function(x){rnorm(1, mean = mean[x], sd = 0.1)})

m1 <- gamlss(ys~pb(xs))
plot(xs, ys)
lines(seq, mean, col="red")
lines(xs[order(xs)], fitted(m1)[order(xs)], col="green")

intercept <- m1$mu.coefficients[1] # 0.5495853
weight <- m1$mu.coefficients[2] # -0.0002851018
coefficients <- c(m1$mu.coefSmo[[1]]$coef) # c(-0.170704842, -0.066451626,  0.026591530,  0.119289203,  0.159657021,  0.149185418,  0.086505094,  0.003904402, -0.100156999, -0.188811997, -0.238717366, -0.237884900, -0.197802945, -0.090559794,  0.012576273,  0.101003289,  0.169210741,  0.181836117,  0.143883546,  0.061281663, -0.038608572, -0.136215586, -0.232871483)
knots <- m1$mu.coefSmo[[1]]$knots # c(-5.039,  0.010,  5.059, 10.108, 15.157, 20.206, 25.255, 30.304, 35.353, 40.402, 45.451, 50.500, 55.549, 60.598, 65.647, 70.696, 75.745, 80.794, 85.843, 90.892, 95.941 100.990 106.039)

enter image description here

How can I obtain the green function, knowing only the intercept, weight, coefficients and knots? I currently plot this function using fitted(m1). However, this is simply a list of $y$ values for the originally inputted list of $x$ values, it is not a function which gives $y$ for any new $x$.


Get this bounty!!!