Models Of Trait Macroevolution On Trees

Description

Functions for fitting models of trait evolution on phylogenies for continuous traits. The majority of functions described in Thomas and Freckleton (2012) and include functions that allow for tests of variation in the rates of trait evolution.

Author(s)

Mark Puttick <marknputtick@gmail.com>.

Gavin Thomas.

Travis Ingram.

Magnus Clarke.

Rob Freckleton.

David Orme.

Emmanuel Paradis.

References

Thomas GH, & Freckleton R. 2012. Body size diversification in Anolis: novel environments and island effects. MOTMOT: models of trait macroevolution on trees 3, 145-151.

Maximum likelihood rate estimation for traits and phylogenies

Description

Full function for maximum likelihood estimation of rate parameters and comparison to a single rate model.

Usage

ML.RatePhylo(
  rateData,
  rate = NULL,
  fixed = NULL,
  pretty = TRUE,
  rateMIN = 0.001,
  rateMAX = 50,
  common.mean = FALSE,
  lambda.est = TRUE,
  est.CI = FALSE,
  meserr = FALSE,
  file = NULL
)

Arguments

Value

If pretty=FALSE, returns a list containing:

• MLRate Maximum likelihood estimates of the rate parameters

• Lambda Maximum likelihood estimate of lambda

• LCI Approximate lower confidence intervals for rate

• UCI Approximate upper confidence intervals for rate parameters

• means Means for each category

• nParam Number of parameters in the model (how many means and rate categories)

• Max.lik Maximum (log) likeihood

• AIC for maximum likelihood model

• AICc for maximum likelihood model

• LambdaSingle Maximum likelihood estimate of lambda for the single rate model

• Lik1 Likelihood of the equivalent single rate model

• Likelihood ratio statistic of "Max.lik" vs "Lik1"

• P P values for the LR statistic

• df Degrees of freedom for the LR statistic

• AIC.rate1 AIC for single rate model

• AICc.rate1 AICc for single rate model

If pretty=TRUE, prints a nice version of the list to screen. If file is specified the pretty output will be sent to file, not the console.

Note

Unlike phyloMean and likRatePhylo (that use treatment contrasts), the means reported here are the actual values

Author(s)

Gavin Thomas, Rob Freckleton

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

## Read in phylogeny and data from Thomas et al. (2009)
data(anolis.tree)
data(anolis.data)
## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph",
data=anolis.data)
anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)
# A model with a different rate in one of the four groups. The 'fixed' command is used to determine
# whether a particular rate is to be constrained or not. Use '1' to fix a group and 'FALSE' to show
# that the parameter is not fixed and should be estimated. The values should be entered in the same 
# order as the ranking of the groups. That is, group 0 (small islands) takes position one in the 
# fixed vector, group 1 (large island trunk crown and trunk ground) takes position 2 and so on. 
# The default is to allow each group to take a different mean. 
ML.RatePhylo(anolis.rateData, fixed=c(1, FALSE, FALSE, FALSE), pretty=TRUE, lambda.est = FALSE)

Confidence intervals for rate parameters

Description

Calculates approximate confidence intervals for all rate parameters. CIs are esimated for one rate parameters while fixing others at a given value (usually the maximum likelihood estimate).

These are reliable (given the asympotic assumptions of the chi-square distribution) if only two groups are being compared but should be regarded only as a rough approximation for =>3 different rate categories. If the rates are correlated the CIs may be underestimated.

Usage

RatePhylo.allCI(
  rateData,
  MLrate = NULL,
  fixed = NULL,
  rateMIN = 0.001,
  rateMAX = 50,
  common.mean = FALSE,
  lambda.est = TRUE
)

Arguments

Value

rateLci Lower confidence interval for rate estimate

rateUci Upper confidence interval for rate estimate

Author(s)

Gavin Thomas, Rob Freckleton

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

data(anolis.data)
data(anolis.tree)
## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)

# A model with a different rate in each of the four groups. The 'fixed' command is used to determine
# whether a particular rate is to be constrained or not. Use '1' to fix a group and 'FALSE' to show 
# that the parameter is not fixed and should be estimated. The values should be entered in the same
# order as the ranking of the groups. That is, group 0 (small islands) takes position one in the 
# fixed vector, group 1 (large island trunk crown and trunk ground) takes position 2 and so on. The 
# default is to allow each group to take a different mean. 

anole.ML <- optim.likRatePhylo(rateData=anolis.rateData, rate=NULL,
fixed=c(FALSE,FALSE,FALSE, FALSE),
common.mean=FALSE, lambda.est=FALSE)

# Confidence intervals for the first two parameters

RatePhylo.allCI(rateData=anolis.rateData, MLrate = anole.ML$MLRate, 
fixed=c(FALSE, TRUE, TRUE, TRUE), rateMIN = 0.001, rateMAX = 50, 
common.mean = FALSE)

add a fossil to an interior branch of a time-scaled phylogeny

Description

the function takes a time-scaled phylogeny and adds a fossil at a branch in the past

Usage

addFossilToPhy(
  phy,
  inGroup,
  fossil,
  fossilAge,
  minLength = 0.1,
  maxLength = NULL
)

Arguments

Value

the the time-scaled phylogeny with the fossil attached

Author(s)

Mark Puttick

References

Puttick, M. N., Kriwet, J., Wen, W., Hu, S., Thomas, G. H., & Benton, M. J. (2017). Body length of bony fishes was not a selective factor during the biggest mass extinction of all time. Palaeontology, 60, 727-741.

Examples

data(anolis.tree)
plot(anolis.tree)
nodelabels(214, 214)
# add fossil to node 214
in.groups <- node.descendents(x=214, phy=anolis.tree, tip.labels=TRUE)[[2]]
fossilPhy <- addFossilToPhy(anolis.tree, in.groups, fossil="fakeFossil", fossilAge=60)
plot(fossilPhy)

Finch allopatric data

Description

Finch data from Clarke et al. 2017

Usage

data(finches)

Format

An object of class "data.frame".

References

Clarke M, Thomas GH, Freckleton RP. 2017. Trait evolution in adaptive radiations: modelling and measuring interspecific competition on phylogenies. The American Naturalist. 189, 121-137.

Examples

data(finches)
head(allopatric.data)

Anolis phenotype data

Description

Data on anolis phenotype data from Thomas et al. 2009

Usage

data(anolis.data)

Format

An object of class "data.frame".

References

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

data(anolis.data)
head(anolis.data)

Anolis phylogeny

Description

Data on anolis phylogeny from Thomas et al. 2009

Usage

data(anolis.tree)

Format

An object of class "phylo".

References

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

data(anolis.tree)
anolis.tree

Conversion among data and phylogeny objects

Description

Function to generate a "rateData" object containing the discrete explanatory variable, continuous response variable and set of variance co-variance matrices. It loads the trait data and removes species with missing data from the data and vcv matrix. as.rateData requires either a set of matrices in rateMatrix format created using as.rateMatrix or, if no rateMatrix object is input then it requires a phylogeny in "phylo" format. If a "phylo" object is used as.rateData will call as.rateMatrix internally. as.rateMatrix calls the "ape" function vcv.phylo multiple times and this can be slow for large phylogenies. It will often be more efficient to use as.rateMatrix first to create a "rateMatrix" object to pass to as.rateData, particularly if there are many response traits of interest to be fitted to the same phylogeny and set of reconstructed ancestral states.

Usage

as.rateData(
  y,
  x,
  rateMatrix = NULL,
  phy = NULL,
  data,
  meserr.col = NULL,
  meserr.propn = NULL,
  log.y = FALSE,
  report_prune = FALSE
)

Arguments

Value

rateData An object of class "rateData" which is a list containing the response (y) and explanatory (x) variable along with a list of variance-covaraince matrices.

Author(s)

Gavin Thomas

References

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

## Read in phylogeny and data from Thomas et al. (2009)
data(anolis.tree)
data(anolis.data)

## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)  

## Convert data to class rateData with a phylo object as input 
anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = NULL, phy=anolis.tree, data=anolis.data, log.y=TRUE)

Conversion among data and phylogeny objects

Description

Function to generate a rateMatrix object containing a set of variance covariance matrices. Note that as.rateMatrix calls the CAIC function vcv.array multiple times and this can be slow for large phylogenies (though faster than using the ape equivalent vcv.phylo).

Usage

as.rateMatrix(phy, x, data)

Arguments

Value

rateMatrix An object of class "rateMatrix" - a list of matrices describing the expected variances and covariances of between species. Each matrix refers to the variances and covariances for a given state of x (see Thomas et al. 2006).

Author(s)

Gavin Thomas

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Examples

## Read in phylogeny and data from Thomas et al. (2009)
data(anolis.tree)
data(anolis.data)

## Convert data to class rateMatrix
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

Blomberg's K Estimate Blomberg's K (Blomberg et al. 2003)

Description

