Trait distribution at equilibrium

An interesting property of the POUMM is that, in the limit t → ∞, it defines a stationary normal distribution for the heritable component (g) at the tips with mean θ and a variance-covariance matrix:

where d_ij denotes the phylogenetic distance between the tips i and j. The above property proves useful when there is a prior knowledge that the observed population is at equilibrium, because one can use the trait variance in the population, σ_z² = σ_∞² + σ_e² as model parameter. The corresponding parametrization is:

With this parametrization, one can specify an informed prior for σ_z² based on empirical estimates on similar data.

Another important aspect of the above property is that it helps to better understand the selection stregnth parameter α. As it turns out, α can have two different biological interpretations. On the one hand, α defines the rate of convergence of the population mean towards the long-term optimum θ. This rate is bigger for bigger values of α and for bigger deviations from θ. Thus, α is considered as selection strength or rate of adaptation under stabilizing selection. On the other hand, it is possible to assume that the majority of the tips and their mrca’s are far enough from the root, so that Σ_ij can be viewed as an exponentially decreasing function of the phylogenetic distance d_ij. Seen from that angle, the parameter α can be interpreted as the rate of phenotypic decorrelation between tips, due to genetic drift. When interpreting the results of a model fit, it is important to be aware of this dual interpretation of α. In many cases (e.g. in ultrametric macro-evolutionary tree), the only source of information for inferring α are the observed differences between the tips in the tree. Thus, in the absence of additional evidence, it can be erronous to assume that the inferred value of α informs stabilizing selection and an adaptation rate towards θ.

A likelihood ratio test between the ML POUMM and PMM fits can be used to test if the inferred parameter α is significantly above 0. As pointed out in the previous paragraph, a significantly positive α does not necessarily imply stabilizing selection towards θ. Further, it is important to note that the value of α can only be interpreted with respect to the time scale of the phylogeny. It can be more intuitive to consider the phylogenetic half-life, $t_{1/2}=\frac{\ln(2)}{\alpha}$, which equals the time it takes for a species entering a new niche to evolve halfway toward its new expected optimum (Hansen 1997).

Phylogenetic heritability

The term phylogenetic heritability, introduced with the phylogenetic mixed model (PMM) (Housworth, Martins, and Lynch 2004), measures the proportion of phenotypic variance in a population attributable to heritable factors, such as genes, as opposed to non-heritable factors, such as environment and measurement error. Although this concept has been applied mostly in the context of the original PMM, i.e. under the assumption of Brownian motion, the same concept applies to any evolutionary model allowing for the estimation of measurement error (ME) (Hansen and Bartoszek 2012). The phylogenetic heritability is defined as the expected proportion of phenotypic variance attributable to g at the tips of the tree, σ²(g)/[σ²(g) + σ_e²] (Housworth, Martins, and Lynch 2004). This definition is a phylogenetic variant of the definition of broad-sense heritability, H², from quantitative genetics (Lynch and Walsh 1998). However, in the case of a trait evolving along a phylogeny, the expected genotypic variance, σ²(g), and, therefore, the phylogenetic heritability, are functions of time. Depending on the applicaiton, the following three types of phylogenetic heritability might all be of interest:

• Expectation at the mean root-tip distance (t̄): $H_{\bar{t}}^2:=\left[\sigma^2\, \frac{\left(1-e^{-2\alpha \bar{t}}\right)}{2\alpha}\right]/\left[\sigma^2\, \frac{\left(1-e^{-2\alpha \bar{t}}\right)}{2\alpha}+\sigma_e^2\right]$. This definition gives rize to three parametrizations where H_t̄² is a free model parameter, while one of the parameters α, σ or σ_e is a function of H_t̄² and the two other parameters. These parametrizations can be expressed as: where $y=\frac{\bar{t}\sigma^2}{\sigma_e^2 H_{\bar{t}}^2} - \frac{\bar{t}\sigma^2}{\sigma_e^2}$ and W is the Lambert-W function.

• Expectation at equilibrium of the OU-process (t → ∞): $H_{\infty}^2:=\lim_{\bar{t}\to\infty}H_{\bar{t}}^2=\frac{\sigma^2}{\sigma^2+2\alpha\sigma_e^2}$. By taking the limit t̄ → ∞ in equations , and , we obtain the corresponding parametrizations:

