# Using RRPP

#### 2025-02-18

# APPENDIX S2

This is a more detailed version of the section, “Analytical Demonstrations,” within the Methods in Ecology and Evolution article, “RRPP: An R package for fitting linear models to high-dimensional data using residual randomization,” by Michael L. Collyer and Dean C. Adams. The purpose of this vignette is to offer an opportunity for people reading that article to see code and results together. Some changes to code have been made, especially to enhance features that are difficult to introduce in the article, but the overall presentation should be basically the same. As the package has been updated, this document has also been slightly edited to make sure any functional changes have been illustrated.

(Note parallel processing is available on Unix systems for `lm.rrpp`; see `lm.rrpp` help page for description on its use. This is helpful for large data sets.)

Three data sets are provided in `RRPP` (with examples also in help pages), and are used here as examples for various analyses:

    Pupfish (Collyer et al., 2015)
    PupfishHeads (Gilbert, 2016)
    PlethMorph (Adams & Collyer, 2018)

The first two datasets contain landmark-based geometric morphometric data collected from museum samples of Pecos pupfish (Cyprinodon pecosensis), representing body shape and cranial morphology, respectively. Within both data sets, the `$coords` objects are matrices of Procrustes residuals obtained from generalized Procrustes analysis (GPA) of configurations of anatomical landmarks. For the purposes of this example, it is sufficient to recognize that Procrustes residuals embody a highly multivariate dataset representing shape (see the R package, geomorph). The third data set contains averaged linear measurements of 37 species of Plethodontid salamanders, plus a covariance matrix based on a Brownian model of evolution, given the phylogenetic relationship among the species (Adams & Collyer, 2018, for details).

Example: Pupfish Cranial Morphology and Mixed-Model ANOVA

In the first example we use a univariate dependent variable (head size, measured as the centroid size of the cranial landmark configuration; Bookstein, 1991) with a mixed-model design. The following code highlights the analytical steps, with results categorized in Fig. 1:

    library(RRPP)
    data("PupfishHeads")
    PupfishHeads$logHeadSize <- log(PupfishHeads$headSize)
    fit <- lm.rrpp(logHeadSize ~ sex + locality/year,
                   SS.type = "I", data = PupfishHeads,
                   print.progress = FALSE,
                   turbo = FALSE, verbose = TRUE)

    ## Warning:
    ## This is not an error!  It is a friendly warning.
    ##
    ## Because variables in the linear model are redundant,
    ## the linear model design has been truncated (via QR decomposition).
    ## Original X columns: 13
    ## Final X columns (rank): 9
    ## Check coefficients or degrees of freedom in ANOVA to see changes.
    ##
    ## Use suppressWarnings() to turn off these warnings.

    summary(fit)

    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 268
    ## Number of dependent variables: 1
    ## Data space dimensions: 1
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Full Model Analysis of Variance
    ##
    ##                     Df Residual Df       SS Residual SS       Rsq        F
    ## sex + locality/year  8         259 2.718175    10.94887 0.1988854 8.037445
    ##                     Z (from F) Pr(>F)
    ## sex + locality/year   5.097397  0.001
    ##
    ##
    ## Redundancy Analysis (PCA on fitted values and residuals)
    ##
    ##                Trace Proportion Rank
    ## Fitted    0.01018043  0.1988854    1
    ## Residuals 0.04100699  0.8011146    1
    ## Total     0.05118742  1.0000000    1
    ##
    ## Eigenvalues
    ##
    ##                  PC1
    ## Fitted    0.01018043
    ## Residuals 0.04100699
    ## Total     0.05118742

    anova(fit, effect.type = "F")

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##                Df      SS      MS     Rsq       F      Z Pr(>F)
    ## sex             1  0.6091 0.60911 0.04457 14.4087 2.7924  0.002 **
    ## locality        1  0.4712 0.47121 0.03448 11.1466 2.6477  0.002 **
    ## locality:year   6  1.6379 0.27298 0.11984  6.4574 4.2051  0.001 ***
    ## Residuals     259 10.9489 0.04227 0.80111
    ## Total         267 13.6670
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, turbo = FALSE,
    ##     SS.type = "I", data = PupfishHeads, print.progress = FALSE,
    ##     verbose = TRUE)

Of important note, we choose to log-transform head size and include it as a separate variable in the RRPP data frame, `PupfishHeads`. We accomplished this via the code above – rather than using `log(headSize)` - because downstream functions like `predict.lm.rrpp` work better without functions in the formula. ANOVA was performed using random distributions of F-statistics to calculate z-scores and P-values (but one could use alternative distributions - see the `RRPP` help page). The S3 Generic functions (`summary`, `anova`) return summaries that remind the user how random data were generated, the type of SS, and how z-scores were calculated. This particular ANOVA summary is a default that fails to consider the year fish were sampled as a random effect. A mixed-model ANOVA update can be performed by changing the expected mean-square (MS) error estimates in each F calculation:

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

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##                Df      SS      MS     Rsq       F      Z Pr(>F)
    ## sex             1  0.6091 0.60911 0.04457 14.4087 2.7924  0.002 **
    ## locality        1  0.4712 0.47121 0.03448  1.7262 0.7158  0.234
    ## locality:year   6  1.6379 0.27298 0.11984  6.4574 4.2051  0.001 ***
    ## Residuals     259 10.9489 0.04227 0.80111
    ## Total         267 13.6670
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, turbo = FALSE,
    ##     SS.type = "I", data = PupfishHeads, print.progress = FALSE,
    ##     verbose = TRUE)