Blomberg's K Estimate Blomberg's K (Blomberg et al. 2003)

Usage

blomberg.k(phy, y)

Arguments

Value

The estimate of the K statistic

References

Blomberg SP, Garland T, & Ives AR. 2003. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 57, 717-745.

See Also

transformPhylo.ML, the Picante package

Calculate multiple-test cut-off

Description

Calculate the log-likelihood, AIC, or AICc cut-off necessary for type-one error to reach acceptable levels

Usage

calcCutOff(
  phy,
  n = 1000,
  mc.cores = 1,
  model,
  measure = "AICc",
  alpha.error = 0.05,
  ...
)

Arguments

Value

The cut-off requred to produce an type-one error rate equal to quantile.cut.off (default = 0.05) when data are simulated under Brownian motion, and these data are analysed under the appropriate model.

Author(s)

Mark Puttick

See Also

transformPhylo.ML, transformPhylo.ll, transformPhylo, transformPhylo.MCMC

Examples

data(anolis.tree)
set.seed(393)
# calculated necessary AICc cut-off to reduce type-one error to 5% 
# for a timeSlice model with a split at 30Ma (only 5 simulations used,
# it's recommend to use 1000 for analyses)
calcCutOff(anolis.tree, n=5, model="timeSlice", splitTime=30)

Character displacement likelihood ratio test

Description

Conducts a likelihood ratio test between empirical data (phylogeny and trait data), and simumlations from the function chr.disp.sim using an approximate Bayesian computation (ABC) approach (Clarke et al. 2017)

Usage

chr.disp.lrt(emp.tree, emp.data, param.out, posteriorSize = 500)

Arguments

Value

List containing element of 'estimates' with the estimates of sigma and a, with the Brownian motion (a = 0) summarised in column one and the character displacement (a > 0) in column two. 'likelihood' contains the likelihood of the Brownian motion model and the character displacement model, and the likelihood ratio test estimate. If used, there is an estimate of Blomberg's K for the empirical and simulated data.

Author(s)

Magnus Clarke and Mark Puttick

References

Clarke M, Thomas GH, Freckleton RP. 2017. Trait evolution in adaptive radiations: modelling and measuring interspecific competition on phylogenies. The American Naturalist. 189, 121-137.

See Also

chr.disp.sim, chr.disp.param

Examples

## import finch data form Clarke et al. (2017)
data(finches)
## simulate small amount of data 
## (example only - many more datasets are required for accuracy)
param.simulation <- chr.disp.param(finch.tree, n.sim = 100, n.steps=100,
max.sigma = 8, max.a = 8, ntraits=1, 
allopatry=as.matrix(allopatric.data), mc.cores = 1)
chr.disp.lrt(finch.tree, finch.data, param.simulation, 50)

Simulate character displacement data wrapper

Description

Simulates phylogenetic trait data under a character displacement model (Clarke et al. 2017) in which traits interact inter-specifically, with competition between sympatric lineages driving trait values apart

Usage

chr.disp.param(
  phy,
  n.sim = 100,
  n.steps = 1000,
  max.sigma = 8,
  max.a = 8,
  est.blomberg.k = FALSE,
  ntraits = 1,
  sympatry = NA,
  allopatry = NA,
  trait.lim = NA,
  mc.cores = 1
)

Arguments

Value

List containing the simulated data 'simulated.param': a matrix with each row represented an iteration, the sigma (Brownian variance) used in the iteration, the 'a' value used in each iteration, the mean and standard deviation between neighbouring trait values. The 'input.arguments' from the model, the 'input.phy' from the model, and the input 'sympatry' and 'allopatry' matrices.

Author(s)

Magnus Clarke and Mark Puttick

References

Clarke M, Thomas GH, & Freckleton RP. 2017. Trait evolution in adaptive radiations: modeling and measuring interspecific competition on phylogenies. The American Naturalist 189, 121-137.

See Also

chr.disp.sim, chr.disp.lrt

Examples

## import finch data form Clarke et al. (2017)
data(finches)
## simulate small amount of data 
## (example only - many more datasets are required for accuracy)
param.simulation <- chr.disp.param(finch.tree, n.sim = 3, n.steps=100,  
max.sigma = 8, max.a = 8, ntraits=1, 
allopatry=as.matrix(allopatric.data), mc.cores = 1)

Simulate character displacement data

Description

Simulates data under a Brownian motion or character displacement model

Usage

chr.disp.sim(
  phy,
  n.steps = 1000,
  sigma = 1,
  a = 0,
  ntraits = 1,
  sympatry = NA,
  allopatry = NA,
  trait.lim = NA
)

Arguments

Value

A list containing the the simulated data (tval) showing the sigma, a, mean gap and gap standard deviation. Additionally, if used, the user input sympatry (symp) and/or allopatry (allo) matrices

Author(s)

Magnus Clarke and Mark Puttick

References

Clarke M, Thomas GH, Freckleton RP. 2017. Trait evolution in adaptive radiations: modelling and measuring interspecific competition on phylogenies. The American Naturalist. 189, 121-137.

See Also

chr.disp.param, chr.disp.lrt

Examples

## import finch data form Clarke et al. (2017)
data(finches)
emp.tree <- finch.tree
emp.data <- finch.data
## simulate small amount of data 
## (example only - many more datasets are required for accuracy)
sim.data <- chr.disp.sim(emp.tree, n.steps=100,
sigma=1, a=2, ntraits=1, sympatry=NA, allopatry=NA, trait.lim=NA)

prune tree and data to lineages present in a time bin in the past

Description

the function takes a full tree and returns a pruned phylogeny with only tips and lineages found within a time bin preserved. If trait data are supplied the function will return tip states based either on the original tips found in the bin, or tip states inferred from ancestral states

Usage

contemporaryPhy(
  phy,
  maxBin,
  minBin,
  reScale = 0,
  allTraits,
  closest.min = TRUE,
  traits.from.tip = TRUE
)

Arguments

Value

the pruned phylogeny. The object descendants refers to the lineages the branch in the time bin gave rise to before it was pruned. If traits are included a vector of trait values representing species at the tips.

Author(s)

Mark Puttick

References

Puttick, M. N., Kriwet, J., Wen, W., Hu, S., Thomas, G. H., & Benton, M. J. (2017). Body length of bony fishes was not a selective factor during the biggest mass extinction of all time. Palaeontology, 60, 727-741.

Examples

## prune a random tree to taxa present between 4 and 2 units before present
# generate tree
set.seed(20)
tree <- rtree(20)
# generate traits
traits <- rnorm(20)
# plot tree and timeframe
plot(tree)
max.age <- nodeTimes(tree)[1,1]
abline(v=max.age - c(4, 2))
# prune tree to timeframe
cont.tree <- contemporaryPhy(phy=tree, maxBin=4, minBin=2, allTraits=traits)
plot(cont.tree$phy)

Drop tips from a phylogenetic tree while preserving deleted nodes

Description

Wrapper for the ape function drop.tip that preserves the number of nodes affecting each branch. For use with the psi and multipsi models.

Usage

dropTipPartial(phy, tip)

Arguments

Value

Phylogenetic tree in phylo format, with an added element Shid, a vector of numbers of observed but "missing" speciation events per branch, in the same order as the branches in the phylo object

Author(s)

Travis Ingram

References

Ingram, T. 2011. Speciation along a depth gradient in a marine adaptive radiation. Proc. R. Soc. B 278: 613-618.

See Also

transformPhylo.ML

Examples

## Read in phylogeny and data from Thomas et al. (2009)
data(anolis.tree)
data(anolis.data)
## identify tips to drop
tips.to.go <- anolis.tree$tip.label[1:30]
dropTipPartial(phy=anolis.tree, tip=tips.to.go)

Calculate fair proportions phylogenetic diversity metric

Description

Calculate fair proportions phylogenetic diversity metric Note that as.rateMatrix calls the CAIC function vcv.array multiple times and this can be slow for large phylogenies (though faster than using the "ape" equivalent vcv.phylo).

Usage

fairProportions(phy, nodeCount = FALSE)

Arguments

Value

Returns a matrix of fair proportion for all tips in phylogeny and node counts if selected.

Author(s)

Gavin Thomas

References

Redding, D.W. and Mooers, A.O. (2006). Incorporating evolutionary measures into conservation prioritisation. Conservation Biology, 20, 1670-1678.

Isaac, N.J.B., Turvey, S.T., Collen, B., Waterman, C. and Baillie, J.E.M. (2007). Mammals on the EDGE: conservation priorities based on threat and phylogeny. PLoS ONE, 2, e296.

Examples

data(anolis.tree)

fp <- fairProportions(anolis.tree)
fpNodes <- fairProportions(anolis.tree, nodeCount=TRUE)

Finch phenotype data

Description

Finch data from Clarke et al. 2017

Usage

data(finches)

Format

An object of class "data.frame".

References

