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

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