Perform (V)ariance (C)omponent (A)nalysis via REML-Estimation
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.
remlVCA(form, Data, by = NULL, VarVC = TRUE, quiet = FALSE, order.data = TRUE)
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 |
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=σ^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.
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).
Andre Schuetzenmeister andre.schuetzenmeister@roche.com
## Not run:
# 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
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.