Clarke M, Thomas GH, Freckleton RP. 2017. Trait evolution in adaptive radiations: modelling and measuring interspecific competition on phylogenies. The American Naturalist. 189, 121-137.

Examples

data(finches)
head(finch.data)

Finch tree

Description

Finch data from Clarke et al. 2017

Usage

data(finches)

Format

An object of class "phylo"

References

Clarke M, Thomas GH, Freckleton RP. 2017. Trait evolution in adaptive radiations: modelling and measuring interspecific competition on phylogenies. The American Naturalist. 189, 121-137.

Examples

data(finches)
head(finch.tree)

Log-likelihood rate estimation for traits and phylogenies

Description

This function calculates the log-likelihood, phylogenetic mean, and Brownian variance for a trait and a phylogeny transformed according to variation in relative rates.

Usage

likRatePhylo(
  rateData,
  rate = NULL,
  common.mean = FALSE,
  lambda.est = TRUE,
  lambda = 1,
  meserr = FALSE,
  sigmaScale = NULL
)

Arguments

Value

ll log-likelihood of the model

mu phylogenetically corrected mean(s)

s2 Brownian variance

Note

The means are output as treatment contrasts.

Author(s)

Gavin Thomas, Rob Freckleton

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

data(anolis.tree)
data(anolis.data)

## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)

## Calculate phylogenetic mean, variance, log likelihood for a model where the first 

# mean only
phyloMean(rateData=anolis.rateData, rate = c(1,2,0.1,1), common.mean = FALSE)

# variance only
phyloVar(rateData=anolis.rateData, rate = c(1,2,0.1,1), common.mean = FALSE) 

# mean, variance and log-likelihood
likRatePhylo(rateData=anolis.rateData, rate = c(1,2,0.1,1), common.mean = FALSE) 

Log-likelihood estimation for traits and phylogenies

Description

This function calculates the log-likelihood and Brownian (co)variance for a trait(s) and a phylogeny using phylogenetically independent contrasts Note that as.rateMatrix calls the CAIC function vcv.array multiple times and this can be slow for large phylogenies (though faster than using the "ape" equivalent vcv.phylo).

Usage

likTraitPhylo(y, phy, covPIC = TRUE, brCov = NULL)

Arguments

Details

The phylo object must be rooted and fully dichotomous

Value

brownianVariance Brownian variance (or covariance for multiple traits) given the data and phylogeny

logLikelihood The log-likelihood of the model and data

Author(s)

Gavin Thomas, Rob Freckleton

References

Felsenstein J. 1973. Maximum-likelihood estimation of evolutionary trees from continuous characters. Am. J. Hum. Genet. 25, 471-492.

Felsenstein J. 1985. Phylogenies and the comparative method. American Naturalist 125, 1-15.

Freckleton RP & Jetz W. 2009. Space versus phylogeny: disentangling phylogenetic and spatial signals in comparative data. Proc. Roy. Soc. B 276, 21-30.

Examples

data(anolis.tree)
data(anolis.data)
## calculate Brownian variance log-likelihood of female SVL
female.svl <- matrix(anolis.data[,"Female_SVL"], 
dimnames=list(rownames(anolis.data)))
input.data <- sortTraitData(phy=anolis.tree, y=female.svl, log.trait=TRUE)
likTraitPhylo(phy = input.data$phy, y=input.data$trait)

Mammal data

Description

Mammal data from Slater 2013

Usage

data(mammals)

Format

An list containing mass and errors of mammal body mass and mammal phylogeny of class "phylo"

References

Phylogenetic evidence for a shift in the mode of mammalian body size evolution at the Cretaceous-Palaeogene boundary. Methods in Ecology and Evolution. 4: 734-744.

Examples

data(mammals)

plot the output from transformPhylo.MCMC

Description

Plots a histogram of the estimated parameter and a trace of the results

Usage

mcmc.plot(
  mcmc.input,
  y.limit = NULL,
  x.limit = NULL,
  label.text = NULL,
  cex.axis = 1,
  cex.labels = 0.7,
  col.hist = "green4",
  col.trace = "navy"
)

Arguments

Value

Two plots showing the histogram of the estimated parameter value and a trace of the MCMC estimation

Author(s)

Mark Puttick

Examples

library(motmot)
data(anolis.tree)
data(anolis.data)
attach(anolis.data)
male.length <- matrix(Male_SVL, dimnames=list(rownames(anolis.data)))
sortedData <- sortTraitData(anolis.tree, male.length)
phy <- sortedData$phy
male.length <- sortedData$trait
phy.clade <- extract.clade(phy, 182)
male.length.clade <- as.matrix(male.length[match(phy.clade$tip.label, rownames(male.length)),])
## not run
# please note, this model will be need to run for longer to achieve convergence
# lambda.mcmc <- transformPhylo.MCMC(y=male.length.clade, phy=phy.clade, 
# model="lambda", mcmc.iteration=100, burn.in=0.1)
# mcmc.plot(lambda.mcmc)

Identify nodes and tips descended from a node

Description

Obtains a vector of the tips and nodes subtending from a node in a phylogeny.

Usage

node.descendents(x, phy, tip.labels = FALSE)

Arguments

Details

This function is stolen from the clade.members function in the CAIC package but returns both node and tip id's.

Value

Returns a vector of node and tip ids descended from the tip(s) "x". If tip.labels=TRUE then returns a list of node ids and tip labels.

Note

as.rateMatrix calls the CAIC function vcv.array multiple times and this can be slow for large phylogenies (though faster than using the ape equivalent vcv.phylo).

Author(s)

Gavin Thomas, David Orme

Get times for nodes and tips

Description

Produces branching and tip times for ultrametric and non-ultrametric trees

Usage

nodeTimes(phy)

Arguments

Value

Returns a matrix corresponging the phy "edge" matrix showning internal and external node times

Note

nodeTimes is essentially a re-written version of the ape branching.times.

Author(s)

Mark Puttick, Emmanuel Paradis

Examples

## Read in phylogeny from Thomas et al. (2009)
data(anolis.tree)
anolis.node.times <- nodeTimes(phy=anolis.tree)

Maximum likelihood rate estimation for traits and phylogenies

Description

Function for the maximum likelihood estimation of rate parameters on a trait and phylogeny.

Usage

optim.likRatePhylo(
  rateData,
  rate = NULL,
  fixed = NULL,
  rateMIN = 0.001,
  rateMAX = 50,
  common.mean = FALSE,
  lambda.est = TRUE,
  meserr = FALSE
)

Arguments

Value

MLRate Maximum likelihood estimates of the rate parameters

Max.lik Maximum (log) likeihood

AIC AIC for maximum likelihood model

AICc AICc for maximum likelihood model

convergence convergence value from optim

n.parameters Number of parameters in the model (how many means and rate categories)

Author(s)

Gavin Thomas

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

data(anolis.tree)
data(anolis.data)

## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)

# A model with a different rate in each of the four groups. The 'fixed' command is used to determine
# whether a particular rate is to be constrained or not. Use '1' to fix a group and 'FALSE' to show
# that the parameter is not fixed and should be estimated. The values should be entered in the same 
# order as the ranking of the groups. That is, group 0 (small islands) takes position one in the 
# fixed vector, group 1 (large island trunk crown and trunk ground) takes position 2 and so on. 
# The default is to allow each group to take a different mean. 

optim.likRatePhylo(anolis.rateData, rate=c(1,1,1,1), common.mean=TRUE, lambda.est=FALSE)

Calculation of Brownian (co)variance using independent contrasts.

Description

Calculates the Brownian variance (single trait) or variance-covariance matrix (mutliple traits) using phylogenetically independent contrasts.

Usage

phyloCovar(x, phy, estimator = "unbiased")

Arguments

Value

brownianVariance Brownian variance (or covariance for multiple traits) given the data and phylogeny

Author(s)

Gavin Thomas, Rob Freckleton

References

Felsenstein J. 1973. Maximum-likelihood estimation of evolutionary trees from continuous characters. Am. J. Hum. Genet. 25, 471-492.

Felsenstein J. 1985. Phylogenies and the comparative method. American Naturalist 125, 1-15.

Freckleton RP & Jetz W. 2009. Space versus phylogeny: disentangling phylogenetic and spatial signals in comparative data. Proc. Roy. Soc. B 276, 21-30.

Examples

data(anolis.tree)
data(anolis.data)
## calculate Brownian variance of female SVL
female.svl <- matrix(anolis.data[,"Female_SVL"], 
dimnames=list(rownames(anolis.data)))
input.data <- sortTraitData(phy=anolis.tree, y=female.svl, log.trait=TRUE)
phyloCovar(x=input.data$trait, phy = input.data$phy)

Calculation of phylogenetically corrected mean.

Description

This function calculates the phylogenetic mean of the data given the tree and model of evolution

Usage

