Package 'RRPP'

Description

Function adds trajectories to a principal component plot

Usage

add.trajectories(
  TP,
  traj.pch = 21,
  traj.col = 1,
  traj.lty = 1,
  traj.lwd = 1,
  traj.cex = 1.5,
  traj.bg = 1,
  start.bg = 3,
  end.bg = 2
)

Arguments

Details

The function adds trajectories to a plot made by plot.trajectory.analysis. This function has a restricted set of plot parameters based on the number of trajectories to be added to the plot.

Author(s)

Michael Collyer

References

Adams, D. C., and M. M. Cerney. 2007. Quantifying biomechanical motion using Procrustes motion analysis. J. Biomech. 40:437-444.

Adams, D. C., and M. L. Collyer. 2007. The analysis of character divergence along environmental gradients and other covariates. Evolution 61:510-515.

Adams, D. C., and M. L. Collyer. 2009. A general framework for the analysis of phenotypic trajectories in evolutionary studies. Evolution 63:1143-1154.

Collyer, M. L., and D. C. Adams. 2007. Analysis of two-state multivariate phenotypic change in ecological studies. Ecology 88:683-692.

Collyer, M. L., and D. C. Adams. 2013. Phenotypic trajectory analysis: comparison of shape change patterns in evolution and ecology. Hystrix 24: 75-83.

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 115:357-365.

See Also

plot.default and par

Description

Function adds a tree based on a description of edges from a class phylo object to an existing plot made from an ordinate object.

Usage

add.tree(
  OP,
  tree,
  edge.col = 1,
  edge.lty = 1,
  edge.lwd = 1,
  anc.pts = FALSE,
  return.ancs = FALSE,
  ...
)

Arguments

Details

With some ordinate plots, it might be desirable to add a tree connecting points in a prescribed way, which would be tedious using points or lines. This function will project a tree from an object of class phylo into a plot with class, plot.ordinate. Using an edges matrix from a phylo object, this function will systematically connect plot points with lines that pass through estimated ancestral character points in the same plot space. Ancestral states are estimated assuming a Brownian motion model of evolutionary divergence.

Author(s)

Michael Collyer

See Also

lines and points

Examples

# Examples use residuals from a regression of salamander morphological 
# traits against body size (snout to vent length, SVL).
# Observations are species means and a phylogenetic covariance matrix
# describes the relatedness among observations.

data("PlethMorph")
Y <- as.data.frame(PlethMorph[c("TailLength", "HeadLength", 
"Snout.eye", "BodyWidth", 
"Forelimb", "Hindlimb")])
Y <- as.matrix(Y)
R <- lm.rrpp(Y ~ SVL, data = PlethMorph, 
iter = 0, print.progress = FALSE)$LM$residuals

PCA <- ordinate(R, scale. = TRUE)
pc.plot <- plot(PCA, pch = 19, col = "blue")

add.tree(pc.plot, tree = PlethMorph$tree, anc.pts = TRUE, 
pch = 19, cex = 0.5, col = "red")

Description

Computes an analysis of variance (ANOVA) table using distributions of random statistics from lm.rrpp. ANOVA can be performed on one model or multiple models. If the latter, the first model is considered a null model for comparison to other models. The ANOVA is functionally similar to a non-parametric likelihood ratio test for all null-full model comparisons Residuals from the null model will be used to generate random pseudo-values via RRPP for evaluation of subsequent models. The permutation schedule from the null model will be used for random permutations. This function does not correct for improper null models. One must assure that the null model is nested within the other models. Illogical results can be generated if this is not the case.

Usage

## S3 method for class 'lm.rrpp'
anova(
  object,
  ...,
  effect.type = c("F", "cohenf", "SS", "MS", "Rsq", "LR"),
  error = NULL,
  print.progress = TRUE
)

Arguments

Author(s)

Michael Collyer

Examples

## Not run: 
# See examples for lm.rrpp to see how anova.lm.rrpp works in conjunction
# with other functions

data(Pupfish)
names(Pupfish)
Pupfish$logSize <- log(Pupfish$CS) # better to not have functions in formulas

# Single-Model ANOVA

fit <- lm.rrpp(coords ~ logSize + Sex*Pop, SS.type = "I", 
data = Pupfish, print.progress = FALSE, iter = 999) 
anova(fit)
anova(fit, effect.type = "MS")
anova(fit, effect.type = "Rsq")
anova(fit, effect.type = "cohenf")

# Multi-Model ANOVA (like a Likelihood Ratio Test)
fit.size <- lm.rrpp(coords ~ logSize, SS.type = "I", data = Pupfish, 
print.progress = FALSE, iter = 999) 
fit.sex <- lm.rrpp(coords ~ logSize + Sex, SS.type = "I", data = Pupfish, 
print.progress = FALSE, iter = 999) 
fit.pop <- lm.rrpp(coords ~ logSize + Pop, SS.type = "I", data = Pupfish, 
print.progress = FALSE, iter = 999) 
anova(fit.size, fit.sex, fit.pop, 
print.progress = FALSE) # compares two models to the first

# see lm.rrpp examples for mixed model ANOVA example and how to vary SS type

## End(Not run)

Description

Computes an analysis of variance (ANOVA) table using distributions of random statistics from lm.rrpp. This function is the same as anova.lm.rrpp but includes statistics specific to measurement.error, and with restrictions on how P-values and effect.sizes are calculated.

Usage

## S3 method for class 'measurement.error'
anova(object, ...)

Arguments

Author(s)

Michael Collyer

Examples

# See measurement.error help file examples for use.

Description

For any lm.rrpp object, a vector of coefficients can be used for a specific test of a vector of betas (specific population parameters). This test follows the form of (b - beta) in the numerator of a t-statistic, where beta can be a value other than 0 (or 0). However, for this test, a vector (Beta) of length, p, is used for the p variables in the lm.rrpp fit. If Beta is a vector of 0s, this test is essentially the same as the test performed for coef.lm.rrpp. However, it is possible to test null hypotheses for beta values other than 0, sensu Cicuéndez et al. (2023).

This function can use either the square-root of the inner-product of vectors of coefficients (distance, d) or generalized inner-product based on the inverse of the residual covariance matrix (Mahalanobis distance, md) as statistics. In most cases, either will likely yield similar (or same) P-values. However, Mahalanobis distance might be preferred for generalized least squares fits, which do not have consistent residual covariance matrices for null (intercept only) models over RRPP permutations (the distances are thus standardized by the residual covariances). If high-dimensional data are analyzed, a generalized inverse of the residual covariance matrix will be used because of singular covariance matrices. Results are less trustworthy with Mahalanobis distances, in these cases.

The coefficient number should be provided for specific tests. One can determine this with, e.g., coef(fit). If it is not provided (NULL), tests will be performed on all possible vectors of coefficients (rows of coefficients matrix). These tests will be performed sequentially. If a null model is not specified, then for each vector of coefficients, the corresponding parameter is dropped from the linear model design matrix to make a null model. This process is analogous in some ways to a leave-one-out cross-validation (LOOCV) analysis, testing each coefficient against models containing parameters for all other coefficients. For example, for a linear model fit, y ~ x1 + x2 + 1, where x1 and x2 are single-parameter covariates, the analysis would first drop the intercept, then x1, then x2, performing three sequential analyses. This option could require large amounts of computation time for large models, high-dimensional data, many RRPP permutations, or any combination of these. The test results previously reported via coef.lm.rrpp can be found using X.null. One would have to be cognizant of the null model used for each coefficient, based on which term it represents. The function, reveal.model.designs could help determine terms to include in a null model. Regardless, such tests have to be performed iteratively now, but do not require verbose results for initial lm.rrpp fits.

Difference between coef.lm.rrpp test and betaTest

The test for coef.lm.rrpp uses the square-root of inner-products of vectors (d) as a test statistic and only tests the null hypothesis that the length of the vector is 0. The significance of the test is based on random values produced by RRPP, based on the matrices of coefficients that are produced in all permutations. The null models for generating RRPP distributions are consistent with those used for ANOVA, as specified in the lm.rrpp fit by choice of SS type. Therefore, the random coefficients are consistent with those produced by RRPP for generating random estimates used in ANOVA.

The betaTest analysis allows different null hypotheses to be used (vector length is not necessarily 0) and unless otherwise specified, uses a null model that lacks one vector of parameters and a full model that contains all vectors of parameters, for the parameter for which coefficients are estimated. This is closest to a type III SS method of estimation, but each parameter is dropped from the model, rather than terms potentially comprising several parameters. Additionally, betaTest calculates Mahalanobis distance, in addition to Euclidean distance, for vectors of coefficients. This statistic is probably better for more types of models (like generalized least squares fits).

High-dimensional data

If data are high-dimensional (more variables than observations), or even highly multivariate, using Mahalanobis distance can require significant computation time and will require using a generalized inverse. One might wish to consider first whether using principal component scores or other ordinate scores could achieve the same goal. (See ordinate.) For example, one could use the first few principal components as a surrogate for a high-dimensional trait, and test whether the surrogate trait is different than Beta. This would require that the PC scores make sense compared to the original variables, but it would be more computationally tractable.

Generalized Least Squares

To the extent that is possible, tests for GLS estimated coefficients should use Mahalanobis distance. The reason is that the covariance matrix for the data (not to be confused with the residual covariance matrix of a linear model) might not be consistent across RRPP permutations. To assure that random distances are comparable in terms of scale, a generalized (Mahalanobis) distance is safer. However, this can impose a computational burden for high-dimensional data (see above).

Usage

betaTest(
  fit,
  X.null = NULL,
  include.md = FALSE,
  coef.no = NULL,
  Beta = NULL,
  print.progress = FALSE
)

Arguments

Value

Function returns a list with the following components:

Author(s)

Michael Collyer

References

Tejero-Cicuéndez, H., I. Menéndez, A. Talavera, G. Riaño, B. Burriel-Carranza, M. Simó-Riudalbas, S. Carranza, and D.C. Adams. 2023. Evolution along allometric lines of least resistance: Morphological differentiation in Pristurus geckos. Evolution. 77:2547–2560.

See Also

coef.lm.rrpp

Examples

## Not run: 
data(PlethMorph)
fit <- lm.rrpp(TailLength ~ SVL, 
data = PlethMorph,
verbose = TRUE)

## Allometry test (Beta = 0)

T1 <- betaTest(fit, coef.no = 2, Beta = 0)
summary(T1)

# Including Mahalanobis distance

T1 <- betaTest(fit, coef.no = 2, 
Beta = 0, include.md = TRUE)
summary(T1)


# compare to
coef(fit, test = TRUE)

# Note that if Beta is not provided

T1 <- betaTest(fit, coef.no = 2)
summary(T1)

# Note that if coef.no is not provided

T1 <- betaTest(fit)
summary(T1)

# Note that if X.null is provided

T1 <- betaTest(fit, X.null = model.matrix(fit)[, 1],
coef.no = 2)
summary(T1)


## Isometry test (Beta = 1)
# Failure to reject H0 suggests isometric-like association.

T2 <- betaTest(fit, coef.no = 2, Beta = 1)
summary(T2)


## More complex tests

# Multiple covariates

fit2 <- lm.rrpp(HeadLength ~ SVL + TailLength, 
data = PlethMorph,
SS.type = "II",
verbose = TRUE)

fit.null1 <- lm.rrpp(HeadLength ~ SVL, 
data = PlethMorph,
verbose = TRUE)

fit.null2 <- lm.rrpp(HeadLength ~ TailLength, 
data = PlethMorph,
verbose = TRUE)

## allometries
T3 <- betaTest(fit2, fit.null2, coef.no = 2, Beta = 0)
T4 <- betaTest(fit2, fit.null1, coef.no = 3, Beta = 0)
summary(T3)
summary(T4)