• Empirical (time-independent) version of the heritability based on the sample phenotypic variance s²(z) : H_e² := 1 − σ_e²/s²(z). This definition is useful when the tree is non-ultrametric but there is sufficient evidence that the empirical distribution of the trait is stationary along the tree. In this case, s²(z) coincides with the sum of the OU variance at equilibrium and σ_e². The corresponding parametrization is:

To learn about different parametrizations supported by the POUMM package, read the documentation page ?specifyPOUMM.

Is the POUMM an appropriate model for the data?

The first step to answering that question is to visualize the data and check for obvious violations of the POUMM assumptions. The POUMM method expects that the trait-values at the tips are a sample from a multivariate normal distribution. With an ultrametric species tree, where all tips are equally distant from the root, this assumption translates in having all trait-values be realizations of identically distributed normal random variables. In the case of a non-ultrametric tree, it is far more useful to look at a sequence of box-whisker or violin plots of the trait-values, gouped by their root-tip distance.

Once visualizing the data has confirmed its normality, we recommend comaparing the POUMM-fit with a fit from a NULL-model such as the phylogenetic mixed model (PMM) (Housworth, Martins, and Lynch 2004). Since the PMM is nested in the POUMM, i.e. in the limit α → 0, the POUMM model is equivalent to a PMM model with the same initial genotypic value g₀ and unit-time variance σ, it is easy to fit a PMM model to the data by fixing the value of the parameter α to 0:

specPMM <- specifyPMM(z[1:N], tree)
fitPMM <- POUMM(z[1:N], tree, spec = specPMM, doMCMC=FALSE)

Now a likelihood-ratio test between the maximum likelihood fits clearly shows that the POUMM significantly outperforms the PMM to the data:

lmtest::lrtest(fitPMM, fitPOUMM2)

## Likelihood ratio test
## 
## Model 1: fitPMM
## Model 2: fitPOUMM2
##   #Df LogLik Df Chisq Pr(>Chisq)    
## 1   3  -12.9                        
## 2   5   16.3  2  58.5      2e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since lrtest only uses the ML-fit, to save time, we disabled the MCMC fit by specifying doMCMC = FALSE. In real situations, though, it is always recommended to enable the MCMC fit, since it can improve the ML-fit if it finds a region of higher likelihood in the parameter space that has not been discovered by the ML-fit. Moreover, comparing posterior samples for the PMM and POUMM fit can, in principle, be more informative and a likelihood-ratio test.

As an exersise, we can generate data under the PMM model and see if a POUMM fit on that data remains significantly better than a PMM fit:

gBM <- rVNodesGivenTreePOUMM(tree, g0, alpha = 0, theta = 0, sigma = sigma)
zBM <- gBM + e

fitPMM_on_zBM <- POUMM(zBM[1:N], tree, spec = specPMM, doMCMC = FALSE)
fitPOUMM_on_zBM <- POUMM(zBM[1:N], tree, doMCMC = FALSE)

lmtest::lrtest(fitPMM_on_zBM, fitPOUMM_on_zBM)

## Likelihood ratio test
## 
## Model 1: fitPMM_on_zBM
## Model 2: fitPOUMM_on_zBM
##   #Df LogLik Df Chisq Pr(>Chisq)
## 1   3   -114                    
## 2   5   -113  2  1.65       0.44

Assuming that the trait undergoes stabilizing selection, what are the long-term optimum value and the rate of convergence towards it?

To answer this question, consider the estimated values of the POUMM-parameters θ and α. Note that the parameter θ is relevant only if the value of the parameter α is significantly positive. One could accept that the ML-estimate for α is significantly positive if a likelihood ratio test between a ML PMM and POUMM fits gives a p-value below a critical level (see the question above for an example). An inisignificant value of α reveals that the hypothesis of neutral drift (Brownian motion) cannot be rejected.

What is the long-term equilibrium normal distribution of the trait values in the population?

The long-term mean is given by the estimated parameter theta. The long-term variance is given by the sum sigmaG2tInf+sigmae^2 taken from the summary table of a POUMM-fit.

To what extent are the observed trait-values determined by heritable (i.e. genetic) versus non-heritable (i.e. environmental) factors?

In other words, what is the proportion of observable phenotypic variance attributable to the phylogeny? To answer this question, the POUMM package allows to estimate the phylogenetic heritability of the trait. Assuming that the tree represents the genetic relationship between individuals in a population, H_t̄² provides an estimate for the broad-sense heritability H² of the trait in the population.

