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.0754 eps = 0.000001     

## 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.4959 eps = 0.204786     
GAMLSS-RS iteration  2: Global Deviance = 370.4036 eps = 0.133550     
GAMLSS-RS iteration  3: Global Deviance = 359.232 eps = 0.030160     
GAMLSS-RS iteration  4: Global Deviance = 334.3919 eps = 0.069147     
GAMLSS-RS iteration  5: Global Deviance = 322.9054 eps = 0.034350     
GAMLSS-RS iteration  6: Global Deviance = 318.0482 eps = 0.015042     
GAMLSS-RS iteration  7: Global Deviance = 316.2218 eps = 0.005742     
GAMLSS-RS iteration  8: Global Deviance = 315.5587 eps = 0.002096     
GAMLSS-RS iteration  9: Global Deviance = 315.3218 eps = 0.000750     
GAMLSS-RS iteration 10: Global Deviance = 315.2367 eps = 0.000269     
GAMLSS-RS iteration 11: Global Deviance = 315.2063 eps = 0.000096     
GAMLSS-RS iteration 12: Global Deviance = 315.1955 eps = 0.000034     
GAMLSS-RS iteration 13: Global Deviance = 315.1917 eps = 0.000012     
GAMLSS-RS iteration 14: Global Deviance = 315.1902 eps = 0.000004     

## 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.1178 eps = 0.165409     
GAMLSS-RS iteration  2: Global Deviance = 1551.3069 eps = 0.005647     
GAMLSS-RS iteration  3: Global Deviance = 1551.1747 eps = 0.000085     
GAMLSS-RS iteration  4: Global Deviance = 1551.1728 eps = 0.000001     

## 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.4785 eps = 0.056256     
GAMLSS-RS iteration  3: Global Deviance = 1426.6643 eps = 0.006832     
GAMLSS-RS iteration  4: Global Deviance = 1425.7728 eps = 0.000624     
GAMLSS-RS iteration  5: Global Deviance = 1425.6942 eps = 0.000055     
GAMLSS-RS iteration  6: Global Deviance = 1425.6873 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)
 }