Random Effects

Description

There are two ways of fitting a random effect within gamlss2. The first, using s(), is for a simple random effect, that is, when only one factor is entered the model as a smoother. This method uses the function s() of the package mgcv with argument bs = “re”. For example, if area is factor with several levels, s(area, bs = “re”) will sringh the levels of area towards their mean level. The second, more general way, allows to fit more complicated random effect models using the function re(). The function re() is an interface connecting gamlss2 with the specialised package for random effects nlme.

Here we document only the re() function only but we also give 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 parameter 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
b <- gamlss2(distance ~ s(age,k=3) + re(random =~ age11 | Subject),
  data = Orthodont)

GAMLSS-RS iteration  1: Global Deviance = 326.6569 eps = 0.392363     
GAMLSS-RS iteration  2: Global Deviance = 326.0931 eps = 0.001726     
GAMLSS-RS iteration  3: Global Deviance = 326.0758 eps = 0.000053     
GAMLSS-RS iteration  4: Global Deviance = 326.0757 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", elements = "model")
summary(st)

                             Length Class  Mode
mu.s(age)                     9     -none- list
mu.re(random=~age11|Subject) 19     lme    list

## 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.3879 eps = 0.204987     
GAMLSS-RS iteration  2: Global Deviance = 370.2269 eps = 0.133744     
GAMLSS-RS iteration  3: Global Deviance = 340.2986 eps = 0.080837     
GAMLSS-RS iteration  4: Global Deviance = 325.4086 eps = 0.043755     
GAMLSS-RS iteration  5: Global Deviance = 319.3603 eps = 0.018586     
GAMLSS-RS iteration  6: Global Deviance = 316.706 eps = 0.008311     
GAMLSS-RS iteration  7: Global Deviance = 315.7332 eps = 0.003071     
GAMLSS-RS iteration  8: Global Deviance = 315.3835 eps = 0.001107     
GAMLSS-RS iteration  9: Global Deviance = 315.2588 eps = 0.000395     
GAMLSS-RS iteration 10: Global Deviance = 315.2156 eps = 0.000137     
GAMLSS-RS iteration 11: Global Deviance = 315.2007 eps = 0.000047     
GAMLSS-RS iteration 12: Global Deviance = 315.1956 eps = 0.000016     
GAMLSS-RS iteration 13: Global Deviance = 315.1933 eps = 0.000007     

## 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 = 1553.7576 eps = 0.168812     
GAMLSS-RS iteration  2: Global Deviance = 1550.0142 eps = 0.002409     
GAMLSS-RS iteration  3: Global Deviance = 1549.9734 eps = 0.000026     
GAMLSS-RS iteration  4: Global Deviance = 1549.9731 eps = 0.000000     

## 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.1007 eps = 0.185747     
GAMLSS-RS iteration  2: Global Deviance = 1436.4746 eps = 0.056255     
GAMLSS-RS iteration  3: Global Deviance = 1426.662 eps = 0.006831     
GAMLSS-RS iteration  4: Global Deviance = 1425.7706 eps = 0.000624     
GAMLSS-RS iteration  5: Global Deviance = 1425.692 eps = 0.000055     
GAMLSS-RS iteration  6: Global Deviance = 1425.6854 eps = 0.000004     

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