This adjustment illustrates that the head size variation does not significantly differ between localities, with respect to the variation among sampling events. The anova function can also be used for multi-model inference.

    fit.sex <- lm.rrpp(logHeadSize ~ sex,
                       data = PupfishHeads,
                       print.progress = FALSE)
    fit.sex.loc<- lm.rrpp(logHeadSize ~ sex + locality,
                          data = PupfishHeads,
                          print.progress = FALSE)
    fit.sex.loc.year<- lm.rrpp(logHeadSize ~ sex + locality/year,
                               data = PupfishHeads,
                               print.progress = FALSE)

    ## Warning:
    ## This is not an error!  It is a friendly warning.
    ##
    ## Because variables in the linear model are redundant,
    ## the linear model design has been truncated (via QR decomposition).
    ## Original X columns: 13
    ## Final X columns (rank): 9
    ## Check coefficients or degrees of freedom in ANOVA to see changes.
    ##
    ## Use suppressWarnings() to turn off these warnings.

    anova(fit.sex, fit.sex.loc, fit.sex.loc.year, print.progress= FALSE)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Effect sizes (Z) based on F distributions
    ##
    ##                                   ResDf Df    RSS      SS      MS      Rsq
    ## logHeadSize ~ sex (Null)            266  1 13.058                 0.000000
    ## logHeadSize ~ sex + locality        265  1 12.587 0.47121 0.47121 0.034478
    ## logHeadSize ~ sex + locality/year   259  7 10.949 2.10907 0.30130 0.154318
    ## Total                               267    13.667
    ##                                        F      Z     P Pr(>F)
    ## logHeadSize ~ sex (Null)
    ## logHeadSize ~ sex + locality      9.9207 2.5237 0.004
    ## logHeadSize ~ sex + locality/year 7.1273 4.5980 0.001
    ## Total

One might wish to also look at individual model coefficients, and ascertain which have the largest effect:

    coef(fit, test = TRUE)

    ##
    ## Please be aware that the test of coefficients will be deprecated.
    ## Use betaTest, which is more flexible;
    ## e.g., betaTest(fit, null.method = 'terms')
    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 268
    ## Number of dependent variables: 1
    ## Data space dimensions: 1
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Statistics (distances) of coefficients with
    ## 95 percent confidence intervals, effect sizes,
    ## and probabilities of exceeding observed values,
    ## based on 1000 random permutations, using RRPP
    ##
    ##                            d.obs  UCL (95%)         Zd Pr(>d)
    ## sexMale               0.09579938 0.05535808  2.6097241  0.003
    ## localitySH            0.09164189 0.05788327  2.5097854  0.004
    ## localitySH:year1998   0.06402169 0.13652635  0.4170950  0.355
    ## localityLake:year1999 0.35681331 0.13286852  3.7362990  0.001
    ## localitySH:year1999   0.01328613 0.13259002 -1.1032264  0.856
    ## localitySH:year2000   0.07165589 0.10621936  0.9403694  0.196
    ## localityLake:year2001 0.23035259 0.13911422  2.5419582  0.003
    ## localityLake:year2002 0.31050534 0.13479886  3.3397163  0.001

This function produces a table much like `summary.lm` output, but with bootstrap-generated confidence intervals of coefficients.

It might be of interest to visualize model predictions for certain effects, holding constant other effects. For example, if we want to look at confidence intervals to compare male and female head sizes, holding constant the effects of locality and sampling period, we could do the following:

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

    ## Warning:
    ## This is not an error!  It is a friendly warning.
    ##
    ## Not all variables in model accounted for in newdata.
    ## Missing variables will be averaged from observed data for prediction.
    ##
    ## Use suppressWarnings to turn off these warnings.

    plot(sizePreds)

The plots are perfectly amenable (e.g., point type and color, line thickness, alternative labels, and additional text can be added or adjusted with typical `par` arguments).

    plot(sizePreds, pch = 21, cex = 3, bg = c(2,4), lwd = 2)

Finally, the SS type can be also toggled easily by refitting the model:

    fit2 <- lm.rrpp(logHeadSize ~ sex + locality/year,
                    SS.type = "II", data = PupfishHeads, print.progress = FALSE)

    ## Warning:
    ## This is not an error!  It is a friendly warning.
    ##
    ## Because variables in the linear model are redundant,
    ## the linear model design has been truncated (via QR decomposition).
    ## Original X columns: 13
    ## Final X columns (rank): 9
    ## Check coefficients or degrees of freedom in ANOVA to see changes.
    ##
    ## Use suppressWarnings() to turn off these warnings.

    fit3 <- lm.rrpp(logHeadSize ~ sex + locality/year,
                    SS.type = "III", data = PupfishHeads, print.progress = FALSE)

    ## Warning:
    ## This is not an error!  It is a friendly warning.
    ##
    ## Because variables in the linear model are redundant,
    ## the linear model design has been truncated (via QR decomposition).
    ## Original X columns: 13
    ## Final X columns (rank): 9
    ## Check coefficients or degrees of freedom in ANOVA to see changes.
    ##
    ## Use suppressWarnings() to turn off these warnings.

    anova(fit)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##                Df      SS      MS     Rsq       F      Z Pr(>F)
    ## sex             1  0.6091 0.60911 0.04457 14.4087 2.7924  0.002 **
    ## locality        1  0.4712 0.47121 0.03448 11.1466 2.6477  0.002 **
    ## locality:year   6  1.6379 0.27298 0.11984  6.4574 4.2051  0.001 ***
    ## Residuals     259 10.9489 0.04227 0.80111
    ## Total         267 13.6670
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, turbo = FALSE,
    ##     SS.type = "I", data = PupfishHeads, print.progress = FALSE,
    ##     verbose = TRUE)

    anova(fit2)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type II
    ## Effect sizes (Z) based on F distributions
    ##
    ##                Df      SS      MS     Rsq       F      Z Pr(>F)
    ## sex             1  0.5554 0.55537 0.04064 13.1376 2.6737  0.002 **
    ## locality:year   6  1.6379 0.27298 0.11984  6.4574 4.2051  0.001 ***
    ## Residuals     259 10.9489 0.04227 0.80111
    ## Total         267 13.6670
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, SS.type = "II",
    ##     data = PupfishHeads, print.progress = FALSE)

    anova(fit3)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type III
    ## Effect sizes (Z) based on F distributions
    ##
    ##                Df      SS      MS     Rsq       F      Z Pr(>F)
    ## sex             1  0.5554 0.55537 0.04064 13.1376 2.6737  0.002 **
    ## locality        1  0.0573 0.05733 0.00419  1.3562 0.7005  0.255
    ## locality:year   6  1.6379 0.27298 0.11984  6.4574 4.2051  0.001 ***
    ## Residuals     259 10.9489 0.04227 0.80111
    ## Total         267 13.6670
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, SS.type = "III",
    ##     data = PupfishHeads, print.progress = FALSE)

