20.12 – Phylogenetically independent contrasts

Introduction

Assumption of independence among the subjects in a study, key assumption. Cross species comparisons as tests of correlated evolution of traits has a long history in biology research (Cornwallis and Griffin 2024). Methods have transitioned from observations to rigorous, phylogenetic tree-based statistical approaches (phylogenetic comparative methods), and continue to provide important analytical tools, eg., host adaptation and cross-species transmission of Covid-19 virus (Attwood et al 2022).

Comparisons among species for the association between body size and metabolic rate, for example Kleiber’s Law, are classic cases (Fig 1, after Table 2 in Kleiber 1947). We’ll use Kleiber’s data basal metabolic rate (BMR) and body mass for a collection of species to explore two phylogenetic comparative methods, PIC and PGLS.

Figure 1. Scatter plot of BMR on Body Mass, log-transformed.

Typical approaches statistically, use of ANOVA or linear regression approaches. A basic assumption of ANOVA is that sampling units are independent 13.1 – ANOVA assumptions.

Prior to the 1980s, it was rarely appreciated in comparative analysis that species are not independent sample units (Harvey and Pagel 1991); evolution produced nested hierarchical relationships among the species. We recognize this with phylogenies (Felsenstein 1985, Harvey and Pagel 1991, Martins 1996, Garland et al 2005).  Mice and rats share more recent common ancestor, and cattle and pigs share more recent common ancestor, than do mice and cattle, for example. Felsenstein (1985, 1988), largely credited for making the argument that Type I error likely if phylogeny ignored, and, importantly, provided an algorithm, Phylogenetic Independent Contrasts, PIC, which provided a simple way to correct for phylogenetic nonindependence. Felsentein’s landmark 1985 paper has been cited more than ten thousand times (Feb 2024).

Logic of PIC

Treating comparative data, eg, species, as a collection of independent samples implies that the evolutionary history was a spontaneous burst, or star-like phylogeny — an unresolved phylogenetic tree with branches radiating from a single ancestral node (Fig 2).

Figure 2. Star tree phylogeny (same image shown Figure 5, 20.11 – Plot a Newick tree).

But what nature provides is nonindependence (Fig 3, for more about star phylogeny in PIC see discussion in Garland et al 2005), Darwin’s “groups under groups,” or a nested hiearchy. Species are not independent units and phylogenetic history should be accounted for during statistical analysis.

Figure 3. A cladogram for same species, showing the hierarchical, nested relationships among taxa, what nature actually provides (same image shown Figure 2, 20.11 – Plot a Newick tree).

Let’s revisit Figure 1. What do the points look like if we line up the known phylogeny for Kleiber’s six species against the plotted data? That gives us Figure 4.

Figure 4. Kleiber’s Law and the confounding effect of shared ancestry.

Dashed connecting lines in Figure 4 trace each species from its tree position to its phenotypic value, demonstrating how phylogenetic clustering violates the independent data assumptions of standard OLS regression.

Begin R code for Figure 4, with remarks — code generated with Gemini AI assistance.

# 1. Reset the data frame 
kleiber_data <- data.frame(
Species = c("Mouse", "Rat", "Rabbit", "Cat", "Dog", "Cow"),
Mass = c(0.021, 0.282, 2.45, 3.1, 13.4, 439.0),
BMR = c(3.38, 32.4, 134.0, 151.0, 506.0, 6135.0),
stringsAsFactors = FALSE
)
kleiber_data$log_Mass <- log10(kleiber_data$Mass)
kleiber_data$log_BMR <- log10(kleiber_data$BMR)

# 2. Build the exact matching standalone 6-species tree
library(ape)
newick_tree_6 <- "((((Mouse:10,Rat:10):10,Rabbit:20):20,(Cat:20,Dog:20):20):15,Cow:55);"
my_tree <- read.tree(text = newick_tree_6)

# 3. Setup a side-by-side 2-panel layout
layout(matrix(c(1, 2), nrow = 1), widths = c(1, 1.3))
par(mar = c(5, 4, 4, 1))

# PANEL 1: Plot the tree hierarchy 
plot(my_tree, main = "1. Nested Tree", 
edge.width = 3, edge.color = "darkgrey", cex = 1.1)

# FIXED COORDINATE MAP: 'ape' plots leaves bottom-to-top exactly matching the tree order
# We map the tree tip labels to simple integers 1 through 6
tree_heights <- 1:6
names(tree_heights) <- my_tree$tip.label

