Skip to content

Example 5: Within Species Variation

Eric W. Goolsby edited this page May 11, 2016 · 6 revisions

Like the models described in Felsenstein (2008) and Bruggeman et al. (2009), Rphylopars can directly incorporate raw individual-level data (multiple within-species observations). To demonstrate this feature, we simulate a 10-species dataset with 3 traits and multiple observations per species, with 50 data points missing completely at random. In this example, each species has 3 intraspecific observations, but note that there is no requirement for the number of within-species observations to be identical across species (i.e., any species may have 0, 1, or more observations).

library(Rphylopars)
library(phytools) # for simulating pure-birth phylogenies
set.seed(1) # Set the seed for reproducible results
trait_cov <- matrix(c(4,2,2.2,2,3,1.5,2.2,1.5,2.5),nrow = 3,ncol = 3)
trait_cov # Phylogenetic trait covariance
within_species_cov <- diag(c(2,6,4))
within_species_cov # Simulated within-species covariance
tree <- pbtree(n = 100)
sim_data_intra <- simtraits(v = trait_cov,intraspecific = within_species_cov,tree = tree,
            nreps = 3,nmissing = 50)

Note that the simulated within-species (phenotypic) covariance matrix is a diagonal matrix, so that within-species observations are uncorrelated. We may fit a model assuming no within-species correlation by setting pheno_correlated=FALSE (by default, pheno_correlated=TRUE). Rphylopars automatically detects multiple within-species observations and estimates within-species covariance (pheno_correlated=TRUE). Here again, AIC or BIC can be used to pick the best-supported model (correlated or uncorrelated within-species observations).

p_BM_uncorrelated_pheno <- phylopars(trait_data = sim_data_intra$trait_data,
                           tree = sim_data_intra$tree,pheno_correlated = FALSE)
p_BM_correlated_pheno <- phylopars(trait_data = sim_data_intra$trait_data,
                           tree = sim_data_intra$tree,pheno_correlated = TRUE)

BIC(p_BM_correlated_pheno)
BIC(p_BM_uncorrelated_pheno)

p_BM_uncorrelated_pheno # best supported model by BIC

The % variance explained by the phylogeny is obtained as 100 - the ratio of the phenotypic variance to the raw variance, where the raw variance is calculated among all the observations in the usual way (without reference to the phylogeny).

Next, we simulate a dataset under the same conditions as above, but with correlated within-species (phenotypic) observations. That is, the simulated within-species (phenotypic) covariance matrix is a full (non-diagonal) matrix.

set.seed(1) # Set the seed for reproducible results
trait_cov <- matrix(c(4,2,2.2,2,3,1.5,2.2,1.5,2.5),nrow = 3,ncol = 3)
trait_cov # Phylogenetic trait covariance
within_species_cov <- matrix(c(2,1.5,1.8,1.5,6,3.5,1.8,3.5,4),nrow = 3,ncol = 3)
within_species_cov # Simulated within-species covariance
tree <- pbtree(n = 100)
sim_data_intra <- simtraits(v = trait_cov,intraspecific = within_species_cov,tree = tree,
            nreps = 3,nmissing = 50)

p_BM_uncorrelated_pheno <- phylopars(trait_data = sim_data_intra$trait_data,
                           tree = sim_data_intra$tree,pheno_correlated = FALSE)
p_BM_correlated_pheno <- phylopars(trait_data = sim_data_intra$trait_data,
                           tree = sim_data_intra$tree,pheno_correlated = TRUE)

BIC(p_BM_correlated_pheno)
BIC(p_BM_uncorrelated_pheno)

p_BM_correlated_pheno  # best supported model by BIC

Models with within-species variation can also incorporate alternative evolutionary models by setting model to the desired evolutionary model (e.g., "OU", "EB", "mvOU").

References

You can’t perform that action at this time.