# compare to
coef(fit2, test = TRUE)

## isometries
T5 <- betaTest(fit2, fit.null2, coef.no = 2, Beta = 1)
T6 <- betaTest(fit2, fit.null1, coef.no = 3, Beta = 1)

summary(T5)
summary(T6) 

# Intercept test
T7 <- betaTest(fit2, fit.null1, coef.no = 1)
summary(T7)

# multivariate data

PlethMorph$Y <- cbind(PlethMorph$HeadLength, PlethMorph$TailLength)
fit3 <- lm.rrpp(Y ~ SVL, 
data = PlethMorph,
verbose = TRUE)

T8 <- betaTest(fit3, coef.no = 2, Beta = c(0, 0))
T9 <- betaTest(fit3, coef.no = 2, Beta = c(1, 1))

summary(T8)
summary(T9)

## GLS example

fit4 <- lm.rrpp(TailLength ~ SVL, 
data = PlethMorph,
Cov = PlethMorph$PhyCov,
verbose = TRUE)

T10 <- betaTest(fit4, include.md = TRUE)

summary(T10)

# compare to
coef(fit4, test = TRUE)

anova(fit4)


## End(Not run)

Description

The following function has been deprecated in RRPP

Usage

classify()

Details

This function has been deprecated. Use prep.lda instead.

Description

Computes ordinary or generalized least squares coefficients over the permutations of an lm.rrpp model fit with predefined random permutations. For each coefficient vector, the Euclidean distance is calculated as an estimate of the amount of change in Y, the n x p matrix of dependent variables; larger distances mean more change in location in the data space associated with a one unit change in the model design, for the parameter described. Random coefficients are based on either RRPP or FRPP, as defined by the lm.rrpp model fit.

The test argument will be deprecated

The following details will also be removed:

This function can be used to test the specific coefficients of an lm.rrpp fit. The test statistics are the distances (d), which are also standardized (Z-scores). The Z-scores might be easier to compare, as the expected values for random distances can vary among coefficient vectors.

If RRPP is used, all distributions of coefficient vector distances are based on appropriate null models, as defined by SS type. Please be aware that this can result in two seemingly strange but reasonable phenomena. First, if type II or type III SS is used, the intercept will not appear in test results (because the function seeks model parameter differences to know for which coefficients to calculate Euclidean distances). Even if it appears for type I SS, this is merely an artifact of sequential model building and there really is no meaningful test of intercept = 0. Second, Euclidean distances might not always be logical, especially when viewing univariate coefficients, in which case the expected d is |b|. Coefficients without a test are based on the full model; tests are based on the estimates of coefficients (b), given a null model. For example, for a model, y ~ b1 + b2 + b3, with type I SS, b2 will be estimated and tested, using a null model, y ~ b1 and a full model, y ~ b1 + b2. The estimate for b2 might not be the same in the test as when estimated from the model, y ~ b1 + b2 + b3. Therefore, the d statistic might not reflect what one would expect from the full model (like when using type III SS).

Difference between coef.lm.rrpp test and betaTest

The test for coef.lm.rrpp uses the square-root of inner-products of vectors (d) as a test statistic and only tests the null hypothesis that the length of the vector is 0. The significance of the test is based on random values produced by RRPP, based on the matrices of coefficients that are produced in all permutations. The null models for generating RRPP distributions are consistent with those used for ANOVA, as specified in the lm.rrpp fit by choice of SS type. Therefore, the random coefficients are consistent with those produced by RRPP for generating random estimates used in ANOVA.

The betaTest analysis allows different null hypotheses to be used (vector length is not necessarily 0) and unless otherwise specified, uses a null model that lacks one vector of parameters and a full model that contains all vectors of parameters, for the parameter for which coefficients are estimated. This is closest to a type III SS method of estimation, but each parameter is dropped from the model, rather than terms potentially comprising several parameters. Additionally, betaTest calculates Mahalanobis distance, in addition to Euclidean distance, for vectors of coefficients. This statistic is probably better for more types of models (like generalized least squares fits).

coef.lmm.rrpp

The coef.lmm.rrpp function computes combined fixed and random effects by subject for lmm.rrpp model fits. This function does not operate the same as coef.lm.rrpp, for which coefficients and fixed effects are the same, and tests can be performed on fixed effects. It is merely a wrapper to combine the results of ranef.lmm.rrpp and fixef.lmm.rrpp.

Usage

## S3 method for class 'lm.rrpp'
coef(object, SE = FALSE, test = FALSE, confidence = 0.95, ...)

Arguments

Author(s)

Michael Collyer

See Also

betaTest

Examples

## Not run: 
# See examples for lm.rrpp to see how anova.lm.rrpp works in conjunction
# with other functions

data(Pupfish)
names(Pupfish)
Pupfish$logSize <- log(Pupfish$CS)

fit <- lm.rrpp(coords ~ logSize + Sex*Pop, 
SS.type = "I", data = Pupfish, verbose = TRUE) 

coef(fit) # Just coefficients
coef(fit, SE = TRUE) # Coefficients with SE
coef(fit, test = TRUE, 
confidence = 0.99) # Test of coefficients

## End(Not run)

Description

The coef.lmm.rrpp function computes combined fixed and random effects by subject for lmm.rrpp model fits. This function does not operate the same as coef.lm.rrpp, for which coefficients and fixed effects are the same, and tests can be performed on fixed effects. It is merely a wrapper to combine the results of ranef.lmm.rrpp and fixef.lmm.rrpp.

Usage

## S3 method for class 'lmm.rrpp'
coef(object, type = c("matrix", "list"), ...)

Arguments

Author(s)

Michael Collyer

Description

Function attempts to coerce plot information from an RRPP plot object to an amenable ggplot object.

Usage

convert2ggplot(object)

Arguments

Details

This function will attempt to use the plot arguments from an RRPP plot object to make a ggplot that can be additionally updated, as desired. Not all plot characteristics can be converted. For example, text arguments are not currently passed to ggplot, as the text function and geom_text arguments do not easily align. However, one can use text arguments produced by a RRPP plot object and geom_text to augment a ggplot object the way they like.

This function assumes no responsibility for arguments made by ggplot. It merely produces a ggplot object that should resemble an RRPP plot default. Any augmentation of ggplot objects can be done either by direct intervention of the ggplot produced or reformatting the initial RRPP plot produced. One should not expect direct correspondence between R base plot parameters and ggplot parameters. For example, error bars will generally appear as different widths, without an easy way to control them, changing from one format to the other.

Author(s)

Michael Collyer

Examples

## Not run: 
### Linear Model Example

data(Pupfish)
fit <- lm.rrpp(coords ~ log(CS) + Sex*Pop, SS.type = "I", 
data = Pupfish, print.progress = FALSE) 

# Predictions (holding alternative effects constant)

shapeDF <- expand.grid(Sex = levels(Pupfish$Sex), 
Pop = levels(Pupfish$Pop))
rownames(shapeDF) <- paste(shapeDF$Sex, shapeDF$Pop, sep = ".")

shapePreds <- predict(fit, shapeDF)
summary(shapePreds, PC = TRUE)

# Plot prediction

P <- plot(shapePreds, PC = TRUE, ellipse = TRUE)
convert2ggplot(P)

### Ordination Example

data("PlethMorph")

Y <- as.data.frame(PlethMorph[c("TailLength", "HeadLength", 
                               "Snout.eye", "BodyWidth", 
                               "Forelimb", "Hindlimb")])
Y <- as.matrix(Y)
R <- lm.rrpp(Y ~ SVL, data = PlethMorph, 
             print.progress = FALSE)$LM$residuals

# PCA (on correlation matrix)


PCA.ols <- ordinate(R, scale. = TRUE)
PCA.ols$rot
prcomp(R, scale. = TRUE)$rotation # should be the same


PCA.gls <- ordinate(R, scale. = TRUE, 
                   transform. = FALSE, 
                   Cov = PlethMorph$PhyCov)
               
P <- plot(PCA.gls)
convert2ggplot(P)                   

## End(Not run)

Description

A function to find the effect size (Z-score) of a target, from a vector of values presumably obtained in random permutations.

Usage

effect.size(x, center = TRUE, ...)

Arguments

Author(s)

Michael Collyer

Description

Simulated Fish Data for Measurement Error Analysis

Details

Data as simulated in Collyer and Adams (2024), resembling a fish shape, comprising Procrustes coordinates for 11 anatomical landmarks. Data represent 120 configurations for 60 subjects, each with two replicates of measurement. The 60 subjects represent 20 subjects each from three groups.

Author(s)

Michael Collyer

References

Collyer, M.L, and D.C. Adams. 2024. Interrogating random and systematic measurement error in morphometric data. Evolutionary Biology. In press.

Description

Extract fitted values

Usage

## S3 method for class 'lm.rrpp'
fitted(object, ...)

Arguments

Author(s)

Michael Collyer

Examples

# See examples for lm.rrpp

Description

Returns fixed effects from a lmm.rrpp object.

Usage

## S3 method for class 'lmm.rrpp'
fixef(object, ...)

Arguments

Author(s)

Michael Collyer

Description

A data set containing variables for flower production and temperature, averaged yearly for 30 years, from a dry tropical forest.

Details

The data in this data set are samples from a larger study that compared flower production between forest types over 30 years, based on climate variables. The data in this set are from a dry tropical forest, include variables for annual flower production, temperature, and year. These data are suitable for considering the temporal non-independence of observations in generalized least squares (GLS) estimation, for the example in the book, Permutational Biometry (Collyer & Adams, 2025).

References

Pau, S., Wolkovich, E. M., Cook, B. I., Nytch, C. J., Regetz, J., Zimmerman, J. K., and Joseph Wright, S. (2013). Clouds and temperature drive dynamic changes in tropical flower production. Nature Climate Change, 3:838–842.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Reduces a plot.measurement.error to a single research subject. This can be for cases when many overlapping subjects in a plot obscure interpretation for specific subjects.

Usage

focusMEonSubjects(x, subjects = NULL, shadow = TRUE, ...)

Arguments

Author(s)

Michael Collyer

Description

A function mostly for internal processing but can be used to extract ANOVA statistics for other uses, such as plotting histograms

Usage

getANOVAStats(fit, stat = c("SS", "MS", "Rsq", "F", "cohenf", "all"))

Arguments

Author(s)

Michael Collyer

Examples

## Not run: 
data(Pupfish)
fit <- lm.rrpp(coords ~ log(CS) + Sex*Pop, SS.type = "I", 
data = Pupfish, print.progress = FALSE, iter = 999) 
anova(fit)
Fstats <- getANOVAStats(fit, stat = "F")
par(mfrow = c(2, 2))
hist(Fstats$Fs[1,], breaks = 50, main = "log(CS)", xlab = "F")
abline(v = Fstats$Fs[1, 1])
hist(Fstats$Fs[2,], breaks = 50, main = "Sex", xlab = "F")
abline(v = Fstats$Fs[2, 1])
hist(Fstats$Fs[3,], breaks = 50, main = "Pop", xlab = "F")
abline(v = Fstats$Fs[3, 1])
hist(Fstats$Fs[4,], breaks = 50, main = "Sex:Pop", xlab = "F")
abline(v = Fstats$Fs[4, 1])

## End(Not run)

Description

A function mostly for internal processing but can be used to extract either the model covariance matrix (Cov) or the projection matrix for transformations made from the covariance matrix (Pcov), which is basically the square-root of the covariance matrix. This matrix is the model covariance used for estimation, not the residual covariance matrix (see getResCov). There are also options for S3 or S4 format versions, or a forcing of symmetry for Pcov.

Usage

