GAMLSS Random Effect Coefficients How do I extract the coefficients of my random effects in a Gamlss model?
Let's take a simple example of a sample of individuals with intercepts which are normally distributed. Additional normally distributed error is added on top of this.
library(data.table)
library(nlme)
library(gamlss)

my.mean <- 4 # I wish to retrieve this value
my.sd <- 3 # I wish to retrieve this value
my.sd.id <- 2 # I wish to retrieve this value

set.seed(1234)
n <- 4 * 10^3
t <- 2

my.data <- CJ(id = seq(n), t = seq(t))
my.data[, mean_id := rnorm(1, my.mean, my.sd), by = id]
my.data[, value := mean_id + rnorm(t, 0, my.sd.id), by = id]

my.data[t == 1, sd(mean_id)] # 2.988226 which is ~my.sd
my.data[t == 1, sd(value)] # 3.57642 which is ~sqrt(my.sd^2 + my.sd.id^2)

l <- lme(
  fixed = value ~ 1,
  random = ~ 1|id,
  data = my.data
)
l
# Fixed(Intercept): 4.033439 (~my.mean)
# Random (Intercept): 2.974774 (~my.sd)
# Residual: 1.987995 (~my.sd.id)

m <- gamlss(
  formula = value ~ 1 + random(factor(id), df = NULL),
  sigma.formula = ~ 1,
  family = NO(),
  data = my.data
)
m
# Mu Coefficients (Intercept): 4.033 (~my.mean)
# Sigma Coefficients (Intercept): 0.4244 (exp(0.4244) = 1.528673 ??)

Although the mu coefficient from the Gamlss model corresponds to the expected value, I cannot find information about the random intercepts nor the residual error term.
The Gamlss Sigma Coefficients (Intercept) also seems to depend on the number of values (t) per id, which I think affects the internal optimising of lambda and the df. It converges to log(my.sd.id) for large t, whereas if I set df = 1, the value is approximately sqrt(my.sd^2 + my.sd.id^2).
 A: According to the documentation:

The function random() can be seen as a smoother for use with factors
in gamlss(). It allows the fitted values for a factor predictor to be
shrunk towards the overall mean, where the amount of shrinking depends
either on lambda, or on the equivalent degrees of freedom or on the
estimated sigma parameter (default).

If I understand it correctly, the equivalent of the lme method is one where GAMLSS does not shrink the parameters. We can force this by setting df = 0 or lambda = 10^9. However, I could only find this argument in the random() notation of the random effect and not in the re() notation. Vice versa, I could only find the VarCorr for the re() method.
So one silly solution is to run the model twice:
m1 <- gamlss(
  formula = value ~ 1 + re(random = ~ 1|id),
  sigma.formula = ~ 1,
  family = NO(),
  data = my.data
)

m2 <- gamlss(
  formula = value ~ 1 + random(factor(id), df = 0),
  sigma.formula = ~ 1,
  family = NO(),
  data = my.data
)

Subsequently, the parameter estimates can be extracted from the combination of these two models.
my.mean.estimated <- m1$mu.coefficients[[1]] # ~4
my.sd.estimated <- as.numeric(VarCorr(getSmo(m1))[[3]]) # ~3
my.sd.id.estimated <- sqrt(exp(m2$sigma.coefficients[[1]])^2 - my.sd.estimated^2) # ~2