Example: Pupfish Body Shape and High-Dimensional Data

In the second example, we highlight the `RRPP` ability to efficiently handle large data computations. For this demonstration, a 54(n) × 112 (p) matrix of Procrustes residuals are the data. In every one of the 1,000 random permutations, RRPP shuffles residual vectors the same way for four different null models, estimates coefficients for four different full models, estimates the SS as the difference between residual SS (RSS) for four null-full model comparisons, and calculates the total SS, before calculating MS, R2, F, Cohen’s f2, and Euclidean distances of coefficient vectors across all 1,000 permutations. This process, plus packaging of results, took approximately 0.5 seconds on a notebook computer, without any parallel processing. In the second example, the initial steps are quite the same as the first example:

    data(Pupfish)
    Pupfish$logSize <- log(Pupfish$CS)
    fit <- lm.rrpp(coords ~ logSize + Sex*Pop, SS.type = "I",
                   data = Pupfish, print.progress = FALSE,
                   turbo = FALSE, verbose = TRUE)
    summary(fit, formula = FALSE)

    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 54
    ## Number of dependent variables: 112
    ## Data space dimensions: 53
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Full Model Analysis of Variance
    ##
    ##     Df Residual Df         SS Residual SS       Rsq        F Z (from F) Pr(>F)
    ## fit  4          49 0.03224913  0.02408373 0.5724746 16.40327   7.783173  0.001
    ##
    ##
    ## Redundancy Analysis (PCA on fitted values and residuals)
    ##
    ##                  Trace Proportion Rank
    ## Fitted    0.0006084742  0.5724746    4
    ## Residuals 0.0004544100  0.4275254   49
    ## Total     0.0010628843  1.0000000   53
    ##
    ## Eigenvalues
    ##
    ##                    PC1          PC2          PC3          PC4          PC5
    ## Fitted    0.0003849870 0.0001764119 0.0000308588 0.0000162165
    ## Residuals 0.0001619045 0.0000683959 0.0000504732 0.0000375901 0.0000204706
    ## Total     0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
    ##                    PC6          PC7          PC8          PC9         PC10
    ## Fitted
    ## Residuals 0.0000158227 0.0000145861 0.0000120066 0.0000083853 0.0000074828
    ## Total     0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
    ##                   PC11         PC12         PC13         PC14         PC15
    ## Fitted
    ## Residuals 0.0000068974 0.0000062331 0.0000054188 0.0000045485 0.0000043871
    ## Total     0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
    ##                   PC16         PC17         PC18         PC19         PC20
    ## Fitted
    ## Residuals 0.0000032928 0.0000031713 0.0000027759 0.0000024382 0.0000023968
    ## Total     0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
    ##                   PC21         PC22         PC23         PC24         PC25
    ## Fitted
    ## Residuals 0.0000020800 0.0000016410 0.0000013892 0.0000012122 0.0000010647
    ## Total     0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
    ##                   PC26         PC27         PC28         PC29         PC30
    ## Fitted
    ## Residuals 0.0000009268 0.0000008974 0.0000008083 0.0000007337 0.0000006036
    ## Total     0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
    ##                   PC31         PC32         PC33         PC34         PC35
    ## Fitted
    ## Residuals 0.0000005856 0.0000004746 0.0000004236 0.0000004042 0.0000003381
    ## Total     0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
    ##                   PC36         PC37         PC38         PC39         PC40
    ## Fitted
    ## Residuals 0.0000003127 0.0000002573 0.0000002486 0.0000002226 0.0000002059
    ## Total     0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
    ##                   PC41         PC42         PC43         PC44         PC45
    ## Fitted
    ## Residuals 0.0000001783 0.0000001509 0.0000001282 0.0000001152 0.0000000896
    ## Total     0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
    ##                   PC46         PC47         PC48         PC49         PC50
    ## Fitted
    ## Residuals 0.0000000737 0.0000000708 0.0000000488 0.0000000468
    ## Total     0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
    ##                   PC51         PC52         PC53
    ## Fitted
    ## Residuals
    ## Total     0.0000000525 0.0000000447 0.0000000395

    anova(fit)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df       SS        MS     Rsq       F      Z Pr(>F)
    ## logSize    1 0.014019 0.0140193 0.24886 28.5232 4.9982  0.001 ***
    ## Sex        1 0.007615 0.0076151 0.13518 15.4933 4.6798  0.001 ***
    ## Pop        1 0.007661 0.0076606 0.13599 15.5860 4.3464  0.001 ***
    ## Sex:Pop    1 0.002954 0.0029542 0.05244  6.0105 3.4054  0.002 **
    ## Residuals 49 0.024084 0.0004915 0.42753
    ## Total     53 0.056333
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = coords ~ logSize + Sex * Pop, turbo = FALSE, SS.type = "I",
    ##     data = Pupfish, print.progress = FALSE, verbose = TRUE)

    coef(fit, test = TRUE)

    ##
    ## Please be aware that the test of coefficients will be deprecated.
    ## Use betaTest, which is more flexible;
    ## e.g., betaTest(fit, null.method = 'terms')
    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 54
    ## Number of dependent variables: 112
    ## Data space dimensions: 53
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Statistics (distances) of coefficients with
    ## 95 percent confidence intervals, effect sizes,
    ## and probabilities of exceeding observed values,
    ## based on 1000 random permutations, using RRPP
    ##
    ##                       d.obs  UCL (95%)       Zd Pr(>d)
    ## logSize          0.11101234 0.04871922 4.148348  0.001
    ## SexM             0.02755503 0.01300605 4.120183  0.001
    ## PopSinkhole      0.02895656 0.01256411 4.121852  0.001
    ## SexM:PopSinkhole 0.03002431 0.01772916 3.469951  0.002

