What is the “Felsenstein problem”?

Not entirely sure if the name “Felsenstein problem” is appropriate but Felsenstein was among the first to clearly point out that interspecific data violate a basic assumption of traditional statistics: the independence of observations. Ignoring the non-independence of observation can lead to false interpretation of patterns, as he convincingly demonstrated in 1985. This paper proved to be very influential for the field of comparative biology.

In this tutorial we will explore how independent contrasts may mitigate the problem of phylogenetic covariance by means of simulations.

#Load libraries.
library(ape)
library(geiger)
library(phytools)

#Simulate trees and data

#The rcoal() function in 'ape' simulates coalescent (clustered), ultrametric trees.
phy <- rcoal(n=100) # n specifies number of tips on the tree 
  
#Liam's fastBM() creates continuous trait under the BM model (plus BM with trend, and bounded BM)
x <- fastBM(phy, sig2=0.2) # sig2 is the the instantaneous rate variance of the BM process
y <- fastBM(phy, sig2=0.2) # y is here simulated completely independent of x
y1 <- 0.7*x + y # y1 is forced to be strongly correlated with x
  
#Let's examine what we actually simulated.
fit <- lm(y~x) #OLS regression
summary(fit)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34731 -0.13497  0.00177  0.14482  0.35196 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33494    0.04589   7.299 7.68e-11 ***
## x           -0.52526    0.06419  -8.182 1.03e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1638 on 98 degrees of freedom
## Multiple R-squared:  0.4059, Adjusted R-squared:  0.3998 
## F-statistic: 66.95 on 1 and 98 DF,  p-value: 1.034e-12
cor.test(x, y) #Correlation test
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = -8.1824, df = 98, p-value = 1.034e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7408066 -0.5037033
## sample estimates:
##      cor 
## -0.63709
#And here the same for x and y1.
fit2 <- lm(y1~x)
summary(fit2)
## 
## Call:
## lm(formula = y1 ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34731 -0.13497  0.00177  0.14482  0.35196 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33494    0.04589   7.299 7.68e-11 ***
## x            0.17474    0.06419   2.722  0.00768 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1638 on 98 degrees of freedom
## Multiple R-squared:  0.07029,    Adjusted R-squared:  0.06081 
## F-statistic: 7.409 on 1 and 98 DF,  p-value: 0.00768
cor.test(x, y1)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y1
## t = 2.722, df = 98, p-value = 0.00768
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07248391 0.43869984
## sample estimates:
##       cor 
## 0.2651274

Are the simulated traits correlated?

#Let's also visualize the data and tree => x and y
par(mfrow=c(1,2)) 
plot(x, y, asp=1)
abline(fit, col="red") #Adding the regression line
plot(phy, show.tip.label=FALSE) #Plotting the tree next to the scatter plot

#And the same for x and y1
par(mfrow=c(1,2)) 
plot(x, y1, asp=1)
abline(fit2, col="red")
plot(phy, show.tip.label=FALSE)

How may tree shape influence patterns in the scatterplot and correlation?

How can we calculate independent contrasts?

There are multiple options for contrast calculation in R; we will use the implementation in the ape package (pic()). For details, please refer to Felsenstein (1985).

x.contrasts <- pic(x, phy, scaled=TRUE, var.contrast=TRUE) #If you want contrasts standardized by branch length set scaled=TRUE
y.contrasts <- pic(y, phy, scaled=TRUE, var.contrast=TRUE)
y1.contrasts <- pic(y1, phy, scaled=TRUE, var.contrast=TRUE)

#We obtain a 2 cloumn matrix, with contrasts and their variances in the 1st and 2nd column, respectively:
head(x.contrasts)
##       contrasts  variance
## 101  0.20579832 4.0738872
## 102  0.20088697 0.5182375
## 103 -0.02606677 0.7292417
## 104  0.38464557 0.4857205
## 105  0.14792834 0.4429421
## 106  0.48631309 0.4252179

Before proceeding we must test if the contrasts are adequately standardized. To do this we are plotting the absolute standardized contrasts against their standard deviations.

What is it that we want to see?

