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.6292 eps = 0.000085     
GAMLSS-RS iteration  3: Global Deviance = 326.0947 eps = 0.001636     
GAMLSS-RS iteration  4: Global Deviance = 326.0795 eps = 0.000046     
GAMLSS-RS iteration  5: Global Deviance = 326.0778 eps = 0.000005     

## 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.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.1178 eps = 0.165409     
GAMLSS-RS iteration  2: Global Deviance = 1551.2684 eps = 0.005672     
GAMLSS-RS iteration  3: Global Deviance = 1551.1322 eps = 0.000087     
GAMLSS-RS iteration  4: Global Deviance = 1551.1214 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)
 }