ANOVA results reveal that after accounting for body size allometry, not only are there significant inter-population differences in body shape and sexual dimorphism in body shape, but sexual dimorphism also significantly varies between the two populations. A fuller evaluation of these results is available in Collyer et al. (2015). It is worth taking a moment to realize why this analysis is a valuable alternative to a parametric M-ANOVA. The following code attempts to perform a parametric M-ANOVA:

    fit$LM$data$coords <- Pupfish$coords
    fit.par <- lm(fit$call$f1, data = fit$LM$data)
    all.equal(fit$LM$coefficients, fit.par$coefficients)

    ## [1] TRUE

    summary(manova(fit.par))

    ## Error in summary.manova(manova(fit.par)): residuals have rank 49 < 112

Although both functions return the same coefficients, the error summarizing the attempted parametric M-ANOVA clearly indicates the limitation of having residual degrees of freedom (rank) lower than the number of variables. Returning to the `lm.rrpp` fit, which does not suffer this problem, we can look at the precision of group mean estimation, accounting for allometric shape variation, by doing the following:

    shapeDF <- expand.grid(Sex = levels(Pupfish$Sex), Pop = levels(Pupfish$Pop))
    rownames(shapeDF) <- paste(shapeDF$Sex, shapeDF$Pop, sep = ".")
    shapePreds <- predict(fit, shapeDF, confidence = 0.95)

    ## Warning:
    ## This is not an error!  It is a friendly warning.
    ##
    ## Not all variables in model accounted for in newdata.
    ## Missing variables will be averaged from observed data for prediction.
    ##
    ## Use suppressWarnings to turn off these warnings.

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

    plot(shapePreds, PC = TRUE, ellipse = TRUE,
         pch = 19, col = 1:NROW(shapeDF)) # with added par arguments

    groups <- interaction(Pupfish$Sex, Pupfish$Pop)
    plot(fit, type = "PC") # generic

    plot(fit, type = "PC", pch = 19, col = groups) # with added par arguments

These plots differ as there is a rotational difference between the covariance matrices estimated with 4 predicted and 54 fitted values. Additionally, the former illustrates prediction precision and the latter sample dispersion. Both functions allow passing `par` arguments to the plot as well as saving plot data for more advanced plotting. The following is a regression-type plot:

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

The function, pairwise, can be used to test pairwise differences between least-squares means with:

    PWT <- pairwise(fit, groups = interaction(Pupfish$Sex, Pupfish$Pop))
    summary(PWT, confidence = 0.95)

    ##
    ## Pairwise comparisons
    ##
    ## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
    ##
    ## RRPP: 1000 permutations
    ##
    ## LS means:
    ## Vectors hidden (use show.vectors = TRUE to view)
    ##
    ## Pairwise distances between means, plus statistics
    ##                                d  UCL (95%)          Z Pr > d
    ## F.Marsh:M.Marsh       0.03579214 0.03278645  2.2202630  0.015
    ## F.Marsh:F.Sinkhole    0.03639758 0.03601481  1.7807393  0.044
    ## F.Marsh:M.Sinkhole    0.03809715 0.04432235 -0.6390591  0.735
    ## M.Marsh:F.Sinkhole    0.03681731 0.04779814 -0.2465676  0.601
    ## M.Marsh:M.Sinkhole    0.02694896 0.03646576 -0.9484513  0.830
    ## F.Sinkhole:M.Sinkhole 0.01939739 0.03295865 -1.7845209  0.962

