| Title: | High-Dimensional Analysis of Variance |
|---|---|
| Description: | Functions and datasets to support Smilde, Marini, Westerhuis and Liland (2025, ISBN: 978-1-394-21121-0) "Analysis of Variance for High-Dimensional Data - Applications in Life, Food and Chemical Sciences". This implements and imports a collection of methods for HD-ANOVA data analysis with common interfaces, result- and plotting functions, multiple real data sets and four vignettes covering a range different applications. |
| Authors: | Kristian Hovde Liland [aut, cre] (ORCID: <https://orcid.org/0000-0001-6468-9423>) |
| Maintainer: | Kristian Hovde Liland <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.8.6 |
| Built: | 2026-06-01 10:41:36 UTC |
| Source: | https://github.com/khliland/hdanova |
APCA function for fitting ANOVA Principal Component Analysis models.
apca( formula, data, add_error = TRUE, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )apca( formula, data, add_error = TRUE, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )
formula |
Model formula accepting a single response (block) and predictors. |
data |
The data set to analyse. |
add_error |
Add error to LS means (default = TRUE). |
contrasts |
Effect coding: "sum" (default = sum-coding), "weighted", "reference", "treatment". |
permute |
Number of permutations to perform (default = 1000). |
perm.type |
Type of permutation to perform, either "approximate" or "exact" (default = "approximate"). |
... |
Additional parameters for the |
An object of class apca, inheriting from the general asca class.
Further arguments and plots can be found in the asca documentation.
Harrington, P.d.B., Vieira, N.E., Espinoza, J., Nien, J.K., Romero, R., and Yergey, A.L. (2005) Analysis of variance–principal component analysis: A soft tool for proteomic discovery. Analytica chimica acta, 544 (1-2), 118–127.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
data(candies) ap <- apca(assessment ~ candy, data=candies) scoreplot(ap) # Numeric effects candies$num <- eff <- 1:165 mod <- apca(assessment ~ candy + assessor + num, data=candies) summary(mod) scoreplot(mod, factor=3, gr.col=rgb(eff/max(eff), 1-eff/max(eff),0), pch.scores="x")data(candies) ap <- apca(assessment ~ candy, data=candies) scoreplot(ap) # Numeric effects candies$num <- eff <- 1:165 mod <- apca(assessment ~ candy + assessor + num, data=candies) summary(mod) scoreplot(mod, factor=3, gr.col=rgb(eff/max(eff), 1-eff/max(eff),0), pch.scores="x")
This is a quite general and flexible implementation of APLS.
apls( formula, data, add_error = TRUE, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )apls( formula, data, add_error = TRUE, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )
formula |
Model formula accepting a single response (block) and predictors. See Details for more information. |
data |
The data set to analyse. |
add_error |
Add error to LS means (default = TRUE). |
contrasts |
Effect coding: "sum" (default = sum-coding), "weighted", "reference", "treatment". |
permute |
Number of permutations to perform (default = 1000). |
perm.type |
Type of permutation to perform, either "approximate" or "exact" (default = "approximate"). |
... |
Additional arguments to |
APLS is a method which decomposes a multivariate response according to one or more design
variables. ANOVA is used to split variation into contributions from factors, and PLS is performed
on the corresponding least squares estimates, i.e., Y = X1 B1 + X2 B2 + ... + E = T1 P1' + T2 P2' + ... + E.
For balanced designs, the PLS components are equivalent to PCA components, i.e., APLS and APCA are equivalent.
This version of APLS encompasses variants of LiMM-PLS, generalized APLS and covariates APLS.
The formula interface is extended with the function r() to indicate random effects and comb() to indicate effects that should be combined. See Examples for use cases.
An apls object containing loadings, scores, explained variances, etc. The object has
associated plotting (asca_plots) and result (asca_results) functions.
Smilde, A., Jansen, J., Hoefsloot, H., Lamers,R., Van Der Greef, J., and Timmerman, M.(2005). ANOVA-Simultaneous Component Analysis (ASCA): A new tool for analyzing designed metabolomics data. Bioinformatics, 21(13), 3043–3048.
Liland, K.H., Smilde, A., Marini, F., and Næs,T. (2018). Confidence ellipsoids for ASCA models based on multivariate regression theory. Journal of Chemometrics, 32(e2990), 1–13.
Martin, M. and Govaerts, B. (2020). LiMM-PCA: Combining ASCA+ and linear mixed models to analyse high-dimensional designed data. Journal of Chemometrics, 34(6), e3232.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# Load candies data data(candies) # Basic APLS model with two factors mod <- apls(assessment ~ candy + assessor, data=candies) print(mod) # APLS model with interaction mod <- apls(assessment ~ candy * assessor, data=candies) print(mod) # Result plotting for first factor loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # No backprojection scoreplot(mod, projections=FALSE) # Spider plot scoreplot(mod, spider=TRUE) # APLS model with compressed response using 5 principal components mod.pca <- apls(assessment ~ candy + assessor, data=candies, pca.in=5) # Mixed Model APLS, random assessor mod.mix <- apls(assessment ~ candy + r(assessor), data=candies) scoreplot(mod.mix) # Mixed Model APLS, REML estimation mod.mix <- apls(assessment ~ candy + r(assessor), data=candies, REML=TRUE) scoreplot(mod.mix) # Load Caldana data data(caldana) # Combining effects in APLS mod.comb <- apls(compounds ~ time + comb(light + time:light), data=caldana) summary(mod.comb) timeplot(mod.comb, factor="light", time="time", comb=2) # Permutation testing mod.perm <- apls(assessment ~ candy * assessor, data=candies, permute=TRUE) summary(mod.perm)# Load candies data data(candies) # Basic APLS model with two factors mod <- apls(assessment ~ candy + assessor, data=candies) print(mod) # APLS model with interaction mod <- apls(assessment ~ candy * assessor, data=candies) print(mod) # Result plotting for first factor loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # No backprojection scoreplot(mod, projections=FALSE) # Spider plot scoreplot(mod, spider=TRUE) # APLS model with compressed response using 5 principal components mod.pca <- apls(assessment ~ candy + assessor, data=candies, pca.in=5) # Mixed Model APLS, random assessor mod.mix <- apls(assessment ~ candy + r(assessor), data=candies) scoreplot(mod.mix) # Mixed Model APLS, REML estimation mod.mix <- apls(assessment ~ candy + r(assessor), data=candies, REML=TRUE) scoreplot(mod.mix) # Load Caldana data data(caldana) # Combining effects in APLS mod.comb <- apls(compounds ~ time + comb(light + time:light), data=caldana) summary(mod.comb) timeplot(mod.comb, factor="light", time="time", comb=2) # Permutation testing mod.perm <- apls(assessment ~ candy * assessor, data=candies, permute=TRUE) summary(mod.perm)
This is a quite general and flexible implementation of ASCA.
asca( formula, data, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )asca( formula, data, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )
formula |
Model formula accepting a single response (block) and predictors. See Details for more information. |
data |
The data set to analyse. |
contrasts |
Effect coding: "sum" (default = sum-coding), "weighted", "reference", "treatment". |
permute |
Number of permutations to perform (default = 1000). |
perm.type |
Type of permutation to perform, either "approximate" or "exact" (default = "approximate"). |
... |
Additional arguments to |
ASCA is a method which decomposes a multivariate response according to one or more design
variables. ANOVA is used to split variation into contributions from factors, and PCA is performed
on the corresponding least squares estimates, i.e., Y = X1 B1 + X2 B2 + ... + E = T1 P1' + T2 P2' + ... + E.
This version of ASCA encompasses variants of LiMM-PCA, generalized ASCA and covariates ASCA. It includes
confidence ellipsoids for the balanced crossed-effect ASCA.
The formula interface is extended with the function r() to indicate random effects and comb() to indicate effects that should be combined. See Examples for use cases.
An asca object containing loadings, scores, explained variances, etc. The object has
associated plotting (asca_plots) and result (asca_results) functions.
Smilde, A., Jansen, J., Hoefsloot, H., Lamers,R., Van Der Greef, J., and Timmerman, M.(2005). ANOVA-Simultaneous Component Analysis (ASCA): A new tool for analyzing designed metabolomics data. Bioinformatics, 21(13), 3043–3048.
Liland, K.H., Smilde, A., Marini, F., and Næs,T. (2018). Confidence ellipsoids for ASCA models based on multivariate regression theory. Journal of Chemometrics, 32(e2990), 1–13.
Martin, M. and Govaerts, B. (2020). LiMM-PCA: Combining ASCA+ and linear mixed models to analyse high-dimensional designed data. Journal of Chemometrics, 34(6), e3232.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# Load candies data data(candies) # Basic ASCA model with two factors mod <- asca(assessment ~ candy + assessor, data=candies) print(mod) # ASCA model with interaction mod <- asca(assessment ~ candy * assessor, data=candies) print(mod) # Result plotting for first factor loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # No backprojection scoreplot(mod, projections=FALSE) # Spider plot scoreplot(mod, spider=TRUE) # ASCA model with compressed response using 5 principal components mod.pca <- asca(assessment ~ candy + assessor, data=candies, pca.in=5) # ASCA model with design-relevant compressed response # using 5 partial least squares components mod.pca <- asca(assessment ~ candy + assessor, data=candies, pls.in=5) # Mixed Model ASCA, random assessor mod.mix <- asca(assessment ~ candy + r(assessor), data=candies) scoreplot(mod.mix) # Mixed Model ASCA, REML estimation mod.mix <- asca(assessment ~ candy + r(assessor), data=candies, REML=TRUE) scoreplot(mod.mix) # Load Caldana data data(caldana) # Combining effects in ASCA mod.comb <- asca(compounds ~ time + comb(light + time:light), data=caldana) summary(mod.comb) timeplot(mod.comb, factor="light", time="time", comb=2) # Permutation testing mod.perm <- asca(assessment ~ candy * assessor, data=candies, permute=TRUE) summary(mod.perm) # Unbalanced data: compare native vs SS-type-aligned effects drop_idx <- which(candies$candy == levels(candies$candy)[1] & candies$assessor == levels(candies$assessor)[1]) candies.unbal <- candies[-drop_idx[seq_len(min(2, length(drop_idx)))], ] mod.native <- asca(assessment ~ candy * assessor, data = candies.unbal) mod.sstype <- asca(assessment ~ candy * assessor, data = candies.unbal, respect_SStype = TRUE) par.old <- par(mfrow = c(2,1), mar = c(4,4,2,1)) scoreplot(mod.native, main = "Top: native (regression-based LS)") scoreplot(mod.sstype, main = "Bottom: respect_SStype = TRUE") par(par.old) perm.native <- permutation(mod.native, permute = 1000, respect_SStype = FALSE) perm.sstype <- permutation(mod.sstype, permute = 1000, respect_SStype = TRUE) summary(perm.native) summary(perm.sstype)# Load candies data data(candies) # Basic ASCA model with two factors mod <- asca(assessment ~ candy + assessor, data=candies) print(mod) # ASCA model with interaction mod <- asca(assessment ~ candy * assessor, data=candies) print(mod) # Result plotting for first factor loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # No backprojection scoreplot(mod, projections=FALSE) # Spider plot scoreplot(mod, spider=TRUE) # ASCA model with compressed response using 5 principal components mod.pca <- asca(assessment ~ candy + assessor, data=candies, pca.in=5) # ASCA model with design-relevant compressed response # using 5 partial least squares components mod.pca <- asca(assessment ~ candy + assessor, data=candies, pls.in=5) # Mixed Model ASCA, random assessor mod.mix <- asca(assessment ~ candy + r(assessor), data=candies) scoreplot(mod.mix) # Mixed Model ASCA, REML estimation mod.mix <- asca(assessment ~ candy + r(assessor), data=candies, REML=TRUE) scoreplot(mod.mix) # Load Caldana data data(caldana) # Combining effects in ASCA mod.comb <- asca(compounds ~ time + comb(light + time:light), data=caldana) summary(mod.comb) timeplot(mod.comb, factor="light", time="time", comb=2) # Permutation testing mod.perm <- asca(assessment ~ candy * assessor, data=candies, permute=TRUE) summary(mod.perm) # Unbalanced data: compare native vs SS-type-aligned effects drop_idx <- which(candies$candy == levels(candies$candy)[1] & candies$assessor == levels(candies$assessor)[1]) candies.unbal <- candies[-drop_idx[seq_len(min(2, length(drop_idx)))], ] mod.native <- asca(assessment ~ candy * assessor, data = candies.unbal) mod.sstype <- asca(assessment ~ candy * assessor, data = candies.unbal, respect_SStype = TRUE) par.old <- par(mfrow = c(2,1), mar = c(4,4,2,1)) scoreplot(mod.native, main = "Top: native (regression-based LS)") scoreplot(mod.sstype, main = "Bottom: respect_SStype = TRUE") par(par.old) perm.native <- permutation(mod.native, permute = 1000, respect_SStype = FALSE) perm.sstype <- permutation(mod.sstype, permute = 1000, respect_SStype = TRUE) summary(perm.native) summary(perm.sstype)
Various plotting procedures for asca objects.
## S3 method for class 'asca' loadingplot(object, factor = 1, comps = 1:2, ...) ## S3 method for class 'asca' scoreplot( object, factor = 1, comps = 1:2, within_level = "all", pch.scores = 19, pch.projections = 1, gr.col = NULL, projections = TRUE, spider = FALSE, ellipsoids, confidence, xlim, ylim, xlab, ylab, legendpos, ... ) permutationplot(object, factor = 1, xlim, xlab = "SSQ", main, ...) rotationplot(object, factor = 1, xlim, xlab = "SSQ", main, ...)## S3 method for class 'asca' loadingplot(object, factor = 1, comps = 1:2, ...) ## S3 method for class 'asca' scoreplot( object, factor = 1, comps = 1:2, within_level = "all", pch.scores = 19, pch.projections = 1, gr.col = NULL, projections = TRUE, spider = FALSE, ellipsoids, confidence, xlim, ylim, xlab, ylab, legendpos, ... ) permutationplot(object, factor = 1, xlim, xlab = "SSQ", main, ...) rotationplot(object, factor = 1, xlim, xlab = "SSQ", main, ...)
object |
|
factor |
|
comps |
|
... |
additional arguments to underlying methods. |
within_level |
MSCA parameter for chosing plot level (default = "all"). |
pch.scores |
|
pch.projections |
|
gr.col |
|
projections |
Include backprojections in score plot (default = TRUE). |
spider |
Draw lines between group centers and backprojections (default = FALSE). |
ellipsoids |
|
confidence |
|
xlim |
|
ylim |
|
xlab |
|
ylab |
|
legendpos |
|
main |
Plot title. |
Usage of the functions are shown using generics in the examples in asca.
Plot routines are available as
scoreplot.asca and loadingplot.asca.
The plotting routines have no return.
Smilde, A., Jansen, J., Hoefsloot, H., Lamers,R., Van Der Greef, J., and Timmerman, M.(2005). ANOVA-Simultaneous Component Analysis (ASCA): A new tool for analyzing designed metabolomics data. Bioinformatics, 21(13), 3043–3048.
Liland, K.H., Smilde, A., Marini, F., and Næs,T. (2018). Confidence ellipsoids for ASCA models based on multivariate regression theory. Journal of Chemometrics, 32(e2990), 1–13.
Martin, M. and Govaerts, B. (2020). LiMM-PCA: Combining ASCA+ and linear mixed models to analyse high-dimensional designed data. Journal of Chemometrics, 34(6), e3232.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# For end-to-end examples using scoreplot()/loadingplot(), see # ?asca, ?apca, ?apls, ?limmpca and ?msca.# For end-to-end examples using scoreplot()/loadingplot(), see # ?asca, ?apca, ?apls, ?limmpca and ?msca.
Standard result computation and extraction functions for ASCA (asca).
## S3 method for class 'hdanova' print(x, ...) ## S3 method for class 'hdanova' summary(object, extended = TRUE, df = FALSE, ...) ## S3 method for class 'summary.hdanova' print(x, digits = 2, ...) ## S3 method for class 'asca' loadings(object, factor = 1, ...) ## S3 method for class 'asca' scores(object, factor = 1, ...) projections(object, ...) ## S3 method for class 'asca' projections(object, factor = 1, ...)## S3 method for class 'hdanova' print(x, ...) ## S3 method for class 'hdanova' summary(object, extended = TRUE, df = FALSE, ...) ## S3 method for class 'summary.hdanova' print(x, digits = 2, ...) ## S3 method for class 'asca' loadings(object, factor = 1, ...) ## S3 method for class 'asca' scores(object, factor = 1, ...) projections(object, ...) ## S3 method for class 'asca' projections(object, factor = 1, ...)
x |
|
... |
additional arguments to underlying methods. |
object |
|
extended |
Extended output in summary (default = TRUE). |
df |
Show degrees of freedom in summary (default = FALSE). |
digits |
|
factor |
|
Usage of the functions are shown using generics in the examples in asca.
Explained variances are available (block-wise and global) through blockexpl and print.rosaexpl.
Object printing and summary are available through:
print.asca and summary.asca.
Scores and loadings have their own extensions of scores() and loadings() through
scores.asca and loadings.asca. Special to ASCA is that scores are on a
factor level basis, while back-projected samples have their own function in projections.asca.
Returns depend on method used, e.g. projections.asca returns projected samples,
scores.asca return scores, while print and summary methods return the object invisibly.
Smilde, A., Jansen, J., Hoefsloot, H., Lamers,R., Van Der Greef, J., and Timmerman, M.(2005). ANOVA-Simultaneous Component Analysis (ASCA): A new tool for analyzing designed metabolomics data. Bioinformatics, 21(13), 3043–3048.
Liland, K.H., Smilde, A., Marini, F., and Næs,T. (2018). Confidence ellipsoids for ASCA models based on multivariate regression theory. Journal of Chemometrics, 32(e2990), 1–13.
Martin, M. and Govaerts, B. (2020). LiMM-PCA: Combining ASCA+ and linear mixed models to analyse high-dimensional designed data. Journal of Chemometrics, 34(6), e3232.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# For end-to-end examples using summary(), scores(), loadings(), and # projections(), see ?asca and ?apca.# For end-to-end examples using summary(), scores(), loadings(), and # projections(), see ?asca and ?apca.
Biplot for ASCA models
## S3 method for class 'asca' biplot( x, factor = 1, comps = 1:2, xlim = NULL, ylim = NULL, col = "darkgray", expand = 1, labels, legendpos, ... )## S3 method for class 'asca' biplot( x, factor = 1, comps = 1:2, xlim = NULL, ylim = NULL, col = "darkgray", expand = 1, labels, legendpos, ... )
x |
|
factor |
Factor number or name. |
comps |
|
xlim |
|
ylim |
|
col |
|
expand |
|
labels |
optional. If |
legendpos |
|
... |
Additional arguments to |
No return, only a plot.
# Load candies data data(candies) # Basic ASCA model with two factors and interaction mod <- asca(assessment ~ candy * assessor, data=candies) # Biplot biplot(mod) # Biplot with named loadings biplot(mod, labels="names")# Load candies data data(candies) # Basic ASCA model with two factors and interaction mod <- asca(assessment ~ candy * assessor, data=candies) # Biplot biplot(mod) # Biplot with named loadings biplot(mod, labels="names")
This is a convenience function for making data.frames that are easily
indexed on a block-wise basis.
block.data.frame(X, block_inds = NULL, to.matrix = TRUE)block.data.frame(X, block_inds = NULL, to.matrix = TRUE)
X |
Either a single |
block_inds |
Named |
to.matrix |
|
A data.frame which can be indexed block-wise.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# Random data M <- matrix(rnorm(200), nrow = 10) # .. with dimnames dimnames(M) <- list(LETTERS[1:10], as.character(1:20)) # A named list for indexing inds <- list(B1 = 1:10, B2 = 11:20) X <- block.data.frame(M, inds) str(X)# Random data M <- matrix(rnorm(200), nrow = 10) # .. with dimnames dimnames(M) <- list(LETTERS[1:10], as.character(1:20)) # A named list for indexing inds <- list(B1 = 1:10, B2 = 11:20) X <- block.data.frame(M, inds) str(X)
A dataset containing 67 metabolites from plants grown under different light and temperature conditions. This subset of the data contains only the light effect and time effect for limited conditions, while the full data also contains gene expressions.
data(caldana)data(caldana)
A data.frame having 140 rows and 3 variables:
Light levels
Time of measurement
Metabolic compounds
Caldana C, Degenkolbe T, Cuadros-Inostroza A, Klie S, Sulpice R, Leisse A, et al. High-density kinetic analysis of the metabolomic and transcriptomic response of Arabidopsis to eight environmental conditions. Plant J. 2011;67(5):869-884.
A dataset containing 9 sensory attributes for 5 candies assessed by 11 trained assessors.
data(candies)data(candies)
A data.frame having 165 rows and 3 variables:
Matrix of sensory attributes
Factor of assessors
Factor of candies
Luciano G, Næs T. Interpreting sensory data by combining principal component analysis and analysis of variance. Food Qual Prefer. 2009;20(3):167-175.
Flexible dummy-coding allowing for all R's built-in types of contrasts and optional dropping of a factor level to reduce rank defficiency probability.
dummycode(Y, contrast = "contr.sum", drop = TRUE)dummycode(Y, contrast = "contr.sum", drop = TRUE)
Y |
|
contrast |
Contrast type, default = "contr.sum". |
drop |
|
matrix made by dummy-coding the input vector.
vec <- c("a","a","b","b","c","c") dummycode(vec)vec <- c("a","a","b","b","c","c") dummycode(vec)
This function attempts to apply model.frame and extend the
result with columns of interactions.
extended.model.frame(formula, data, ..., sep = ".")extended.model.frame(formula, data, ..., sep = ".")
formula |
a model formula or terms object or an R object. |
data |
a data.frame, list or environment (see |
... |
further arguments to pass to |
sep |
separator in contraction of names for interactions (default = "."). |
A data.frame that includes everything a model.frame
does plus interaction terms.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
dat <- data.frame(Y = c(1,2,3,4,5,6), X = factor(LETTERS[c(1,1,2,2,3,3)]), W = factor(letters[c(1,2,1,2,1,2)])) extended.model.frame(Y ~ X*W, dat)dat <- data.frame(Y = c(1,2,3,4,5,6), X = factor(LETTERS[c(1,1,2,2,3,3)]), W = factor(letters[c(1,2,1,2,1,2)])) extended.model.frame(Y ~ X*W, dat)
Extracts and sums the LS estimates for a given factor combination
from an object of class hdanova. If add_residuals is TRUE,
the residuals are added to the LS estimates. If substract is TRUE,
the returned matrix is the data with chosen estimates subtracted.
extract_estimates(object, factors, subtract = FALSE, add_residuals = FALSE)extract_estimates(object, factors, subtract = FALSE, add_residuals = FALSE)
object |
|
factors |
|
subtract |
|
add_residuals |
|
A matrix of the extracted estimates.
Model fitting and related outputs: hdanova,
asca_results and predict.hdanova.
# Load candies data data(candies) # Basic HDANOVA model with two factors and interaction mod <- hdanova(assessment ~ candy * assessor, data=candies) # Extract estimates for the interaction inter <- extract_estimates(mod, c("candy:assessor")) # Visualize the interaction effect image(t(inter), main="Interaction effect", xlab="Attribute", ylab="Sample")# Load candies data data(candies) # Basic HDANOVA model with two factors and interaction mod <- hdanova(assessment ~ candy * assessor, data=candies) # Extract estimates for the interaction inter <- extract_estimates(mod, c("candy:assessor")) # Visualize the interaction effect image(t(inter), main="Interaction effect", xlab="Attribute", ylab="Sample")
This function provides a high-dimensional analysis of variance (HDANOVA) method
which can be used alone or as part of a larger analysis, e.g., ASCA, APCA, LiMM-PCA, MSCA or PC-ANOVA. It
can be called directly or through the convenience functions asca, apca,
limmpca, msca and pcanova.
hdanova( formula, data, subset, weights, na.action, family, scale = FALSE, add_error = FALSE, aug_error = "denominator", pca.in = FALSE, pls.in = FALSE, contrasts = "contr.sum", unrestricted = FALSE, SStype = "II", respect_SStype = FALSE, REML = NULL, REML_ssq_method = c("exact_refit", "wald", "ls"), equal_baseline = FALSE, use_ED = FALSE )hdanova( formula, data, subset, weights, na.action, family, scale = FALSE, add_error = FALSE, aug_error = "denominator", pca.in = FALSE, pls.in = FALSE, contrasts = "contr.sum", unrestricted = FALSE, SStype = "II", respect_SStype = FALSE, REML = NULL, REML_ssq_method = c("exact_refit", "wald", "ls"), equal_baseline = FALSE, use_ED = FALSE )
formula |
Model formula accepting a single response (block) and predictors. See Details for more information. |
data |
The data set to analyse. |
subset |
Expression for subsetting the data before modelling. |
weights |
Optional object weights. |
na.action |
How to handle NAs (no action implemented). |
family |
Error distributions and link function for Generalized Linear Models. |
scale |
Scaling of the response matrix. Defaults to |
add_error |
Add error to LS means, e.g., for APCA. |
aug_error |
Augment score matrices in backprojection. Default = "denominator" (of F test), "residual" (force error term), numeric value (alpha-value in LiMM-PCA). |
pca.in |
Compress response before ASCA (number of components). |
pls.in |
Compress response before ASCA using PLS (number of components). |
contrasts |
Effect coding: "contr.sum" (default), "contr.weighted" (not for lme4 models), "contr.reference", "contr.treatment". Can also be specified shorthand without "contr.", e.g., "sum", or as a named vector for different coding per factor. |
unrestricted |
Use unrestricted ANOVA decomposition (default = FALSE). |
SStype |
Type of sum-of-squares: "I" = sequential, "II" (default) = last term, obeying marginality, "III" = last term, not obeying marginality. |
respect_SStype |
Logical; if |
REML |
Parameter to mixlm: NULL (default) = sum-of-squares, TRUE = REML, FALSE = ML.
If supplied without any |
REML_ssq_method |
Method for REML mixed-model SSQ decomposition:
|
equal_baseline |
Experimental: Set to |
use_ED |
Use "effective dimensions" for score rescaling in LiMM-PCA. |
Scaling of the response matrix can be done by setting the scale parameter. If scale=TRUE,
each column is scaled by its standard deviation (autoscaling). A numeric value can be provided to scale
the columns by specific quantities. If scale is a character string, the first element
is interpreted as a factor name and the second element is interpreted as a factor level, whose samples
the standard deviations are calculated from (reference group scaling).
SStype and respect_SStype (fixed and MoM mixed models):
SStype controls how each effect's contribution is isolated from the others.
Type I (sequential) assigns variance in the order terms appear in the formula;
Type II (default) tests each term against all others that do not contain it, respecting
marginality; Type III tests each term against the full model, ignoring marginality.
For balanced designs the three types give identical results; differences arise in
unbalanced or non-orthogonal designs.
By default (respect_SStype = FALSE) the effect matrices in LS are regression
projections — each effect's columns are projected independently from the full model
coefficient matrix. This is consistent with the legacy ASCA workflow. When
respect_SStype = TRUE the LS matrices are re-derived from QR contrasts
that match the chosen SStype, so the Frobenius norm of LS[[eff]] equals
the corresponding ANOVA sum of squares. The regression matrices are still available in
LS_regression. For balanced, fully-crossed designs the two sets are numerically
identical; for unbalanced designs they will differ, particularly for Type II and III.
REML SSQ strategies (mixed models with r() terms):
When REML = TRUE or REML = FALSE, the model is fitted with lme4 and the
classical ANOVA table is replaced by one of three decomposition strategies controlled by
REML_ssq_method:
"exact_refit" (default): For each effect a reduced model is re-fitted
(REML or ML as specified) and the SSQ is the difference in log-likelihoods scaled
to the sum-of-squares metric. Concretely, if and are the
log-likelihoods of the full and reduced models, the statistic is
This is the most principled approach but requires one additional model fit per effect, making it slower for large models.
"wald": Uses the Wald chi-square statistic from the full model divided
by the residual variance to obtain an approximate SSQ. For a fixed effect with
coefficient vector and covariance
,
No additional model fits are required; fast but less accurate for small samples or near-singular random effects.
"ls": Projects the REML/ML coefficient estimates through the fixed-effect
design matrix to recover a least-squares-style SSQ,
where contains the columns of the design matrix for effect .
Computationally cheap and
numerically stable but ignores the mixed-model covariance structure; best treated
as an approximation.
Permutation statistics vs. fitted-model SSQ:
Permutation testing (permutation) always uses QR-based projection
statistics computed on the fixed-effect design matrix, regardless of REML_ssq_method.
When respect_SStype = FALSE (default) the permutation statistic is a
regression-projection norm; when respect_SStype = TRUE it is an SS-type-aligned
QR contrast norm. Neither is identical to the fitted-model REML SSQ values reported in
object$ssq, because REML/ML decompositions account for the random-effect
covariance structure whereas permutation statistics do not. The p-values are still
valid under their respective null hypotheses; the SSQ values should not be compared
across the two sources.
An hdanova object containing loadings, scores, explained variances, etc. The object has
associated plotting (asca_plots) and result (asca_results) functions.
# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) summary(mod)# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) summary(mod)
This function mimics parts of the LiMM-PCA framework, combining ASCA+ and linear mixed models to analyse high-dimensional designed data. The default is to use REML estimation and scaling of the backprojected errors. See examples for alternatives.
limmpca( formula, data, pca.in = 5, aug_error = 0.05, use_ED = FALSE, REML = TRUE, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), SStype = "III", ... )limmpca( formula, data, pca.in = 5, aug_error = 0.05, use_ED = FALSE, REML = TRUE, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), SStype = "III", ... )
formula |
Model formula accepting a single response (block) and predictors. See Details for more information. |
data |
The data set to analyse. |
pca.in |
Compress response before ASCA (number of components), default = 5. |
aug_error |
Error term of model ("denominator", "residual", numeric alpha-value). The latter implies the first with a scaling factor. |
use_ED |
Use Effective Dimensions instead of degrees of freedom when scaling. |
REML |
Use restricted maximum likelihood estimation. Alternatives: TRUE (default), FALSE (ML), NULL (least squares). |
contrasts |
Effect coding: "sum" (default = sum-coding), "weighted", "reference", "treatment". |
permute |
Number of permutations to perform (default = 1000). |
perm.type |
Type of permutation to perform, either "approximate" or "exact" (default = "approximate"). |
SStype |
Type of sum-of-squares: "I" = sequential, "II" = last term, obeying marginality, "III" (default) = last term, not obeying marginality. |
... |
Additional arguments to |
The Sum of Squares for the model is dependent on the SStype of the model. For SStype = "I" and SStype = "II" the SSQ is based on LLR (possibly inflating large contributions), while it is directly estimated from the model for SStype = "III". SStype = "III" is the default for LiMM-PCA and should be combined with sum coding. Sum of Squares for the random effects are based on the variance components.
An object of class limmpca, inheriting from the general asca class.
Martin, M. and Govaerts, B. (2020). LiMM-PCA: Combining ASCA+ and linear mixed models to analyse high-dimensional designed data. Journal of Chemometrics, 34(6), e3232.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
Prediction on new data: predict_asca_family.
# Load candies data data(candies) # Default LiMM-PCA model with two factors and interaction, 5 PCA components mod <- limmpca(assessment ~ candy*r(assessor), data=candies) summary(mod) scoreplot(mod, factor = "candy") # LiMM-PCA with least squares estimation and 8 PCA components modLS <- limmpca(assessment ~ candy*r(assessor), data=candies, REML=NULL, pca.in=8) summary(modLS) scoreplot(modLS, factor = "candy") # Load Caldana data data(caldana) # Combining effects in LiMM-PCA (assuming light is a random factor) mod.comb <- limmpca(compounds ~ time + comb(r(light) + r(time:light)), data=caldana, pca.in=8) summary(mod.comb)# Load candies data data(candies) # Default LiMM-PCA model with two factors and interaction, 5 PCA components mod <- limmpca(assessment ~ candy*r(assessor), data=candies) summary(mod) scoreplot(mod, factor = "candy") # LiMM-PCA with least squares estimation and 8 PCA components modLS <- limmpca(assessment ~ candy*r(assessor), data=candies, REML=NULL, pca.in=8) summary(modLS) scoreplot(modLS, factor = "candy") # Load Caldana data data(caldana) # Combining effects in LiMM-PCA (assuming light is a random factor) mod.comb <- limmpca(compounds ~ time + comb(r(light) + r(time:light)), data=caldana, pca.in=8) summary(mod.comb)
Extract log-likelihood from fitted HDANOVA objects. For fixed-effects
and MoM ('r()' with REML = NULL) models, the log-likelihood is computed
from the residual sum of squares. For REML/ML mixed models ('r()' with
REML = TRUE or FALSE), the log-likelihood is extracted directly
from the stored lme4 model objects.
## S3 method for class 'hdanova' logLik(object, ...)## S3 method for class 'hdanova' logLik(object, ...)
object |
A fitted |
... |
Reserved for generic compatibility. |
The log-likelihood is computed as follows:
Fixed-effects and MoM models: The multivariate Gaussian log-likelihood is
where is the residual for observation and response ,
and is the ML variance estimate for response
column . The total log-likelihood is .
REML/ML mixed models: The log-likelihood is the sum of the individual response-specific log-likelihoods from the fitted lme4 models,
where is the log-likelihood of the lme4::lmerMod fit for response .
An object of class logLik with the computed log-likelihood value,
degrees of freedom ('df'), and the number of observations ('nobs').
Model fitting: hdanova.
Information criteria: AIC and BIC.
data(candies) mod <- hdanova(assessment ~ candy + assessor, data = candies) ll <- logLik(mod) print(ll) ## Not run: # For mixed models: mod_reml <- hdanova(assessment ~ candy + r(assessor), data = candies, REML = TRUE) ll_reml <- logLik(mod_reml) print(ll_reml) ## End(Not run)data(candies) mod <- hdanova(assessment ~ candy + assessor, data = candies) ll <- logLik(mod) print(ll) ## Not run: # For mixed models: mod_reml <- hdanova(assessment ~ candy + r(assessor), data = candies, REML = TRUE) ll_reml <- logLik(mod_reml) print(ll_reml) ## End(Not run)
Extraction functions to retrieve the model.frame and model.matrix of an asca object.
## S3 method for class 'asca' model.frame(formula, ...) ## S3 method for class 'asca' model.matrix(object, ...)## S3 method for class 'asca' model.frame(formula, ...) ## S3 method for class 'asca' model.matrix(object, ...)
formula |
The |
... |
Not implemented |
object |
The |
A data.frame or matrix object.
Core model constructors: asca, apca,
apls, limmpca, msca and
pcanova.
# Load candies data data(candies) # Basic ASCA model with two factors mod <- asca(assessment ~ candy + assessor, data=candies) # Extract model frame and model matrix mf <- model.frame(mod) head(mf) mm <- model.matrix(mod) par.old <- par(mar=c(3,3,3,1), mgp=c(1,0.7,0)) image(t(mm[seq(165,1,-1),]), main="Model Matrix", xlab="dummy values", ylab="samples", axes=FALSE) par(par.old)# Load candies data data(candies) # Basic ASCA model with two factors mod <- asca(assessment ~ candy + assessor, data=candies) # Extract model frame and model matrix mf <- model.frame(mod) head(mf) mm <- model.matrix(mod) par.old <- par(mar=c(3,3,3,1), mgp=c(1,0.7,0)) image(t(mm[seq(165,1,-1),]), main="Model Matrix", xlab="dummy values", ylab="samples", axes=FALSE) par(par.old)
This MSCA implementation assumes a single factor to be used as between-individuals factor.
msca( formula, data, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )msca( formula, data, contrasts = "contr.sum", permute = FALSE, perm.type = c("approximate", "exact"), ... )
formula |
Model formula accepting a single response (block) and predictors. See Details for more information. |
data |
The data set to analyse. |
contrasts |
Effect coding: "sum" (default = sum-coding), "weighted", "reference", "treatment". |
permute |
Number of permutations to perform (default = 1000). |
perm.type |
Type of permutation to perform, either "approximate" or "exact" (default = "approximate"). |
... |
Additional arguments to |
An asca object containing loadings, scores, explained variances, etc. The object has
associated plotting (asca_plots) and result (asca_results) functions.
Smilde, A., Jansen, J., Hoefsloot, H., Lamers,R., Van Der Greef, J., and Timmerman, M.(2005). ANOVA-Simultaneous Component Analysis (ASCA): A new tool for analyzing designed metabolomics data. Bioinformatics, 21(13), 3043–3048.
Liland, K.H., Smilde, A., Marini, F., and Næs,T. (2018). Confidence ellipsoids for ASCA models based on multivariate regression theory. Journal of Chemometrics, 32(e2990), 1–13.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# Load candies data data(candies) # Basic MSCA model with a single factor mod <- msca(assessment ~ candy, data=candies) print(mod) summary(mod) # Result plotting for first factor loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # Within scores scoreplot(mod, factor="within") # Within scores per factor level par.old <- par(mfrow=c(3,2), mar=c(4,4,2,1), mgp=c(2,0.7,0)) for(i in 1:length(mod$scores.within)) scoreplot(mod, factor="within", within_level=i, main=paste0("Level: ", names(mod$scores.within)[i]), panel.first=abline(v=0,h=0,col="gray",lty=2)) par(par.old) # Permutation testing mod.perm <- asca(assessment ~ candy * assessor, data=candies, permute=TRUE) summary(mod.perm)# Load candies data data(candies) # Basic MSCA model with a single factor mod <- msca(assessment ~ candy, data=candies) print(mod) summary(mod) # Result plotting for first factor loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # Within scores scoreplot(mod, factor="within") # Within scores per factor level par.old <- par(mfrow=c(3,2), mar=c(4,4,2,1), mgp=c(2,0.7,0)) for(i in 1:length(mod$scores.within)) scoreplot(mod, factor="within", within_level=i, main=paste0("Level: ", names(mod$scores.within)[i]), panel.first=abline(v=0,h=0,col="gray",lty=2)) par(par.old) # Permutation testing mod.perm <- asca(assessment ~ candy * assessor, data=candies, permute=TRUE) summary(mod.perm)
This is a quite general and flexible implementation of PC-ANOVA.
pcanova(formula, data, ncomp = 0.9, contrasts = "contr.sum", ...)pcanova(formula, data, ncomp = 0.9, contrasts = "contr.sum", ...)
formula |
Model formula accepting a single response (block) and predictor names separated by + signs. |
data |
The data set to analyse. |
ncomp |
The number of components to retain, proportion of variation or default = minimum cross-validation error. |
contrasts |
Effect coding: "sum" (default = sum-coding), "weighted", "reference", "treatment". |
... |
Additional parameters for the |
PC-ANOVA works in the opposite order of ASCA. First the response matrix is decomposed using ANOVA. Then the components are analysed using ANOVA with respect to a design or grouping in the data. The latter can be ordinary fixed effects modelling or mixed models.
A pcanova object containing loadings, scores, explained variances, etc. The object has
associated plotting (pcanova_plots) and result (pcanova_results) functions.
Luciano G, Næs T. Interpreting sensory data by combining principal component analysis and analysis of variance. Food Qual Prefer. 2009;20(3):167-175.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# Load candies data data(candies) # Basic PC-ANOVA model with two factors, cross-validated opt. of #components mod <- pcanova(assessment ~ candy + assessor, data = candies) print(mod) # PC-ANOVA model with interaction, minimum 90% explained variance mod <- pcanova(assessment ~ candy * assessor, data = candies, ncomp = 0.9) print(mod) summary(mod) # Tukey group letters for 'candy' per component lapply(mod$models, function(x) mixlm::cld(mixlm::simple.glht(x, effect = "candy"))) # Result plotting loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # Mixed Model PC-ANOVA, random assessor mod.mix <- pcanova(assessment ~ candy + r(assessor), data=candies, ncomp = 0.9) scoreplot(mod.mix) # Fixed effects summary(mod.mix)# Load candies data data(candies) # Basic PC-ANOVA model with two factors, cross-validated opt. of #components mod <- pcanova(assessment ~ candy + assessor, data = candies) print(mod) # PC-ANOVA model with interaction, minimum 90% explained variance mod <- pcanova(assessment ~ candy * assessor, data = candies, ncomp = 0.9) print(mod) summary(mod) # Tukey group letters for 'candy' per component lapply(mod$models, function(x) mixlm::cld(mixlm::simple.glht(x, effect = "candy"))) # Result plotting loadingplot(mod, scatter=TRUE, labels="names") scoreplot(mod) # Mixed Model PC-ANOVA, random assessor mod.mix <- pcanova(assessment ~ candy + r(assessor), data=candies, ncomp = 0.9) scoreplot(mod.mix) # Fixed effects summary(mod.mix)
Various plotting procedures for pcanova objects.
## S3 method for class 'pcanova' scoreplot(object, factor = 1, comps = 1:2, col = "factor", ...)## S3 method for class 'pcanova' scoreplot(object, factor = 1, comps = 1:2, col = "factor", ...)
object |
|
factor |
|
comps |
|
col |
|
... |
additional arguments to underlying methods. |
Usage of the functions are shown using generics in the examples in pcanova.
Plot routines are available as
scoreplot.pcanova and loadingplot.pcanova.
The plotting routines have no return.
Luciano G, Næs T. Interpreting sensory data by combining principal component analysis and analysis of variance. Food Qual Prefer. 2009;20(3):167-175.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# For end-to-end examples using scoreplot()/loadingplot() on PC-ANOVA # objects, see ?pcanova.# For end-to-end examples using scoreplot()/loadingplot() on PC-ANOVA # objects, see ?pcanova.
Standard result computation and extraction functions for ASCA (pcanova).
## S3 method for class 'pcanova' summary(object, ...) ## S3 method for class 'summary.pcanova' print(x, digits = 2, ...) ## S3 method for class 'pcanova' print(x, ...) ## S3 method for class 'pcanova' summary(object, ...)## S3 method for class 'pcanova' summary(object, ...) ## S3 method for class 'summary.pcanova' print(x, digits = 2, ...) ## S3 method for class 'pcanova' print(x, ...) ## S3 method for class 'pcanova' summary(object, ...)
object |
|
... |
additional arguments to underlying methods. |
x |
|
digits |
|
Usage of the functions are shown using generics in the examples in pcanova.
Explained variances are available (block-wise and global) through blockexpl and print.rosaexpl.
Object printing and summary are available through:
print.pcanova and summary.pcanova.
Scores and loadings have their own extensions of scores() and loadings() through
scores.pcanova and loadings.pcanova. Special to ASCA is that scores are on a
factor level basis, while back-projected samples have their own function in projections.pcanova.
Returns depend on method used, e.g. projections.pcanova returns projected samples,
scores.pcanova return scores, while print and summary methods return the object invisibly.
Luciano G, Næs T. Interpreting sensory data by combining principal component analysis and analysis of variance. Food Qual Prefer. 2009;20(3):167-175.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
# For examples using print()/summary() for PC-ANOVA objects, see # ?pcanova.# For examples using print()/summary() for PC-ANOVA objects, see # ?pcanova.
Wrapper for the adonis2 function to allow ordinary formula input.
permanova(formula, data, ...)permanova(formula, data, ...)
formula |
Model formula accepting a single response matrix and predictors.
See details in |
data |
The data set to analyse. |
... |
Additional arguments to |
An ANOVA table with permutation-based p-values.
Related wrappers and model families: prc,
asca, apca, apls,
limmpca, msca and pcanova.
data(caldana) (pr <- permanova(compounds ~ light * time, caldana))data(caldana) (pr <- permanova(compounds ~ light * time, caldana))
Permutation testing for hdanova() objects using SS-type-aligned QR
effect matrices. This function is intended as an opt-in alternative to the
legacy regression-based permutation() workflow.
permutation( object, permute = 1000, perm.type = c("approximate", "exact"), respect_SStype = NULL, unique.digits = 12, unique.frac = 0.95, exhaustive.warn = TRUE )permutation( object, permute = 1000, perm.type = c("approximate", "exact"), respect_SStype = NULL, unique.digits = 12, unique.frac = 0.95, exhaustive.warn = TRUE )
object |
A |
permute |
Number of permutations to perform (default = 1000). |
perm.type |
Type of permutation to perform, either |
respect_SStype |
Logical or |
unique.digits |
Number of digits used when rounding permutation SSQ values before
checking uniqueness (default = 12). Set to |
unique.frac |
Minimum fraction of unique rounded SSQ values required to avoid
warning (default = 0.95). Set to |
exhaustive.warn |
Logical; if |
The permutation statistics are computed from SS-type-aligned QR reduced/full
model contrasts rather than the legacy regression-based LS matrices. Fixed
models and mixed MoM models are supported. Approximate permutation uses a
relaxed global shuffle of observations; exact permutation uses permissible
block-restricted shuffles. For REML/ML mixed models with
respect_SStype = TRUE, a warning is issued to highlight that
permutation statistics and fitted-model REML SSQ decompositions are based on
different computational definitions.
An updated hdanova object with permute results.
This function performs Partial Least Squares (PLS) on a hdanova.
pls(object, ...) ## Default S3 method: pls(object, ...)pls(object, ...) ## Default S3 method: pls(object, ...)
object |
A |
... |
Additional arguments (not used). |
For residuals, PCA is performed instead of PLS as there is no natural response.
An updated hdanova object with PLS results.
Main wrapper: apls.
Related decomposition: sca.
Plotting and summaries: asca_plots and asca_results.
# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) mod <- pls(mod) scoreplot(mod)# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) mod <- pls(mod) scoreplot(mod)
Wrapper for the prc function to allow for formula input.
prc(formula, data, ...)prc(formula, data, ...)
formula |
Model formula accepting a single response (block) and predictors. If no predictor is called 'time', time is assumed to be the second predictor. |
data |
The data set to analyse. |
... |
Additional arguments to |
An object of class prc.
Main methods: asca, apca, limmpca, msca, pcanova, prc and permanova.
Workhorse function underpinning most methods: hdanova.
Extraction of results and plotting: asca_results, asca_plots, pcanova_results and pcanova_plots
data(caldana) (pr <- prc(compounds ~ light * time, caldana)) summary(pr)data(caldana) (pr <- prc(compounds ~ light * time, caldana)) summary(pr)
Reconstructs the HDANOVA base on newdata via predict.hdanova()
and then applies the class-specific decomposition step:
sca() for ASCA/APCA/MSCA and pls() for APLS.
By default, decomposition is done by projection onto training component
spaces. A refit mode is also available.
## S3 method for class 'asca' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'apca' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'msca' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'apls' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'limmpca' predict(object, newdata, decomposition = c("project", "refit"), ...)## S3 method for class 'asca' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'apca' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'msca' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'apls' predict(object, newdata, decomposition = c("project", "refit"), ...) ## S3 method for class 'limmpca' predict(object, newdata, decomposition = c("project", "refit"), ...)
object |
A fitted |
newdata |
A data frame containing variables from the original model formula. |
decomposition |
Decomposition mode: |
... |
Reserved for compatibility; forwarded to |
A predicted object of the same high-level class as object.
Base prediction engine: predict.hdanova.
Related model constructors: asca, apca,
apls, msca and limmpca.
data(candies) test_idx <- seq(3, nrow(candies), by = 3) train_idx <- setdiff(seq_len(nrow(candies)), test_idx) candies_train <- candies[train_idx, ] candies_test <- candies[test_idx, ] mod_asca <- asca(assessment ~ candy * assessor, data = candies_train) pred_asca <- predict(mod_asca, newdata = candies_test) scoreplot(mod_asca, factor="candy", legend=TRUE) with(pred_asca$projected, points(candy[,1], candy[,2], pch="x", cex=0.8, col=as.numeric(pred_asca$model.frame$candy))) pred_asca_refit <- predict(mod_asca, newdata = candies_test, decomposition = "refit") mod_apca <- apca(assessment ~ candy + assessor, data = candies_train) pred_apca <- predict(mod_apca, newdata = candies_test) mod_msca <- msca(assessment ~ candy, data = candies_train) pred_msca <- predict(mod_msca, newdata = candies_test) mod_apls <- apls(assessment ~ candy + assessor, data = candies_train) pred_apls <- predict(mod_apls, newdata = candies_test) mod_limmpca <- limmpca(assessment ~ candy + r(assessor), data = candies_train, pca.in = 3) pred_limmpca <- predict(mod_limmpca, newdata = candies_test)data(candies) test_idx <- seq(3, nrow(candies), by = 3) train_idx <- setdiff(seq_len(nrow(candies)), test_idx) candies_train <- candies[train_idx, ] candies_test <- candies[test_idx, ] mod_asca <- asca(assessment ~ candy * assessor, data = candies_train) pred_asca <- predict(mod_asca, newdata = candies_test) scoreplot(mod_asca, factor="candy", legend=TRUE) with(pred_asca$projected, points(candy[,1], candy[,2], pch="x", cex=0.8, col=as.numeric(pred_asca$model.frame$candy))) pred_asca_refit <- predict(mod_asca, newdata = candies_test, decomposition = "refit") mod_apca <- apca(assessment ~ candy + assessor, data = candies_train) pred_apca <- predict(mod_apca, newdata = candies_test) mod_msca <- msca(assessment ~ candy, data = candies_train) pred_msca <- predict(mod_msca, newdata = candies_test) mod_apls <- apls(assessment ~ candy + assessor, data = candies_train) pred_apls <- predict(mod_apls, newdata = candies_test) mod_limmpca <- limmpca(assessment ~ candy + r(assessor), data = candies_train, pca.in = 3) pred_limmpca <- predict(mod_limmpca, newdata = candies_test)
Reconstructs an HDANOVA-style object on newdata without refitting by
reusing stored coefficients and projection objects from the fitted model.
This implementation supports fixed-effects, mixed MoM (r() with
REML = NULL), and REML/ML mixed workflows (r() with
REML = TRUE or FALSE).
## S3 method for class 'hdanova' predict(object, newdata, ...)## S3 method for class 'hdanova' predict(object, newdata, ...)
object |
A fitted |
newdata |
A data frame containing the variables used in the original model formula. |
... |
Reserved for generic compatibility; runtime overrides are not supported. |
An hdanova-family object computed on newdata.
ASCA-family prediction wrappers: predict_asca_family.
Model constructors: hdanova, asca,
apca, apls, limmpca,
msca and pcanova.
data(candies) # Train/test split (every third sample to test) test_idx <- seq(3, nrow(candies), by = 3) train_idx <- setdiff(1:nrow(candies),test_idx) candies_train <- candies[train_idx, ] candies_test <- candies[test_idx, ] # Fixed-effects model prediction mod <- hdanova(assessment ~ candy + assessor, data = candies_train) pred <- predict(mod, newdata = candies_test) var_idx <- seq_len(ncol(mod$LS$candy)) old.par <- par(mfrow = c(1,2), mar = c(4,4,2,1), mgp = c(2,0.7,0)) image(x = var_idx, y = seq_along(train_idx), z = t(mod$LS$candy), xaxt = "n", yaxt = "n", main = "Original candy LS", xlab = "Variable index", ylab = "Train sample index") axis(1, at = var_idx, labels = var_idx) axis(2, at = seq_along(train_idx), labels = train_idx) image(x = var_idx, y = seq_along(test_idx), z = t(pred$LS$candy), xaxt = "n", yaxt = "n", main = "Predicted candy LS", xlab = "Variable index", ylab = "Test sample index") axis(1, at = var_idx, labels = var_idx) axis(2, at = seq_along(test_idx), labels = test_idx) par(old.par) # Mixed MoM model prediction (r() with REML = NULL) mod_mom <- hdanova(assessment ~ candy + r(assessor), data = candies_train) pred_mom <- predict(mod_mom, newdata = candies_test) cat("Mixed MoM model prediction successful.\n") cat("SSQ names:", paste(names(pred_mom$ssq), collapse = "|"), "\n") cat("dfDenom:", paste(pred_mom$dfDenom, collapse = "|"), "\n") # REML mixed model prediction (r() with REML = TRUE) mod_reml <- hdanova(assessment ~ candy + r(assessor), data = candies_train, REML = TRUE) pred_reml <- predict(mod_reml, newdata = candies_test) cat("REML mixed model prediction successful.\n") cat("SSQ names:", paste(names(pred_reml$ssq), collapse = "|"), "\n") cat("dfDenom:", paste(pred_reml$dfDenom, collapse = "|"), "\n")data(candies) # Train/test split (every third sample to test) test_idx <- seq(3, nrow(candies), by = 3) train_idx <- setdiff(1:nrow(candies),test_idx) candies_train <- candies[train_idx, ] candies_test <- candies[test_idx, ] # Fixed-effects model prediction mod <- hdanova(assessment ~ candy + assessor, data = candies_train) pred <- predict(mod, newdata = candies_test) var_idx <- seq_len(ncol(mod$LS$candy)) old.par <- par(mfrow = c(1,2), mar = c(4,4,2,1), mgp = c(2,0.7,0)) image(x = var_idx, y = seq_along(train_idx), z = t(mod$LS$candy), xaxt = "n", yaxt = "n", main = "Original candy LS", xlab = "Variable index", ylab = "Train sample index") axis(1, at = var_idx, labels = var_idx) axis(2, at = seq_along(train_idx), labels = train_idx) image(x = var_idx, y = seq_along(test_idx), z = t(pred$LS$candy), xaxt = "n", yaxt = "n", main = "Predicted candy LS", xlab = "Variable index", ylab = "Test sample index") axis(1, at = var_idx, labels = var_idx) axis(2, at = seq_along(test_idx), labels = test_idx) par(old.par) # Mixed MoM model prediction (r() with REML = NULL) mod_mom <- hdanova(assessment ~ candy + r(assessor), data = candies_train) pred_mom <- predict(mod_mom, newdata = candies_test) cat("Mixed MoM model prediction successful.\n") cat("SSQ names:", paste(names(pred_mom$ssq), collapse = "|"), "\n") cat("dfDenom:", paste(pred_mom$dfDenom, collapse = "|"), "\n") # REML mixed model prediction (r() with REML = TRUE) mod_reml <- hdanova(assessment ~ candy + r(assessor), data = candies_train, REML = TRUE) pred_reml <- predict(mod_reml, newdata = candies_test) cat("REML mixed model prediction successful.\n") cat("SSQ names:", paste(names(pred_reml$ssq), collapse = "|"), "\n") cat("dfDenom:", paste(pred_reml$dfDenom, collapse = "|"), "\n")
Rotation testing for HDANOVA. This function performs
random orthogonal rotations of exchangeable residual units for each approved
effect and adds the resulting null distributions to the hdanova object.
rotation( object, rotate = 1000, unique.digits = 12, unique.frac = 0.95, block.type = c("denominator", "global") )rotation( object, rotate = 1000, unique.digits = 12, unique.frac = 0.95, block.type = c("denominator", "global") )
object |
A |
rotate |
Number of random rotations to perform (default = 1000). |
unique.digits |
Number of digits used when rounding rotation SSQ values before
checking uniqueness (default = 12). Set to |
unique.frac |
Minimum fraction of unique rounded SSQ values required to avoid
warning (default = 0.95). Set to |
block.type |
Rotation blocking strategy. |
An updated hdanova object with rotation-test results stored in
object$permute for compatibility with existing summary and plotting tools.
Base fitting: hdanova.
Permutation alternative: permutation.
Plot helper: rotationplot.
# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) # Rotation test modRot <- rotation(mod) summary(modRot) # Plot null distribution for "candy" effect rotationplot(modRot, factor="candy")# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) # Rotation test modRot <- rotation(mod) summary(modRot) # Plot null distribution for "candy" effect rotationplot(modRot, factor="candy")
This function performs Simultaneous Component Analysis (SCA) on a hdanova object.
sca(object)sca(object)
object |
A |
An updated hdanova object with SCA results.
Model constructors and wrappers: hdanova,
asca, apca, limmpca and
msca.
Plotting and summaries: asca_plots and asca_results.
# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) mod <- sca(mod) scoreplot(mod)# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) mod <- sca(mod) scoreplot(mod)
This function flips the sign of a selected component in a selected factor of an asca object.
This affects both scores, loadings and projected data.
signflip(object, factor, comp)signflip(object, factor, comp)
object |
|
factor |
|
comp |
|
An asca object with the sign of the selected component flipped.
# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) mod <- sca(mod) old.par <- par(mfrow=c(1,2), mar=c(4,4,1,1)) scoreplot(mod, factor="candy") loadingplot(mod, factor="candy") par(old.par) # Flip the sign of the first component of the candy factor mod <- signflip(mod, factor="candy", comp=1) old.par <- par(mfrow=c(1,2), mar=c(4,4,1,1)) scoreplot(mod, factor="candy") loadingplot(mod, factor="candy") par(old.par)# Load candies data data(candies) # Basic HDANOVA model with two factors mod <- hdanova(assessment ~ candy + assessor, data=candies) mod <- sca(mod) old.par <- par(mfrow=c(1,2), mar=c(4,4,1,1)) scoreplot(mod, factor="candy") loadingplot(mod, factor="candy") par(old.par) # Flip the sign of the first component of the candy factor mod <- signflip(mod, factor="candy", comp=1) old.par <- par(mfrow=c(1,2), mar=c(4,4,1,1)) scoreplot(mod, factor="candy") loadingplot(mod, factor="candy") par(old.par)
Timeplot for Combined Effects
timeplot( object, factor, time, comb, comp = 1, ylim, x_time = FALSE, xlab = time, ylab = paste0("Score ", comp), lwd = 2, ... )timeplot( object, factor, time, comb, comp = 1, ylim, x_time = FALSE, xlab = time, ylab = paste0("Score ", comp), lwd = 2, ... )
object |
|
factor |
|
time |
|
comb |
|
comp |
|
ylim |
|
x_time |
|
xlab |
|
ylab |
|
lwd |
|
... |
additional arguments to |
Nothing
data("caldana") mod.comb <- asca(compounds ~ time + comb(light + time:light), data=caldana) # Default time axis timeplot(mod.comb, factor="light", time="time", comb=2) # Non-equispaced time axis (using time levels) timeplot(mod.comb, factor="light", time="time", comb=2, x_time=TRUE) # Second component timeplot(mod.comb, factor="light", time="time", comb=2, comp=2, x_time=TRUE)data("caldana") mod.comb <- asca(compounds ~ time + comb(light + time:light), data=caldana) # Default time axis timeplot(mod.comb, factor="light", time="time", comb=2) # Non-equispaced time axis (using time levels) timeplot(mod.comb, factor="light", time="time", comb=2, x_time=TRUE) # Second component timeplot(mod.comb, factor="light", time="time", comb=2, comp=2, x_time=TRUE)
Perform a model update while removing a chosen factor. Hierarchical corresponds to type "II" sum-of-squares, i.e., obeying marginality, while non-hierarchical corresponds to type "III" sum-of-squares.
update_without_factor(model, fac, hierarchical = TRUE)update_without_factor(model, fac, hierarchical = TRUE)
model |
|
fac |
|
hierarchical |
|
An updated model object is returned. If the supplied model is of
type lmerMod and no random effects are left, the model is
automatically converted to a linear model before updating.