fitLMM: Fit Linear Mixed Model by ANOVA or REML

In VCA: Variance Component Analysis

Fit Linear Mixed Model by ANOVA or REML

Description

Function serves as interface to functions anovaMM and remlMM for fitting a linear mixed model (LMM) either by ANOVA or REML. All arguments applicable to either one of these functions can be specified (see anovaMM or remlMM for details).

Usage

fitLMM(
  form,
  Data,
  method = c("anova", "reml"),
  scale = TRUE,
  VarVC = TRUE,
  ...
)

Arguments

form	(formula) specifiying the linear mixed model, random effects need to be identified by enclosing them in round brackets, i.e. ~a/(b) will model factor 'a' as fixed and 'b' as random
Data	(data.frame) containing all variables referenced in 'form', note that variables can only be of type "numeric", "factor" or "character". The latter will be automatically converted to "factor"
method	(character) either "anova" to use ANOVA Type-I estimation of variance components or "reml" to use restricted maximum likelihood (REML) estimation of variance component
scale	(logical) TRUE = scale values of the response aiming to avoid numerical problems when numbers are either very small or very large, FALSE = use original scale
VarVC	(logical) TRUE = variance-covariance matrix of variance components will be computed, FALSE = it will not be computed
...	additional arguments to be passed to function anovaMM or function remlMM.

Details

Besides offering a convenient interface to both functions for fitting a LMM, this function also provides all elements required for standard task of fitted models, e.g. prediction, testing general linear hypotheses via R-package multcomp, etc. (see examples).

Author(s)

Andre Schuetzenmeister [e-mail]

See Also

fitVCA, anovaMM, remlMM

Examples

## Not run: 
data(dataEP05A2_2)

# assuming 'day' as fixed, 'run' as random
# Note: default method is "anova"
fitLMM(y~day/(run), dataEP05A2_2)

# explicitly request "reml"
fitLMM(y~day/(run), dataEP05A2_2, method="reml")

# assuming both as random leads to same results as
# calling anovaVCA (ANOVA is the default)
fitLMM(y~(day)/(run), dataEP05A2_2)
anovaVCA(y~day/run, dataEP05A2_2)

# now using REML-estimation
fitLMM(y~(day)/(run), dataEP05A2_2, "reml")
remlVCA(y~day/run, dataEP05A2_2)

# use different approaches to estimating the covariance of 
# variance components (covariance parameters)
# create unbalanced data
dat.ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]
m1.ub <- fitLMM(y~day/(run), dat.ub, VarVC.method="scm")
# VarVC.method="gb" is an approximation not relying on quadratic forms
m2.ub <- fitLMM(y~day/(run), dat.ub, VarVC.method="gb")		
# REML-estimated variance components usually differ from ANOVA-estimates
# and so do the variance-covariance matrices
m3.ub <- fitLMM(y~day/(run), dat.ub, "reml", VarVC=TRUE)		
V1.ub <- round(vcovVC(m1.ub), 12)
V2.ub <- round(vcovVC(m2.ub), 12)
V3.ub <- round(vcovVC(m3.ub), 12)

# fit a larger random model
data(VCAdata1)
fitMM1 <- fitLMM(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 <- fitLMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])
# use REML on this (balanced) data
fitMM2.2 <- fitLMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,], "reml")

# 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 <- fitLMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])

# fit same model for each sample using by-processing
lst <- fitLMM(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
Ortho <- Orthodont
Ortho$age2 <- Ortho$age - 11
Ortho$Subject <- factor(as.character(Ortho$Subject))
fit.anovaMM1 <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho)

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

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

# compare results
fit.lme
fit.anovaMM1
fit.anovaMM2
fit.anovaMM3

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

test.fixef(fit.anovaMM3, L=cmat)

# fit LMM with fixed lot and device effects and test for lot-differences
data(VCAdata1)
fitS5 <- fitLMM(y~(lot+device)/(day)/(run), subset(VCAdata1, sample==5), "reml")
fitS5

# apply Tukey-HSD test to screen for lot differences
library(multcomp)
res.tuk <- glht(fitS5, linfct=mcp(lot="Tukey"))
summary(res.tuk)

# compact letter display
res.tuk.cld <- cld(res.tuk, col=paste0("gray", c(90,60,75)))
plot(res.tuk.cld)

## End(Not run)

VCA documentation built on May 29, 2024, 1:48 a.m.