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.
(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 pecosnesis), 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)
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
## logHeadSize ~ sex + locality/year 8 259 2.718175 10.94887
## Rsq F Z (from F) Pr(>F)
## logHeadSize ~ sex + locality/year 0.1988854 8.037445 12.6872 0.001##
## 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 1.7127 0.001 **
## locality 1 0.4712 0.47121 0.03448 11.1466 1.6665 0.001 **
## locality:year 6 1.6379 0.27298 0.11984 6.4574 3.2300 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 = "I",
## data = PupfishHeads, print.progress = FALSE)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:
##
## 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 1.7127 0.001 **
## locality 1 0.4712 0.47121 0.03448 1.7262 0.6967 0.248
## locality:year 6 1.6379 0.27298 0.11984 6.4574 3.2300 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 = "I",
## data = PupfishHeads, print.progress = FALSE)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)
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 13.058 0.044568
## logHeadSize ~ sex + locality 265 1 12.587 0.47121 0.47121 0.079045
## logHeadSize ~ sex + locality/year 259 7 10.949 2.10907 0.30130 0.198885
## Total 267 13.667
## F Z Pr(>F)
## logHeadSize ~ sex (Null)
## logHeadSize ~ sex + locality 9.9207 1.6155 0.001 **
## logHeadSize ~ sex + locality/year 7.1273 3.5277 0.001 **
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1One might wish to also look at individual model coefficients, and ascertain which have the largest effect:
##
## 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)
## (Intercept) 2.96477681 3.04519769 -3.3910750 0.999
## sexMale 0.09579938 0.05816582 4.2020830 0.002
## localitySH 0.09164189 0.05771827 3.7900752 0.001
## localitySH:year1998 0.06402169 0.13525163 0.2373692 0.344
## localityLake:year1999 0.35681331 0.13628574 6.9373198 0.001
## localitySH:year1999 0.01328613 0.13859748 -0.9637889 0.846
## localitySH:year2000 0.07165589 0.10185529 0.9919253 0.167
## localityLake:year2001 0.23035259 0.13141374 4.1803649 0.002
## localityLake:year2002 0.31050534 0.13751014 5.9007583 0.001This 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: Not all variables in model accounted for in newdata.
## Missing variables will be averaged from observed data for prediction.
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).
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)
fit3 <- lm.rrpp(logHeadSize ~ sex + locality/year,
SS.type = "III", data = PupfishHeads, 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)
## sex 1 0.6091 0.60911 0.04457 14.4087 1.7127 0.001 **
## locality 1 0.4712 0.47121 0.03448 11.1466 1.6665 0.001 **
## locality:year 6 1.6379 0.27298 0.11984 6.4574 3.2300 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 = "I",
## data = PupfishHeads, 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
## 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 1.7486 0.001 **
## locality 1 0.4712 0.47121 0.03448 11.1466 1.5918 0.001 **
## locality:year 6 1.6379 0.27298 0.11984 6.4574 2.9795 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)##
## 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 1.74864 0.001 **
## locality 1 0.0573 0.05733 0.00419 1.3562 0.73543 0.228
## locality:year 6 1.6379 0.27298 0.11984 6.4574 2.97954 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)
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)
## fit 4 49 0.03224913 0.02408373 0.5724746 16.40327 24.93101
## Pr(>F)
## fit 0.001##
## 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 5.1789 0.001 **
## Sex 1 0.007615 0.0076151 0.13518 15.4933 5.5329 0.001 **
## Pop 1 0.007661 0.0076606 0.13599 15.5860 6.1375 0.001 **
## Sex:Pop 1 0.002954 0.0029542 0.05244 6.0105 5.2777 0.001 **
## 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, SS.type = "I", data = Pupfish,
## print.progress = 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
##
## 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)
## (Intercept) 1.10733195 1.02055707 14.274555 0.001
## logSize 0.11101234 0.04757252 7.978520 0.001
## SexM 0.02755503 0.01314346 7.517624 0.001
## PopSinkhole 0.02895656 0.01257359 8.923784 0.001
## SexM:PopSinkhole 0.03002431 0.01755391 5.716307 0.001ANOVA 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)
identical(fit$LM$coefficients, fit.par$coefficients)## [1] TRUE## Error in summary.manova(manova(fit.par)): residuals have rank 49 < 112Although 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: Not all variables in model accounted for in newdata.
## Missing variables will be averaged from observed data for prediction.
plot(shapePreds, PC = TRUE, ellipse = TRUE,
pch = 19, col = 1:NROW(shapeDF)) # 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.03289829 2.4842340 0.011
## F.Marsh:F.Sinkhole 0.03639758 0.03554347 1.9212320 0.026
## F.Marsh:M.Sinkhole 0.03809715 0.04432661 -0.6238091 0.723
## M.Marsh:F.Sinkhole 0.03681731 0.04763945 -0.2487848 0.579
## M.Marsh:M.Sinkhole 0.02694896 0.03670302 -1.0029384 0.842
## F.Sinkhole:M.Sinkhole 0.01939739 0.03262909 -1.5591950 0.952Much 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)
## fit2 7 46 0.03425362 0.02207925 0.6080574 10.19488 22.05171
## Pr(>F)
## fit2 0.001##
## 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 0.024084
## coords ~ logSize * Sex * Pop 46 3 0.022079 0.0020045 0.00066816
## Total 53 0.056333
## Rsq F Z Pr(>F)
## coords ~ logSize + Sex * Pop (Null) 0.57247
## coords ~ logSize * Sex * Pop 0.60806 1.3921 3.6279 0.001 **
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PW2 <- 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##
## 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)
##
## Slope vector lengths
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## 0.09846857 0.08899635 0.08048951 0.13005765
##
## Pairwise absolute difference (d) between vector lengths, plus statistics
## d UCL (95%) Z Pr > d
## F.Marsh:M.Marsh 0.009472222 0.08066795 -0.9266208 0.805
## F.Marsh:F.Sinkhole 0.017979061 0.08701897 -0.6684921 0.691
## F.Marsh:M.Sinkhole 0.031589079 0.08573929 -0.1404949 0.485
## M.Marsh:F.Sinkhole 0.008506839 0.04838718 -0.7095093 0.706
## M.Marsh:M.Sinkhole 0.041061301 0.04776048 1.3887438 0.093
## F.Sinkhole:M.Sinkhole 0.049568140 0.04244181 2.4759957 0.023##
## 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)
##
## Slope vector lengths
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## 0.09846857 0.08899635 0.08048951 0.13005765
##
## Pairwise absolute differences (d) between slope lengths
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 0.000000000 0.009472222 0.017979061 0.03158908
## M.Marsh 0.009472222 0.000000000 0.008506839 0.04106130
## F.Sinkhole 0.017979061 0.008506839 0.000000000 0.04956814
## M.Sinkhole 0.031589079 0.041061301 0.049568140 0.00000000
##
## Pairwise 95% upper confidence limits between slope lengths
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 0.00000000 0.08066795 0.08701897 0.08573929
## M.Marsh 0.08066795 0.00000000 0.04838718 0.04776048
## F.Sinkhole 0.08701897 0.04838718 0.00000000 0.04244181
## M.Sinkhole 0.08573929 0.04776048 0.04244181 0.00000000
##
## Pairwise effect sizes (Z) between slope lengths
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 0.0000000 -0.9266208 -0.6684921 -0.1404949
## M.Marsh -0.9266208 0.0000000 -0.7095093 1.3887438
## F.Sinkhole -0.6684921 -0.7095093 0.0000000 2.4759957
## M.Sinkhole -0.1404949 1.3887438 2.4759957 0.0000000
##
## Pairwise P-values between slope lengths
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 1.000 0.805 0.691 0.485
## M.Marsh 0.805 1.000 0.706 0.093
## F.Sinkhole 0.691 0.706 1.000 0.023
## M.Sinkhole 0.485 0.093 0.023 1.000summary(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 85.99703 -0.2996071 0.582
## F.Marsh:F.Sinkhole 0.6318267 50.81498 82.18877 -0.2229753 0.538
## F.Marsh:M.Sinkhole 0.4461018 63.50614 82.21052 0.5581365 0.260
## M.Marsh:F.Sinkhole 0.7175129 44.15048 65.21340 0.1341814 0.395
## M.Marsh:M.Sinkhole 0.5629975 55.73665 66.72786 1.0332053 0.146
## F.Sinkhole:M.Sinkhole 0.4627734 62.43378 59.87174 2.1810996 0.036Because 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 (dimensions of data after PCoA of distance matrix)
## Data space dimensions: 53 (dimensions of data after PCoA of distance matrix)
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Full Model Analysis of Variance
##
## Df Residual Df SS Residual SS Rsq
## D ~ logSize + Sex * Pop 4 49 0.03224913 0.02408373 0.5724746
## F Z (from F) Pr(>F)
## D ~ logSize + Sex * Pop 16.40327 24.93101 0.001##
## 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
## coords ~ logSize + Sex * Pop 4 49 0.03224913 0.02408373
## Rsq F Z (from F) Pr(>F)
## coords ~ logSize + Sex * Pop 0.5724746 16.40327 24.93101 0.001##
## 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 5.1789 0.001 **
## Sex 1 0.007615 0.0076151 0.13518 15.4933 5.5329 0.001 **
## Pop 1 0.007661 0.0076606 0.13599 15.5860 6.1375 0.001 **
## Sex:Pop 1 0.002954 0.0029542 0.05244 6.0105 5.2777 0.001 **
## 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)##
## 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 5.1789 0.001 **
## Sex 1 0.007615 0.0076151 0.13518 15.4933 5.5329 0.001 **
## Pop 1 0.007661 0.0076606 0.13599 15.5860 6.1375 0.001 **
## Sex:Pop 1 0.002954 0.0029542 0.05244 6.0105 5.2777 0.001 **
## 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, SS.type = "I", data = Pupfish,
## print.progress = FALSE)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)
fitGLS <- lm.rrpp(TailLength ~ SVL,
data = PlethMorph,
Cov = PlethMorph$PhyCov,
print.progress = FALSE)
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 2.6789 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, data = PlethMorph, print.progress = FALSE)##
## 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.35 189.346 0.3426 18.24 1.8545 0.002 **
## Residuals 35 363.33 10.381 0.6574
## Total 36 552.68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = TailLength ~ SVL, data = PlethMorph, Cov = PlethMorph$PhyCov,
## print.progress = FALSE)##
## 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)
## (Intercept) 1.2992281 78.5454347 -5.143547 1.000
## SVL 0.9681511 0.3648632 7.290075 0.001##
## 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)
## (Intercept) 1.2992281 78.5454347 -5.143547 1.000
## SVL 0.9681511 0.3648632 7.290075 0.001Although 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 <- rrpp.data.frame(PlethMorph, Y=Y)
fitOLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph,
print.progress = FALSE)
fitGLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph,
Cov = PlethMorph$PhyCov,
print.progress = FALSE)
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 3.5609 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, data = PlethMorph, print.progress = FALSE)##
## 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.36 395.36 0.34715 18.611 2.3489 0.002 **
## Residuals 35 743.49 21.24 0.65285
## Total 36 1138.85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ SVL, data = PlethMorph, Cov = PlethMorph$PhyCov,
## print.progress = FALSE)sizeDF <- data.frame(SVL = sort(PlethMorph$SVL))
plot(predict(fitOLSm, sizeDF), PC= TRUE) # Correlated 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 useable results for adaptive programming. The output from a lm.rrpp fit is arranged hierarchically, i.e.:
## $names
## [1] "call" "LM" "ANOVA" "PermInfo"
##
## $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, R2, F, and Cohen’s f2 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: Evalutions, 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). Western Kentucky University.