phylopars: Estimation of phylogenetic and phenotypic covariance...

Description Usage Arguments Value Author(s) References Examples

View source: R/phylopars_update.R

Description

This function estimates parameters for the phylogenetic and phenotypic variance-covariance matrices for datasets with missing observations and multiple within-species observations. This function can also be used to fit altnerative evolutionary models, including Ornstein-Uhlenbeck, Early-Burst, star phylogeny, or Pagel's lambda, kappa, or delta. Reconstructed ancestral states and predicted species means (i.e., for missing data), along with prediction variances, are also provided.

Usage

1
2
3
4
5
6
7
8
phylopars(trait_data, tree, model = "BM", pheno_error, phylo_correlated = TRUE,
pheno_correlated = TRUE, REML = TRUE, full_alpha = TRUE, phylocov_start,
phenocov_start, model_par_start, phylocov_fixed, phenocov_fixed, model_par_fixed,
skip_optim = FALSE, skip_EM = FALSE, EM_Fels_limit = 1000, repeat_optim_limit = 1,
EM_missing_limit = 50, repeat_optim_tol = 0.01, model_par_evals = 10, max_delta = 10000,
EM_verbose = FALSE, optim_verbose = FALSE, npd = FALSE,
nested_optim = FALSE, usezscores = TRUE, phenocov_list = list(), ret_args = FALSE,
ret_level = 1)

Arguments

trait_data

A data frame with the first column labeled "species" (with species names matching tips on the phylogeny) and one column per trait. Each row corresponds to a single observation, and multiple observations for species are allowed. Missing data should be represented with NA.

tree

An object of class phylo

model

Model of evolution. Default is "BM". Alternative evolutionary models include "mvOU" (for the multivariate Ornstein-Uhlenbeck), or univariate tree transformations: "OU" "lambda", "kappa", "delta", "EB", "star".

pheno_error

If TRUE (default, unless <=1 observation per species is provided), parameters are estimated assuming within-species variation.

phylo_correlated

If TRUE (default), parameters are estimated assuming traits are correlated.

pheno_correlated

If TRUE (default), parameters are estimated assuming within-species observations traits are correlated.

REML

If TRUE (default), the algorithm will return REML estimates. If FALSE, maximum likelihood estimates will be returned.

full_alpha

Only applicable for the multivariate OU (model="mvOU"). If TRUE (default), a fully parametrized alpha matrix is fit. If FALSE, a diagonal alpha matrix is fit.

phylocov_start

Optional starting value for phylogenetic trait variance-covariance matrix. Must be of dimension n_traits by n_traits.

phenocov_start

Optional starting value for phenotypic trait variance-covariance matrix. Must be of dimension n_traits by n_traits.

model_par_start

Optional starting parameters for the evolutionary model. For model="mvOU", must be of dimension n_traits by n_traits. Otherwise, must be a single value.

phylocov_fixed

Optional fixed value for phylogenetic trait variance-covariance matrix. Must be of dimension n_traits by n_traits.

phenocov_fixed

Optional starting value for phenotypic trait variance-covariance matrix. Must be of dimension n_traits by n_traits.

model_par_fixed

Optional fixed parameter for the evolutionary model. For model="mvOU", must be of dimension n_traits by n_traits. Otherwise, must be a single value.

skip_optim

Whether to skip BFGS optimization (not recommended unless all parameters are fixed).

skip_EM

Whether to skip Expectation-Maximiation prior to generating starting parameters for BFGS optimization (not recommended unless providing fixed parameters).

EM_Fels_limit

Whether to skip Expectation-Maximiation prior to generating starting parameters for BFGS optimization (not recommended unless providing fixed parameters).

repeat_optim_limit

The number of times to repeat numerical optimization (default is 1).

EM_missing_limit

Maximum number of iterations for EM.

repeat_optim_tol

Maximum tolerance for repeated numerical optimization (only relevant if repeat_optim_limit>1).

model_par_evals

Number of times to evaluate univariate tree transformation models along the range of possible parameter values. Used to generate informed starting values for alternative evolutionary models if nested_optim=TRUE.

max_delta

Maximum allowed difference between the log-likelihood for EM-generated starting parameters and new parameters tried under numerical optimization. Extremely large deltas are likely to be numerical artifiacts. Prevents artificial convergence.

EM_verbose

Whether to print the log-likelihood during Expectation-Maximization.

optim_verbose

Whether to print log-likelihooods during numerical optimization.

npd

