Performs Phylogenetic linear regression evaluating intraspecific variability in response and/or predictor variables and uncertainty in trees topology.

tree_intra_phylm(
  formula,
  data,
  phy,
  Vy = NULL,
  Vx = NULL,
  y.transf = NULL,
  x.transf = NULL,
  n.intra = 10,
  n.tree = 2,
  distrib = "normal",
  model = "lambda",
  track = TRUE,
  ...
)

Arguments

formula

The model formula: response~predictor.

data

Data frame containing species traits and species names as row names.

phy

A phylogeny (class 'phylo', see ?ape).

Vy

Name of the column containing the standard deviation or the standard error of the response variable. When information is not available for one taxon, the value can be 0 or NA.

Vx

Name of the column containing the standard deviation or the standard error of the predictor variable. When information is not available for one taxon, the value can be 0 or NA

y.transf

Transformation for the response variable (e.g. log or sqrt). Please use this argument instead of transforming data in the formula directly (see also details below).

x.transf

Transformation for the predictor variable (e.g. log or sqrt). Please use this argument instead of transforming data in the formula directly (see also details below).

n.intra

Number of times to repeat the analysis generating a random value for response and/or predictor variables. If NULL, n.intra = 30

n.tree

Number of times to repeat the analysis with n different trees picked randomly in the multiPhylo file. If NULL, n.tree = 2

distrib

A character string indicating which distribution to use to generate a random value for the response and/or predictor variables. Default is normal distribution: "normal" (function rnorm). Uniform distribution: "uniform" (runif) Warning: we recommend to use normal distribution with Vx or Vy = standard deviation of the mean.

model

The phylogenetic model to use (see Details). Default is lambda.

track

Print a report tracking function progress (default = TRUE)

...

Further arguments to be passed to phylolm

Value

The function tree_intra_phylm returns a list with the following components:

formula: The formula

data: Original full dataset

sensi.estimates: Coefficients, aic and the optimised value of the phylogenetic parameter (e.g. lambda) for each regression using a value in the interval of variation and a different phylogenetic tree.

N.obs: Size of the dataset after matching it with tree tips and removing NA's.

stats: Main statistics for model parameters.CI_low and CI_high are the lower and upper limits of the 95

all.stats: Complete statistics for model parameters. Fields coded using all describe statistics due to both intraspecific variation and phylogenetic uncertainty. Fields coded using intra describe statistics due to intraspecific variation only. Fields coded using tree describe statistics due to phylogenetic uncertainty only. sd is the standard deviation. CI_low and CI_high are the lower and upper limits of the 95

sp.pb: Species that caused problems with data transformation (see details above).

Details

This function fits a phylogenetic linear regression model using phylolm to n trees (n.tree), randomly picked in a multiPhylo file. The regression is also repeated n.intra times. At each iteration the function generates a random value for each row in the dataset using the standard deviation or errors supplied and assuming a normal or uniform distribution. To calculate means and se for your raw data, you can use the summarySE function from the package Rmisc.

#' All phylogenetic models from phylolm can be used, i.e. BM, OUfixedRoot, OUrandomRoot, lambda, kappa, delta, EB and trend. See ?phylolm for details.

Currently, this function can only implement simple linear models (i.e. \(trait~ predictor\)). In the future we will implement more complex models.

Output can be visualised using sensi_plot.

Warning

When Vy or Vx exceed Y or X, respectively, negative (or null) values can be generated, this might cause problems for data transformation (e.g. log-transformation). In these cases, the function will skip the simulation. This problem can be solved by increasing times, changing the transformation type and/or checking the target species in output$sp.pb.

References

Paterno, G. B., Penone, C. Werner, G. D. A. sensiPhy: An r-package for sensitivity analysis in phylogenetic comparative methods. Methods in Ecology and Evolution 2018, 9(6):1461-1467

Martinez, P. a., Zurano, J.P., Amado, T.F., Penone, C., Betancur-R, R., Bidau, C.J. & Jacobina, U.P. (2015). Chromosomal diversity in tropical reef fishes is related to body size and depth range. Molecular Phylogenetics and Evolution, 93, 1-4

Ho, L. S. T. and Ane, C. 2014. "A linear-time algorithm for Gaussian and non-Gaussian trait evolution models". Systematic Biology 63(3):397-408.

See also

Examples

# Load data: data(alien) # run PGLS accounting for intraspecific and phylogenetic variation: intra.tree <- tree_intra_phylm(gestaLen ~ adultMass, data = alien$data, phy = alien$phy, Vy = "SD_gesta", n.intra = 3, n.tree = 3, y.transf = log, x.transf = log)
#> Warning: distrib = normal: make sure that standard deviation is provided for Vx and/or Vy
#> Warning: NA's in response or predictor, rows with NA's were removed
#> Warning: Some phylo tips do not match species in data (this can be due to NA removal) species were dropped from phylogeny or data
#> Used dataset has 84 species that match data and phylogeny
#> | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%
# To check summary results: summary(intra.tree)
#> mean.all CI_low_all CI_high_all mean.intra CI_low_intra #> intercept 2.367 2.316 2.418 2.367 2.366 #> se.intercept 0.340 0.339 0.341 0.340 0.339 #> pval.intercept 0.000 0.000 0.000 0.000 0.000 #> estimate 0.147 0.140 0.153 0.147 0.147 #> se.estimate 0.021 0.020 0.022 0.021 0.021 #> pval.estimate 0.000 0.000 0.000 0.000 0.000 #> CI_high_intra mean.tree CI_low_tree CI_high_tree #> intercept 2.368 2.367 2.366 2.368 #> se.intercept 0.340 0.340 0.339 0.340 #> pval.intercept 0.000 0.000 0.000 0.000 #> estimate 0.147 0.147 0.147 0.147 #> se.estimate 0.021 0.021 0.021 0.021 #> pval.estimate 0.000 0.000 0.000 0.000
# Visual diagnostics sensi_plot(intra.tree, uncer.type = "all") #or uncer.type = "tree", uncer.type = "intra"