getModelCov(
  fit,
  type = c("Cov", "Pcov"),
  format = c("S3", "S4"),
  forceSym = TRUE
)

Arguments

Author(s)

Michael Collyer

Description

A function mostly for internal processing but can be used to obtain terms, design matrices, or QR decompositions used for each reduced or full model that is fitted in an lm.rrpp fit.

Usage

getModels(fit, attribute = c("terms", "X", "qr", "all"))

Arguments

Author(s)

Michael Collyer

Description

A function mostly for internal processing but can be used to extract RRPP permutation information for other reasons.

Usage

getPermInfo(
  fit,
  attribute = c("perms", "perm.method", "block", "perm.seed", "perm.schedule", "all")
)

Arguments

Author(s)

Michael Collyer

Description

A function mostly for internal processing but can be used to extract the residual covariance matrix. This matrix is the residual covariance matrix, not the model covariance matrix used for estimation (see getModelCov). Options for averaging over degrees of freedom or number of observations, plus standardization, are also available.

Usage

getResCov(fit, useDf = TRUE, standardize = FALSE)

Arguments

Author(s)

Michael Collyer

Description

A function mostly for internal processing but can be used to extract The terms for each reduced and full model used in an lm.rrpp fit.

Usage

getTerms(fit)

Arguments

Author(s)

Michael Collyer

Description

Function performs analyses concerned with the repeatability (reliability) of multivariate data (measurements) collected from the same research subjects. Although there is no requirement for repeated measurements on all research subjects, the analysis assumes that multiple observations are made.

Usage

ICCstats(
  fit,
  subjects = NULL,
  with_in = NULL,
  groups = NULL,
  multivariate = FALSE,
  print.AOV = TRUE
)

Arguments

Details

Function uses ANOVA statistics or SSCP matrices to find the ratio of among-subject to within-subject variance. The former is a dispersion-based approach and the latter is a multivariate generalization of the ICC statistic (as a matrix product). The multivariate generalizations of the statistics described by Liljequist et al. (2019) are used to find matrix products, from which eigenanalysis is performed, providing ICC statistics by eigenvectors.

Three statistics describe the ICC for the population, agreement of measurements among subjects, and consistency between measurements. The last statistic does not necessarily measure the sameness between measurements but the consistency of change between measurements, which might be indicative of a systematic measurement error. If groups are used, these three statistics are repeated, using the SSCP for groups-adjusted data. This approach accounts for group differences, which would avoid large subject variation compared to measurement error inflating ICC values. If there are inherently disparate groups from which subjects are sampled, this approach can elucidate better agreement and consistency in light of group differences.

This function is most useful for analyses performed with measurement.error, but any lm.rrpp fit can be used, so long as research subjects can be defined.

It is essential that all arguments are terms that can be found in the model frame of the model fit, as provoke by ANOVA. Using anova(fit) will elucidate the row names of the ANOVA that could be used.

Value

Objects of class "ICCstats" return the following:

Author(s)

Michael Collyer

References

Liljequist, D., Elfving, B., & Skavberg Roaldsen, K. (2019). Intraclass correlation–A discussion and demonstration of basic features. PloS one, 14(7), e0219854.

Examples

## Not run: 
# Measurement error analysis on simulated data of fish shapes

data(fishy)

# Analysis unconcerned with groups 

ME1 <- measurement.error(
  Y = "coords",
  subjects = "subj",
  replicates = "reps",
  data = fishy)

anova(ME1)
ICCstats(ME1, subjects = "Subjects", with_in = "Systematic ME")

# Analysis concerned with groups 

ME2 <- measurement.error(
  Y = "coords",
  subjects = "subj",
  replicates = "reps",
  groups = "groups",
  data = fishy)

anova(ME2)
ICCstats(ME2, subjects = "Subjects", 
  with_in = "Systematic ME", groups = "groups")
ICCstats(ME2, subjects = "Subjects", 
  with_in = c("Systematic ME", "Systematic ME:Groups"), 
  groups = "groups")
  
## End(Not run)

Description

Function produces both a list of inter-subject Euclidean distance matrices, based on replicate measurements of the same subjects, and one matrix that summarizes the variability among the inter-subject distances, across subjects. This function can be considered a tool for the evaluation of subject estimate precision. The function, plot.interSubVar can produce a heat map of the inter-subject variability.

Usage

interSubVar(ME, type = c("range", "sd", "var", "cv"))

Arguments

Value

An object of class interSubVar is a list containing the following

Author(s)

Michael Collyer

Examples

## Not run: 
# Measurement error analysis on simulated data of fish shapes

data(fishy)

# Analysis unconcerned with groups 

ME1 <- measurement.error(
  Y = "coords",
  subjects = "subj",
  replicates = "reps",
  data = fishy)

anova(ME1)
ICCstats(ME1, subjects = "Subjects", with_in = "Systematic ME")
plot(ME1)

# Analysis concerned with groups 

ME2 <- measurement.error(
  Y = "coords",
  subjects = "subj",
  replicates = "reps",
  groups = "groups",
  data = fishy)

anova(ME2)
ICCstats(ME2, subjects = "Subjects", 
  with_in = "Systematic ME", groups = "groups")
P <- plot(ME2)
focusMEonSubjects(P, subjects = 18:20, shadow = TRUE)
 
## End(Not run)

Description

Function performs a K-comonent analysis, which is an eigen decomposition of a ratio of ordinary least-squares (OLS) to generalized least squares (GLS) covariance matrices.

Usage

kcomp(Y, Cov, transform. = TRUE, scale. = FALSE, tol = NULL, rank. = NULL)

Arguments

Details

The function performs a K-component analysis (KCA), as described in Mitteroecker et al. (2025). The analysis is similar to many that can be performed with the ordinate function, but whereas that function finds ordinate scores for a rotation or shear of a data space, KCA performs a relative eigen decomposition of a matrix product that represents a ratio between two covariance matrices (Mitteroecker and Bookstein, 2014). The eigenvalues are helpful for understanding the importance of the GLS estimation of a covariance matrix. For example, if the GLS estimation is made with respect to a matrix of phylogenetic covariances, the eigenvalues express the distribution of phylogenetic signal across components of the data space.

The matrix product (K matrix) is the inverse of the OLS covariance matrix times the GLS covariance matrix. Redundancies between these matrices will mean the K matrix is likely not full rank. Therefore, an algorithm is used to compute scores in appropriate dimensions. First, data are aligned to the square-root matrix of the GLS covariance matrix (Collyer and Adams, 2021), using the ordinate function. Second, relative eigen decomposition is performed, and the number of real eigenvectors that can be computed are retained. Lastly, the same number of aligned components are used for projection of data onto relative eigenvectors (the K components).

Value

An object of class kcomp is a list containing the following

Author(s)

Michael Collyer

References

Collyer, M.L. and D.C. Adams. 2021. Phylogenetically-aligned Component Analysis. Methods in Ecology and evolution. In press.

Bookstein, F. L., & Mitteroecker, P. (2014). Comparing covariance matrices by relative eigenanalysis, with applications to organismal biology. Evolutionary biology, 41, 336-350.

Mitteroecker, P., Collyer, M. L., & Adams, D. C. (2024). Exploring Phylogenetic Signal in Multivariate Phenotypes by Maximizing Blomberg’s K. Systematic Biology, syae035.

See Also

plot.kcomp, ordinate, plot.default, rrpp.data.frame

Examples

data(PlethMorph)
 PCA <- ordinate(as.data.frame(PlethMorph[c("TailLength",
 "HeadLength", "Snout.eye", "BodyWidth","Forelimb", "Hindlimb")],
 scale. = TRUE))
 Y <- PCA$x
 KCA <- kcomp(Y, Cov = PlethMorph$PhyCov, tol = 0.001)

Description

Function performs a linear model fit over many random permutations of data, using a randomized residual permutation procedure.

Usage

