Function fits Linear Mixed Models (LMM) using Restricted Maximum Likelihood (REML).
(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.frame) containing all variables referenced in 'form'
(factor, character) variable specifying groups for which the analysis should be performed individually, i.e. by-processing
(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
(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
(logical) TRUE = will suppress any warning, which will be issued otherwise
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
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=σ^2/se(σ^2), and σ^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
As for objects returned by function
anovaMM linear hypotheses of fixed effects or LS Means can be
tested with functions
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).
Andre Schuetzenmeister email@example.com
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## 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) 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+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)