Much like the tukeyHSD function in the R stats package, `pairwise` will generate tables with confidence intervals and P-values for the pairwise statistic, Euclidean distance between least-squares means. This function could also be used for pairwise comparison of slopes in analysis of covariance (ANCOVA) designs.

    fit2 <- lm.rrpp(coords ~ logSize * Sex * Pop, SS.type = "I",
                    data = Pupfish, print.progress = FALSE, iter = 999)
    summary(fit2, formula = FALSE)

    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 54
    ## Number of dependent variables: 112
    ## Data space dimensions: 53
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Full Model Analysis of Variance
    ##
    ##      Df Residual Df         SS Residual SS       Rsq        F Z (from F) Pr(>F)
    ## fit2  7          46 0.03425362  0.02207925 0.6080574 10.19488   9.289796  0.001
    ##
    ##
    ## Redundancy Analysis (PCA on fitted values and residuals)
    ##
    ##                  Trace Proportion Rank
    ## Fitted    0.0006462947  0.6080575    7
    ## Residuals 0.0004165896  0.3919426   46
    ## Total     0.0010628843  1.0000000   53
    ##
    ## Eigenvalues
    ##
    ##                    PC1          PC2          PC3          PC4          PC5
    ## Fitted    0.0003957338 0.0001850501 0.0000318297 0.0000187072 0.0000081143
    ## Residuals 0.0001395372 0.0000679693 0.0000501429 0.0000336509 0.0000189840
    ## Total     0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
    ##                    PC6          PC7          PC8          PC9         PC10
    ## Fitted    0.0000038627 0.0000029970
    ## Residuals 0.0000153088 0.0000133860 0.0000109149 0.0000081753 0.0000072355
    ## Total     0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
    ##                   PC11         PC12         PC13         PC14         PC15
    ## Fitted
    ## Residuals 0.0000066482 0.0000053575 0.0000051028 0.0000040219 0.0000034070
    ## Total     0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
    ##                   PC16         PC17         PC18         PC19         PC20
    ## Fitted
    ## Residuals 0.0000030820 0.0000029693 0.0000025181 0.0000024262 0.0000019439
    ## Total     0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
    ##                   PC21         PC22         PC23         PC24         PC25
    ## Fitted
    ## Residuals 0.0000018143 0.0000016024 0.0000013239 0.0000011070 0.0000009227
    ## Total     0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
    ##                   PC26         PC27         PC28         PC29         PC30
    ## Fitted
    ## Residuals 0.0000008477 0.0000008014 0.0000007148 0.0000006218 0.0000005587
    ## Total     0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
    ##                   PC31         PC32         PC33         PC34         PC35
    ## Fitted
    ## Residuals 0.0000005125 0.0000004234 0.0000003470 0.0000003184 0.0000002807
    ## Total     0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
    ##                   PC36         PC37         PC38         PC39         PC40
    ## Fitted
    ## Residuals 0.0000002721 0.0000002379 0.0000002230 0.0000001699 0.0000001561
    ## Total     0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
    ##                   PC41         PC42         PC43         PC44         PC45
    ## Fitted
    ## Residuals 0.0000001439 0.0000001208 0.0000000942 0.0000000763 0.0000000673
    ## Total     0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
    ##                   PC46         PC47         PC48         PC49         PC50
    ## Fitted
    ## Residuals 0.0000000496
    ## Total     0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
    ##                   PC51         PC52         PC53
    ## Fitted
    ## Residuals
    ## Total     0.0000000525 0.0000000447 0.0000000395

    anova(fit, fit2, print.progress = FALSE)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Effect sizes (Z) based on F distributions
    ##
    ##                                     ResDf Df      RSS        SS         MS
    ## coords ~ logSize + Sex * Pop (Null)    49  1 0.024084
    ## coords ~ logSize * Sex * Pop           46  3 0.022079 0.0020045 0.00066816
    ## Total                                  53    0.056333
    ##                                          Rsq      F      Z     P Pr(>F)
    ## coords ~ logSize + Sex * Pop (Null) 0.000000
    ## coords ~ logSize * Sex * Pop        0.035583 1.3921 1.1451 0.136
    ## Total

    PW2 <- pairwise(fit2, fit.null = fit, groups = groups,
                    covariate = Pupfish$logSize, print.progress = FALSE)
    PW2

    ##
    ## Pairwise comparisons
    ##
    ## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
    ##
    ## RRPP: 1000 permutations

    summary(PW2, confidence = 0.95,
            test.type = "dist") # distances between slope vector lengths

    ##
    ## Pairwise comparisons
    ##
    ## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
    ##
    ## RRPP: 1000 permutations
    ##
    ## Slopes (vectors of variate change per one unit of covariate
    ##         change, by group):
    ## Vectors hidden (use show.vectors = TRUE to view)
    ##
    ## Pairwise distances between slope vector
    ##             (end-points), plus statistics
    ##                                d UCL (95%)          Z Pr > d
    ## F.Marsh:M.Marsh       0.08279940 0.1640452 -0.9689342  0.829
    ## F.Marsh:F.Sinkhole    0.07848121 0.1621176 -0.9763904  0.838
    ## F.Marsh:M.Sinkhole    0.12322729 0.1610484  0.7139190  0.236
    ## M.Marsh:F.Sinkhole    0.06418279 0.1095190 -0.3981426  0.652
    ## M.Marsh:M.Sinkhole    0.10863856 0.1132783  1.5445168  0.063
    ## F.Sinkhole:M.Sinkhole 0.11706696 0.1001323  2.2066562  0.015

    summary(PW2, confidence = 0.95,
            test.type = "dist", stat.table = FALSE)

    ##
    ## Pairwise comparisons
    ##
    ## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
    ##
    ## RRPP: 1000 permutations
    ##
    ## Slopes (vectors of variate change per one unit of covariate
    ##         change, by group):
    ## Vectors hidden (use show.vectors = TRUE to view)
    ##
    ## Pairwise distances between slope vector (end-points
    ##               F.Marsh    M.Marsh F.Sinkhole M.Sinkhole
    ## F.Marsh    0.00000000 0.08279940 0.07848121  0.1232273
    ## M.Marsh    0.08279940 0.00000000 0.06418279  0.1086386
    ## F.Sinkhole 0.07848121 0.06418279 0.00000000  0.1170670
    ## M.Sinkhole 0.12322729 0.10863856 0.11706696  0.0000000
    ##
    ## Pairwise 95% Upper confidence limits between slopes
    ##              F.Marsh   M.Marsh F.Sinkhole M.Sinkhole
    ## F.Marsh    0.0000000 0.1640452  0.1621176  0.1610484
    ## M.Marsh    0.1640452 0.0000000  0.1095190  0.1132783
    ## F.Sinkhole 0.1621176 0.1095190  0.0000000  0.1001323
    ## M.Sinkhole 0.1610484 0.1132783  0.1001323  0.0000000
    ##
    ## Pairwise effect sizes (Z) between slopes
    ##               F.Marsh    M.Marsh F.Sinkhole M.Sinkhole
    ## F.Marsh     0.0000000 -0.9689342 -0.9763904   0.713919
    ## M.Marsh    -0.9689342  0.0000000 -0.3981426   1.544517
    ## F.Sinkhole -0.9763904 -0.3981426  0.0000000   2.206656
    ## M.Sinkhole  0.7139190  1.5445168  2.2066562   0.000000
    ##
    ## Pairwise P-values between slopes
    ##            F.Marsh M.Marsh F.Sinkhole M.Sinkhole
    ## F.Marsh      1.000   0.829      0.838      0.236
    ## M.Marsh      0.829   1.000      0.652      0.063
    ## F.Sinkhole   0.838   0.652      1.000      0.015
    ## M.Sinkhole   0.236   0.063      0.015      1.000

    summary(PW2, confidence = 0.95,
            test.type = "VC",
            angle.type = "deg") # correlation between slope vectors (and angles)

    ##
    ## Pairwise comparisons
    ##
    ## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
    ##
    ## RRPP: 1000 permutations
    ##
    ## Slopes (vectors of variate change per one unit of covariate
    ##         change, by group):
    ## Vectors hidden (use show.vectors = TRUE to view)
    ##
    ## Pairwise statistics based on slopes vector correlations (r)
    ##           and angles, acos(r)
    ## The null hypothesis is that r = 1 (parallel vectors).
    ## This null hypothesis is better treated as the angle
    ##           between vectors = 0
    ##                               r    angle UCL (95%)           Z Pr > angle
    ## F.Marsh:M.Marsh       0.6139591 52.12367  83.58445 -0.13196270      0.549
    ## F.Marsh:F.Sinkhole    0.6318267 50.81498  83.57546 -0.07310658      0.527
    ## F.Marsh:M.Sinkhole    0.4461018 63.50614  81.89039  0.66791084      0.273
    ## M.Marsh:F.Sinkhole    0.7175129 44.15048  63.60905  0.24275646      0.391
    ## M.Marsh:M.Sinkhole    0.5629975 55.73665  65.19951  1.04510764      0.154
    ## F.Sinkhole:M.Sinkhole 0.4627734 62.43378  60.22153  1.77232308      0.037