lm.rrpp(
  f1,
  iter = 999,
  turbo = FALSE,
  seed = NULL,
  int.first = FALSE,
  RRPP = TRUE,
  full.resid = FALSE,
  block = NULL,
  SS.type = c("I", "II", "III"),
  data = NULL,
  Cov = NULL,
  print.progress = FALSE,
  Parallel = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

The function fits a linear model using ordinary least squares (OLS) or generalized least squares (GLS) estimation of coefficients over any number of random permutations of the data. A permutation procedure that randomizes vectors of residuals is employed. This procedure can randomize two types of residuals: residuals from null models or residuals from an intercept model. The latter is the same as randomizing full values, and is referred to as as a full randomization permutation procedure (FRPP); the former uses the residuals from null models, which are defined by the type of sums of squares and cross-products (SSCP) sought in an analysis of variance (ANOVA), and is referred to as a randomized residual permutation procedure (RRPP). Types I, II, and III SSCPs are supported.

Users define the SSCP type, the permutation procedure type, whether a covariance matrix is included (GLS estimation), and a few arguments related to computations. Results comprise observed linear model results (coefficients, fitted values, residuals, etc.), random sums of squares (SS) across permutation iterations, and other parameters for performing ANOVA and other hypothesis tests, using empirically-derived probability distributions.

lm.rrpp emphasizes estimation of standard deviates of observed statistics as effect sizes from distributions of random outcomes. When performing ANOVA, using the anova function, the effect type (statistic choice) can be varied. See anova.lm.rrpp for more details. Please recognize that the type of SS must be chosen prior to running lm.rrpp and not when applying anova to the lm.rrpp fit, as design matrices for the linear model must be created first. Therefore, SS.type is an argument for lm.rrpp and effect.type is an argument for anova.lm.rrpp. If MANOVA statistics are preferred, eigenvalues can be added with manova.update and statistics summarized with summary.manova.lm.rrpp. See manova.update for examples.

The coef.lm.rrpp function can be used to test the specific coefficients of an lm.rrpp fit. The test statistics are the distances (d), which are also standardized (Z-scores). The Z-scores might be easier to compare, as the expected values for random distances can vary among coefficient vectors (Adams and Collyer 2016).

ANOVA vs. MANOVA

Two SSCP matrices are calculated for each linear model effect, for every random permutation: R (Residuals or Random effects) and H, the difference between SSCPs for "full" and "reduced" models. (Full models contain and reduced models lack the effect tested; SSCPs are hypothesized to be the same under a null hypothesis, if there is no effect. The difference, H, would have a trace of 0 if the null hypothesis were true.) In RRPP, ANOVA and MANOVA correspond to two different ways to calculate statistics from R and H matrices.

ANOVA statistics are those that find the trace of R and H SSCP matrices before calculating subsequent statistics, including sums of squares (SS), mean squares (MS), and F-values. These statistics can be calculated with univariate data and provide univariate-like statistics for multivariate data. These statistics are dispersion measures only (covariances among variables do not contribute) and are the same as "distance-based" stats proposed by Goodall (1991) and Anderson (2001). MANOVA stats require multivariate data and are implicitly affected by variable covariances. For MANOVA, the inverse of R times H (invR.H) is first calculated for each effect, then eigenanalysis is performed on these matrix products. Multivariate statistics are calculated from the positive, real eigenvalues. In general, inferential conclusions will be similar with either approach, but effect sizes might differ.

ANOVA tables are generated by anova.lm.rrpp on lm.rrpp fits and MANOVA tables are generated by summary.manova.lm.rrpp, after running manova.update on lm.rrpp fits.

Currently, mixed model effects are only possible with $ANOVA statistics, not $MANOVA.

More detail is found in the vignette, ANOVA versus MANOVA.

Notes for RRPP 0.5.0 and subsequent versions

The output from lm.rrpp has changed, compared to previous versions. First, the $LM component of output no longer includes both OLS and GLS statistics, when GLS fits are performed. Only GLS statistics (coefficients, residuals, fitted values) are provided and noted with a "gls." tag. GLS statistics can include those calculated when weights are input (similar to the lm argument). Unlike previous versions, GLS and weighted LS statistics are not labeled differently, as weighted LS is one form of generalized LS estimation. Second, a new object, $Models, is included in output, which contains the linear model fits (lm attributes ) for all reduced and full models that are possible to estimate fits.

Notes for RRPP 0.3.1 and subsequent versions

F-values via RRPP are calculated with residual SS (RSS) found uniquely for any model terms, as per Anderson and ter Braak (2003). This method uses the random pseudo-data generated by each term's null (reduced) model, meaning RSS can vary across terms. Previous versions used an intercept-only model for generating random pseudo-data. This generally has appropriate type I error rates but can have elevated type I error rates if the observed RSS is small relative to total SS. Allowing term by term unique RSS alleviates this concern.

Value

An object of class lm.rrpp is a list containing the following

Author(s)

Michael Collyer

References

Anderson MJ. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology 26: 32-46.

Anderson MJ. and C.J.F. ter Braak. 2003. Permutation tests for multi-factorial analysis of variance. Journal of Statistical Computation and Simulation 73: 85-113.

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 115:357-365.

Adams, D.C. and M.L. Collyer. 2016. On the comparison of the strength of morphological integration across morphometric datasets. Evolution. 70:2623-2631.

Adams, D.C and M.L. Collyer. 2018. Multivariate phylogenetic anova: group-clade aggregation, biological challenges, and a refined permutation procedure. Evolution. 72:1204-1215.

ter Braak, C.J.F. 1992. Permutation versus bootstrap significance tests in multiple regression and ANOVA. pp .79–86 In Bootstrapping and Related Techniques. eds K-H. Jockel, G. Rothe & W. Sendler.Springer-Verlag, Berlin. lm for more on linear model fits.

See Also

procD.lm and procD.pgls within geomorph;

Examples

## Not run: 

# Examples use geometric morphometric data
# See the package, geomorph, for details about obtaining such data

data("PupfishHeads")
names(PupfishHeads)

# Head Size Analysis (Univariate)-------------------------------------------------------

fit <- lm.rrpp(log(headSize) ~ sex + locality/year, SS.type = "I", 
data = PupfishHeads, print.progress = FALSE, iter = 999)
summary(fit)
anova(fit, effect.type = "F") # Maybe not most appropriate
anova(fit, effect.type = "Rsq") # Change effect type, but still not 
# most appropriate

# Mixed-model approach (most appropriate, as year sampled is a random 
# effect:

anova(fit, effect.type = "F", error = c("Residuals", "locality:year", 
"Residuals"))

# Change to Type III SS

fit <- lm.rrpp(log(headSize) ~ sex + locality/year, SS.type = "III", 
data = PupfishHeads, print.progress = FALSE, iter = 999,
verbose = TRUE)
summary(fit)
anova(fit, effect.type = "F", error = c("Residuals", "locality:year", 
"Residuals"))

# Coefficients Test

coef(fit, test = TRUE)

# Predictions (holding alternative effects constant)

sizeDF <- data.frame(sex = c("Female", "Male"))
rownames(sizeDF) <- c("Female", "Male")
sizePreds <- predict(fit, sizeDF)
summary(sizePreds)
plot(sizePreds)

# Diagnostics plots of residuals

plot(fit)

# Body Shape Analysis (Multivariate) -----------

data(Pupfish)
names(Pupfish)

# Note:

dim(Pupfish$coords) # highly multivariate!
fit <- lm.rrpp(coords ~ log(CS) + Sex*Pop, SS.type = "I", 
data = Pupfish, print.progress = FALSE, iter = 999,
verbose = TRUE) 
summary(fit, formula = FALSE)
anova(fit) 
coef(fit, test = TRUE)

# Predictions (holding alternative effects constant)

shapeDF <- expand.grid(Sex = levels(Pupfish$Sex), 
Pop = levels(Pupfish$Pop))
rownames(shapeDF) <- paste(shapeDF$Sex, shapeDF$Pop, sep = ".")
shapeDF

shapePreds <- predict(fit, shapeDF)
summary(shapePreds)
summary(shapePreds, PC = TRUE)

# Plot prediction

plot(shapePreds, PC = TRUE)
plot(shapePreds, PC = TRUE, ellipse = TRUE)

# Diagnostics plots of residuals

plot(fit)

# PC-plot of fitted values

groups <- interaction(Pupfish$Sex, Pupfish$Pop)
plot(fit, type = "PC", pch = 19, col = as.numeric(groups))

# Regression-like plot

plot(fit, type = "regression", reg.type = "PredLine", 
    predictor = log(Pupfish$CS), pch=19,
    col = as.numeric(groups))

# Body Shape Analysis (Distances) ----------

D <- dist(Pupfish$coords) # inter-observation distances
length(D)
Pupfish$D <- D

fitD <- lm.rrpp(D ~ log(CS) + Sex*Pop, SS.type = "I", 
data = Pupfish, print.progress = FALSE, iter = 999) 

# These should be the same:
summary(fitD, formula = FALSE)
summary(fit, formula = FALSE) 

# GLS Example (Univariate) -----------

data(PlethMorph)
fitOLS <- lm.rrpp(TailLength ~ SVL, data = PlethMorph, 
print.progress = FALSE, iter = 999)
fitGLS <- lm.rrpp(TailLength ~ SVL, data = PlethMorph, Cov = PlethMorph$PhyCov, 
print.progress = FALSE, iter = 999)

anova(fitOLS)
anova(fitGLS)

sizeDF <- data.frame(SVL = sort(PlethMorph$SVL))

# Prediction plots

# By specimen
plot(predict(fitOLS, sizeDF)) # Correlated error
plot(predict(fitGLS, sizeDF)) # Independent error

# With respect to independent variable (using abscissa)
plot(predict(fitOLS, sizeDF), abscissa = sizeDF) # Correlated error
plot(predict(fitGLS, sizeDF), abscissa = sizeDF) # Independent error


# GLS Example (Multivariate) -----------

Y <- as.matrix(cbind(PlethMorph$TailLength,
PlethMorph$HeadLength,
PlethMorph$Snout.eye,
PlethMorph$BodyWidth,
PlethMorph$Forelimb,
PlethMorph$Hindlimb))
PlethMorph$Y <- Y
fitOLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph, 
print.progress = FALSE, iter = 999)
fitGLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph, 
Cov = PlethMorph$PhyCov,
print.progress = FALSE, iter = 999)

anova(fitOLSm)
anova(fitGLSm)

# Prediction plots

# By specimen
plot(predict(fitOLSm, sizeDF)) # Correlated error
plot(predict(fitGLSm, sizeDF)) # Independent error

# With respect to independent variable (using abscissa)
plot(predict(fitOLSm, sizeDF), abscissa = sizeDF) # Correlated error
plot(predict(fitGLSm, sizeDF), abscissa = sizeDF) # Independent error

## End(Not run)

Description

Function performs a linear model fit over many random permutations of data, using a randomized residual permutation procedure restricted to subjects.

Usage

lm.rrpp.ws(
  f1,
  subjects,
  iter = 999,
  turbo = FALSE,
  seed = NULL,
  int.first = FALSE,
  RRPP = TRUE,
  data,
  Cov = NULL,
  delta = 0.001,
  gamma = c("sample", "equal"),
  print.progress = FALSE,
  verbose = FALSE,
  Parallel = FALSE,
  ...
)

Arguments

Details

The function fits a linear model using ordinary least squares (OLS) or generalized least squares (GLS) estimation of coefficients over any number of random permutations of the data, but the permutations are mostly restricted to occur with subject blocks for any model terms other than subjects. All functionality should resemble that of lm.rrpp. However, an argument for research subjects is also required. The purpose of this function is to account for the non-independence among observations of research subjects (due to sampling within subjects), while also allowing for the non-independence among subjects to be considered (Adams and Collyer, submitted).

By comparison, the covariance matrix option in lm.rrpp must have a one-to-one match to observations, which can be matched by the row names of the data. In this function, the covariance matrix can be the same one used in lm.rrpp but the number of observations can be greater. For example, if subjects are species or some other level of taxonomic organization, data can comprise measurements on individuals. Users have the option to expand the covariance matrix for subjects or input one they have generated.

Irrespective of covariance matrix type, the row names of the data matrix must match the subjects. This step assures that the analysis can proceed in lm.rrpp. It is also best to make sure to use an rrpp.data.frame, so that the subjects can be a name in that data frame. For example, if research subjects are species and data (observations) are collected from individuals within species, then a procedure like the following should produce results:

rownames(Y) <- species

rdf <- rrpp.data.frame(Y = Y, subjects = species, x = x)

fit <- lm.rrpp.ws(Y ~ species * x, subject = species, data = rdf, Cov = myCov, ...)

where ... means other arguments. The covariances in the the Covariance matrix can be sorted by the subjects factor but data will not be sorted. Therefore, names matching the subjects is essential. Additionally, subjects must be a factor in the data frame or a factor in the global environment. It cannot be part of a list. Something like subjects <- mylist$species will not work. Assuring that data and subjects are in the same rrpp.data.frame object as data is the best way to avoid errors.

Most attributes for this analysis are explained with lm.rrpp. The notable different attributes for this function are that: (1) a covariance matrix for the non-independence of subjects can be either a symmetric matrix that matches in dimensions the number of subjects or the number of observations; (2) a parameter (delta) that can range between near 0 and 1 to calibrate the covariances between observations of different subjects; and (3) a parameter (gamma) that is either 1 (equal) or the square-root of the subject sample size (sample) to calibrate the covariances among observations within subjects. If delta = 0, it is expected that the covariance between individual observations, between subjects, is the same as expected from the covariance matrix, as if observations were the single observations made on subjects. As delta approaches 1, the observations become more independent, as if it is expected that the many observations would not be expected to be as correlated as if from one observation. Increasing delta might be useful, if, for example, many individuals are sampled within species, from different locations, different age groups, etc. Alternatively, the sample size (n_i) for subject i can also influence the trend of inter-subject covariances. If more individual observations are sampled, the correlation between subjects might be favored to be smaller compared to fewer observations. The covariances can be adjusted to allow for greater independence among observations to be assumed for larger samples.

A design matrix, X, is constructed with 0s and 1s to indicate subjects association, and it is used to expand the covariance matrix (C) by XCt(X), where t(X) is the matrix transpose. The parameters in X are multiplied by exp(-delta * gamma) to scale the covariances. (If delta = 0 and gamma = 1, they are unscaled.)

These options for scaling covariances could be important for data with hierarchical organization. For example, data sampled from multiple species with expected covariances among species based on phylogenetic distances, might be expected to not covary as strongly if sampling encounters other strata like population, sex, and age. An a priori expectation is that covariances among observations would be expected to be smaller than between species, if only one observation per species were made.

If one wishes to have better control over between-subject and within-subject covariances, based on either a model or empirical knowledge, a covariance matrix should be generated prior to analysis. One can input a covariance matrix with dimensions the same as XCt(X), if they prefer to define covariances in an alternative way. A function to generate such matrices based on separate inter-subject and intra-subject covariance matrices is forthcoming.

IMPORTANT. It is assumed that either the levels of the covariance matrix (if subject by subject) match the subject levels in the subject argument, or that the order of the covariance matrix (if observation by observation) matches the order of the observations in the data. No attempt is made to reorder a covariance matrix by observations and row-names of data are not used to re-order the covariance matrix. If the covariance matrix is small (same in dimension as the number of subject levels), the function will compile a large covariance matrix that is correct in terms of order, but this is based on the subjects argument, only.