Whether to find the nearest positive-definite matrix for all covariance matrices during numerical optimization (slow – only set to TRUE if converging to singular matrices).

nested_optim

Only relevant if fitting a univariate alternative evolutionary model. Tries multiple tree transformation parameter values along the range of possible values to make informed starting parameters. Slower than the default (nested_optim=FALSE), in which all parameters are estimated simultaneously.

usezscores

Whether or not ot use centered and standardized data during numerical optimization (recommended).

phenocov_list

An optional named list of species-specific within-species covariance matrices to be held fixed, as in Ives et al (2007). This option forces pheno_error and pheno_correlated to be FALSE, and uses mean species values instead of raw data. Raw variance should be divided by the number of observations per species (i.e., squared standard errors). See Ives et al (2007) for more details.

ret_args

For internal use only.

ret_level

For internal use only.

Value

An object of class phylopars. For models with phenotypic (within-species) covariance, the estimated percentage of variance explained by the phylogeny is provided as 100*(1 - phenotypic_variance/raw_variance), where raw_variance is the variance of all observations for a given trait across species (var(PPE$trait_data[,2:ncol(PPE$trait_data)],na.rm=TRUE)).

logLik

The log-likelihood of the model

pars

A list composed of phylogenetic trait covariance and phenotypic (within-species) trait covariance, if estimated

model

The model of evolution (e.g., BM, OU, lambda, etc.), and any additional evolutionary model parameters estimated. For OU models, stationary covariance is calculated from both phylogenetic covariance (Sigma) and alpha (see Supplement 1 of Clavel et al. 2015).

mu

The estimate ancestral state at the root of the tree.

npars

The total number of parameters estimated by optimization (used for AIC and BIC).

anc_recon

Reconstructed ancestral states and species means. Row names correspond to species names (for the first 1:nspecies rows), and the remaining row names correspond to node numbers on a tree with edges in postorder: reorder(tree,"postorder"). Or, if node labels were included on the original tree, row names correspond to node labels.

anc_var

Variance of reconstructed ancestral estimates and imputed species means.

anc_cov

Covariance of estimates among variables.

tree

The phylogenetic tree supplied to phylopars

trait_data

The trait data supplied to phylopars

REML

TRUE if REML, FALSE if ML.

Author(s)

Eric W. Goolsby eric.goolsby.evolution@gmail.com, Cecile Ane, Jorn Bruggeman

References

Bruggeman J, Heringa J and Brandt BW. (2009) PhyloPars: estimation of missing parameter values using phylogeny. Nucleic Acids Research 37: W179-W184. Clavel, J., Escarguel, G. & Merceron, G. (2015) mvmorph: an r package for fitting multivariate 261 evolutionary models to morphometric data. Methods in Ecology and Evolution, 6, 131-1319. Felsenstein, J. (2008) Comparative methods with sampling error and within-species variation: contrasts revisited and revised. American Naturalist, 171, 713-725. Ho L.S.T., Ane C. 2014. A linear-time algorithm for Gaussian and non-Gaussian trait evolution models. Syst. Biol. 63:397-408.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
  # simulate data
  sim_data <- simtraits(ntaxa = 15,ntraits = 4,nreps = 3,nmissing = 10)

  # estimate parameters under Brownian motion
  # pheno_error = TRUE assumes intraspecific variation
  # pheno_correlated = FALSE assumes intraspecific variation is not correlated
  # phylo_correlated = TRUE assumed traits are correlated

  PPE <- phylopars(trait_data = sim_data$trait_data,tree = sim_data$tree,
    pheno_error = TRUE,phylo_correlated = TRUE,pheno_correlated = TRUE)

  PPE

  PPE$anc_recon # Ancestral state reconstruction and species mean prediction
  PPE$anc_var # Prediction variance


  ###NOT RUN
  # estimate parameters under multivariate OU
  # PPE_OU <- phylopars(trait_data = sim_data$trait_data,tree = sim_data$tree,
  #    model="mvOU",pheno_error = TRUE,phylo_correlated = TRUE,
  #    pheno_correlated = TRUE)
  #
  # PPE

Example output

Loading required package: ape
Phylogenetic trait variance-covariance
           V1         V2        V3        V4
V1 0.10379427 0.05966994 0.1357337 0.1112332
V2 0.05966994 0.16409266 0.2665444 0.2531963
V3 0.13573368 0.26654445 0.4513081 0.4203385
V4 0.11123322 0.25319631 0.4203385 0.3951567

Phenotypic trait variance-covariance
             V1           V2          V3         V4
