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(fixed =~ 1, random = NULL, ...)

Arguments

fixed	A formula that specifies the fixed effects of the nlme{lme} model. In most cases, this can also be included in the gamlss2 formula.
random	A formula specifying the random effect part of the model, as in the nlme{lme()} function.
…	For the re() function, the dots argument is used to specify additional control arguments for the nlme{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 ~ s(age,k=3) + re(random =~ age11 | Subject)

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

GAMLSS-RS iteration  1: Global Deviance = 326.6566 eps = 0.392364     
GAMLSS-RS iteration  2: Global Deviance = 326.0942 eps = 0.001721     
GAMLSS-RS iteration  3: Global Deviance = 326.0782 eps = 0.000049     
GAMLSS-RS iteration  4: Global Deviance = 326.0769 eps = 0.000004     

## 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: x$data 
       AIC     BIC    logLik
  448.9844 462.395 -219.4922

Random effects:
 Formula: ~age11 | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 2.0910678 (Intr)
age11       0.2157607 0.518 
Residual    1.1939543       

Variance function:
 Structure: fixed weights
 Formula: ~weights_w 
Fixed effects:  x$fixed 
                    Value Std.Error DF       t-value p-value
(Intercept) -2.074577e-07  0.404892 81 -5.123779e-07       1

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.31636025 -0.47633680  0.02042094  0.46926855  3.94312339 

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 ~ s(Time) + re(random =~ sin1 | Mare,
  correlation = corARMA(q = 2), control = lmeControl(maxIter = 100))

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

GAMLSS-RS iteration  1: Global Deviance = 1560.1177 eps = 0.165409     
GAMLSS-RS iteration  2: Global Deviance = 1551.2682 eps = 0.005672     
GAMLSS-RS iteration  3: Global Deviance = 1551.132 eps = 0.000087     
GAMLSS-RS iteration  4: Global Deviance = 1551.1212 eps = 0.000006     

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