The covariance matrix is important for describing the expected covariances among observations, especially knowing observations between and within subjects are not independent. However, the randomization of residuals in a permutation procedure (RRPP) is also important for testing inter-subject and intra-subject effects. There are two RRPP philosophies used. If the variable for subjects is part of the formula, the subject effect is evaluated with type III sums of squares and cross-products (estimates SSCPs between a model with all terms and a model lacking subject term), and RRPP performed for all residuals of the reduced model. Effects for all other terms are evaluated with type II SSCPs and RRPP restricted to randomization of reduced model residuals, within subject blocks. This assures that subject effects are held constant across permutations, so that intra-subject effects are not confounded by inter-subject effects.

More details will be made and examples provided after publication of articles introducing the novel RRPP approach.

The lm.rrpp arguments not available for this function include: full.resid, block, and SS.type. These arguments are fixed because of the within-subject blocking for tests, plus the requirement for type II SS for within-subject effects.

Value

An object of class lm.rrpp.ws is a list containing the following

Author(s)

Michael Collyer

References

Adams, D.C and M.L Collyer. (submitted) Extended phylogenetic regression models for comparing within-species patterns across the Tree of Life. Methods in Ecology and Evolution

ter Braak, C.J.F. 1992. Permutation versus bootstrap significance tests in multiple regression and ANOVA. pp .79–86 In Bootstrapping and Related Techniques. eds K-H. Jockel, G. Rothe & W. Sendler.Springer-Verlag, Berlin. lm for more on linear model fits.

See Also

lm.rrpp; measurement.error

Examples

## Not run: 
data(fishy)

suppressWarnings(fit <- lm.rrpp.ws(coords ~ subj + groups * reps,
  subjects = "subj", 
  data = fishy))

anova(fit)

## End(Not run)

Description

Function performs a mixed linear model fit over many random permutations of data, using a randomized residual permutation procedure restricted to subjects. This function is likely to evolve greater flexibility in the future.

Usage

lmm.rrpp(
  fixed,
  subjects,
  type = c("intercepts", "slopes"),
  slopeTerm = NULL,
  estimation = c("REML", "ML", "LS"),
  iter = 999,
  data,
  ...
)

Arguments

Details

The function fits a mixed linear model using ordinary least squares (OLS) or generalized least squares (GLS) estimation of coefficients for fixed effects, and maximum likelihood (ML) or restricted ML (REML) for random effects, over any number of random permutations of the data, but the permutations are mostly restricted to occur with subject blocks for any model terms other than subjects. Most functionality should resemble that of lm.rrpp and lm.rrpp.ws, with some caveats. First, estimation of random effects uses the methodology of the lmer function (Bates et al., 2015), for only one random term based on a subjects designation. Second, this means that a "random slopes" model involves an interaction between subjects and a numeric variable (not a categorical variable). Third, although any number of fixed effects can be used, only one fixed effect can interaction with subjects. Finally, the function is intended to allow analysis of variance (ANOVA), like with lm.rrpp.ws, but not treated random effects as fixed. There is not a coefficients test associated with this function, as with other lm.rrpp that include only fixed effects.

Although this function allows covariance matrices to be used for generalized least squares (GLS) estimation of fixed effects, it does not currently expand the covariance mnatrix based only on subjects, as does lm.rrpp.ws. However, that function can be used to obtain a covariance matrix for this function. The covariance structure (Lambda matrix) among subjects is estimated via ML or REML (Bates et al., 2015) and is not impacted by the expected covariances among fixed effects, however, the expected covariances will influence interactions between fixed and random effects. The Lambda matrix is estimated as a single matrix, even for multivariate data. This matrix is obtained by projecting data onto the first component of a singular value decomposition of the product between the data and a hat matrix of the the fixed effects (basically projection following a principal component analysis of model fitted values). This is a computationally proficient step that is similar to averaging covariance parameters (theta) over all variables. If independent Lambda matrices are desired, then this function can be used on, e.g., each principal component of data, followed by combining the coefficients. For ANOVA, the precise estimation of Lambda is not as important for tests of fixed effects, as it means a linear transformation of coefficients that would be consistent across random permutations of RRPP. It might be more important for model comparison and selection, in which case whether to use LS, ML, REML, and whether to use separate Lambda matrices is a much deeper consideration.

A data frame (preferably an rrpp.data.frame object) must be used for this function. Not all "downstream" functions (e.g., logLik, predict) are guaranteed to work as expected, as this function is still very much under development.

The lm.rrpp arguments not available for this function include: full.resid, block, and SS.type. These arguments are fixed because of the within-subject blocking for tests, plus the requirement for type II SS for within-subject effects.

Value

An object of class lmm.rrpp is a list containing the following

Author(s)

Michael Collyer

References

Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.

Limitations

This function currently only works for "Random Intercepts" and "Random slopes" models, where the intercepts or slopes of a single subject term (random effect) are free to vary. If a Random Slopes model is used, a single covariate (numeric variable) must be specified for slopes. Generally, this is term that is also among the fixed effects. These limitations are imposed to allow ML or REML estimation. More complex models can be entertained with lmm.rrpp.ws, but using only "LS" estimation. If "LS" estimation is used in this function, results should be the same with a comparable formula, using lmm.rrpp.ws

If a ML or REML fit is near singular (one or more covariance parameters are 0), the analysis will be arrested rather than continued with some random effects rendered as 0. This should be an indication that the model is over-fit with the type of estimation used.

These limitations might be addressed with future versions of this function, as theoretical advances are made.

See Also

lm.rrpp; lm.rrpp.ws; measurement.error

Examples

## Not run: 
# Random intercepts models

data(ReptAbund)
fit_ls <- lmm.rrpp(fixed = reptile_abundance ~ Year, 
subjects = "Site",      
type = "intercepts", 
estimation = "LS", 
data = ReptAbund)

fit_ml <- lmm.rrpp(fixed = reptile_abundance ~ Year, 
subjects = "Site",      
type = "intercepts", 
estimation = "ML", 
data = ReptAbund)

fit_reml <- lmm.rrpp(fixed = reptile_abundance ~ Year, 
subjects = "Site",      
type = "intercepts", 
estimation = "REML", 
data = ReptAbund) 

anova(fit_ls, error = c("Site", "Residuals"))
anova(fit_ml, error = c("Site", "Residuals"))
anova(fit_reml, error = c("Site", "Residuals"))

coef(fit_ls, type = "list")
coef(fit_ml, type = "list")
coef(fit_reml, type = "list")

# Random slopes models
# none should not work!

fit_ls <- lmm.rrpp(fixed = reptile_abundance ~ Year, 
subjects = "Site",      
type = "slopes",
slopeTerm = "Year", 
estimation = "LS", 
data = ReptAbund)

fit_ml <- lmm.rrpp(fixed = reptile_abundance ~ Year, 
subjects = "Site",      
type = "slopes",
slopeTerm = "Year", 
estimation = "ML", 
data = ReptAbund)


fit_reml <- lmm.rrpp(fixed = reptile_abundance ~ Year, 
subjects = "Site",      
type = "slopes",
slopeTerm = "Year",  
estimation = "REML", 
data = ReptAbund) 

anova(fit_ls, error = c("Year:Site", "Year:Site", "Residuals"))
anova(fit_ml, error = c("Year:Site", "Year:Site", "Residuals"))

coef(fit_ls, type = "list")
coef(fit_ml, type = "list")
 
## End(Not run)

Description

logLik.lm.rrpp returns the log-likelihood of an lm.rrpp object or an lmm.rrpp object. Ridge regularization will be performed for ill-conditioned or singular residual covariance matrices, but dimension reduction could be augmented via projection, using the arguments, tol and pc.no. See ordinate for details.

Usage

## S3 method for class 'lm.rrpp'
logLik(
  object,
  tol = NULL,
  pc.no = NULL,
  Z = TRUE,
  verbose = FALSE,
  gls.null = FALSE,
  ...
)

Arguments

Details

Note that options are limited for lmm.rrpp objects. It is currently not possible to use RRPP to find Z scores.

Author(s)

Michael Collyer

Description

Function performs a leave-one-out cross-validation estimate (method = "fit") or ordination scores from estimates (method = "PC"). The latter can be helpful for determining if apparent group differences in ordination plots arise merely from data dimensionality. The choice of method essentially chooses one of two functions.

Usage

looCV(fit, method = c("fit", "PC"), ...)

Arguments

Details

The fit method is standard, using n - 1 observations as training data, and estimating the response for 1 datum, not included in the training data, based on the model design of an original model fit. This method provides alternative fitted values than the original fit.

The PC method uses the strategy of Thioulouse et al. (2021) to perform N ordinations for N observations, in which each of the N observations are left out of the estimation of linear model coefficients, but the vector of data for the left-out observation is projected on the eigenvectors of the fitted values obtained from the leave-one-out cross-validation (jackknife) strategy. The purpose of this diagnostic tool is to determine whether apparent "group differences" in an ordination plot (using the function, ordinate) are because of high-dimensional data (number of variables exceed number of observations) rather than real differences. An apparent group difference is common for high-dimensional data, when variables are far greater in number than observations (Cardini et al., 2019). However, leave-one-out cross-validation can help elucidate whether an observed visual difference is spurious.

This function differs from the strategy of Thioulouse et al. (2021) in two important ways. First, this function uses the linear model design from a lm.rrpp fit, and can contain any number of independent variables, rather than a single factor for groups. Second, after obtaining leave-one-out cross-validated scores, a Procrustes alignment between cross-validated scores and "observed" (real) scores is performed, which minimizes summed squared distances between the alternative ordinations. This latter step assures comparisons are appropriate.

The type = "PC" plot from plot.lm.rrpp has the same scores as obtained from ordinate(Y, A = H), using the ordinate function, where H is a hat matrix (that can be calculated from plot.lm.rrpp output), and Y is a matrix of data. This function updates H for every possible case that one row of Y is left out (meaning the rotation matrix from ordinate is updated N times). If the H matrix is robust in spite of dropped data and design matrix parameters, the result will be similar to the original ordination. If apparent group differences are spurious, H will tend to change, as will data projections.

The functions summary.looCV and plot.looCV are essential for evaluating results. These support functions compare eigenvalues and projected scores, between observed and cross-validated cases.

This function should be viewed as a diagnostic tool and not as a data transformation tool! The cross-validated scores will not retain Euclidean distances among observations. This could cause problems in analyses that substitute cross-validated scores as data.

Value

An object of class looCV is a list containing the following

Author(s)

Michael Collyer

References

Thioulouse, J., Renaud, S., Dufour, A. B., & Dray, S. (2021). Overcoming the Spurious Groups Problem in Between-Group PCA. Evolutionary Biology, In press.

Cardini, A., O’Higgins, P., & Rohlf, F. J. (2019). Seeing distinct groups where there are none: spurious patterns from between-group PCA. Evolutionary Biology, 46(4), 303-316.

See Also

summary.looCV, plot.looCV, looCV_table

Examples

# Example with real group differences

data(Pupfish)
fit <- lm.rrpp(coords ~ Pop * Sex, 
data = Pupfish, iter = 0)
CV1f <- looCV(fit)

summary(CV1f)

CV1p <- looCV(fit, method = "PC")

summary(CV1p)

group <- interaction(Pupfish$Pop, Pupfish$Sex)
plot(CV1p, flip = 1, pch = 19, col = group)

# Example with apparent but not real group differences

n <- NROW(Pupfish$coords)
p <- NCOL(Pupfish$coords)
set.seed(1001)
Yr <- matrix(rnorm(n * p), n, p) # random noise

