Random Effects

Description

Random effects can be included in gamlss2 models in two ways.
The first uses s() with bs = “re” for simple random effects, i.e., when a single factor is entered into the model as a smoother. This approach relies on the s() function from the mgcv package. For example, if area is a factor with several levels, then s(area, bs = “re”) shrinks the level-specific effects of area towards their overall mean.

The second, more general approach uses the model term constructor re(), which provides an interface to the specialised random-effects functionality in the nlme package. This allows fitting more complex random-effects structures.

This documentation focuses on the re() function, but we also provide examples using s(…, bs = “re”).

Usage

re(random, correlation = NULL, ...)

Arguments

random	A formula specifying the random effect part of the model, as in the lme function.
correlation	An optional correlation object, see lme.
…	For the re() function, the dots argument is used to specify additional control arguments for the lme function, such as the method and correlation arguments.

Details

Both functions set up model terms that can be estimated using a backfitting algorithm, e.g., the default RS algorithm.

Value

Function s with bs = “re” returns a smooth specification object of class “re.smooth.spec”, see also smooth.construct.re.smooth.spec.

The re() function returns a special model term specification object, see specials for details.

See Also

gamlss2, smooth.construct.re.smooth.spec, s, lme

Examples

library("gamlss2")



## orthdontic measurement data
data("Orthodont", package = "nlme")

## model using lme()
m <- lme(distance ~ I(age-11), data = Orthodont,
  random =~ I(age-11) | Subject, method = "ML")

## using re(), function I() is not supported,
## please transform all variables in advance
Orthodont$age11  <- Orthodont$age - 11

## estimation using the re() constructor,
## setup formula first
f <- distance ~ age11 + re(~age11|Subject)

## estimate model
b <- gamlss2(f, data = Orthodont)

GAMLSS-RS iteration  1: Global Deviance = 326.6623 eps = 0.392353     
GAMLSS-RS iteration  2: Global Deviance = 326.6623 eps = 0.000000     

## compare fitted values
plot(fitted(b, model = "mu"), fitted(m))
abline(0, 1, col = 4)

## extract summary for re() model term
st <- specials(b, model = "mu", term = "age11", elements = "model")
summary(st)