When the goal is to estimate H_t̄² (H2tMean), it is imortant to specify an uninformed prior for it. Looking at the densities for chain 1 (red) on the previous figures, it becomes clear that the default prior favors values of H2tMean, which are either close to 0 or close to 1. Since by definition H_t̄² ∈ [0, 1], a reasonable uninformed prior for it is the standard uniform distribution. We set this prior by using the specifyPOUMM_ATH2tMeanSeG0 function. This specifies that the POUMM fit should be done on a parametrization  < α, θ, H_t̄², σ_e, g₀> rather than the standard parametrization  < α, θ, σ, σ_e, g₀>. It also specifies a uniform prior for H_t̄². You can explore the members of the specification list to see the different settings:

specH2tMean <- specifyPOUMM_ATH2tMeanSeG0(z[1:N], tree, nSamplesMCMC = 4e5)
# Mapping from the sampled parameters to the standard POUMM parameters:
specH2tMean$parMapping

## function (par) 
## {
##     if (is.matrix(par)) {
##         par[, 3] <- sigmaOU(par[, 3], par[, 1], par[, 4], tMean)
##         colnames(par) <- c("alpha", "theta", "sigma", "sigmae", 
##             "g0")
##     }
##     else {
##         par[3] <- sigmaOU(par[3], par[1], par[4], tMean)
##         names(par) <- c("alpha", "theta", "sigma", "sigmae", 
##             "g0")
##     }
##     par
## }
## <environment: 0x55ba00a6a278>

# Prior for the MCMC sampling
specH2tMean$parPriorMCMC

## function (par) 
## {
##     dexp(par[1], rate = tMean/6.931, log = TRUE) + dnorm(par[2], 
##         zMean, 2 * zSD, TRUE) + dunif(par[3], min = 0, max = 1, 
##         log = TRUE) + dexp(par[4], rate = 2/zSD, log = TRUE) + 
##         dnorm(par[5], zMean, 2 * zSD, log = TRUE)
## }
## <environment: 0x55ba00a6a278>

# Bounds for the maximum likelihood search
specH2tMean$parLower

##   alpha   theta H2tMean  sigmae      g0 
##  0.0000 -2.2073  0.0000  0.0000  0.5176

specH2tMean$parUpper

##   alpha   theta H2tMean  sigmae      g0 
## 10.8379  5.9389  0.9900  0.5335  3.2139

Then we fit the model:

fitH2tMean <- POUMM(z[1:N], tree, spec = specH2tMean)

plot(fitH2tMean, stat = c("H2tMean", "H2e", "H2tInf", "sigmae"), 
     doZoomIn = TRUE, doPlot = TRUE)

summary(fitH2tMean)[stat %in% c("H2tMean", "H2e", "H2tInf", "sigmae")]

##       stat     N    MLE PostMean           HPD   ESS  G.R.
##     <char> <int>  <num>    <num>        <list> <num> <num>
## 1: H2tMean   500 0.3057   0.3155 0.1788,0.4708   720     1
## 2:  sigmae   500 0.2013   0.2032 0.1830,0.2266   720     1
## 3:     H2e   500 0.4307   0.4183 0.2783,0.5296   720     1
## 4:  H2tInf   500 0.3073   0.3213 0.1781,0.4743   720     1

Now we see that the prior density for H2tMean is nearly uniform. It becomes clear that the process has converged to its long-term heritability since the intervals for H2tMean and H2tInf are nearly the same. Notice, though, that the estimate for the empirical heritability H2e is shifted towards 1 compared to H2tMean and H2tInf. This shows an important difference between H2e and the time-dependent formulae for phylogenetic heritability: H2e takes into account all values of z including those at the very beginning when the process was far away from equilibrium. Thus the estimated phenotypic variance over all trait-values at all times can be substantially bigger compared to the current trait-variance in the population:

# Compare global empirical heritability
H2eGlobal <- H2e(z[1:N], sigmae = coef(fitH2tMean)['sigmae'])
# versus recent empirical heritability
H2eRecent <- H2e(z[1:N], tree, sigmae = coef(fitH2tMean)['sigmae'], tFrom = 5)
print(c(H2eGlobal, H2eRecent))

## [1] 0.4307 0.3409

To learn more about different ways to specify the POUMM fit, read the documentation page ?specifyPOUMM.