fit2 <-lm.rrpp(Yr ~ Pop * Sex, 
data = Pupfish, iter = 0)
CV2p <- looCV(fit2, method = "PC")
summary(CV2p)
group <- interaction(Pupfish$Pop, Pupfish$Sex)
plot(CV2p, pch = 19, col = group)

Description

Function simply produces a table of MSE and cross-validated MSE for any number of compared models. The table is sorted from lowest to highest cross-validated MSE. The looCV objects must use method = 'fit' or a table cannot be generated.

Usage

looCV_table(...)

Arguments

Author(s)

Michael Collyer

See Also

looCV, model.comparison

Examples

# Example with snakeHS data

data(snakeHS)
fitS <- lm.rrpp(HS ~ SVL, data = snakeHS,
                iter = 0)
fitSR <- lm.rrpp(HS ~ SVL + Region, data = snakeHS,
                 iter = 0)
fitSxR <- lm.rrpp(HS ~ SVL * Region, data = snakeHS,
                  iter = 0)

CVS <- looCV(fitS)
CVSR <- looCV(fitSR)
CVSxR <- looCV(fitSxR)

summary(CVS)
summary(CVSR)
summary(CVSxR)

looCV_table(CVS, CVSR, CVSxR)

Description

Function performs likelihood ratio tests on an lm.rrpp fit, using RRPP or FRPP. Likelihood ratio statistics are calculated for every random permutation, and the effect size is estimated from the distribution of random statistics. The likelihood ratio tests has some resemblance to MANOVA, especially using Wilks' lambda. Sums of squares and cross-products (SSCP) matrices are calculated over the random permutations of a lm.rrpp fit. SSCP matrices are computed, as are the inverse of R times H (invR.H), where R is a SSCP for the residuals or random effects and H is the difference between SSCP matrices of full and reduced models (see manova.update). From invR.H, Wilks lambda is first estimated, and the likelihood ratio stat is then estimated as -n * log(Wilks).

This function does one of two things. It either performs an update using manova.update, using Wilks' lambda as the test statistic, converting Wilks' lambda to likelihood ratio statistics or it uses the results from a previously performed update to calculate new statistics.

Usage

lr_test(fit, verbose = FALSE, ...)

Arguments

Author(s)

Michael Collyer

References

Adams, D. C., and M. L. Collyer. 2024. Extended phylogenetic regression models for comparing within-species patterns across the tree of life. Methods in Ecology and Evolution. In review.

Examples

# Body Shape Analysis (Multivariate) ----------------

## Not run: 
data(Pupfish)

# Although not recommended as a practice, this example will use only
# three principal components of body shape for demonstration.  
# A larger number of random permutations should also be used.

Pupfish$shape <- ordinate(Pupfish$coords)$x[, 1:3]


fit <- lm.rrpp(shape ~ log(CS) + Sex, SS.type = "I", 
data = Pupfish, print.progress = FALSE, iter = 499) 
summary(fit, formula = FALSE)
anova(fit) # ANOVA table

# MANOVA

fit.m <- manova.update(fit, print.progress = FALSE, tol = 0.001)
summary(fit.m, test = "Roy")
summary(fit.m, test = "Wilks")


# Likelihood Ratio Test

LRT <- lr_test(fit.m)
summary(LRT)


## End(Not run)

Description

Home range and body mass data for 49 species of mammals.

Details

The data set is from a study of mammal species home range size (squared kilometers) and its covariation with body mass or whether mammals are carnivores or ungulates (Garland et al., 1993). the data set also includes a covariance matrix based on phylogenetic relatedness among species, to account for the non-independence of observations. These data are suitable data for demonstrating generalized least squares (GLS) estimation, used in the book, Permutational Biometry (Collyer & Adams, 2025).

References

Garland, T., Dickerman, A. W., Janis, C. M., and Jones, J. A. (1993). Phylogenetic analysis of covariance by computer simulation. Systematic Biology, 42:265–292.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Phylogenetic tree data that can be used along with mammalHR.

Details

A class phylo object containing data for the phylogenetic tree that pertains to the mammalHR data set. Although these data are not required for the generalized least squares (GLS) estimation example, used in the book, Permutational Biometry (Collyer & Adams, 2025), they provide additional graphical support for the example.

References

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Function updates a lm.rrpp fit to add $MANOVA, which like $ANOVA, provides statistics or matrices typically associated with multivariate analysis of variance (MANOVA).

MANOVA statistics or sums of squares and cross-products (SSCP) matrices are calculated over the random permutations of a lm.rrpp fit. SSCP matrices are computed, as are the inverse of R times H (invR.H), where R is a SSCP for the residuals or random effects and H is the difference between SSCP matrices of full and reduced models (see below). From invR.H, MANOVA statistics are calculated, including Roy's maximum root (eigenvalue), Pillai trace, Hotelling-Lawley trace, and Wilks lambda (via summary.manova.lm.rrpp).

The manova.update to add $MANOVA to lm.rrpp fits requires more computation time than the $ANOVA statistics that are computed automatically in lm.rrpp. Generally, the same inferential conclusions will be found with either approach, when observations outnumber response variables. For high-dimensional data (more variables than observations) data are projected into a Euclidean space of appropriate dimensions (rank of residual covariance matrix). One can vary the tolerance for eigenvalue decay or specify the number of PCs, if a smaller set of PCs than the maximum is desired. This is advised if there is strong correlation among variables (the data space could be simplified to fewer dimensions), as spurious results are possible. Because distributions of MANOVA stats can be generated from the random permutations, there is no need to approximate F-values, like with parametric MANOVA. By restricting analysis to the real, positive eigenvalues calculated, all statistics can be calculated (but Wilks lambda, as a product but not a trace, might be less reliable as variable number approaches the number of observations).

ANOVA vs. MANOVA

Two SSCP matrices are calculated for each linear model effect, for every random permutation: R (Residuals or Random effects) and H, the difference between SSCPs for "full" and "reduced" models. (Full models contain and reduced models lack the effect tested; SSCPs are hypothesized to be the same under a null hypothesis, if there is no effect. The difference, H, would have a trace of 0 if the null hypothesis were true.) In RRPP, ANOVA and MANOVA correspond to two different ways to calculate statistics from R and H matrices.

ANOVA statistics are those that find the trace of R and H SSCP matrices before calculating subsequent statistics, including sums of squares (SS), mean squares (MS), and F-values. These statistics can be calculated with univariate data and provide univariate-like statistics for multivariate data. These statistics are dispersion measures only (covariances among variables do not contribute) and are the same as "distance-based" stats proposed by Goodall (1991) and Anderson (2001). MANOVA stats require multivariate data and are implicitly affected by variable covariances. For MANOVA, the inverse of R times H (invR.H) is first calculated for each effect, then eigen-analysis is performed on these matrix products. Multivariate statistics are calculated from the positive, real eigenvalues. In general, inferential conclusions will be similar with either approach, but effect sizes might differ.

Two important differences between manova.update and summary.manova (for lm objects) are that manova.update does not attempt to normalize residual SSCP matrices (unneeded for non-parametric statistical solutions) and (2) uses a generalized inverse of the residual SSCP, if needed, when the number of variables could render eigen-analysis problematic. This approach is consistent with covariance regularization methods that attempt to make covariance matrices positive-definite for calculating model likelihoods or multivariate statistics. If the number of observations far exceeds the number of response variables, observed statistics from manova.update and summary.manova will be quite similar. If the number of response variables approaches or exceeds the number of observations, manova.update statistics will be much more reliable.

ANOVA tables are generated by anova.lm.rrpp on lm.rrpp fits and MANOVA tables are generated by summary.manova.lm.rrpp, after running manova.update on lm.rrpp fits.

Currently, mixed model effects are only possible with $ANOVA statistics, not $MANOVA.

More detail is found in the vignette, ANOVA versus MANOVA.

Usage

manova.update(
  fit,
  error = NULL,
  tol = 1e-07,
  PC.no = NULL,
  print.progress = TRUE,
  verbose = NULL
)

Arguments

Value

An object of class lm.rrpp is updated to include class manova.lm.rrpp, and the object, $MANOVA, which includes

Author(s)

Michael Collyer

References

Goodall, C.R. 1991. Procrustes methods in the statistical analysis of shape. Journal of the Royal Statistical Society B 53:285-339.

Anderson MJ. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology 26: 32-46.

Examples

## Not run: 
# Body Shape Analysis (Multivariate) ----------------

data(Pupfish)

# Although not recommended as a practice, this example will use only
# three principal components of body shape for demonstration.  
# A larger number of random permutations should also be used.

Pupfish$shape <- ordinate(Pupfish$coords)$x[, 1:3]

Pupfish$logSize <- log(Pupfish$CS) 

fit <- lm.rrpp(shape ~ logSize + Sex, SS.type = "I", 
data = Pupfish, print.progress = FALSE, iter = 499) 
summary(fit, formula = FALSE)
anova(fit) # ANOVA table

# MANOVA

fit.m <- manova.update(fit, print.progress = FALSE, tol = 0.001)
summary(fit.m, test = "Roy")
summary(fit.m, test = "Pillai")

fit.m$MANOVA$eigs$obs # observed eigenvalues
fit.m$MANOVA$SSCP$obs # observed SSCP
fit.m$MANOVA$invR.H$obs # observed invR.H 

# Distributions of test statistics

summ.roy <- summary(fit.m, test = "Roy")
dens <- apply(summ.roy$rand.stats, 1, density)
par(mfcol = c(1, length(dens)))
for(i in 1:length(dens)) {
     plot(dens[[i]], xlab = "Roy max root", ylab = "Density",
     type = "l", main = names(dens)[[i]])
     abline(v = summ.roy$rand.stats[1, i], col = "red")
}
par(mfcol = c(1,1))

## End(Not run)

Description

model.frame.lm.rrpp returns the model frame constructed for an lm.rrpp object.

Usage

## S3 method for class 'lm.rrpp'
model.frame(formula, ...)

Arguments

Author(s)

Michael Collyer

Description

model.matrix.lm.rrpp returns the design matrix constructed for an lm.rrpp object.

Usage