Linear mixed-effects model fit by maximum likelihood
  Data: structure(list(age11 = c(-3, -1, 1, 3, -3, -1, 1, 3, -3, -1,  1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3,  -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1,  1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3,  -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1,  1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3, -3, -1, 1, 3,  -3, -1, 1, 3, -3, -1, 1, 3), Subject = structure(c(15L, 15L,  15L, 15L, 3L, 3L, 3L, 3L, 7L, 7L, 7L, 7L, 14L, 14L, 14L, 14L,  2L, 2L, 2L, 2L, 13L, 13L, 13L, 13L, 5L, 5L, 5L, 5L, 6L, 6L, 6L,  6L, 11L, 11L, 11L, 11L, 16L, 16L, 16L, 16L, 4L, 4L, 4L, 4L, 8L,  8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 12L, 12L, 12L,  12L, 1L, 1L, 1L, 1L, 20L, 20L, 20L, 20L, 23L, 23L, 23L, 23L,  25L, 25L, 25L, 25L, 26L, 26L, 26L, 26L, 21L, 21L, 21L, 21L, 19L,  19L, 19L, 19L, 22L, 22L, 22L, 22L, 24L, 24L, 24L, 24L, 18L, 18L,  18L, 18L, 17L, 17L, 17L, 17L, 27L, 27L, 27L, 27L), levels = c("M16",  "M05", "M02", "M11", "M07", "M08", "M03", "M12", "M13", "M14",  "M09", "M15", "M06", "M04", "M01", "M10", "F10", "F09", "F06",  "F01", "F05", "F07", "F02", "F08", "F03", "F04", "F11"), class = c("ordered",  "factor")), response_z = c(3.95740740740771, 1.63703703703704,  4.31666666666667, 4.9962962962963, -0.542592592592602, -0.862962962962968,  -1.68333333333333, 0.496296296296297, 0.957407407407398, -0.862962962962968,  -0.683333333333334, 1.4962962962963, 3.4574074074074, 4.13703703703703,  1.81666666666666, 0.996296296296297, -2.0425925925926, 0.137037037037032,  -2.18333333333333, -0.00370370370369955, 2.4574074074074, 2.13703703703704,  2.31666666666666, 2.4962962962963, -0.0425925925926016, -1.36296296296297,  -0.183333333333334, 0.496296296296297, 1.9574074074074, -1.86296296296297,  -0.183333333333334, -0.5037037037037, 0.957407407407398, -2.86296296296297,  6.31666666666666, -0.00370370370369955, 5.45740740740739, 4.63703703703703,  6.31666666666666, 5.49629629629629, 0.957407407407398, -0.362962962962968,  -1.18333333333333, -1.0037037037037, -0.542592592592602, 0.137037037037032,  -0.683333333333334, 1.9962962962963, -5.0425925925926, 1.13703703703703,  1.31666666666667, 3.4962962962963, 0.457407407407402, 2.13703703703704,  0.81666666666667, -0.00370370370369955, 0.957407407407402, 1.13703703703703,  1.31666666666667, 3.9962962962963, -0.0425925925926016, -1.86296296296297,  -1.18333333333333, -1.0037037037037, -1.0425925925926, -3.36296296296297,  -3.18333333333333, -3.0037037037037, -1.0425925925926, -1.86296296296297,  -0.683333333333334, -0.5037037037037, -1.5425925925926, 0.637037037037032,  -0.183333333333334, -0.00370370370369955, 1.4574074074074, 1.13703703703703,  0.316666666666666, 0.496296296296297, -0.542592592592602, -0.362962962962968,  -2.18333333333333, -2.5037037037037, -2.0425925925926, -2.36296296296297,  -3.68333333333334, -3.5037037037037, -0.542592592592602, -0.862962962962968,  -1.68333333333333, -1.0037037037037, 0.957407407407398, -0.362962962962968,  -1.18333333333333, -2.0037037037037, -2.04259259259261, -2.36296296296297,  -2.68333333333333, -4.5037037037037, -5.5425925925926, -4.36296296296297,  -5.68333333333334, -6.5037037037037, 2.4574074074074, 1.63703703703704,  3.31666666666666, 1.9962962962963), weights_w = c(0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369, 0.116596837376369,  0.116596837376369, 0.116596837376369, 0.116596837376369)), terms = ~age11 +      Subject, row.names = c(NA, 108L), class = "data.frame") 
       AIC      BIC    logLik
  449.2116 462.6223 -219.6058

Random effects:
 Formula: ~age11 | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 2.0906352 (Intr)
age11       0.2149235 0.521 
Residual    3.8365552       

Variance function:
 Structure: fixed weights
 Formula: ~weights_w 
Fixed effects:  response_z ~ 1 
                   Value Std.Error DF      t-value p-value
(Intercept) 6.064444e-16 0.4048917 81 1.497794e-15       1

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.305976369 -0.487428882  0.007598099  0.482237876  3.922787577 

Number of Observations: 108
Number of Groups: 27 

## random intercepts and slopes with s() using AIC
a <- gamlss2(distance ~ s(age,k=3) + s(Subject, bs = "re") + s(Subject, age11, bs = "re"),
  data = Orthodont)

GAMLSS-RS iteration  1: Global Deviance = 427.4959 eps = 0.204786     
GAMLSS-RS iteration  2: Global Deviance = 370.4092 eps = 0.133537     
GAMLSS-RS iteration  3: Global Deviance = 340.4089 eps = 0.080992     
GAMLSS-RS iteration  4: Global Deviance = 325.4676 eps = 0.043892     
GAMLSS-RS iteration  5: Global Deviance = 319.3778 eps = 0.018710     
GAMLSS-RS iteration  6: Global Deviance = 317.1286 eps = 0.007042     
GAMLSS-RS iteration  7: Global Deviance = 316.3204 eps = 0.002548     
GAMLSS-RS iteration  8: Global Deviance = 316.0276 eps = 0.000925     
GAMLSS-RS iteration  9: Global Deviance = 315.913 eps = 0.000362     
GAMLSS-RS iteration 10: Global Deviance = 315.865 eps = 0.000151     
GAMLSS-RS iteration 11: Global Deviance = 315.8413 eps = 0.000075     
GAMLSS-RS iteration 12: Global Deviance = 315.828 eps = 0.000042     
GAMLSS-RS iteration 13: Global Deviance = 315.8178 eps = 0.000032     
GAMLSS-RS iteration 14: Global Deviance = 315.8112 eps = 0.000020     
GAMLSS-RS iteration 15: Global Deviance = 315.8054 eps = 0.000018     
GAMLSS-RS iteration 16: Global Deviance = 315.8023 eps = 0.000009     

