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 = 359.232 eps = 0.331768     
GAMLSS-RS iteration  2: Global Deviance = 318.4956 eps = 0.113398     
GAMLSS-RS iteration  3: Global Deviance = 315.5943 eps = 0.009109     
GAMLSS-RS iteration  4: Global Deviance = 315.3577 eps = 0.000749     
GAMLSS-RS iteration  5: Global Deviance = 315.3255 eps = 0.000102     
GAMLSS-RS iteration  6: Global Deviance = 315.3171 eps = 0.000026     
GAMLSS-RS iteration  7: Global Deviance = 315.3099 eps = 0.000022     
GAMLSS-RS iteration  8: Global Deviance = 315.3036 eps = 0.000020     
GAMLSS-RS iteration  9: Global Deviance = 315.2977 eps = 0.000018     
GAMLSS-RS iteration 10: Global Deviance = 315.2923 eps = 0.000017     
GAMLSS-RS iteration 11: Global Deviance = 315.2873 eps = 0.000015     
GAMLSS-RS iteration 12: Global Deviance = 315.2826 eps = 0.000014     
GAMLSS-RS iteration 13: Global Deviance = 315.2782 eps = 0.000013     
GAMLSS-RS iteration 14: Global Deviance = 315.2741 eps = 0.000013     
GAMLSS-RS iteration 15: Global Deviance = 315.2701 eps = 0.000012     
GAMLSS-RS iteration 16: Global Deviance = 315.2664 eps = 0.000011     
GAMLSS-RS iteration 17: Global Deviance = 315.2628 eps = 0.000011     
GAMLSS-RS iteration 18: Global Deviance = 315.2595 eps = 0.000010     
GAMLSS-RS iteration 19: Global Deviance = 315.2563 eps = 0.000010     
GAMLSS-RS iteration 20: Global Deviance = 315.2533 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.3816 eps = 0.164733     
GAMLSS-RS iteration  2: Global Deviance = 1561.3704 eps = 0.000007     

## 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 = 1426.6646 eps = 0.236800     
GAMLSS-RS iteration  2: Global Deviance = 1425.7016 eps = 0.000674     
GAMLSS-RS iteration  3: Global Deviance = 1425.6912 eps = 0.000007     

## 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 = 53766.4377 eps = 0.449136     
GAMLSS-RS iteration  2: Global Deviance = 53756.9844 eps = 0.000175     
GAMLSS-RS iteration  3: Global Deviance = 53755.6652 eps = 0.000024     
GAMLSS-RS iteration  4: Global Deviance = 53754.912 eps = 0.000014     
GAMLSS-RS iteration  5: Global Deviance = 53754.4264 eps = 0.000009     

## 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)