## S3 method for class 'lm.rrpp'
model.matrix(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Simulated motion paths

Author(s)

Dean Adams

References

Adams, D. C., and M. L. Collyer. 2009. A general framework for the analysis of phenotypic trajectories in evolutionary studies. Evolution 63:1143-1154.

Description

Handle missing values in rrpp.data.frame objects

Usage

## S3 method for class 'rrpp.data.frame'
na.omit(object, ...)

Arguments

Author(s)

Michael Collyer

Examples

y <- matrix(rnorm(15), 5, 3)
x <- rnorm(5)
rdf <- rrpp.data.frame(x = x, y = y, d = dist(y))
rdf$x[1] <- NA # create missing data
rdf

ndf <- na.omit(rdf)
ndf

Description

A data set containing variables for phytoplanton biomass, phosphorous levels, trophic status, and locations (coordinates) of 70 temperate lakes.

Details

The data in this data set are samples from a larger study that compared phytoplankton biomass from temperate lakes in Alberta, Canada, with different levels of eutrophication (Loewen et al., 2020). These data are suitable for considering the spatial non-independence of observations in generalized least squares (GLS) estimation, for the example in the book, Permutational Biometry (Collyer & Adams, 2025).

References

Loewen, C. J. G., Wyatt, F. R., Mortimer, C. A., Vinebrooke, R. D., and Zurawell, R. W. (2020). Multiscale drivers of phytoplankton communities in north‐temperate lakes. Ecological Applications, 30(5):e02102.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Snout-to-vent-length (SVL) data for two species of salamanders

Details

The data set is a subset of data from two species of salamanders (Plethodon cinereus and P. hoffmani) in eastern North America (based on Adams and Rohlf, 2000; Adams, 2000), collected from the Smithsonian Museum (Washington DC, USA). Data include species, locality (whether in sympatry or allopatry), and SVL in mm. These data are suitable data for demonstrating single factor and factorial analysis of variance (ANOVA). (ANCOVA), used in the book, Permutational Biometry (Collyer & Adams, 2025).

Author(s)

Dean Adams

References

Adams, D. C. (2000). Divergence of trophic morphology and resource use among populations of Plethodon cinereus and P. hoffmani in Pennsylvania. In Bruce, R. C., Jaeger, R. G., and Houck, L. D., editors, The Biology of Plethodontid Salamanders, pages 383–394. Springer US.

Adams, D. C. and Rohlf, F. J. (2000). Ecological character displacement in Plethodon: Biomechanical differences found from a geometric morphometric study. Proceedings of the National Academy of Sciences, 97:4106–4111.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Data for 37 species of plethodontid salamanders. Variables include snout to vent length (SVL) as species size, tail length, head length, snout to eye length, body width, forelimb length, and hind limb length, all measured in mm. A grouping variable is also included for functional guild size. A variable for species names is also included. The data set also includes a phylogenetic covariance matrix based on a Brownian model of evolution, to assist in generalized least squares (GLS) estimation.

Details

The covariance matrix was estimated with the vcv.phylo function of the R package, ape, based on the tree described in Adams and Collyer (2018).

Author(s)

Michael Collyer and Dean Adams

References

Adams, D.C and Collyer, M.L. 2018. Multivariate phylogenetic anova: group-clade aggregation, biological challenges, and a refined permutation procedure. Evolution, 72: 1204-1215.

Description

This function produces a heat map for inter-subject variability, based on results from a measurement.error object. The function, interSubVar, must first be used on the measurement.error object to obtain variability statistics. This function use the image function to produce plots. It does little to manipulate such plots, but any argument for image can be manipulated here, as well as the graphical parameters that can be adjusted within image.

Usage

## S3 method for class 'interSubVar'
plot(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Plot Function for RRPP

Usage

## S3 method for class 'kcomp'
plot(x, axis1 = 1, axis2 = 2, flip = NULL, include.axes = TRUE, ...)

Arguments

Value

An object of class "plot.kcomp" is a list with components that can be used in other plot functions, such as the type of plot, points, a group factor, and other information depending on the plot parameters used.

Author(s)

Michael Collyer

Description

Plot Function for RRPP

Usage

## S3 method for class 'lm.rrpp'
plot(
  x,
  type = c("diagnostics", "regression", "PC"),
  resid.type = c("p", "n"),
  fitted.type = c("o", "t"),
  predictor = NULL,
  reg.type = c("PredLine", "RegScore"),
  ...
)

Arguments

Author(s)

Michael Collyer

References

Drake, A. G., and C. P. Klingenberg. 2008. The pace of morphological change: Historical transformation of skull shape in St Bernard dogs. Proc. R. Soc. B. 275:71-76.

Adams, D. C., and A. Nistri. 2010. Ontogenetic convergence and evolution of foot morphology in European cave salamanders (Family: Plethodontidae). BMC Evol. Biol. 10:1-10.

Examples

## Not run: 
# Univariate example
data(PlethMorph)
fitGLS <- lm.rrpp(TailLength ~ SVL, data = PlethMorph, Cov = PlethMorph$PhyCov, 
 print.progress = FALSE, iter = 0)
 
 par(mfrow = c(2, 2))
 plot(fitGLS)
 plot(fitGLS, resid.type = "n") # use normalized (transformed) residuals
 plot(fitGLS, resid.type = "n", fitted.type = "t") # use also transformed fitted values
 
 # Multivariate example
 
Y <- as.matrix(cbind(PlethMorph$TailLength,
PlethMorph$HeadLength,
PlethMorph$Snout.eye,
PlethMorph$BodyWidth,
PlethMorph$Forelimb,
PlethMorph$Hindlimb))
PlethMorph$Y <- Y
fitGLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph, 
Cov = PlethMorph$PhyCov,
print.progress = FALSE, iter = 0)

par(mfrow = c(2, 2))
 plot(fitGLSm)
 plot(fitGLSm, resid.type = "n") # use normalized (transformed) residuals
 plot(fitGLSm, resid.type = "n", fitted.type = "t") # use also transformed fitted values
 par(mfrow = c(1, 1))
 
## End(Not run)

Description

Plot Function for RRPP

Usage

## S3 method for class 'looCV'
plot(x, axis1 = 1, axis2 = 2, flip = NULL, ...)

Arguments

Author(s)

Michael Collyer

Description

This function produces multivariate signal-to-noise ratio plots for measurement.error objects. See the function, plot.interSubVar for plotting the inter-subject variability from a measurement.error object, after applying the function, interSubVar.

Usage

## S3 method for class 'measurement.error'
plot(
  x,
  separate.by.groups = TRUE,
  add.connectors = TRUE,
  add.labels = FALSE,
  use.std.vectors = FALSE,
  titles = NULL,
  add.legend = TRUE,
  ...
)

Arguments

Author(s)

Michael Collyer

Description

Plot Function for RRPP

Usage

## S3 method for class 'model.comparison'
plot(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Plot Function for RRPP

Usage

## S3 method for class 'ordinate'
plot(x, axis1 = 1, axis2 = 2, flip = NULL, include.axes = TRUE, ...)

Arguments

Value

An object of class "plot.ordinate" is a list with components that can be used in other plot functions, such as the type of plot, points, a group factor, and other information depending on the plot parameters used.

Author(s)

Michael Collyer

Description

Plot Function for RRPP

Usage

## S3 method for class 'predict.lm.rrpp'
plot(x, PC = FALSE, ellipse = FALSE, abscissa = NULL, label = TRUE, ...)

Arguments

Author(s)

Michael Collyer

Examples

# See lm.rrpp help file for examples.

Description

Function generates a principal component plot for trajectories

Usage

## S3 method for class 'trajectory.analysis'
plot(x, ...)

Arguments

Details

The function calculates and plots principal components of fitted values from lm.rrpp that are passed onto trajectory.analysis, and projects data onto them. This function is a set.up, and add.trajectories is needed to add trajectories to the plot. By having two stages of control, the plotting functions are more flexible. This function also returns plotting information that can be valuable for making individualized plots, if add.trajectories is not preferred.

Value

If an object is assigned, it will return:

Author(s)

Michael Collyer

References

Adams, D. C., and M. M. Cerney. 2007. Quantifying biomechanical motion using Procrustes motion analysis. J. Biomech. 40:437-444.

Adams, D. C., and M. L. Collyer. 2007. The analysis of character divergence along environmental gradients and other covariates. Evolution 61:510-515.

Adams, D. C., and M. L. Collyer. 2009. A general framework for the analysis of phenotypic trajectories in evolutionary studies. Evolution 63:1143-1154.

Collyer, M. L., and D. C. Adams. 2007. Analysis of two-state multivariate phenotypic change in ecological studies. Ecology 88:683-692.

Collyer, M. L., and D. C. Adams. 2013. Phenotypic trajectory analysis: comparison of shape change patterns in evolution and ecology. Hystrix 24: 75-83.

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 115:357-365.

See Also

plot.default and par

Examples

# See trajectory.analysis help file for examples

Description

A data set of body size for males and females form different lizard species in different geographic locations.

Details

These data are a subset of the data collected originally by Kaliontzopoulou et al.(2015). Variables include snout-to-vent-length (SVL), sex, species, and habitat preference. These data are suitable for comparing effect sizes in the book, Permutational Biometry (Collyer & Adams, 2025).

References

Kaliontzopoulou, A., Carretero, M. A., and Adams, D. C. (2015). Ecomorphological variation in male and female wall lizards and the macroevolution of sexual dimorphism in relation to habitat use. Journal of Evolutionary Biology, 28:80–94.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

A function that transforms a vector of statistics or data to be approximately normally distributed.

Usage

powerTrans(
  y,
  eps = 0.001,
  target = NULL,
  interp = FALSE,
  forceYJ = FALSE,
  useStDev = FALSE
)

Arguments

Details

This function uses a Box-Cox power transformation for a vector of values, (presumably statistics from RRPP) with the following conditions:

(1) If any values are 0 or negative, either a Yeo-Johnson transformation or a linear transformation of y - min(y) + 1e-7 is performed to optimize lambda with only positive values. There might be cases where one works better than the other. This function will determine which produces more normally distributed data, indicating which transformation was used.

(2) The lambda parameter is optimized over a range of -2 to +2, as is typical for most applications.

(3) The first (or targeted) value (assumed to be the observed value in a distribution of random values) is removed from the vector to perform optimization of lambda, but is returned to the vector for transformation. This assures that large, aberrant values, such as might be found via RRPP for large effects, do not have leverage for the optimization.

A Yeo-Johnson transformation can be forced, which for positive values merely shifts the distribution by adding 1. For most RRPP statistics (SS. MS, Rsq, F, d), this is not needed for effect sizes, as values are positive. However, for certain distributions, this might be a good idea for comparisons, because values might be strongly positive or bimodal in nature.

Using standard deviates (useStDev) is a good idea for statistics with aberrant ranges (like 0 to 1, or 1,000 to 2,000). If used, a Yeo-Johnson or linear transformation is implicit, since standardizing introduces negative values.

Value

A list containing the following:

Author(s)

Michael Collyer

References

Box, G. E. P. and Cox, D. R. (1964) An analysis of transformations (with discussion). Journal of the Royal Statistical Society B, 26, 211–252.

Yeo, I. and Johnson, R. (2000) A new family of power transformations to improve normality or symmetry. Biometrika, 87, 954-959.

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'anova.lm.rrpp'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'betaTest'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'coef.lm.rrpp'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'ICCstats'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'kcomp'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'lm.rrpp'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'looCV'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'lr_test'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'measurement.error'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Examples

## Not run: 

# Measurement error analysis on simulated data of fish shapes

data(fishy)

# Analysis unconcerned with groups 

ME1 <- measurement.error(
  Y = "coords",
  subjects = "subj",
  replicates = "reps",
  data = fishy)

anova(ME1)
ICCstats(ME1, subjects = "Subjects", with_in = "Systematic ME")
plot(ME1)

# Analysis concerned with groups 

ME2 <- measurement.error(
  Y = "coords",
  subjects = "subj",
  replicates = "reps",
  groups = "groups",
  data = fishy)

anova(ME2)
ICCstats(ME2, subjects = "Subjects", 
  with_in = "Systematic ME", groups = "groups")
P <- plot(ME2)
focusMEonSubjects(P, subjects = 18:20, shadow = TRUE)

## End(Not run)

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'model.comparison'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'ordinate'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'pairwise'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'pairwise.model.Z'
print(x, stats.table = TRUE, ...)

Arguments

Author(s)

Dean Adams and Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'predict.lm.rrpp'
print(x, PC = FALSE, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'summary.betaTest'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'summary.lm.rrpp'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'summary.manova.lm.rrpp'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'summary.ordinate'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'summary.pairwise'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'summary.trajectory.analysis'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'trajectory.analysis'
print(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Landmark data from Cyprinodon pecosensis body shapes, with indication of Sex and Population from which fish were sampled (Marsh or Sinkhole).

Details

These data were previously aligned with GPA. Centroid size (CS) is also provided. See the geomorph package for details.

Author(s)

Michael Collyer

References

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 113: doi:10.1038/hdy.2014.75.

Description

Landmark data from Cyprinodon pecosensis head shapes, with variables for sex, month and year sampled, locality, head size, and coordinates of landmarks for head shape, per specimen. These data are a subset of a larger data set.

Details

The variable, "coords", are data that were previously aligned with GPA. The variable, "headSize", is the Centroid size of each vector of coordinates. See the geomorph package for details.

Author(s)

Michael Collyer

References

Gilbert, M.C. 2016. Impacts of habitat fragmentation on the cranial morphology of a threatened desert fish (Cyprinodon pecosensis). Masters Thesis, Western Kentucky University.

Description

Simulated metabolic rate data from experimental pupfish

Details

These data are simulated data meant to mimic a subset of the data from Collyer & Stockwell (2004). Variables include experimental treatment (Control versus artificially infected), mass in grams, standard length in mm, and metabolic rate (oxygen consumption, mg O2 per L per hour). These data do not have a temporal component like the original data, and do not include aberrant or missing observations. These data are suitable data for demonstrating linear regression and analysis of covariance (ANCOVA), used in the book, Permutational Biometry (Collyer & Adams, 2025).

Author(s)

Michael Collyer

References

Collyer, M. L., & Stockwell, C. A. (2004). Experimental evidence for costs of parasitism for a threatened species, White Sands pupfish (Cyprinodon tularosa). Journal of Animal Ecology, 73(5), 821-830.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

A function to find the probability of values greater or lesser than target, from a vector of statistics presumably obtained in random permutations.

Usage

pval(s, target = NULL, greater = TRUE)

Arguments

Author(s)

Michael Collyer

Description

This function performs a QR decomposition (factorization) on a linear model design matrix (X) and returns useful results for subsequent analysis. This is intended as an internal function but can be used externally. Because base::qr and Matrix::qr have different options for QR algorithms, this function assures that results are consistent for other RRPP function use, whether X is a dense or sparse matrix.

Usage

QRforX(X, returnQ = TRUE, ...)

Arguments

Value

An object of class QR is a list containing the following:

Author(s)

Michael Collyer

Examples

## Simple Example
data(Pupfish)
fit <- lm.rrpp(coords ~ Pop, data = Pupfish, print.progress = FALSE)
QR <- QRforX(model.matrix(fit))
QR$Q
QR$R
QR$rank
QR$S4

## Not run, but one could get base::qr and Matrix::qr results as

# base::qr(as.matrix(QR$X))
# Matrix::qr(QR$X)

## Complex example

data("PupfishHeads")
fit <- suppressWarnings(lm.rrpp(headSize ~ sex + 
locality/year, data = PupfishHeads))
X <- model.matrix(fit)
dim(X) # Already reduced
colnames(X)
X <- model.matrix(terms(fit), fit$LM$data)
dim(X) # Retains redundant parameters
colnames(X)
QR <- QRforX(X)
dim(QR$X)
dim(QR$Q)
dim(QR$R)
QR$rank

Description

Returns random effects from a lmm.rrpp object.

Usage

## S3 method for class 'lmm.rrpp'
ranef(object, type = c("matrix", "list"), ...)

Arguments

Author(s)

Michael Collyer

Description

Reptile abundance data from several sites in Western Australia, originally collected to examine recovery after a large fire.

Details

The data set is from a study examining reptile abundance from several sites in Western Australia (Davis and Doherty, 2015). The original study examined the recovery of reptile communities in ten remnant localities following a large fire in 2009. Half of localities were in the burn zone and the remaining half were unaffected by the fire. Data were sampled over five years. These data are suitable data for demonstrating mixed linear models, used in the book, Permutational Biometry (Collyer & Adams, 2025).

References

Davis, R. A. and Doherty, T. S. (2015). Rapid recovery of an urban remnant reptile community following summer wildfire. PLOS ONE, 10(5):e0127925.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Extract residuals

Usage

## S3 method for class 'lm.rrpp'
residuals(object, ...)

Arguments

Author(s)

Michael Collyer

Examples

# See examples for lm.rrpp

Description

Create a data frame for lm.rrpp analysis, when covariance or distance matrices are used

Usage

rrpp.data.frame(...)

Arguments

Details

This function is not much different than data.frame but is more flexible to allow distance matrices and covariance matrices to be included. Essentially, this function creates a list, much like an object of class data.frame is also a list. However, rrpp.data.frame is less concerned with coercing the list into a matrix and more concerned with matching the number of observations (n). It is wise to use this function with any lm.rrpp analysis so that lm.rrpp does not have to search the global environment for needed data.

It is assumed that multiple data sets for the same subjects are in the same order.

See lm.rrpp for examples.

Author(s)

Michael Collyer

Examples

# Why use a rrpp.data.frame?
y <- matrix(rnorm(30), 10, 3)
x <- rnorm(10)
df <- data.frame(x = x, y = y)
df
rdf <- rrpp.data.frame(x = x, y = y)
rdf # looks more like a list

is.list(df)
is.list(rdf)

d <- dist(y) # distance matrix as data

# One can try this but it will result in an error
# df <- data.frame(df, d = d) 
rdf <- rrpp.data.frame(rdf, d = d) # works

fit <- lm.rrpp(d ~ x, data = rdf, iter = 99)
summary(fit)

Description

A data set containing variables for snake head size, snout-to-vent-length, sex, and geographic information.

Details

These data were collected from prairie rattlesnakes (Crotalus viridis viridis), from several geographic locations in the upper Midwest of the United States (Smith and Collyer, 2008). Variables include a spatial measure of head size (HS), snout-to-vent-length (SVL, cm), sex, and the site and region where 107 snakes were captured. These data are suitable for demonstrating model selection in the book, Permutational Biometry (Collyer & Adams, 2025).

References

Smith, M. T. and Collyer, M. L. (2008). The biology of rattlesnakes, chapter Regional variation and sexual dimorphism in head form of the prairie rattlesnake (Crotalus viridid viridis): Comparisons using new analytical techniques and collection methods, pages 79–90. Loma Linda University Press.

Collyer, M. L., & Adams D. C. (2025). Permutational Biometry: Volume 1, Univariate Data. Iowa State University Digital Press. NOT YET PUBLISHED.

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'anova.lm.rrpp'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'betaTest'
summary(object, confidence = 0.95, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'coef.lm.rrpp'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'ICCstats'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'kcomp'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'lm.rrpp'
summary(object, formula = TRUE, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'looCV'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'lr_test'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'manova.lm.rrpp'
summary(object, test = c("Roy", "Pillai", "Hotelling-Lawley", "Wilks"), ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'measurement.error'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'model.comparison'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'ordinate'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

See pairwise for further description.

Usage

## S3 method for class 'pairwise'
summary(
  object,
  stat.table = TRUE,
  test.type = c("dist", "stdist", "mdist", "VC", "DL", "var"),
  angle.type = c("rad", "deg"),
  confidence = 0.95,
  show.vectors = FALSE,
  p.adjust.method = "none",
  ...
)

Arguments

Details

The following summarize the test that can be performed:

Distance between vectors, "dist"

Vectors for LS means or slopes originate at the origin and point to some location, having both a magnitude and direction. A distance between two vectors is the inner-product of of the vector difference, i.e., the distance between their endpoints. For LS means, this distance is the difference between means. For multivariate slope vectors, this is the difference in location between estimated change for the dependent variables, per one-unit change of the covariate considered. For univariate slopes, this is the absolute difference between slopes.

Standardized distance between vectors, "stdist"

Same as the distance between vectors, but distances are divided by the model standard error (square-root of the trace of the residual covariance matrix). Pairwise tests with this statistic should be consistent with ANOVA results.

Mahalanobis distance between vectors, "mdist"

Similar to the standardized distance between vectors but the inverse of the residual covariance matrix is used in calculation of the distance, rather than dividing the Euclidean distance between means and dividing by the standard error. If the residual covariance matrix is singular, Mahalanobis distances will not be calculated. Pairwise tests with this statistic should be consistent with MANOVA results.

Vector correlation, "VC"

If LS mean or slope vectors are scaled to unit size, the vector correlation is the inner-product of the scaled vectors. The arccosine (acos) of this value is the angle between vectors, which can be expressed in radians or degrees. Vector correlation indicates the similarity of vector orientation, independent of vector length.

Difference between vector lengths, "DL"

If the length of a vector is an important attribute – e.g., the amount of multivariate change per one-unit change in a covariate – then the absolute value of the difference in vector lengths is a practical statistic to compare vector lengths, rather than the estimates the vectors make. Let d1 and d2 be the distances (length) of vectors. Then |d1 - d2| is a statistic that compares their lengths. For slope vectors, this is a comparison of rates. For comparison, if vectors are rates, "dist" finds the difference between estimates per unit change of, e.g., time, size, etc., which could be large, even for small rates of change, if vectors point in dissimilar directions. "DL" is a comparison of rates, irrespective of direction.

Variance, "var"

Vectors of residuals from a linear model indicate can express the distances of observed values from fitted values. Mean squared distances of values (variance), by group, can be used to measure the amount of dispersion around estimated values for groups. Absolute differences between variances are used as test statistics to compare mean dispersion of values among groups. Variance degrees of freedom equal n, the group size, rather than n-1, as the purpose is to compare mean dispersion in the sample. (Additionally, tests with one subject in a group are possible, or at least not a hindrance to the analysis.)

The argument, test.type is used to select one of the tests above. See pairwise for examples.

Notes for RRPP 0.6.2 and subsequent versions

In previous versions of pairwise, summary.pairwise had three test types: "dist", "VC", and "var". When one chose "dist", for LS mean vectors, the statistic was the inner-product of the vector difference. For slope vectors, "dist" returned the absolute value of the difference between vector lengths, which is "DL" in 0.6.2 and subsequent versions. This update uses the same calculation, irrespective of vector types. Generally, "DL" is the same as a contrast in rates for slope vectors, but might not have much meaning for LS means. Likewise, "dist" is the distance between vector endpoints, which might make more sense for LS means than slope vectors. Nevertheless, the user has more control over these decisions with version 0.6.2 and subsequent versions.

Notes for RRPP 2.0.4 and subsequent versions

The test types, standardized distance between vectors, "stdist", and Mahalanobis distances between vectors were added. The former simply divides the distance between vectors by model standard error (square-root of the trace of the residual covariance matrix). This is a multivariate generalization of a t-statistic for the difference between means. It is not the same as Hotelling squared-T-statistic, which requires incorporating the inverse of the residual covariance matrix in the calculation of the distance, a calculation that also requires a non-singular covariance matrix. However, the Mahalanobis distances are similar (and proportional) to the Hotelling squared-T-statistic. Pairwise tests using Mahalanobis distances represent a non-parametric analog to the parametric Hotelling squared-T test. Both tests should be better for GLS model fits compared to Euclidean distances between means, as the total sums of squares are more likely to vary across random permutations. In general, if ANOVA is performed a pairwise test with "stdist" is a good association; if MANOVA is performed, a pairwise test with "mdist" is a good association.

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'pairwise.model.Z'
summary(object, ...)

Arguments

Author(s)

Dean Adams and Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'predict.lm.rrpp'
summary(object, ...)

Arguments

Author(s)

Michael Collyer

Description

Print/Summary Function for RRPP

Usage

## S3 method for class 'trajectory.analysis'
summary(
  object,
  stat.table = TRUE,
  attribute = c("MD", "TC", "SD"),
  angle.type = c("rad", "deg"),
  confidence = 0.95,
  p.adjust.method = "none",
  show.trajectories = FALSE,
  ...
)

Arguments

Author(s)

Michael Collyer

Description

terms.lm.rrpp returns the terms constructed for an lm.rrpp object.

Usage

## S3 method for class 'lm.rrpp'
terms(x, ...)

Arguments

Author(s)

Michael Collyer

Description

Calculate vector correlations for a matrix (by rows). Used for pairwise comparisons.

Usage

vec.cor.matrix(M)

Arguments

Author(s)

Michael Collyer