phyloMean(
  rateData,
  rate = NULL,
  common.mean = FALSE,
  lambda.est = TRUE,
  lambda = 1,
  meserr = FALSE
)

Arguments

Value

mu phylogenetically corrected mean

Note

The means are output as treatment contrasts.

Author(s)

Gavin Thomas, Rob Freckleton

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

 ## Read in phylogeny and data from Thomas et al. (2009)
data(anolis.tree)
data(anolis.data)

## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)

# A model with a different rate in each of the four groups. The 'fixed' command is used to determine
# whether a particular rate is to be constrained or not. Use '1' to fix a group and 'FALSE' to show
# that the parameter is not fixed and should be estimated. The values should be entered in the same 
# order as the ranking of the groups. That is, group 0 (small islands) takes position one in the 
# fixed vector, group 1 (large island trunk crown and trunk ground) takes position 2 and so on. 
# The default is to allow each group to take a different mean. 

phyloMean(anolis.rateData, rate=c(1,1,1,1), common.mean=FALSE)
# common mean for all groups
phyloMean(anolis.rateData, rate=c(1,1,1,1), common.mean=TRUE)

Calculation of Brownian variance.

Description

This function calculates the phylogenetic variance (Brownian variance, or rate) of the data given the tree and model of evolution

Usage

phyloVar(
  rateData,
  rate = NULL,
  common.mean = FALSE,
  lambda.est = TRUE,
  lambda = 1,
  meserr = FALSE
)

Arguments

Value

phylo.var phylogenetic variance (Brownian variance)

Author(s)

Gavin Thomas, Rob Freckleton

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624. @references Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

 ## Read in phylogeny and data from Thomas et al. (2009)
data(anolis.tree)
data(anolis.data)

## Convert data to class rateData with a rateMatrix object as input
anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)

anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)

# A model with a different rate in each of the four groups. The 'fixed' command is used to determine
# whether a particular rate is to be constrained or not. Use '1' to fix a group and 'FALSE' to show
# that the parameter is not fixed and should be estimated. The values should be entered in the same 
# order as the ranking of the groups. That is, group 0 (small islands) takes position one in the 
# fixed vector, group 1 (large island trunk crown and trunk ground) takes position 2 and so on. 
# The default is to allow each group to take a different mean. 

phyloVar(anolis.rateData, rate=c(1,2,1,1), common.mean=FALSE)
# common mean for all groups
phyloVar(anolis.rateData, rate=c(1,2,1,1), common.mean=TRUE)

Fast PLGS estimation based on contrasts

Description

Estimates regression parameters for a phylogenetic generalised least-squares analysis using the fast constrasts method (Felsenstein 1973; 1985; Freckleton 2012). This implementation is applicable for continuous traits only and not factors

Usage

pic.pgls(
  formula,
  phy,
  y,
  lambda = "ML",
  return.intercept.stat = FALSE,
  meserr = NULL
)

Arguments

Value

A list containing the model, model summary, intercept, estimate of Lambda, model log-Likelihood, model AIC

Author(s)

Mark N Puttick and Rob Freckleton

See Also

pgls

Examples

# Data and phylogeny
data(anolis.tree)
anolis.tree$node.label <- NULL
lm.data <- transformPhylo.sim(phy=anolis.tree, n=2, model="bm")
dat <- data.frame(x = lm.data[,1], y = lm.data[,2],
names = anolis.tree$tip, row.names = anolis.tree$tip)
picModel <- pic.pgls(formula=y ~  x,
phy=anolis.tree, y = dat, lambda=1,
return.intercept.stat=FALSE)

Identify shifts in the rate of trait diversification through time

Description

Summarises phenotypic rate variation on phylogenies through

Usage

## S3 method for class 'timeSlice.ML'
plot(
  x,
  ...,
  cutoff = 4,
  AICc = TRUE,
  lowerBound = 1e-08,
  upperBound = 1000,
  phylo.plot = TRUE,
  colour.ramp = c("blue", "red"),
  cex.plot = 1,
  model.average = FALSE
)

Arguments

Details

This functions summarises the output of a "timeSlice" model in transformPhylo.ML (see below). The best overall model is chosen based on AIC (or AICc if AICc=TRUE). The cut-off point for improvement in AIC score between successively more complex models can be defined using cutoff. The default cutoff is 4 but this is somewhat arbitrary and a "good" cut-off may well vary between data sets so it may well be worth exploring different cutoffs.

Value

ModelFit Summary of the best optimal rate shift model or the model average of each split time (if model averaging was used).

Rates Summary of the rate parameters from the best rate shift model or the model averaged rates through time.

optimalTree A phylo object with branch lengths scaled relative to rate and a plot of estimated rates through time with their associated CIs.

Author(s)

Mark Puttick

References

To Add

See Also

transformPhylo.ML

Examples

data(anolis.tree)
data(anolis.data)
attach(anolis.data)
male.length <- matrix(Male_SVL, dimnames=list(rownames(anolis.data)))
sortedData <- sortTraitData(anolis.tree, male.length)
phy <- sortedData$phy
male.length <- sortedData$trait
phy.clade <- extract.clade(phy, 182)
male.length.clade <- as.matrix(male.length[match(phy.clade$tip.label, 
rownames(male.length)),])
timeSlice.10.ml <- transformPhylo.ML(y=male.length.clade, phy=phy.clade, model="timeSlice", 
splitTime=c(10))
outputSummary <- plot(timeSlice.10.ml, cutoff=0.001, cex=0.5, 
colour.ramp=c("blue", "red"))

Tree plotting for rates

Description

Plots trees with colours based on rates of trait evolution. Also provides simple coloured plotting for trait values using the ace function in the ape library.

Usage

## S3 method for class 'traitMedusa.model'
plot(x, y = NULL, ..., reconType = "rates", palette = "hotspot.colors")

Arguments

Value

Returns a data frame of colours used in plot along with rate (or ancestral state) range for each colour.

Author(s)

Gavin Thomas, Mark Puttick

See Also

transformPhylo.ML, summary.traitMedusa.

Examples

# Data and phylogeny
data(anolis.tree)
data(anolis.data)

# female SVL data
female.svl <- matrix(anolis.data[,"Female_SVL"], 
dimnames=list(rownames(anolis.data)))
input.data <- sortTraitData(phy=anolis.tree, y=female.svl, log.trait=TRUE)

# arbitarily reduce data size for speed in this example
phy.clade <- extract.clade(input.data[[1]], 182)
male.length.clade <- as.matrix(input.data[[2]][match(input.data[[1]]$tip.label, 
rownames(input.data[[2]])),])
# Identify rate shifts and print and plot results with up to one rate shifts 
# and minimum clade size of 10.
anolisSVL_MEDUSA <- transformPhylo.ML(male.length.clade, phy=phy.clade, 
model="tm1",minCladeSize=10, nSplits=1)
anolisSVL_MEDUSA_out <- summary(anolisSVL_MEDUSA, cutoff=1, AICc=FALSE)
colours <- plot(x = anolisSVL_MEDUSA_out,
reconType = "rates", type = "fan", cex=0.6, edge.width=3)

Sample hidden speciation events along branches of a tree

Description

Uses estimated speciation and extinction rates to sample the number of speciation events 'hidden' by subsequent extinction on each branch of a tree following Bokma (2008). For use with the psi and multipsi models.

Usage

sampleShid(phy, la = NULL, mu = NULL, useMean = FALSE)

Arguments

Details

The expected number of hidden speciation events are calculated for each branch given its start and end times, and estimates of lambda and mu which are assumed to be constant across the tree. To properly account for uncertainty in the effect of extinction on the number of nodes affecting each branch of a tree, it may be appropriate to repeat model-fitting on many realizations of Sobs on the tree of interest (similar to evaluating phylogenetic uncertainty)

Value

Phylogenetic tree in phylo format, with an added element Sobs, a vector of numbers of hidden speciation events per branch, in the same order as the branches in the phylo object

Author(s)

Travis Ingram

References

Bokma, F. 2008. Detection of "punctuated equilibrium" by Bayesian estimation of speciation and extinction rates, ancestral character states, and rates of anagenetic and cladogenetic evolution on a molecular phylogeny. Evolution 62: 2718-2726.

Ingram, T. 2011. Speciation along a depth gradient in a marine adaptive radiation. Proc. R. Soc. B 278: 613-618.

See Also

transformPhylo.ML

Sort data and remove missing entries for tree and trait data

Description

Plots a phylogeny with lines representing the value of a continuous trait

Usage

sortTraitData(
  phy,
  y,
  data.name = NULL,
  log.trait = TRUE,
  pass.ultrametric = FALSE
)

Arguments

Value

phy Tree with missing data pruned

trait Rearranged data with missing species removed

Author(s)

Mark Puttick

Examples

