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.4591 eps = 0.080856     
GAMLSS-RS iteration  4: Global Deviance = 325.4834 eps = 0.043986     
GAMLSS-RS iteration  5: Global Deviance = 319.1487 eps = 0.019462     
GAMLSS-RS iteration  6: Global Deviance = 316.6824 eps = 0.007727     
GAMLSS-RS iteration  7: Global Deviance = 315.7553 eps = 0.002927     
GAMLSS-RS iteration  8: Global Deviance = 315.4122 eps = 0.001086     
GAMLSS-RS iteration  9: Global Deviance = 315.2857 eps = 0.000401     
GAMLSS-RS iteration 10: Global Deviance = 315.2397 eps = 0.000145     
GAMLSS-RS iteration 11: Global Deviance = 315.2211 eps = 0.000059     
GAMLSS-RS iteration 12: Global Deviance = 315.2141 eps = 0.000022     
GAMLSS-RS iteration 13: Global Deviance = 315.2112 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)