#Making absolute contrasts:
abs.x <- abs(x.contrasts[,1])
abs.y <- abs(y.contrasts[,1])
abs.y1 <- abs(y1.contrasts[,1])

#Getting SD from contrast variance:
sd.x <-  sqrt(x.contrasts[,2])
sd.y <-  sqrt(y.contrasts[,2])
sd.y1 <-  sqrt(y1.contrasts[,2])

#Checking if there are any patterns for X contrasts.
fit.diagn <- lm(abs.x~sd.x)
summary(fit.diagn)
## 
## Call:
## lm(formula = abs.x ~ sd.x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36484 -0.21669 -0.06282  0.14759  0.91020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.38278    0.03818  10.026   <2e-16 ***
## sd.x        -0.12549    0.11193  -1.121    0.265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.273 on 97 degrees of freedom
## Multiple R-squared:  0.01279,    Adjusted R-squared:  0.002616 
## F-statistic: 1.257 on 1 and 97 DF,  p-value: 0.265
cor.test(sd.x, abs.x)
## 
##  Pearson's product-moment correlation
## 
## data:  sd.x and abs.x
## t = -1.1212, df = 97, p-value = 0.265
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.30373704  0.08622996
## sample estimates:
##        cor 
## -0.1131073
par(mfrow=c(1,1))
plot(sd.x, abs.x, asp=1)
abline(fit.diagn, col="red")

#Checking if there are any patterns for Y contrasts.
fit.diagn <- lm(abs.y~sd.y)
summary(fit.diagn)
## 
## Call:
## lm(formula = abs.y ~ sd.y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32866 -0.20750 -0.06535  0.13000  0.88783 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.35666    0.03833   9.305 4.25e-15 ***
## sd.y        -0.17209    0.11237  -1.531    0.129    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.274 on 97 degrees of freedom
## Multiple R-squared:  0.02361,    Adjusted R-squared:  0.01354 
## F-statistic: 2.345 on 1 and 97 DF,  p-value: 0.1289
cor.test(sd.y, abs.y)
## 
##  Pearson's product-moment correlation
## 
## data:  sd.y and abs.y
## t = -1.5315, df = 97, p-value = 0.1289
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.34072836  0.04512898
## sample estimates:
##        cor 
## -0.1536517
plot(sd.y, abs.y, asp=1)
abline(fit.diagn, col="red")

#Checking if there are any patterns for Y1 contrasts.
fit.diagn <- lm(abs.y1~sd.y1)
summary(fit.diagn)
## 
## Call:
## lm(formula = abs.y1 ~ sd.y1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44072 -0.23438 -0.08965  0.16587  0.88257 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.49500    0.04375  11.315   <2e-16 ***
## sd.y1       -0.31901    0.12825  -2.487   0.0146 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3128 on 97 degrees of freedom
## Multiple R-squared:  0.05996,    Adjusted R-squared:  0.05027 
## F-statistic: 6.187 on 1 and 97 DF,  p-value: 0.01458
cor.test(sd.y1, abs.y1)
## 
##  Pearson's product-moment correlation
## 
## data:  sd.y1 and abs.y1
## t = -2.4873, df = 97, p-value = 0.01458
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.42187943 -0.04985884
## sample estimates:
##        cor 
## -0.2448606
plot(sd.y1, abs.y1, asp=1)
abline(fit.diagn, col="red")

After assuring that the contrasts are adequately standardized, we can now assess correlated evolution of X and Y. Remember that X and Y should be uncorrelated!

contrast.fit <- lm(y.contrasts[,1] ~ x.contrasts[,1] -1) # OLS regression forced through the origin.