data(anolis.tree)
data(anolis.data)
attach(anolis.data)
male.length <- matrix(Male_SVL, dimnames=list(rownames(anolis.data)))
any(is.na(male.length[,1]))
data.sorted <- sortTraitData(anolis.tree, male.length)
phy <- data.sorted[[1]]
male.length <- data.sorted[[2]]

Identify shifts in the rate of trait diversification

Description

Summarises phenotypic rate variation on phylogenies.

Usage

## S3 method for class 'traitMedusa'
summary(
  object,
  ...,
  cutoff = 4,
  AICc = TRUE,
  lowerBound = 1e-08,
  upperBound = 200,
  print.warnings = FALSE
)

Arguments

Details

This functions summarises the output of a "medusa" model in transformPhylo.ML (see below). The best overall model is chosen based on AIC (or AICc if AICc=TRUE). The cut-off point for improvement in AIC score between successively more complex models can be defined using cutoff. The default cutoff is 4 but this is somewhat arbitrary and a "good" cut-off may well vary between data sets so it may well be worth exploring different cutoffs.

Value

ModelFit Summary of the best optimal rate shift model.

Rates Summary of the rate parameters from the best rate shift model.

optimalTree A phylo object with branch lengths scaled relative to rate.

Author(s)

Gavin Thomas

References

Alfaro ME, Santini F, Brock CD, Alamillo H, Dornburg A, Carnevale G, Rabosky D & Harmon LJ. 2009. Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates. PNAS 106, 13410-13414.

O'Meara BC, Ane C, Sanderson MJ & Wainwright PC. 2006. Testing for different rates of continuous trait evolution using likelihood. Evolution 60, 922-933

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

See Also

transformPhylo.ML

Examples

data(anolis.tree)
data(anolis.data)
attach(anolis.data)
male.length <- matrix(Male_SVL, dimnames=list(rownames(anolis.data)))
sortedData <- sortTraitData(anolis.tree, male.length)
phy <- sortedData$phy
male.length <- sortedData$trait
phy.clade <- extract.clade(phy, 182)
male.length.clade <- as.matrix(male.length[match(phy.clade$tip.label, rownames(male.length)),])
tm1 <- transformPhylo.ML(male.length.clade, phy=phy.clade, model="tm1", minCladeSize=10, nSplits=1)
tm1_out <- summary(tm1, cutoff=1)

plot a univariate continuous trait data on a phylogeny

Description

Plots a phylogeny with lines representing the value of a continuous trait

Usage

traitData.plot(
  y,
  phy,
  col.label = "red",
  col.tree = "black",
  col.hist = "navy",
  cex.plot = 0.7,
  cex.tips = 0.7,
  show.tips = FALSE,
  include.hist = FALSE,
  n.split = 5,
  lwd.traits = 1,
  show.axis = TRUE,
  axis.text = NULL,
  transform.axis.label = NULL,
  at = NULL,
  labels = NULL,
  axis.text.line = 1,
  offset.bars = 1,
  ...
)

Arguments

Value

A plot with the trait values shown at the tips, and a histrogram of the trait values

Author(s)

Mark Puttick

Examples

data(anolis.tree)
data(anolis.data)
attach(anolis.data)
male.length <- matrix(Male_SVL, dimnames=list(rownames(anolis.data)))
sortedData <- sortTraitData(anolis.tree, male.length)
phy <- sortedData$phy
male.length <- sortedData$trait
traitData.plot(y=male.length, phy)

Phylogenetic tree transformations

Description

Transforms the branch lengths of a phylo object according to a model of trait evolution (see details).

Usage

transformPhylo(
  phy,
  model = NULL,
  y = NULL,
  meserr = NULL,
  kappa = NULL,
  lambda = NULL,
  delta = NULL,
  alpha = NULL,
  psi = NULL,
  lambda.sp = NULL,
  nodeIDs = NULL,
  rateType = NULL,
  branchRates = NULL,
  cladeRates = NULL,
  splitTime = NULL,
  timeRates = NULL,
  acdcRate = NULL,
  branchLabels = NULL,
  cophenetic.dist = NULL,
  vcv.matrix = NULL,
  mode.order = NULL,
  mode.param = NULL,
  rate.var = NULL
)

Arguments

Details

Transforms the branch lengths of a phylo object according to one of the following models:

• model="bm"- Brownian motion (constant rates random walk)

• model="kappa" - fits Pagel's kappa by raising all branch lengths to the power kappa. As kappa approaches zero, trait change becomes focused at branching events. For complete phylogenies, if kappa approaches zero this infers speciational trait change.

• model="lambda" - fits Pagel's lambda to estimate phylogenetic signal by multiplying all internal branches of the tree by lambda, leaving tip branches as their original length (root to tip distances are unchanged);

• model="delta" - fits Pagel's delta by raising all node depths to the power delta. If delta <1, trait evolution is concentrated early in the tree whereas if delta >1 trait evolution is concentrated towards the tips. Values of delta above one can be difficult to fit reliably.

• model="free" - fits Mooer's et al's (1999) free model where each branch has its own rate of trait evolution. This can be a useful exploratory analysis but it is slow due to the number of parameters, particularly for large trees.

