VCA (version 1.3.1)

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)

Arguments

form

(formula) specifying the model to be fit, a response variable left of the '~' is mandatory

Data

(data.frame) storing 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 warning, which will be issued otherwise

Details

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).

See Also

remlVCA, VCAinference, ranef.VCA, residuals.VCA, anovaVCA, anovaMM, plotRandVar, test.fixef, test.lsmeans, lmer

Examples

## Not run: 
# 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)
# 
# # use different approaches to estimating the covariance of
# # variance components (covariance parameters)
# dat.ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]			# get unbalanced data
# m1.ub <- remlMM(y~day/(run), dat.ub, SSQ.method="qf", VarVC.method="scm")
# m2.ub <- remlMM(y~day/(run), dat.ub, SSQ.method="qf", VarVC.method="gb")		# is faster
# V1.ub <- round(vcovVC(m1.ub), 12)
# V2.ub <- round(vcovVC(m2.ub), 12)
# all(V1.ub == V2.ub)
# 
# # make it explicit that "gb" is faster than "scm"
# # compute variance-covariance matrix of VCs 10-times
# 
# system.time(for(i in 1:500) vcovVC(m1.ub))	# "scm"
# system.time(for(i in 1:500) vcovVC(m2.ub))	# "gb"
# 
# 
# # 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+Sex:age2+(Subject)+(Subject):age2, 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)
# ## End(Not run)