# 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. ```r 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 ## 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 ```r ## 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) ``` ## See Also `remlMM`, `VCAinference`, `ranef.VCA`, `residuals.VCA`, `anovaVCA`, `anovaMM`, `plotRandVar`, `lmer` ## Author(s) Andre Schuetzenmeister andre.schuetzenmeister@roche.com

remlVCA function