remlMM: Fit Linear Mixed Models via REML

Description

Function fits Linear Mixed Models (LMM) using Restricted Maximum Likelihood (REML).

Usage

remlMM(
  form,
  Data,
  by = NULL,
  VarVC = TRUE,
  cov = TRUE,
  quiet = FALSE,
  order.data = TRUE
)

Arguments

form: (formula) specifying the model to be fit, a response variable left of the '~' is mandatory, random terms have to be enclosed in brackets (see details for definition of valid model terms)
Data: (data.frame) containing all variables referenced in 'form'
by: (factor, character) variable specifying groups for which the analysis should be performed individually, i.e. by-processing
VarVC: (logical) TRUE = the variance-covariance matrix of variance components will be approximated using the method found in Giesbrecht & Burns (1985), which also serves as basis for applying a Satterthwaite approximation of the degrees of freedom for each variance component, FALSE = leaves out this step, no confidence intervals for VC will be available
cov: (logical) TRUE = in case of non-zero covariances a block diagonal matrix will be constructed, FALSE = a diagonal matrix with all off-diagonal element being equal to zero will be contructed
quiet: (logical) TRUE = will suppress any messages or warning, which will be issued otherwise
order.data: (logical) TRUE = class-variables will be ordered increasingly, FALSE = ordering of class-variables will remain as is

Author

Andre Schuetzenmeister andre.schuetzenmeister@roche.com

Details

The model is formulated exactly as in function anovaMM, i.e. random terms need be enclosed by round brackets. All terms appearing in the model (fixed or random) need to be compliant with the regular expression "^[^[\.]]?[[:alnum:]_\.]*$", i.e. they may not start with a dot and may then only consist of alpha-numeric characters, dot and underscore. Otherwise, an error will be issued.

Here, a LMM is fitted by REML using the lmer function of the lme4-package. For all models the Giesbrechnt & Burns (1985) approximation of the variance-covariance matrix of variance components (VC) can be applied ('VarVC=TRUE'). A Satterthwaite approximation of the degrees of freedom for all VC and total variance is based on this approximated matrix using $df=2Z^2$, where $Z$ is the Wald statistic $Z=\sigma^2/se(\sigma^2)$, and $\sigma^2$ is here used for an estimated variance. The variance of total variability, i.e. the sum of all VC is computed via summing up all elements of the variance-covariance matrix of the VC. One can constrain the variance-covariance matrix of random effects $G$ to be either diagonal ('cov=FALSE'), i.e. all random effects are indpendent of each other (covariance is 0). If 'cov=TRUE' (the default) matrix $G$ will be constructed as implied by the model returned by function lmer.

As for objects returned by function anovaMM linear hypotheses of fixed effects or LS Means can be tested with functions test.fixef and test.lsmeans. Note, that option "contain" does not work for LMM fitted via REML.

Note, that for large datasets approximating the variance-covariance matrix of VC is computationally expensive and may take very long. There is no Fisher-information matrix available for 'merMod' objects, which can serve as approximation. To avoid this time-consuming step, use argument 'VarVC=FALSE' but remember, that no confidence intervals for any VC will be available. If you use Microsoft's R Open, formerly known as Revolution-R, which comes with Intel's Math Kernel Library (MKL), this will be automatically detected and an environment-optimized version will be used, reducing the computational time considerably (see examples).

Examples

Run this code

if (FALSE) {
data(dataEP05A2_2)

# assuming 'day' as fixed, 'run' as random
remlMM(y~day/(run), dataEP05A2_2)

# assuming both as random leads to same results as
# calling anovaVCA
remlMM(y~(day)/(run), dataEP05A2_2)
anovaVCA(y~day/run, dataEP05A2_2)
remlVCA(y~day/run, dataEP05A2_2)

# fit a larger random model
data(VCAdata1)
fitMM1 <- remlMM(y~((lot)+(device))/(day)/(run), VCAdata1[VCAdata1$sample==1,])
fitMM1
# now use function tailored for random models
fitRM1 <- anovaVCA(y~(lot+device)/day/run, VCAdata1[VCAdata1$sample==1,])
fitRM1

# there are only 3 lots, take 'lot' as fixed 
fitMM2 <- remlMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])

# the following model definition is equivalent to the one above,
# since a single random term in an interaction makes the interaction
# random (see the 3rd reference for details on this topic)
fitMM3 <- remlMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])

# fit same model for each sample using by-processing
lst <- remlMM(y~(lot+(device))/day/run, VCAdata1, by="sample")
lst

# fit mixed model originally from 'nlme' package

library(nlme)
data(Orthodont)
fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 

# re-organize data for using 'remlMM'
Ortho <- Orthodont
Ortho$age2 <- Ortho$age - 11
Ortho$Subject <- factor(as.character(Ortho$Subject))
fit.remlMM1 <- remlMM(distance~Sex*age2+(Subject)*age2, Ortho)

# use simplified formula avoiding unnecessary terms
fit.remlMM2 <- remlMM(distance~Sex+age2+Sex:age2+(Subject)+age2:(Subject), Ortho)

# and exclude intercept
fit.remlMM3 <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)

# now use exclude covariance of per-subject intercept and slope
# as for models fitted by function 'anovaMM'
fit.remlMM4 <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho, cov=FALSE)

# compare results
fit.lme
fit.remlMM1
fit.remlMM2
fit.remlMM3
fit.remlMM4

# are there a sex-specific differences?
cmat <- getL(fit.remlMM3, c("SexMale-SexFemale", "SexMale:age2-SexFemale:age2")) 
cmat

test.fixef(fit.remlMM3, L=cmat)
}