Because the Procrustes residuals are projected into a Euclidean tangent space (see geomorph function, gpagen; Adams et al., 2017), this analysis could be performed with an object of class `dist` (values from lower half of a distance matrix) representing the inter-specimen shape (Euclidean) distances, using the following code:

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

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

    summary(fitD)

    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 54
    ## Number of dependent variables: 53
    ## Data space dimensions: 53
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Full Model Analysis of Variance
    ##
    ##                     Df Residual Df         SS Residual SS       Rsq        F
    ## logSize + Sex * Pop  4          49 0.03224913  0.02408373 0.5724746 16.40327
    ##                     Z (from F) Pr(>F)
    ## logSize + Sex * Pop   7.783173  0.001
    ##
    ##
    ## Redundancy Analysis (PCA on fitted values and residuals)
    ##
    ##                  Trace Proportion Rank
    ## Fitted    0.0006084742  0.5724746    4
    ## Residuals 0.0004544100  0.4275254   49
    ## Total     0.0010628843  1.0000000   53
    ##
    ## Eigenvalues
    ##
    ##                    PC1          PC2          PC3          PC4          PC5
    ## Fitted    0.0003849870 0.0001764119 0.0000308588 0.0000162165
    ## Residuals 0.0001619045 0.0000683959 0.0000504732 0.0000375901 0.0000204706
    ## Total     0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
    ##                    PC6          PC7          PC8          PC9         PC10
    ## Fitted
    ## Residuals 0.0000158227 0.0000145861 0.0000120066 0.0000083853 0.0000074828
    ## Total     0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
    ##                   PC11         PC12         PC13         PC14         PC15
    ## Fitted
    ## Residuals 0.0000068974 0.0000062331 0.0000054188 0.0000045485 0.0000043871
    ## Total     0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
    ##                   PC16         PC17         PC18         PC19         PC20
    ## Fitted
    ## Residuals 0.0000032928 0.0000031713 0.0000027759 0.0000024382 0.0000023968
    ## Total     0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
    ##                   PC21         PC22         PC23         PC24         PC25
    ## Fitted
    ## Residuals 0.0000020800 0.0000016410 0.0000013892 0.0000012122 0.0000010647
    ## Total     0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
    ##                   PC26         PC27         PC28         PC29         PC30
    ## Fitted
    ## Residuals 0.0000009268 0.0000008974 0.0000008083 0.0000007337 0.0000006036
    ## Total     0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
    ##                   PC31         PC32         PC33         PC34         PC35
    ## Fitted
    ## Residuals 0.0000005856 0.0000004746 0.0000004236 0.0000004042 0.0000003381
    ## Total     0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
    ##                   PC36         PC37         PC38         PC39         PC40
    ## Fitted
    ## Residuals 0.0000003127 0.0000002573 0.0000002486 0.0000002226 0.0000002059
    ## Total     0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
    ##                   PC41         PC42         PC43         PC44         PC45
    ## Fitted
    ## Residuals 0.0000001783 0.0000001509 0.0000001282 0.0000001152 0.0000000896
    ## Total     0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
    ##                   PC46         PC47         PC48         PC49         PC50
    ## Fitted
    ## Residuals 0.0000000737 0.0000000708 0.0000000488 0.0000000468
    ## Total     0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
    ##                   PC51         PC52         PC53
    ## Fitted
    ## Residuals
    ## Total     0.0000000525 0.0000000447 0.0000000395

    summary(fit)

    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 54
    ## Number of dependent variables: 112
    ## Data space dimensions: 53
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Full Model Analysis of Variance
    ##
    ##                     Df Residual Df         SS Residual SS       Rsq        F
    ## logSize + Sex * Pop  4          49 0.03224913  0.02408373 0.5724746 16.40327
    ##                     Z (from F) Pr(>F)
    ## logSize + Sex * Pop   7.783173  0.001
    ##
    ##
    ## Redundancy Analysis (PCA on fitted values and residuals)
    ##
    ##                  Trace Proportion Rank
    ## Fitted    0.0006084742  0.5724746    4
    ## Residuals 0.0004544100  0.4275254   49
    ## Total     0.0010628843  1.0000000   53
    ##
    ## Eigenvalues
    ##
    ##                    PC1          PC2          PC3          PC4          PC5
    ## Fitted    0.0003849870 0.0001764119 0.0000308588 0.0000162165
    ## Residuals 0.0001619045 0.0000683959 0.0000504732 0.0000375901 0.0000204706
    ## Total     0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
    ##                    PC6          PC7          PC8          PC9         PC10
    ## Fitted
    ## Residuals 0.0000158227 0.0000145861 0.0000120066 0.0000083853 0.0000074828
    ## Total     0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
    ##                   PC11         PC12         PC13         PC14         PC15
    ## Fitted
    ## Residuals 0.0000068974 0.0000062331 0.0000054188 0.0000045485 0.0000043871
    ## Total     0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
    ##                   PC16         PC17         PC18         PC19         PC20
    ## Fitted
    ## Residuals 0.0000032928 0.0000031713 0.0000027759 0.0000024382 0.0000023968
    ## Total     0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
    ##                   PC21         PC22         PC23         PC24         PC25
    ## Fitted
    ## Residuals 0.0000020800 0.0000016410 0.0000013892 0.0000012122 0.0000010647
    ## Total     0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
    ##                   PC26         PC27         PC28         PC29         PC30
    ## Fitted
    ## Residuals 0.0000009268 0.0000008974 0.0000008083 0.0000007337 0.0000006036
    ## Total     0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
    ##                   PC31         PC32         PC33         PC34         PC35
    ## Fitted
    ## Residuals 0.0000005856 0.0000004746 0.0000004236 0.0000004042 0.0000003381
    ## Total     0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
    ##                   PC36         PC37         PC38         PC39         PC40
    ## Fitted
    ## Residuals 0.0000003127 0.0000002573 0.0000002486 0.0000002226 0.0000002059
    ## Total     0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
    ##                   PC41         PC42         PC43         PC44         PC45
    ## Fitted
    ## Residuals 0.0000001783 0.0000001509 0.0000001282 0.0000001152 0.0000000896
    ## Total     0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
    ##                   PC46         PC47         PC48         PC49         PC50
    ## Fitted
    ## Residuals 0.0000000737 0.0000000708 0.0000000488 0.0000000468
    ## Total     0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
    ##                   PC51         PC52         PC53
    ## Fitted
    ## Residuals
    ## Total     0.0000000525 0.0000000447 0.0000000395

    anova(fitD)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df       SS        MS     Rsq       F      Z Pr(>F)
    ## logSize    1 0.014019 0.0140193 0.24886 28.5232 4.9982  0.001 ***
    ## Sex        1 0.007615 0.0076151 0.13518 15.4933 4.6798  0.001 ***
    ## Pop        1 0.007661 0.0076606 0.13599 15.5860 4.3464  0.001 ***
    ## Sex:Pop    1 0.002954 0.0029542 0.05244  6.0105 3.4054  0.002 **
    ## Residuals 49 0.024084 0.0004915 0.42753
    ## Total     53 0.056333
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = D ~ logSize + Sex * Pop, SS.type = "I", data = Pupfish,
    ##     print.progress = FALSE)

    anova(fit)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df       SS        MS     Rsq       F      Z Pr(>F)
    ## logSize    1 0.014019 0.0140193 0.24886 28.5232 4.9982  0.001 ***
    ## Sex        1 0.007615 0.0076151 0.13518 15.4933 4.6798  0.001 ***
    ## Pop        1 0.007661 0.0076606 0.13599 15.5860 4.3464  0.001 ***
    ## Sex:Pop    1 0.002954 0.0029542 0.05244  6.0105 3.4054  0.002 **
    ## Residuals 49 0.024084 0.0004915 0.42753
    ## Total     53 0.056333
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = coords ~ logSize + Sex * Pop, turbo = FALSE, SS.type = "I",
    ##     data = Pupfish, print.progress = FALSE, verbose = TRUE)