#The slope of the fitted line and the sumamry stats are extracted in the next 2 lines:
contrast.fit
## 
## Call:
## lm(formula = y.contrasts[, 1] ~ x.contrasts[, 1] - 1)
## 
## Coefficients:
## x.contrasts[, 1]  
##          0.01939
summary(contrast.fit)
## 
## Call:
## lm(formula = y.contrasts[, 1] ~ x.contrasts[, 1] - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23380 -0.24267 -0.02833  0.24431  1.08803 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## x.contrasts[, 1]  0.01939    0.09484   0.204    0.838
## 
## Residual standard error: 0.4205 on 98 degrees of freedom
## Multiple R-squared:  0.0004263,  Adjusted R-squared:  -0.009773 
## F-statistic: 0.04179 on 1 and 98 DF,  p-value: 0.8384
cor.test(x.contrasts[,1], y.contrasts[,1])
## 
##  Pearson's product-moment correlation
## 
## data:  x.contrasts[, 1] and y.contrasts[, 1]
## t = 0.28177, df = 97, p-value = 0.7787
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1697727  0.2247406
## sample estimates:
##        cor 
## 0.02859755
#Let's compare results to a regression performed on the "raw" data.
fit <- lm(y~x)
summary(fit)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34731 -0.13497  0.00177  0.14482  0.35196 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33494    0.04589   7.299 7.68e-11 ***
## x           -0.52526    0.06419  -8.182 1.03e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1638 on 98 degrees of freedom
## Multiple R-squared:  0.4059, Adjusted R-squared:  0.3998 
## F-statistic: 66.95 on 1 and 98 DF,  p-value: 1.034e-12
cor.test(x, y)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = -8.1824, df = 98, p-value = 1.034e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7408066 -0.5037033
## sample estimates:
##      cor 
## -0.63709
#Visualization is always helpful:
par(mfrow=c(1,2)) 
plot(x, y, asp=1, main="raw data") # Plotting raw data, aspect ratio of plot is 1.
abline(fit, col="red") # # Adding the regression line for raw data
plot(x.contrasts[,1], y.contrasts[,1], asp=1, main="contrasts") # Plotting contrasts
abline(contrast.fit, col="red") # Adding contrasts regression

#Now let's assess correlated evolution of X and Y1.
#Remember that X and Y1 should indeed be correlated!
contrast.fit2 <- lm(y1.contrasts[,1] ~ x.contrasts[,1] -1) # OLS regression forced through the origin

#The slope of the fitted line and the sumamry stats are extracted in the next 2 lines:
contrast.fit2
## 
## Call:
## lm(formula = y1.contrasts[, 1] ~ x.contrasts[, 1] - 1)
## 
## Coefficients:
## x.contrasts[, 1]  
##           0.7194
summary(contrast.fit2)
## 
## Call:
## lm(formula = y1.contrasts[, 1] ~ x.contrasts[, 1] - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23380 -0.24267 -0.02833  0.24431  1.08803 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## x.contrasts[, 1]  0.71939    0.09484   7.585 1.93e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4205 on 98 degrees of freedom
## Multiple R-squared:  0.3699, Adjusted R-squared:  0.3635 
## F-statistic: 57.53 on 1 and 98 DF,  p-value: 1.926e-11
cor.test(x.contrasts[,1], y1.contrasts[,1])
## 
##  Pearson's product-moment correlation
## 
## data:  x.contrasts[, 1] and y1.contrasts[, 1]
## t = 7.5936, df = 97, p-value = 1.94e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4698191 0.7210910
## sample estimates:
##       cor 
## 0.6105992
#Let's compare results to a regression performed on the "raw" data.
fit2 <- lm(y1~x)
summary(fit2)
## 
## Call:
## lm(formula = y1 ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34731 -0.13497  0.00177  0.14482  0.35196 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33494    0.04589   7.299 7.68e-11 ***
## x            0.17474    0.06419   2.722  0.00768 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1638 on 98 degrees of freedom
## Multiple R-squared:  0.07029,    Adjusted R-squared:  0.06081 
## F-statistic: 7.409 on 1 and 98 DF,  p-value: 0.00768
cor.test(x, y1)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y1
## t = 2.722, df = 98, p-value = 0.00768
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07248391 0.43869984
## sample estimates:
##       cor 
## 0.2651274
#Visualization is always helpful: ;)
par(mfrow=c(1,2)) 
plot(x, y1, asp=1, main="raw data") # Plotting raw data, aspect ratio of plot is 1.
abline(fit2, col="red") # # Adding the regression line for raw data
plot(x.contrasts[,1], y1.contrasts[,1], asp=1, main="contrasts") # Plotting contrasts
abline(contrast.fit2, col="red") # Adding contrasts regression