V1  0.104495542  0.004192296 -0.09677586  0.1060050
V2  0.004192296  0.314348742  0.13756073 -0.1197072
V3 -0.096775863  0.137560725  0.87295034 -0.8562792
V4  0.106004998 -0.119707185 -0.85627920  3.2660553

% variance explained by phylogeny
      V1       V2       V3       V4 
59.52468 37.63952 35.64670 11.87310 
Brownian motion model
             V1        V2         V3          V4
t5  -0.03517166 0.3268175 -0.6965103  0.14653408
t6  -0.38614568 0.1187140 -1.1646840 -0.23882900
t1   0.76332696 1.1865019  0.9296104  1.58644454
t2   0.12548190 0.3022739 -0.6562137  0.14824441
t4   0.51414585 0.9465803  0.4633422  1.17844612
t14  0.77567630 1.0286782  0.7062168  1.35919931
t15  0.77818946 1.0282062  0.7067193  1.35909771
t7   0.59678740 1.5138613  1.3263588  2.02490624
t8   0.46607999 1.4094018  1.1128478  1.84208207
t3   0.75416993 0.6919872  0.2070215  0.86323973
t12  0.82750661 0.7853073  0.3772326  1.01643025
t13  0.85169474 0.8573505  0.4933064  1.12712450
t11  0.88398245 0.8491062  0.4965950  1.12263942
t9   0.64577357 1.1394214  0.8056583  1.49035710
t10  0.92645948 1.3048322  1.1785964  1.79706265
16   0.38831450 0.7543641  0.1246734  0.86879931
17  -0.15346024 0.2735303 -0.8298250  0.04122459
18   0.49263459 0.8469499  0.3084645  1.02815089
19   0.56698452 0.8866756  0.4013111  1.10343018
20   0.47296319 0.8808391  0.3483881  1.07297433
21   0.54675521 1.0609938  0.6449378  1.35288749
22   0.56112180 1.0986599  0.7064375  1.41116289
23   0.76997708 1.0307054  0.7064671  1.36082498
24   0.54596227 1.2840028  0.9684735  1.67787903
25   0.69564756 0.9249983  0.5177949  1.18934091
26   0.75400129 0.8452204  0.4295062  1.08663467
27   0.84233985 0.8327384  0.4531361  1.08905322
28   0.84202843 0.8314412  0.4511048  1.08708909
29   0.77437274 1.1835565  0.9305542  1.58472821
            V1         V2         V3        V4
t5  0.02209032 0.04164304 0.09948191 0.1541468
t6  0.02202042 0.04041334 0.09790343 0.1485468
t1  0.02593652 0.04741085 0.11095631 0.1610100
t2  0.02841900 0.04077052 0.09443886 0.1547051
t4  0.01941861 0.03426139 0.07987398 0.1339546
t14 0.01328146 0.02740495 0.06566090 0.1172173
t15 0.01309577 0.02599836 0.06170628 0.1144529
t7  0.01794783 0.03194664 0.07549261 0.1292823
t8  0.01794783 0.03194664 0.07549261 0.1292823
t3  0.02103666 0.03677128 0.08420351 0.1390258
t12 0.01561955 0.02974116 0.07130258 0.1300162
t13 0.01559632 0.02822002 0.06772646 0.1230019
t11 0.01579407 0.02852024 0.06839364 0.1235626
t9  0.02319306 0.03253153 0.07710173 0.1305820
t10 0.01894228 0.03260428 0.07576673 0.1302076
16  0.07418290 0.12071590 0.33312634 0.3502643
17  0.03204499 0.05574134 0.14585396 0.1900152
18  0.05602086 0.09191988 0.25386521 0.2808070
19  0.03787354 0.06293471 0.17406615 0.2115244
20  0.02755752 0.04625716 0.12576386 0.1708683
21  0.01860569 0.03290448 0.08821847 0.1369269
22  0.01777569 0.03169581 0.08477574 0.1337525
23  0.01289326 0.02610984 0.06265231 0.1147760
24  0.01848721 0.03278492 0.08588603 0.1353731
25  0.03137829 0.05256282 0.14407020 0.1858091
26  0.02740788 0.04664802 0.12613465 0.1710584
27  0.01428010 0.02633484 0.06660430 0.1214690
28  0.01408934 0.02606453 0.06572364 0.1208437
29  0.02323859 0.03727594 0.09404785 0.1442024

Rphylopars documentation built on Feb. 6, 2021, 9:14 a.m.