The ANOVA results with either method are exactly the same.

Example: Plethodontid Morphology, Phylogenetics, and GLS Estimation

In the third example, we highlight GLS estimation. The following code creates two `lm.rrpp` fits using OLS and GLS, respectively, and evaluates them as in previous examples:

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

    anova(fitOLS)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df     SS     MS     Rsq      F      Z Pr(>F)
    ## SVL        1 3707.9 3707.9 0.76849 116.18 5.7522  0.001 ***
    ## Residuals 35 1117.0   31.9 0.23151
    ## Total     36 4824.9
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = TailLength ~ SVL, turbo = FALSE, data = PlethMorph,
    ##     print.progress = FALSE, verbose = TRUE)

    anova(fitGLS)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Generalized Least-Squares (via OLS projection)
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df     SS      MS     Rsq      F      Z Pr(>F)
    ## SVL        1 189.34 189.343 0.34259 18.239 3.2047  0.001 ***
    ## Residuals 35 363.33  10.381 0.65741
    ## Total     36 552.68
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = TailLength ~ SVL, turbo = FALSE, data = PlethMorph,
    ##     Cov = PlethMorph$PhyCov, print.progress = FALSE, verbose = TRUE)

    coef(fitOLS, test = TRUE)

    ##
    ## Please be aware that the test of coefficients will be deprecated.
    ## Use betaTest, which is more flexible;
    ## e.g., betaTest(fit, null.method = 'terms')
    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 37
    ## Number of dependent variables: 1
    ## Data space dimensions: 1
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Statistics (distances) of coefficients with
    ## 95 percent confidence intervals, effect sizes,
    ## and probabilities of exceeding observed values,
    ## based on 1000 random permutations, using RRPP
    ##
    ##         d.obs UCL (95%)       Zd Pr(>d)
    ## SVL 0.9681511 0.3597148 3.907365  0.001

    coef(fitGLS, test = TRUE)

    ##
    ## Please be aware that the test of coefficients will be deprecated.
    ## Use betaTest, which is more flexible;
    ## e.g., betaTest(fit, null.method = 'terms')
    ##
    ## Linear Model fit with lm.rrpp
    ##
    ## Number of observations: 37
    ## Number of dependent variables: 1
    ## Data space dimensions: 1
    ## Sums of Squares and Cross-products: Type I
    ## Number of permutations: 1000
    ##
    ## Statistics (distances) of coefficients with
    ## 95 percent confidence intervals, effect sizes,
    ## and probabilities of exceeding observed values,
    ## based on 1000 random permutations, using RRPP
    ##
    ##         d.obs UCL (95%)       Zd Pr(>d)
    ## SVL 0.8307223 0.4503532 2.894442  0.001

