remlVCA: Perform (V)ariance (C)omponent (A)nalysis via REML-Estimation

Description

Function performs a Variance Component Analysis (VCA) using Restricted Maximum Likelihood (REML) to fit the random model, i.e. a linear mixed model (LMM) where the intercept is the only fixed effect.

Usage

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

Arguments

form: (formula) specifying the model to be fit, a response variable left of the '~' is mandatory
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
quiet: (logical) TRUE = will suppress any messages or warnings, 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

Here, a variance component model 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) is applied. 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. 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 very much (see examples).

Examples

Run this code

if (FALSE) {

# a VCA standard example
data(dataEP05A2_3)

# fit it by ANOVA first, then by REML
fit0 <- anovaVCA(y~day/run, dataEP05A2_3) 
fit1 <- remlVCA(y~day/run, dataEP05A2_3)
fit0
fit1

# make example unbalanced
set.seed(107)
dat.ub <- dataEP05A2_3[-sample(1:80, 7),]
fit0ub <- anovaVCA(y~day/run, dat.ub) 
fit1ub <- remlVCA(y~day/run, dat.ub) 

# not that ANOVA- and REML-results now differ
fit0ub
fit1ub

### Use the six sample reproducibility data from CLSI EP5-A3
### and fit per sample reproducibility model
data(CA19_9)
fit.all <- remlVCA(result~site/day, CA19_9, by="sample")

reproMat <- data.frame(
Sample=c("P1", "P2", "Q3", "Q4", "P5", "Q6"),
Mean= c(fit.all[[1]]$Mean, fit.all[[2]]$Mean, fit.all[[3]]$Mean, 
fit.all[[4]]$Mean, fit.all[[5]]$Mean, fit.all[[6]]$Mean),
Rep_SD=c(fit.all[[1]]$aov.tab["error","SD"], fit.all[[2]]$aov.tab["error","SD"],
fit.all[[3]]$aov.tab["error","SD"], fit.all[[4]]$aov.tab["error","SD"],
fit.all[[5]]$aov.tab["error","SD"], fit.all[[6]]$aov.tab["error","SD"]),
Rep_CV=c(fit.all[[1]]$aov.tab["error","CV[%]"],fit.all[[2]]$aov.tab["error","CV[%]"],
fit.all[[3]]$aov.tab["error","CV[%]"],fit.all[[4]]$aov.tab["error","CV[%]"],
fit.all[[5]]$aov.tab["error","CV[%]"],fit.all[[6]]$aov.tab["error","CV[%]"]),
WLP_SD=c(sqrt(sum(fit.all[[1]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[2]]$aov.tab[3:4, "VC"])),
sqrt(sum(fit.all[[3]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[4]]$aov.tab[3:4, "VC"])),
sqrt(sum(fit.all[[5]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[6]]$aov.tab[3:4, "VC"]))),
WLP_CV=c(sqrt(sum(fit.all[[1]]$aov.tab[3:4,"VC"]))/fit.all[[1]]$Mean*100,
sqrt(sum(fit.all[[2]]$aov.tab[3:4,"VC"]))/fit.all[[2]]$Mean*100,
sqrt(sum(fit.all[[3]]$aov.tab[3:4,"VC"]))/fit.all[[3]]$Mean*100,
sqrt(sum(fit.all[[4]]$aov.tab[3:4,"VC"]))/fit.all[[4]]$Mean*100,
sqrt(sum(fit.all[[5]]$aov.tab[3:4,"VC"]))/fit.all[[5]]$Mean*100,
sqrt(sum(fit.all[[6]]$aov.tab[3:4,"VC"]))/fit.all[[6]]$Mean*100),
Repro_SD=c(fit.all[[1]]$aov.tab["total","SD"],fit.all[[2]]$aov.tab["total","SD"],
fit.all[[3]]$aov.tab["total","SD"],fit.all[[4]]$aov.tab["total","SD"],
fit.all[[5]]$aov.tab["total","SD"],fit.all[[6]]$aov.tab["total","SD"]),
Repro_CV=c(fit.all[[1]]$aov.tab["total","CV[%]"],fit.all[[2]]$aov.tab["total","CV[%]"],
fit.all[[3]]$aov.tab["total","CV[%]"],fit.all[[4]]$aov.tab["total","CV[%]"],
fit.all[[5]]$aov.tab["total","CV[%]"],fit.all[[6]]$aov.tab["total","CV[%]"]))

for(i in 3:8) reproMat[,i] <- round(reproMat[,i],digits=ifelse(i%%2==0,1,3))
reproMat

# now plot the precision profile over all samples
plot(reproMat[,"Mean"], reproMat[,"Rep_CV"], type="l", main="Precision Profile CA19-9",
xlab="Mean CA19-9 Value", ylab="CV[%]")
grid()
points(reproMat[,"Mean"], reproMat[,"Rep_CV"], pch=16)

# REML-estimation not yes optimzed to the same degree as
# ANOVA-estimation. Note, that no variance-covariance matrix
# for the REML-fit is computed (VarVC=FALSE)!
# Note: A correct analysis would be done per-sample, this is just
#       for illustration.
data(VCAdata1)
# with complete sweeping implemented as FORTRAN-routine fit 
system.time(fit0 <- anovaVCA(y~sample+(device+lot)/day/run, VCAdata1))
system.time(fit1 <- remlVCA(y~sample+(device+lot)/day/run, VCAdata1, VarVC=FALSE))

# The previous example will also be interesting for environments using MKL.
# Run it once in a GNU-R environment and once in a MKL-environment
# and compare computational time of both. Note, that 'VarVC' is now set to TRUE
# and variable "sample" is put into the brackets increasing the number of random
# effects by factor 10. On my Intel Xeon E5-2687W 3.1 GHz workstation it takes
# ~ 400s with GNU-R and ~25s with MKL support (MRO) both run under Windows.
system.time(fit2 <- remlVCA(y~(sample+device+lot)/day/run, VCAdata1, VarVC=TRUE))

# using the SWEEP-Operator is even faster 
system.time(fit3 <- anovaVCA(y~(sample+device+lot)/day/run, VCAdata1))
fit2
fit3
}