anovaVCA: ANOVA-Type Estimation of Variance Components for Random...
In VCA: Variance Component Analysis
anovaVCA R Documentation
ANOVA-Type Estimation of Variance Components for Random Models
Description
This function equates observed ANOVA Type-I sums of squares (SS) to their expected values and solves the resulting system of linear equations
for variance components.
Usage
anovaVCA(
form,
Data,
by = NULL,
NegVC = FALSE,
VarVC.method = c("scm", "gb"),
MME = FALSE,
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
NegVC
(logical) FALSE = negative variance component estimates (VC) will be set to 0 and they will not contribute to the total variance
(as done in SAS PROC NESTED, conservative estimate of total variance). The original ANOVA estimates can be found in element 'VCoriginal'.
The degrees of freedom of the total variance are based on adapted mean squares (MS), i.e. adapted MS are computed as D * VC, where VC is
the column vector with negative VCs set to 0.
TRUE = negative variance component estimates will not be set to 0 and they will contribute to the total variance (original definition of the total variance).
VarVC.method
(character) string specifying whether to use the algorithm given in Searle et al. (1992) which corresponds to VarVC.method="scm" or in
Giesbrecht and Burns (1985) which can be specified via "gb". Method "scm" (Searle, Casella, McCulloch)
is the exact algorithm, "gb" (Giesbrecht, Burns) is termed "rough approximation"
by the authors, but sufficiently exact compared to e.g. SAS PROC MIXED (method=type1) which
uses the inverse of the Fisher-Information matrix as approximation. For balanced designs all
methods give identical results, only in unbalanced designs differences occur.
MME
(logical) TRUE = (M)ixed (M)odel (E)quations will be solved, i.e. 'VCA' object will have additional elements
"RandomEffects", "FixedEffects", "VarFixed" (variance-covariance matrix of fixed effects) and the "Matrices"
element has addional elements corresponding to intermediate results of solving MMEs.
FALSE = do not solve MMEs, which reduces the computation time for very complex models significantly.
quiet
(logical) TRUE = will suppress any warning, 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
For diagnostics, a key parameter is "precision", i.e. the accuracy of a quantification method influenced by varying sources of random error.
This type of experiments is requested by regulatory authorities to proof the quality of diagnostic tests, e.g. quantifying intermediate
precision according to CLSI guideline EP5-A2/A3. No, fixed effects are allowed besides the intercept.
Whenever fixed effects are part of the model to be analyzed, use function anovaMM instead.
Function anovaVCA is tailored for performing Variance Component Analyses (VCA) for random models, assuming all VCs as factor variables, i.e. their levels
correspond to distinct columns in the design matrix (dummy variables). Any predictor variables are automatically converted to factor variables, since continuous
variables may not be used on the right side of the formula 'form'.
ANOVA SS are computed employing the SWEEP-operator (Goodnight 1979, default).
according to Searle et al. (1992) which corresponds to VarVC.method="scm".
Function anovaVCA represents a special form of the "method of moments" approach applicable to arbitrary random models either balanced or unbalanced.
The system of linear equations, which is built from the ANOVA Type-I sums of squares, is closely related to the method used
by SAS PROC VARCOMP, where ANOVA mean squares (MS) are used. The former can be written as ss = C * s
and the latter as ms = D * s, where C and D denote the respective coefficient matrices, s the column-vector
of variance components (VC) to be estimated/predicted, and ss and ms the column vector of ANOVA sum of squares, respectively, mean squares.
Mutliplying element d_ij of matrix D by element c_in of matrix C (i,j = 1,...,n), results in
matrix C. Thus, C can easily be converted to D by the inverse operation. Matrix D is used to estimate
total degrees of freedom (DF) according to Satterthwaite (1946).
The method for computing ANOVA Type-I SS is much faster than fitting the linear model via lm and calling function anova on the 'lm' object
for complex models, where complex refers to the number of columns of the design matrix and the degree of unbalancedness. DF are directly derived from the SWEEP-operator as the number of linearly independent
columns of the partial design matrix corresponding to a specific VC.
Value
(object) of class 'VCA'
Author(s)
Andre Schuetzenmeister andre.schuetzenmeister@roche.com
References
Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York
Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158
Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample
Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486
Satterthwaite, F.E. (1946), An Approximate Distribution of Estimates of Variance Components.,
Biometrics Bulletin 2, 110-114
Gaylor,D.W., Lucas,H.L., Anderson,R.L. (1970), Calculation of Expected Mean Squares by the Abbreviated Doolittle and Square Root Methods.,
Biometrics 26 (4): 641-655
SAS Help and Documentation PROC MIXED, SAS Institute Inc., Cary, NC, USA
SAS Help and Documentation PROC VARCOMP, SAS Institute Inc., Cary, NC, USA
See Also
anovaMM, remlVCA, remlMM, print.VCA, VCAinference,
ranef, plotRandVar, stepwiseVCA
Examples
## Not run:
# load data (CLSI EP05-A2 Within-Lab Precision Experiment)
data(dataEP05A2_2)
# perform ANOVA-estimation of variance components
res <- anovaVCA(y~day/run, dataEP05A2_2)
res
# design with two main effects (ignoring the hierarchical structure of the design)
anovaVCA(y~day+run, dataEP05A2_2)
# compute confidence intervals, perform F- and Chi-Squared tests
INF <- VCAinference(res, total.claim=3.5, error.claim=2)
INF
### load data from package
data(VCAdata1)
data_sample1 <- VCAdata1[VCAdata1$sample==1,]
### plot data for visual inspection
varPlot(y~lot/day/run, data_sample1)
### estimate VCs for 4-level hierarchical design (error counted) for sample_1 data
anovaVCA(y~lot/day/run, data_sample1)
### using different model (ignoring the hierarchical structure of the design)
anovaVCA(y~lot+day+lot:day:run, data_sample1)
### same model with unbalanced data
anovaVCA(y~lot+day+lot:day:run, data_sample1[-c(1,11,15),])
### use the numerical example from the CLSI EP05-A2 guideline (p.25)
data(Glucose,package="VCA")
res.ex <- anovaVCA(result~day/run, Glucose)
### also perform Chi-Squared tests
### Note: in guideline claimed SD-values are used, here, claimed variances are used
VCAinference(res.ex, total.claim=3.4^2, error.claim=2.5^2)
### now use the six sample reproducibility data from CLSI EP5-A3
### and fit per sample reproducibility model
data(CA19_9)
fit.all <- anovaVCA(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)
# load another example dataset and extract the "sample==1" subset
data(VCAdata1)
sample1 <- VCAdata1[which(VCAdata1$sample==1),]
# generate an additional factor variable and random errors according to its levels
sample1$device <- gl(3,28,252)
set.seed(505)
sample1$y <- sample1$y + rep(rep(rnorm(3,,.25), c(28,28,28)),3)
# fit a crossed-nested design with main factors 'lot' and 'device'
# and nested factors 'day' and 'run' nested below
res1 <- anovaVCA(y~(lot+device)/day/run, sample1)
res1
# fit same model for each sample using by-processing
lst <- anovaVCA(y~(lot+device)/day/run, VCAdata1, by="sample")
lst
# now fitting a nonsense model on the complete dataset "VCAdata1"
# where the SWEEP-operator is the new default since package version 1.2
# takes ~5s
system.time(res.sw <- anovaVCA(y~(sample+lot+device)/day/run, VCAdata1))
# applying functions 'anova' and 'lm' in the same manner takes ~ 265s
system.time(res.lm <- anova(lm(y~(sample+lot+device)/day/run, VCAdata1)))
res.sw
res.lm
## End(Not run)
VCA documentation built on May 29, 2024, 1:48 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| anovaVCA | R Documentation |
ANOVA-Type Estimation of Variance Components for Random Models
Description
This function equates observed ANOVA Type-I sums of squares (SS) to their expected values and solves the resulting system of linear equations
for variance components.
Usage
anovaVCA(
form,
Data,
by = NULL,
NegVC = FALSE,
VarVC.method = c("scm", "gb"),
MME = FALSE,
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 |
NegVC |
(logical) FALSE = negative variance component estimates (VC) will be set to 0 and they will not contribute to the total variance
(as done in SAS PROC NESTED, conservative estimate of total variance). The original ANOVA estimates can be found in element 'VCoriginal'.
The degrees of freedom of the total variance are based on adapted mean squares (MS), i.e. adapted MS are computed as |
VarVC.method |
(character) string specifying whether to use the algorithm given in Searle et al. (1992) which corresponds to |
MME |
(logical) TRUE = (M)ixed (M)odel (E)quations will be solved, i.e. 'VCA' object will have additional elements "RandomEffects", "FixedEffects", "VarFixed" (variance-covariance matrix of fixed effects) and the "Matrices" element has addional elements corresponding to intermediate results of solving MMEs. FALSE = do not solve MMEs, which reduces the computation time for very complex models significantly. |
quiet |
(logical) TRUE = will suppress any warning, 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
For diagnostics, a key parameter is "precision", i.e. the accuracy of a quantification method influenced by varying sources of random error.
This type of experiments is requested by regulatory authorities to proof the quality of diagnostic tests, e.g. quantifying intermediate
precision according to CLSI guideline EP5-A2/A3. No, fixed effects are allowed besides the intercept.
Whenever fixed effects are part of the model to be analyzed, use function anovaMM instead.
Function anovaVCA is tailored for performing Variance Component Analyses (VCA) for random models, assuming all VCs as factor variables, i.e. their levels
correspond to distinct columns in the design matrix (dummy variables). Any predictor variables are automatically converted to factor variables, since continuous
variables may not be used on the right side of the formula 'form'.
ANOVA SS are computed employing the SWEEP-operator (Goodnight 1979, default).
according to Searle et al. (1992) which corresponds to VarVC.method="scm".
Function anovaVCA represents a special form of the "method of moments" approach applicable to arbitrary random models either balanced or unbalanced.
The system of linear equations, which is built from the ANOVA Type-I sums of squares, is closely related to the method used
by SAS PROC VARCOMP, where ANOVA mean squares (MS) are used. The former can be written as ss = C * s
and the latter as ms = D * s, where C and D denote the respective coefficient matrices, s the column-vector
of variance components (VC) to be estimated/predicted, and ss and ms the column vector of ANOVA sum of squares, respectively, mean squares.
Mutliplying element d_ij of matrix D by element c_in of matrix C (i,j = 1,...,n), results in
matrix C. Thus, C can easily be converted to D by the inverse operation. Matrix D is used to estimate
total degrees of freedom (DF) according to Satterthwaite (1946).
The method for computing ANOVA Type-I SS is much faster than fitting the linear model via lm and calling function anova on the 'lm' object
for complex models, where complex refers to the number of columns of the design matrix and the degree of unbalancedness. DF are directly derived from the SWEEP-operator as the number of linearly independent
columns of the partial design matrix corresponding to a specific VC.
Value
(object) of class 'VCA'
Author(s)
Andre Schuetzenmeister andre.schuetzenmeister@roche.com
References
Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York
Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158
Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486
Satterthwaite, F.E. (1946), An Approximate Distribution of Estimates of Variance Components., Biometrics Bulletin 2, 110-114
Gaylor,D.W., Lucas,H.L., Anderson,R.L. (1970), Calculation of Expected Mean Squares by the Abbreviated Doolittle and Square Root Methods., Biometrics 26 (4): 641-655
SAS Help and Documentation PROC MIXED, SAS Institute Inc., Cary, NC, USA
SAS Help and Documentation PROC VARCOMP, SAS Institute Inc., Cary, NC, USA
See Also
anovaMM, remlVCA, remlMM, print.VCA, VCAinference,
ranef, plotRandVar, stepwiseVCA
Examples
## Not run:
# load data (CLSI EP05-A2 Within-Lab Precision Experiment)
data(dataEP05A2_2)
# perform ANOVA-estimation of variance components
res <- anovaVCA(y~day/run, dataEP05A2_2)
res
# design with two main effects (ignoring the hierarchical structure of the design)
anovaVCA(y~day+run, dataEP05A2_2)
# compute confidence intervals, perform F- and Chi-Squared tests
INF <- VCAinference(res, total.claim=3.5, error.claim=2)
INF
### load data from package
data(VCAdata1)
data_sample1 <- VCAdata1[VCAdata1$sample==1,]
### plot data for visual inspection
varPlot(y~lot/day/run, data_sample1)
### estimate VCs for 4-level hierarchical design (error counted) for sample_1 data
anovaVCA(y~lot/day/run, data_sample1)
### using different model (ignoring the hierarchical structure of the design)
anovaVCA(y~lot+day+lot:day:run, data_sample1)
### same model with unbalanced data
anovaVCA(y~lot+day+lot:day:run, data_sample1[-c(1,11,15),])
### use the numerical example from the CLSI EP05-A2 guideline (p.25)
data(Glucose,package="VCA")
res.ex <- anovaVCA(result~day/run, Glucose)
### also perform Chi-Squared tests
### Note: in guideline claimed SD-values are used, here, claimed variances are used
VCAinference(res.ex, total.claim=3.4^2, error.claim=2.5^2)
### now use the six sample reproducibility data from CLSI EP5-A3
### and fit per sample reproducibility model
data(CA19_9)
fit.all <- anovaVCA(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)
# load another example dataset and extract the "sample==1" subset
data(VCAdata1)
sample1 <- VCAdata1[which(VCAdata1$sample==1),]
# generate an additional factor variable and random errors according to its levels
sample1$device <- gl(3,28,252)
set.seed(505)
sample1$y <- sample1$y + rep(rep(rnorm(3,,.25), c(28,28,28)),3)
# fit a crossed-nested design with main factors 'lot' and 'device'
# and nested factors 'day' and 'run' nested below
res1 <- anovaVCA(y~(lot+device)/day/run, sample1)
res1
# fit same model for each sample using by-processing
lst <- anovaVCA(y~(lot+device)/day/run, VCAdata1, by="sample")
lst
# now fitting a nonsense model on the complete dataset "VCAdata1"
# where the SWEEP-operator is the new default since package version 1.2
# takes ~5s
system.time(res.sw <- anovaVCA(y~(sample+lot+device)/day/run, VCAdata1))
# applying functions 'anova' and 'lm' in the same manner takes ~ 265s
system.time(res.lm <- anova(lm(y~(sample+lot+device)/day/run, VCAdata1)))
res.sw
res.lm
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.