Although analyses on either model fit indicate a significant relationship between tail length and snout-to-vent length (SVL), the GLS coefficients test and ANOVA show how phylogenetic auto-correlation among species augments the OLS-estimated relationship. The following is a multivariate example:

    Y <- as.matrix(cbind(PlethMorph$TailLength,
    PlethMorph$HeadLength,
    PlethMorph$TailLength,
    PlethMorph$Snout.eye,
    PlethMorph$BodyWidth,
    PlethMorph$Forelimb,
    PlethMorph$Hindlimb))
    PlethMorph$Y <- Y

    fitOLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph,
                       print.progress = FALSE,
                      turbo = FALSE, verbose = TRUE)
    fitGLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph,
                       Cov = PlethMorph$PhyCov,
                       print.progress = FALSE,
                      turbo = FALSE, verbose = TRUE)

    anova(fitOLSm)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Ordinary Least Squares
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df      SS     MS     Rsq     F      Z Pr(>F)
    ## SVL        1  8135.3 8135.3 0.76867 116.3 5.0898  0.001 ***
    ## Residuals 35  2448.3   70.0 0.23133
    ## Total     36 10583.7
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = Y ~ SVL, turbo = FALSE, data = PlethMorph, print.progress = FALSE,
    ##     verbose = TRUE)

    anova(fitGLSm)

    ##
    ## Analysis of Variance, using Residual Randomization
    ## Permutation procedure: Randomization of null model residuals
    ## Number of permutations: 1000
    ## Estimation method: Generalized Least-Squares (via OLS projection)
    ## Sums of Squares and Cross-products: Type I
    ## Effect sizes (Z) based on F distributions
    ##
    ##           Df      SS     MS     Rsq      F      Z Pr(>F)
    ## SVL        1  395.35 395.35 0.34715 18.611 3.0419  0.001 ***
    ## Residuals 35  743.49  21.24 0.65285
    ## Total     36 1138.84
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Call: lm.rrpp(f1 = Y ~ SVL, turbo = FALSE, data = PlethMorph, Cov = PlethMorph$PhyCov,
    ##     print.progress = FALSE, verbose = TRUE)

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

    plot(predict(fitOLSm, sizeDF), PC= TRUE) # Correlated error

    plot(predict(fitGLSm, sizeDF), PC= TRUE) # Independent error

Analytical Summary

On the surface, these three examples and their analyses should seem intuitive to any user of R who has used the `lm` function plus its associated S3 generics (`coef, predict, resid, fitted, summary`, and `anova`), all of which can be used on `lm.rrpp` model fits. The functions, `pairwise` (not an S3 generic) and `anova`, also allow pairwise comparisons of least-squares means or slopes and multi-model inferences, respectively. Advanced users will recognize, however, much more extensive usable results for adaptive programming. The output from a `lm.rrpp` fit is arranged hierarchically, i.e.:

    attributes(fit)

    ## $names
    ## [1] "call"     "LM"       "ANOVA"    "PermInfo" "turbo"    "Models"   "verbose"
    ## [8] "subTest"
    ##
    ## $class
    ## [1] "lm.rrpp"

Within the `$LM` partition, all attributes of the lm function are found, in addition to coefficients for every random permutation. Within the `$ANOVA` partition, the SS type, plus `SS`, `MS`, `R`2, `F`, and Cohen’s `f`2 for all permutations, as well as effect sizes estimated for each of these are provided. Within the `$PermInfo` partition, the number of permutations, type (RRPP or randomization of “full” data values, FRPP), and sampling frame in every permutation (schedule) are provided. Thus, lm.rrpp is the workhorse that makes all downstream analysis efficient.

References

Adams, D. C., & Collyer, M. L. (2018). Multivariate phylogenetic comparative methods: Evaluations, comparisons, and recommendations. Systematic Biology, 67, 14–31.

Adams, D. C., Collyer, M. L., Kaliontzopoulou, A., & Sherratt, E. (2017). Geomorph: Software for geometric morphometric analyses. R package version 3.0.6.

Adams, D. C., Rohlf, F. J., & Slice, D. E. 2013. A field comes of age: Geometric morphometrics in the 21st century. Hystrix, 24, 7–14.

Bookstein, F. L. (1991). Morphometric tools for landmark data: geometry and biology. Cambridge: Cambridge University Press.

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

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