## compare fitted values
plot(fitted(b, model = "mu"), fitted(m))
points(fitted(a, model = "mu"), fitted(m), col = 2)
abline(0, 1, col = 4)

## more complicated correlation structures.
data("Ovary", package = "nlme")

## ARMA model
m <- lme(follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time), data = Ovary, 
  random = pdDiag(~sin(2*pi*Time)), correlation = corARMA(q = 2))

## now with gamlss2(), transform in advance
Ovary$sin1 <- sin(2 * pi * Ovary$Time)
Ovary$cos1 <- cos(2 * pi * Ovary$Time)

## model formula
f <- follicles ~ sin1 + cos1 +
  re(~ 1 | Mare) +
  re(~ sin1 - 1 | Mare, correlation = corARMA(q = 2))

## estimate model
b <- gamlss2(f, data = Ovary)

GAMLSS-RS iteration  1: Global Deviance = 1561.4397 eps = 0.164702     
GAMLSS-RS iteration  2: Global Deviance = 1561.3891 eps = 0.000032     
GAMLSS-RS iteration  3: Global Deviance = 1561.3816 eps = 0.000004     

## smooth random effects
f <- follicles ~ ti(Time) + ti(Mare, bs = "re") + 
  ti(Mare, Time, bs = c("re", "cr"), k = c(11, 5))

g <- gamlss2(f, data = Ovary)

GAMLSS-RS iteration  1: Global Deviance = 1522.1073 eps = 0.185743     
GAMLSS-RS iteration  2: Global Deviance = 1436.5015 eps = 0.056241     
GAMLSS-RS iteration  3: Global Deviance = 1426.6605 eps = 0.006850     
GAMLSS-RS iteration  4: Global Deviance = 1425.7658 eps = 0.000627     
GAMLSS-RS iteration  5: Global Deviance = 1425.6839 eps = 0.000057     
GAMLSS-RS iteration  6: Global Deviance = 1425.6744 eps = 0.000006     

## compare fitted values
par(mfrow = n2mfrow(nlevels(Ovary$Mare)), mar = c(4, 4, 1, 1))

for(j in levels(Ovary$Mare)) {
  ds <- subset(Ovary, Mare == j)

  plot(follicles ~ Time, data = ds)

  f <- fitted(b, model = "mu")[Ovary$Mare == j]
  lines(f ~ ds$Time, col = 4, lwd = 2)

  f <- fitted(g, model = "mu")[Ovary$Mare == j]
  lines(f ~ ds$Time, col = 3, lwd = 2)

  f <- fitted(m)[Ovary$Mare == j]
  lines(f ~ ds$Time, col = 2, lwd = 2)
}

## simulated data
set.seed(1328)

n <- 10000
k <- 500

## generate random effects
f <- as.factor(sample(1:k, size = n, replace = TRUE))
ref <- rnorm(k, sd = 0.6)

## random effects only for sigma
y <- rBCT(n, mu = 10, sigma = exp(-1 + ref[f]), tau = 10)

## estimate model
b <- gamlss2(y ~ 1 | re(~ 1 | f), family = BCT)

GAMLSS-RS iteration  1: Global Deviance = 54284.8778 eps = 0.443825     
GAMLSS-RS iteration  2: Global Deviance = 53798.1163 eps = 0.008966     
GAMLSS-RS iteration  3: Global Deviance = 53767.8702 eps = 0.000562     
GAMLSS-RS iteration  4: Global Deviance = 53761.3751 eps = 0.000120     
GAMLSS-RS iteration  5: Global Deviance = 53758.2934 eps = 0.000057     
GAMLSS-RS iteration  6: Global Deviance = 53757.5692 eps = 0.000013     
GAMLSS-RS iteration  7: Global Deviance = 53756.9982 eps = 0.000010     
GAMLSS-RS iteration  8: Global Deviance = 53756.5246 eps = 0.000008     

## extract fitted random effects
cb <- coef(b, full = TRUE)
cb <- cb[grepl("re", names(cb), fixed = TRUE)]

## compare
plot(cb, ref, main = round(mean((cb - ref)^2), 2))
abline(0, 1, lwd = 2, col = 4)