• model="clade" - fits a model where particular clades are a priori hypothesised to have different rates of trait evolution (see O'Meara et al. 2006; Thomas et al. 2006, 2009). Clades are specified using nodeIDs and are defined as the mrca node. Unique rates for each clade are specified using cladeRates. rateType specifies whether the rate shift occurs in the stem clade or on the single branch leading to the clade.

• model="OU" - fits an Ornstein-Uhlenbeck model - a random walk with a central tendency proportional to alpha. High values of alpha can be interpreted as evidence of evolutionary constraints, stabilising selection or weak phylogenetic signal.

• model="psi" - fits a model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate (Ingram 2010). Note that 'original nodes' from the full phylogeny can be included as an element on the phylogeny (e.g., phy$orig.node) as well as estimates of 'hidden' speciation (e.g., phy$hidden.speciation) if estimates of extinction (mu) are > 0.

• model="multiPsi" fits a model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate but allows seperate values of psi fitted to seperate branches (Ingram 2010; Ingram et al. 2016). Note that 'original nodes' from the full phylogeny can be included as an element on the phylogeny (e.g., phy$orig.node) as well as estimates of 'hidden' speciation (e.g., phy$hidden.speciation) if estimates of extinction (mu) are > 0.

• model="ACDC" fits a model to in which rates can exponentially increased or decrease through time (Blomberg et al. 2003). If the upper bound is < 0, the model is equivalent to the 'Early Burst' model of Harmon et al. 2010. If a nodeIDs is supplied, the model will fit a ACDC model nested within a clade, with a BM fit to the rest of the tree.

• model="timeSlice" A model in which all branch rates change at time(s) in the past.

• model="modeSlice" A model in which all branch modes change at a time or times set a priori by the user.

Value

phy A phylo object with branch lengths scaled according to the given model and parameters

Author(s)

Gavin Thomas, Mark Puttick

References

Ingram T. 2011. Speciation along a depth gradient in a marine adaptive radiation. Proc. Roy. Soc. B. 278, 613-618.

Ingram T, Harrison AD, Mahler L, Castaneda MdR, Glor RE, Herrel A, Stuart YE, and Losos JB. Comparative tests of the role of dewlap size in Anolis lizard speciation. Proc. Roy. Soc. B. 283, 20162199.

Mooers AO, Vamosi S, & Schluter D. 1999. Using phylogenies to test macroevolutionary models of trait evolution: sexual selection and speciation in Cranes (Gruinae). American Naturalist 154, 249-259.

O'Meara BC, Ane C, Sanderson MJ & Wainwright PC. 2006. Testing for different rates of continuous trait evolution using likelihood. Evolution 60, 922-933

Pagel M. 1997. Inferring evolutionary processes from phylogenies. Zoologica Scripta 26, 331-348.

Pagel M. 1999 Inferring the historical patterns of biological evolution. Nature 401, 877-884.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

See Also

transformPhylo.ML, transformPhylo.ll, transformPhylo.MCMC

Examples

data(anolis.tree)
anolis.treeDelta <- transformPhylo(phy=anolis.tree, model="delta", delta=0.5)

Bayesian MCMC for models of trait evolution

Description

Fits Bayesian models for various models of continuous character evolution using a Metropolis-Hastings Markov Chain Monte Carlo (MCMC) approach

Usage

transformPhylo.MCMC(
  y,
  phy,
  model,
  mcmc.iteration = 1000,
  burn.in = 0.1,
  hiddenSpeciation = FALSE,
  full.phy = NULL,
  lowerBound = NULL,
  upperBound = NULL,
  useMean = FALSE,
  random.start = TRUE,
  meserr = NULL,
  covPIC = TRUE,
  lambdaEst = FALSE,
  nodeIDs = NULL,
  branchLabels = NULL,
  acdcScalar = FALSE,
  sample.every = 10
)

Arguments

Details

The method estimates posterior probabilities using a Metropolis-Hastings MCMC approach that places a prior bounded uniform distribution on all parameters with an independence sampler. These prior distributions can be altered by changing the upperBound and lowerBound arguments. The MCMC model will estimate the posterior probability for the following models:

• model="kappa" fits Pagel's kappa by raising all branch lengths to the power kappa. As kappa approaches zero, trait change becomes focused at branching events. For complete phylogenies, if kappa approaches zero this infers speciational trait change. Default bounds from ~0 - 1.

• model="lambda" fits Pagel's lambda to estimate phylogenetic signal by multiplying all internal branches of the tree by lambda, leaving tip branches as their original length (root to tip distances are unchanged). Default bounds from ~0 - 1.

• model="delta" fits Pagel's delta by raising all node depths to the power delta. If delta <1, trait evolution is concentrated early in the tree whereas if delta >1 trait evolution is concentrated towards the tips. Values of delta above one can be difficult to fit reliably. Default bounds from ~0 - 5.

• model="OU" fits an Ornstein-Uhlenbeck model - a random walk with a central tendency proportional to alpha. High values of alpha can be interpreted as evidence of evolutionary constraints, stabilising selection or weak phylogenetic signal. It is often difficult to distinguish among these possibilities. Default bounds from ~0 - 10.

• model="psi" fits a acceleration-deacceleration model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate (Ingram 2010).

• model="ACDC" fits a model to in which rates can exponentially increased or decrease through time (Blomberg et al. 2003). If the upper bound is < 0, the model is equivalent to the 'Early Burst' model of Harmon et al. 2010. Default rate parameter bounds from ln(1e-10) ~ ln(20) divided by the root age.

Value

median The median estimate of the posterior for the parameter

95.HPD The 95 percent Highest Posterior Density for the parameter

ESS Effective Sample Size for the posterior

acceptance.rate The ratio for which new proposals were accepted during the MCMC chain

mcmc.chain Full MCMC chain containing all iterations (including burn-in)

Author(s)

Mark Puttick, Gavin Thomas

See Also

transformPhylo.ML, transformPhylo.ll, transformPhylo

Examples

data(anolis.tree)
data(anolis.data)
attach(anolis.data)
male.length <- matrix(Male_SVL, dimnames=list(rownames(anolis.data)))
sortedData <- sortTraitData(anolis.tree, male.length)
phy <- sortedData$phy
male.length <- sortedData$trait
phy.clade <- extract.clade(phy, 182)
male.length.clade <- as.matrix(male.length[match(phy.clade$tip.label, rownames(male.length)),])
## please note, this model will be need to run for longer to achieve convergence
lambda.mcmc <- transformPhylo.MCMC(y=male.length.clade, phy=phy.clade, 
model="lambda", mcmc.iteration=100, burn.in=0.1)

Maximum likelihood for models of trait evoluion

Description

Fits likelihood models for various models of continuous character evolution. Model fitting is based on maximum-likelihood evaluation using phylogenetically independent contrasts. This is exactly equivalent to, but substantially faster than, GLS approaches.

Usage

transformPhylo.ML(
  y,
  phy,
  model = NULL,
  modelCIs = TRUE,
  nodeIDs = NULL,
  rateType = NULL,
  minCladeSize = 1,
  nSplits = 2,
  splitTime = NULL,
  boundaryAge = 10,
  testAge = 1,
  restrictNode = NULL,
  lambdaEst = FALSE,
  acdcScalar = FALSE,
  branchLabels = NULL,
  hiddenSpeciation = FALSE,
  full.phy = NULL,
  useMean = FALSE,
  profilePlot = FALSE,
  lowerBound = NULL,
  upperBound = NULL,
  covPIC = TRUE,
  n.cores = 1,
  tol = NULL,
  meserr = NULL,
  controlList = c(fnscale = -1, maxit = 100, factr = 1e-07, pgtol = 0, type = 2, lmm =
    5),
  returnPhy = FALSE,
  print.warnings = FALSE,
  mode.order = NULL,
  rate.var = FALSE,
  testShiftTimes = NULL,
  saveAll = TRUE
)

Arguments

Details

This function finds the maximum likelihood parameter values for continuous character evolution. For "kappa", "delta", "OU", "multipsi", and "ACDC" it is possible to fit a 'nested' model of evolution in which the ancestral rate of BM switches to a different node, as specified by nodeIDs or branchLabels for multipsi. The function returns the maximum-likelihood parameter estimates for the following models.

• model="bm" Brownian motion (constant rates random walk).

• model="kappa" fits Pagel's kappa by raising all branch lengths to the power kappa. As kappa approaches zero, trait change becomes focused at branching events. For complete phylogenies, if kappa approaches zero this infers speciational trait change. Default bounds from ~0 - 1.

• model="lambda" fits Pagel's lambda to estimate phylogenetic signal by multiplying all internal branches of the tree by lambda, leaving tip branches as their original length (root to tip distances are unchanged). Default bounds from ~0 - 1.

• model="delta" fits Pagel's delta by raising all node depths to the power delta. If delta <1, trait evolution is concentrated early in the tree whereas if delta >1 trait evolution is concentrated towards the tips. Values of delta above one can be difficult to fit reliably. If a nodeIDs is supplied, the model will fit a delta model nested within a clade, with a BM fit to the rest of the tree. Default bounds from ~0 - 5.

• model="OU" fits an Ornstein-Uhlenbeck model - a random walk with a central tendency proportional to alpha. High values of alpha can be interpreted as evidence of evolutionary constraints, stabilising selection or weak phylogenetic signal. It is often difficult to distinguish among these possibilities. If a nodeIDs is supplied, the model will fit a OU model nested within a clade, with a BM fit to the rest of the tree. For OU models, alternative optimisation are performed with different starting values (1e-8, 0.01, 0.1, 1, 5). Default bounds from ~0 - 10.

• model="ACDC" fits a model to in which rates can exponentially increased or decrease through time (Blomberg et al. 2003). If the upper bound is < 0, the model is equivalent to the 'Early Burst' model of Harmon et al. 2010. If a nodeIDs is supplied, the model will fit a ACDC model nested within a clade, with a BM fit to the rest of the tree. Default rate parameter bounds from ln(1e-10) ~ ln(20) divided by the root age. Note this process starts on the stem branch leading to the MRCA of the common node, unlike the other methods that start at the common node.

• model="trend" fits a model in which the expectated mean change through time is non-zero, signifying a directional evolution to a larger or smaller trait value. This model is only appliacble to non-ultrametric trees.

• model="psi" fits a model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate (Ingram 2010). Note that the algorithm will automatically estimate speciation and extinction estimates, and will incorporate estimates of 'hidden' speciation if death estimates are greater than 0.

• model="multipsi" fits a model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate but allows seperate values of psi fitted to seperate branches (Ingram 2010; Ingram et al. 2016). Note that the algorithm will automatically estimate speciation and extinction estimates, and will incorporate estimates of 'hidden' speciation if death estimates are greater than 0.

• model="free" fits Mooers et al's free model where each branch has its own rate of trait evolution. This can be a useful exploratory analysis but it is slow due to the number of parameters, particularly for large trees. Default rate parameter bounds from ~0 - 200.

• model="clade" fits a model where particular clades are a priori hypothesised to have different rates of trait evolution (see O'Meara et al. 2006; Thomas et al. 2006, 2009). Clades are specified using nodeIDs and are defined as the mrca node. Default rate parameter bounds from ~0 - 200.

• model="tm1" fits "clade" models without any a priori assertion of the location of phenotypic diversification rate shifts. It uses the same AIC approach as the runMedusa function in the geiger package (runMedusa tests for shifts in the rate of lineage diversification). The algorithm first fits a constant-rate Brownian model to the data, it then works iteratively through the phylogeny fitting a two-rate model at each node in turn. Each two-rate model is compared to the constant rate model and the best two-rate model is retained. Keeping the location of this rate shift intact, it then repeats the procedure for a three-rate model and so on. The maximum number of rate shifts can be specified a priori using nSplits. Limits can be applied to the size (species richness) of clades on which to infer new rate shifts using minCladeSize. This can be useful to enable large trees to be handled but should be used cautiously since specifiying a large minimum clade size may result in biologically interesting nested rate shifts being missed. Equally, very small clade sizes may provide poor estimates of rate that may not be informative. Limits on the search can also be placed using restrictNode. This requires a list where each element of the list is a vector of tip names that define monophyletic groups. Rate shifts will not be searched for within any of the defined groups. Default rate parameter bounds from ~0 - 1000.

• model="tm2" this model is similar to "tm1", however, at each node it assesses the fit of two models. The first model is exactly as per "tm1". The second model infers a rate shift on the single branch descending directly from a node but not on any of the descending branches thereafter. Only the best fitting single-branch or whole clade model is retained for the next iteration. If a single-branch shift is favoured, this infers either that there was a rapid shift in trait value along the stem leading to the crown group, or that the members of the clade have undergone parallel shifts. In either case, this can be considered as a change in mean, though separating a single early shift from a clade-parallel shift is not possible with this method.

• model="timeSlice" A model in which all branch rates change at a time or times set a priori by the user. IfDefault rate parameter bounds from ~0 - 1000. If splitTime=NULL, all 1 Ma (as defined by test Age) intervals from the root of the tree - 10 and the youngest tip + 10 will be included in the search. The +/- 10 Ma age can be modified using the argument boundaryAge. At each stage the best fitting model will be stored, and the search will continue until n shifts, with n shifts defined by nSplits. If a single value or vector is used for splitTime, only these ages are included in the search.

• model="modeslice" A model in which all branch modes change at a time or times set a priori by the user.

Value

Returns the maximum log-likelihood and parameter estimates (with 95 percent confidence intervals if specified). If model="bm" instead returns the Brownian (co)variance and log-likelihood. Also returned are the root estimate, the AIC, and AICc.

traitMedusaObject A list in which the first element contains a matrix summarising the parameter estimates and node ids, log-likelihoods, number of parameters (k), AIC and AICc for the best one-rate model, two-rate model, three rate model and so on. The second element is a sub-list where the first element contains all two-rate models, the second element contains all three-rate models and so on. This can be summarised using traitMedusaSummary. The third element is the input trait data. The fourth element is the input phylogeny.

Note

Confidence intervals are based on the assumption of an asymptotic Chi-square distribution. For multi-parameter models (e.g. rate shift models with more than two rates) the confidence intervals are approximate and are calculated for each parameter in turn while holding all other parameters at their maximum likelihood value.

Author(s)

Gavin Thomas, Mark Puttick

References

Alfaro ME, Santini F, Brock CD, Alamillo H, Dornburg A, Carnevale G, Rabosky D & Harmon LJ. 2009. Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates. PNAS 106, 13410-13414.

Blomberg SP, Garland T & Ives AR 2003. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 57, 717-745.

Felsenstein J. 1973. Maximum-likelihood estimation of evolutionary trees from continuous characters. Am. J. Hum. Genet. 25, 471-492.

Felsenstein J. 1985. Phylogenies and the comparative method. American Naturalist 125, 1-15.

Freckleton RP & Jetz W. 2009. Space versus phylogeny: disentangling phylogenetic and spatial signals in comparative data. Proc. Roy. Soc. B 276, 21-30.

Harmon LJ et al. 2010. Early bursts of body size and shape evolution are rare in comparative data. Evolution 57, 717-745.

Ingram T. 2011. Speciation along a depth gradient in a marine adaptive radiation. Proc. Roy. Soc. B. 278, 613-618.

Ingram T,Harrison AD, Mahler L, Castaneda MdR, Glor RE, Herrel A, Stuart YE, and Losos JB. 2016. Comparative tests of the role of dewlap size in Anolis lizard speciation. Proc. Roy. Soc. B. 283, 20162199.

Mooers AO, Vamosi S, & Schluter D. 1999. Using phylogenies to test macroevolutionary models of trait evolution: sexual selection and speciation in Cranes (Gruinae). American Naturalist 154, 249-259.

O'Meara BC, Ane C, Sanderson MJ & Wainwright PC. 2006. Testing for different rates of continuous trait evolution using likelihood. Evolution 60, 922-933

Pagel M. 1997. Inferring evolutionary processes from phylogenies. Zoologica Scripta 26, 331-348.

Pagel M. 1999 Inferring the historical patterns of biological evolution. Nature 401, 877-884.

Pagel M. 1999 Inferring the historical patterns of biological evolution. Nature 401, 877-884.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

See Also

transformPhylo.MCMC, transformPhylo, transformPhylo.ll, blomberg.k

Examples

# Data and phylogeny
data(anolis.tree)
data(anolis.data)
sortedData <- sortTraitData(anolis.tree, anolis.data,
data.name="Male_SVL", log.trait=TRUE, pass.ultrametric=TRUE)
phy <- sortedData$phy
male.length <- sortedData$trait
phy.clade <- extract.clade(phy, 182)
male.length.clade <- as.matrix(male.length[match(phy.clade$tip.label, rownames(male.length)),])

# Brownian motion model
transformPhylo.ML(male.length.clade , phy=phy.clade, model="bm")

# Delta
transformPhylo.ML(male.length.clade , phy=phy.clade, model="delta", upperBound=2)

# The upper confidence interval for kappa is outside the bounds so try increasing
# the upper bound

transformPhylo.ML(male.length.clade , phy=phy.clade, model="delta", upperBound=5)

# Test for different rates in different clades - here with 2 hypothesised
# unusual rates compared to the background

# This fits the non-censored model of O'Meara et al. (2006)
phy.clade$node.label[which(phy.clade$node.label == "3")] <- 2
transformPhylo.ML(male.length.clade, phy=phy.clade, model="clade", nodeIDs=c(49, 54))

# Identify rate shifts and print and plot results with upto three rate shifts
# and minimum clade size of 20.
anolisSVL_MEDUSA <- transformPhylo.ML(male.length.clade, phy=phy.clade, model="tm1",
minCladeSize=10, nSplits=2)

Log-likelhood for models of trait evolution.

Description

Fits likelihood models for various models of continuous character evolution.

Usage

transformPhylo.ll(
  y = NULL,
  phy,
  model = NULL,
  meserr = NULL,
  kappa = NULL,
  lambda = NULL,
  delta = NULL,
  alpha = NULL,
  psi = NULL,
  lambda.sp = NULL,
  nodeIDs = NULL,
  rateType = NULL,
  branchRates = NULL,
  cladeRates = NULL,
  timeRates = NULL,
  splitTime = NULL,
  branchLabels = NULL,
  acdcRate = NULL,
  covPIC = TRUE,
  cophenetic.dist = NULL,
  vcv.matrix = NULL,
  mode.order = NULL,
  mode.param = NULL,
  rate.var = NULL,
  mu = NULL,
  sigma.sq = NULL
)

Arguments

Details

This function fits likelihood models (see below) for continuous character evolution where the parameter values are set a priori. The function returns the log-likihood and the Brownian variance (or variance covariance matrix).

• model="bm" Brownian motion (constant rates random walk).

• model="kappa" fits Pagel's kappa by raising all branch lengths to the power kappa. As kappa approaches zero, trait change becomes focused at branching events. For complete phylogenies, if kappa approaches zero this infers speciational trait change. Default bounds from ~0 - 1.

• model="lambda" fits Pagel's lambda to estimate phylogenetic signal by multiplying all internal branches of the tree by lambda, leaving tip branches as their original length (root to tip distances are unchanged). Default bounds from ~0 - 1.

• model="delta" fits Pagel's delta by raising all node depths to the power delta. If delta <1, trait evolution is concentrated early in the tree whereas if delta >1 trait evolution is concentrated towards the tips. Values of delta above one can be difficult to fit reliably. If a nodeIDs is supplied, the model will fit a delta model nested within a clade, with a BM fit to the rest of the tree. Default bounds from ~0 - 5.

• model="OU" fits an Ornstein-Uhlenbeck model - a random walk with a central tendency proportional to alpha. High values of alpha can be interpreted as evidence of evolutionary constraints, stabilising selection or weak phylogenetic signal. It is often difficult to distinguish among these possibilities. If a nodeIDs is supplied, the model will fit a OU model nested within a clade, with a BM fit to the rest of the tree. For OU models, alternative optimisation are performed with different starting values (1e-8, 0.01, 0.1, 1, 5). Default bounds from ~0 - 10.

• model="ACDC" fits a model to in which rates can exponentially increased or decrease through time (Blomberg et al. 2003). If the upper bound is < 0, the model is equivalent to the 'Early Burst' model of Harmon et al. 2010. If a nodeIDs is supplied, the model will fit a ACDC model nested within a clade, with a BM fit to the rest of the tree. Default rate parameter bounds from ln(1e-10) ~ ln(20) divided by the root age. Note this process starts on the stem branch leading to the MRCA of the common node, unlike the other methods that start at the common node.

• model="trend" fits a model in which the expectated mean change through time is non-zero, signifying a directional evolution to a larger or smaller trait value. This model is only appliacble to non-ultrametric trees.

• model="psi" fits a model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate (Ingram 2010). Note that the algorithm will automatically estimate speciation and extinction estimates, and will incorporate estimates of 'hidden' speciation if death estimates are greater than 0.

• model="multipsi" fits a model to assess to the relative contributions of speciation and gradual evolution to a trait's evolutionary rate but allows seperate values of psi fitted to seperate branches (Ingram 2010; Ingram et al. 2016). Note that the algorithm will automatically estimate speciation and extinction estimates, and will incorporate estimates of 'hidden' speciation if death estimates are greater than 0.

• model="free" fits Mooers et al's free model where each branch has its own rate of trait evolution. This can be a useful exploratory analysis but it is slow due to the number of parameters, particularly for large trees. Default rate parameter bounds from ~0 - 200.

• model="clade" fits a model where particular clades are a priori hypothesised to have different rates of trait evolution (see O'Meara et al. 2006; Thomas et al. 2006, 2009). Clades are specified using nodeIDs and are defined as the mrca node. Default rate parameter bounds from ~0 - 200.

• model="tm1" fits "clade" models without any a priori assertion of the location of phenotypic diversification rate shifts. It uses the same AIC approach as the runMedusa function in the geiger package (runMedusa tests for shifts in the rate of lineage diversification). The algorithm first fits a constant-rate Brownian model to the data, it then works iteratively through the phylogeny fitting a two-rate model at each node in turn. Each two-rate model is compared to the constant rate model and the best two-rate model is retained. Keeping the location of this rate shift intact, it then repeats the procedure for a three-rate model and so on. The maximum number of rate shifts can be specified a priori using nSplits. Limits can be applied to the size (species richness) of clades on which to infer new rate shifts using minCladeSize. This can be useful to enable large trees to be handled but should be used cautiously since specifiying a large minimum clade size may result in biologically interesting nested rate shifts being missed. Equally, very small clade sizes may provide poor estimates of rate that may not be informative. Limits on the search can also be placed using restrictNode. This requires a list where each element of the list is a vector of tip names that define monophyletic groups. Rate shifts will not be searched for within any of the defined groups. Default rate parameter bounds from ~0 - 1000.

• model="tm2" this model is similar to "tm1", however, at each node it assesses the fit of two models. The first model is exactly as per "tm1". The second model infers a rate shift on the single branch descending directly from a node but not on any of the descending branches thereafter. Only the best fitting single-branch or whole clade model is retained for the next iteration. If a single-branch shift is favoured, this infers either that there was a rapid shift in trait value along the stem leading to the crown group, or that the members of the clade have undergone parallel shifts. In either case, this can be considered as a change in mean, though separating a single early shift from a clade-parallel shift is not possible with this method.

• model="timeSlice" A model in which all branch rates change at a time or times set a priori by the user. If Default rate parameter bounds from ~0 - 1000. If splitTime=NULL, all 1 Ma (as defined by test Age) intervals from the root of the tree - 10 and the youngest tip + 10 will be included in the search. The +/- 10 Ma age can be modified using the argument boundaryAge. At each stage the best fitting model will be stored, and the search will continue until n shifts, with n shifts defined by nSplits. If a single value or vector is used for splitTime, only these ages are included in the search.

• model="modeslice" A model in which all branch modes change at a time or times set a priori by the user.

Value

brownianVariance Brownian variance (or covariance for multiple traits) given the data and phylogeny

logLikelihood The log-likelihood of the model and data

Author(s)

Gavin Thomas, Mark Puttick

References

Felsenstein J. 1973. Maximum-likelihood estimation of evolutionary trees from continuous characters. Am. J. Hum. Genet. 25, 471-492. Felsenstein J. 1985. Phylogenies and the comparative method. American Naturalist 125, 1-15. Freckleton RP & Jetz W. 2009. Space versus phylogeny: disentangling phylogenetic and spatial signals in comparative data. Proc. Roy. Soc. B 276, 21-30.

Ingram T. 2011. Speciation along a depth gradient in a marine adaptive radiation. Proc. Roy. Soc. B. 278, 613-618.

Ingram T, Harrison AD, Mahler L, Castaneda MdR, Glor RE, Herrel A, Stuart YE, and Losos JB. 2016. Comparative tests of the role of dewlap size in Anolis lizard speciation. Proc. Roy. Soc. B. 283, 20162199. #' Mooers AO, Vamosi S, & Schluter D. 1999. Using phylogenies to test macroevolutionary models of trait evolution: sexual selection and speciation in Cranes (Gruinae). American Naturalist 154, 249-259. O'Meara BC, Ane C, Sanderson MJ & Wainwright PC. 2006. Testing for different rates of continuous trait evolution using likelihood. Evolution 60, 922-933 Pagel M. 1997. Inferring evolutionary processes from phylogenies. Zoologica Scripta 26, 331-348. Pagel M. 1999 Inferring the historical patterns of biological evolution. Nature 401, 877-884. Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

See Also

transformPhylo.ML,transformPhylo.MCMC, transformPhylo.sim

Examples

# Data and phylogeny
data(anolis.tree)
data(anolis.data)

# anolis.data is not matrix and contains missing data so put together matrix of
# relevant traits (here female and male snout vent lengths) and remove species 
# with missing data from the matrix and phylogeny
sorted.traits <- sortTraitData(anolis.tree, anolis.data,
c("Female_SVL", "Male_SVL"), log.trait=TRUE, pass.ultrametric=TRUE)

tree <- sorted.traits$phy
traits <- sorted.traits$trait

# log likelihood of kappa = 0.1 or 1
transformPhylo.ll(traits, phy=tree, model="kappa", kappa=0.1)
transformPhylo.ll(traits, phy=tree, model="kappa", kappa=1)

# log likelihood of lambda = 0.01 or 1
transformPhylo.ll(traits, phy=tree, model="lambda", lambda=0.01)
transformPhylo.ll(traits, phy=tree, model="lambda", lambda=1)

# log likelihood of delta = 1.5 or 1
transformPhylo.ll(traits, phy=tree, model="delta", delta=1.5)
transformPhylo.ll(traits, phy=tree, model="delta", delta=1)

# log likelihood of alpha = 0.001 or 2
transformPhylo.ll(traits, phy=tree, model="OU", alpha=0.001)
transformPhylo.ll(traits, phy=tree, model="OU", alpha=2)

Phylogenetic tree transformations

Description

Simulates trait data on a tree using a specified model of evolution (see details).

Usage

transformPhylo.sim(
  phy,
  n = 1,
  x = NULL,
  model = NULL,
  returnNodes = FALSE,
  kappa = NULL,
  lambda = NULL,
  delta = NULL,
  alpha = NULL,
  psi = NULL,
  acdcRate = NULL,
  lambda.sp = NULL,
  trend = NULL,
  trend.anc.state = 0,
  nodeIDs = NULL,
  rateType = NULL,
  cladeRates = NULL,
  branchRates = NULL,
  rate = NULL,
  group.means = NULL,
  splitTime = NULL,
  timeRates = NULL,
  branchLabels = NULL,
  rate.var = NULL,
  mode.order = NULL
)

Arguments

Value

Returns a matrix of simulated dated with taxon names as rownames (number of columns=n).

Author(s)

Gavin Thomas, Mark Puttick

References

Ricklefs RE. 2006. Time, species, and the generation of trait variation in clades. Systematic Biology 55, 151-159.

Ricklefs RE. 2006. Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030

See Also

transformPhylo.ML, transformPhylo.ll, transformPhylo, transformPhylo.MCMC

Examples

data(anolis.tree)
data(anolis.data)

# Simulate 10 sets of data with kappa=0.1 using the anolis tree
sim.dat1 <- transformPhylo.sim(phy=anolis.tree, n=10, model="kappa", kappa=0.1)

# Simulate 10 sets of data where rates and means differ between to the categories defined by "x"
x <- anolis.data$geo_ecomorph
names(x) <-  rownames(anolis.data)
sim.dat2 <- transformPhylo.sim(phy=anolis.tree, n=10, x=x, model="mixedRate", rate=c(1,1,2,4),
group.means=c(0,5,0,0))

Conversion among data and phylogeny objects

Description

Transforms the expected variance and covariances among species according to hypotheses of rate variation between lineages.

Usage

transformRateMatrix(rateData, rate = NULL)

Arguments

Value

retMat Rate-transformed variance covariance matrix

Author(s)

Gavin Thomas

References

Thomas GH, Freckleton RP, & Szekely T. 2006. Comparative analyses of the influence of developmental mode on phenotypic diversification rates in shorebirds. Proceedings of the Royal Society B 273, 1619-1624.

Thomas GH, Meiri S, & Phillimore AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63, 2017-2030.

Examples

 data(anolis.tree)
 data(anolis.data)
 ## Convert data to class rateData with a rateMatrix object as input
 anolis.rateMatrix <- as.rateMatrix(phy=anolis.tree, x="geo_ecomorph", data=anolis.data)
 anolis.rateData <- as.rateData(y="Female_SVL", x="geo_ecomorph", 
 rateMatrix = anolis.rateMatrix, phy=NULL, data=anolis.data, log.y=TRUE)
 
 # Tranform the expected variance covariance matrix so that the rates in the first and last 
 # categories are equal (both 1) whereas the rate in the second category is twice as fast (2) and 
 # the rate in the third category is ten times slower.
 
 trans.anolis.rateData <- transformRateMatrix(rateData=anolis.rateData, rate = c(1,2,0.1,1))