# PANEL 2: Plot the scatterplot
plot(kleiber_data$log_BMR ~ kleiber_data$log_Mass, 
pch = 19, cex = 2, col = "blue", 
xlab = "Body mass (log10 kg)", ylab = "BMR (log10 kcal/day)",
xlim = c(min(kleiber_data$log_Mass) - 0.2, max(kleiber_data$log_Mass) + 0.6),
ylim = c(min(kleiber_data$log_BMR) - 0.3, max(kleiber_data$log_BMR) + 0.3),
main = "2. Allometric Trait Space")

# Re-apply your custom 310-degree slanted text labels
text(kleiber_data$log_Mass + 0.05, kleiber_data$log_BMR - 0.03, 
labels = kleiber_data$Species, srt = 310, adj = 0, cex = 0.9, col = "black")

# 4. CROSS-PANEL LINE CORRECTION: Draw lines across the margins safely
par(xpd = NA) # Tells R to ignore the borders of the active scatterplot window

# Calculate the left-hand margin boundaries based on the current plot limits
left_x_boundary <- min(kleiber_data$log_Mass) - 2.8

for(i in 1:nrow(kleiber_data)) {
sp <- kleiber_data$Species[i]

# A. Use a standardized geometric multiplier to align with Panel 1's vertical axis
# This formula maps the tree heights (1-6) into the scatterplot's Y coordinate scale
y_scale_factor <- (tree_heights[sp] - 1) / 5
y_span_range <- max(kleiber_data$log_BMR) - min(kleiber_data$log_BMR)
target_tree_y <- min(kleiber_data$log_BMR) + (y_scale_factor * y_span_range)

# B. Pull the actual coordinates of your blue dot
right_x <- kleiber_data$log_Mass[i]
right_y <- kleiber_data$log_BMR[i]

# C. Draw a clean, connecting line directly out to the dot
segments(x0 = left_x_boundary,
y0 = target_tree_y,
x1 = right_x, 
y1 = right_y, 
col = adjustcolor("darkblue", alpha.f = 0.3), 
lty = 3, lwd = 2)
}

# Reset layout parameters back to normal
layout(1)

End R code.

library(phytools)
phylomorphospace(my_tree, trait_matrix,
label = "off", # Suppresses default horizontal text labels safely
xlab = "Body mass (log10 kg)",
ylab = "BMR (log10 kcal/day)",
xlim = c(min(kleiber_data$log_Mass) - 0.2, max(kleiber_data$log_Mass) + 0.6),
ylim = c(min(kleiber_data$log_BMR) - 0.3, max(kleiber_dat$log_BMR) + 0.3),
control = list(col.edge = "darkgrey", lwd = 3))

Phylogenetic Generalized Least Squares (PGLS).

However, like most innovations, PIC should not be blindly applied in all comparative analysis (eg, unreplicated evolutionary events, Uyeda et al 2018). Fundamentally, PIC is mathematically limited: it primarily handled simple two-variable correlations and required fully resolved, bifurcating evolutionary trees with known branch lengths. At about the same time, Grafen (1989) proposed that phylogenetic non-independence could be treated as a form of spatial autocorrelation — that is, near things are more related than distant things. Instead of computing data contrasts at individual tree nodes, Grafen introduced the phylogenetic regression, which used Generalized Least Squares (GLS) to alter how the regression algorithm treats the residual errors of species data. Improvements to phylogenetic regression are now collectively refered to as the Phylogenetic Generalized Least Squares (PGLS) (Rohlf 2001, Housworth et al 2004, Mazel et al 2016, Adams et al 2019).

Logic of PGLS.

In standard Ordinary Least Squares (OLS) regression, the residual error terms (ε) are assumed to be independent and identically distributed. In matrices, this error structure is represented by an identity matrix (I), where the off-diagonals are zero, implying that knowing about the trait value of a mouse reveals no information about the value for a rat or other rodent species.

kleiber_full <- data.frame(
    Species = c("Mouse", "Rat", "Rabbit", "Macaque", "Cat", "Dog", 
              "Goat", "Sheep", "Chimpanzee", "Woman", "Cow"),
# Exact body mass (W) in kilograms
    Mass = c(0.021, 0.282, 2.45, 4.20, 3.10, 13.40, 
           36.00, 46.40, 38.00, 56.10, 439.00), 
# Exact basal metabolic rate in kcal per 24 hours
    BMR = c(3.38, 32.7, 134.0, 207.0, 151.0, 506.0, 
          826.0, 1145.0, 1090.0, 1398.0, 6135.0), 
stringsAsFactors = FALSE
)

"(((((Mouse:10,Rat:10):10,Rabbit:20):20,((Macaque:25,(Chimpanzee:6,Woman:6):19):10):10):10,(Cat:20,Dog:20):30):10,(Goat:15,Sheep:15):40,Cow:55);"
my_tree <- read.tree(text = expanded_newick)