Package 'coevolve'

Title: Fit Bayesian Generalized Dynamic Phylogenetic Models using 'Stan' or 'JAX'
Description: Fit Bayesian generalized dynamic phylogenetic models using 'Stan' or 'JAX'. The package allows an arbitrary number of variables with different response distributions to coevolve on a phylogenetic tree via an Ornstein-Uhlenbeck evolutionary process, allowing users to assess directionality between coevolving variables. The optional JAX backend enables high-performance sampling with NumPyro and potential GPU acceleration via 'reticulate'.
Authors: Scott Claessens [aut, cre] (ORCID: <https://orcid.org/0000-0002-3562-6981>), Erik Ringen [aut] (ORCID: <https://orcid.org/0000-0002-3565-6961>), Rob Moss [rev] (ORCID: <https://orcid.org/0000-0002-4568-2012>), Krzysztof Bartoszek [rev] (ORCID: <https://orcid.org/0000-0002-5816-4345>)
Maintainer: Scott Claessens <[email protected]>
License: GPL (>= 3)
Version: 1.2.1
Built: 2026-10-07 08:29:37 UTC
Source: https://github.com/ScottClaessens/coevolve

Help Index

Political and religious authority in Austronesian societies

Description

This dataset and associated phylogeny describes the states of political and religious authority in 97 Austronesian societies. These data were originally compiled in Sheehan et al. (2023). Political authority is defined as a right to manage interactions between living human beings, whereas religious authority is defined as a right to manage interactions between living human beings and supernatural agents or powers. Authority is coded as being absent, sublocal (smaller than the local community), local (coextensive with the local community), or supralocal (consisting of more than one local community).

Usage

authority

Format

A list containing a dataset, an associated language phylogeny, and associated geographic coordinates. The phylogeny is a pruned maximum clade credibility tree for 97 Austronesian languages. The coordinates are in a data frame, with longitude and latitude values measured in decimal degrees. The dataset is a data frame with 97 observations and three variables:

language

The name of the Austronesian language linked to each society

political_authority

An ordered factor indicating whether political authority is absent, sublocal, local, or supralocal in each society

religious_authority

An ordered factor indicating whether religious authority is absent, sublocal, local, or supralocal in each society

Source

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45.

Examples

## Not run: 
# fit model to authority data
m <-
  coev_fit(
    data = authority$data,
    variables = list(
      political_authority = "ordered_logistic",
      religious_authority = "ordered_logistic"
    ),
    id = "language",
    tree = authority$phylogeny,
    # arguments for cmdstanr::sample()
    parallel_chains = 4,
    seed = 1,
    # set prior manually
    prior = list(A_offdiag = "normal(0, 2)")
  )
# print model summary
summary(m)
# plot delta theta
coev_plot_delta_theta(m)

## End(Not run)

Extract ancestral state estimates from a fitted coevolutionary model

Description

Extracts posterior estimates of trait values at internal nodes of the phylogenetic tree from a fitted coevfit object. By default, estimates are returned on the latent (eta) scale; use scale = "response" to apply the inverse link for each variable's response distribution. The latent eta parameters in the model correspond directly to nodes in the phylogeny, making ancestral state reconstruction a natural byproduct of model fitting.

Usage

coev_ancestral_states(
  object,
  variables = NULL,
  nodes = "internal",
  tree_id = NULL,
  scale = "latent",
  summary = TRUE,
  prob = 0.95
)

Arguments

object

An object of class coevfit
variables

(optional) A character vector of variable names to include. If NULL (default), all variables from the fitted model are included.
nodes

Which nodes to include. Either "internal" (default) for internal nodes only, "all" for both tips and internal nodes, or an integer vector of specific node IDs (using ape's node numbering: tips = 1:N_tips, internal = (N_tips+1):total).
tree_id

(optional) An integer indicating which tree to use when the model was fit with a multiPhylo object. If NULL (default), uses tree 1. (Integration over trees with topological uncertainty is not yet implemented.)
scale

Character string, either "latent" (default) to return values on the latent eta scale, or "response" to apply the inverse link function for each variable's response distribution.
summary

Logical. If TRUE (default), returns a long-format data frame with posterior point estimates and credible intervals. If FALSE, returns the raw posterior draws.
prob

A single numeric value between 0 and 1 indicating the desired probability mass for the credible interval. Default is 0.95.

Value

If summary = TRUE, a long-format data frame (tibble) with columns:

node

Integer node ID in the reference tree

variable

Character variable name

category

Only present when scale = "response". For ordinal variables, the category label ("cat_1", "cat_2", ...). NA for non-ordinal variables.

estimate

Posterior median

lower

Lower bound of the credible interval

upper

Upper bound of the credible interval

The data frame has attributes ref_tree (the phylo object used), prob, and scale.

If summary = FALSE, a list with elements:

draws

A 3D array (draws x nodes x variables)

ref_tree

The reference phylo object

node_ids

Integer vector of node IDs (matching dim 2 of draws)

variable_names

Character vector of variable names

See Also

coev_fit, coev_plot_trait_values

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  chains = 4,
  parallel_chains = 4,
  seed = 1
)

# extract ancestral states (latent scale, internal nodes, summary)
asr <- coev_ancestral_states(fit)

# extract all nodes including tips
asr_all <- coev_ancestral_states(fit, nodes = "all")

# extract on response scale
asr_resp <- coev_ancestral_states(fit, scale = "response")

# get raw posterior draws
asr_draws <- coev_ancestral_states(fit, summary = FALSE)

## End(Not run)

Calculate delta theta from a fitted coevfit object

Description

Calculate Δθ\Delta\theta from a fitted coevfit object. Δθ\Delta\theta is defined as the difference in the equilibrium value for a "response" trait that results from a standardised increase in the value of another "predictor" trait. This value can be used to assess contingencies and directionality between variables in the coevolutionary process.

Usage

coev_calculate_delta_theta(object, response, predictor)

Arguments

object

An object of class coevfit
response

A character string of length one equal to one of the coevolving variables in the model (strictly case-sensitive)
predictor

A character string of length one equal to one of the coevolving variables in the model (strictly case-sensitive)

Details

This function calculates Δθ\Delta\theta, which is defined as the difference in the equilibrium value for a "response" trait that results from a standardised increase in the value of another "predictor" trait. The function first calculates the equilibrium trait value for the response trait when the predictor trait is held at its empirical median value (the coev_calculate_theta function is used for this purpose, see help(coev_calculate_theta) for further details). The function then calculates the equilibrium trait value for the response variable after increasing the predictor trait from its median by one median absolute deviation. The function then returns the posterior difference between these values. The resulting Δθ\Delta\theta samples can be used to infer whether increases in one trait have a positive or negative selective effect on another trait in the model.

Value

Posterior samples in the draws_array format

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_calculate_delta_theta, coev_plot_delta_theta

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# calculate delta theta
coev_calculate_delta_theta(
  object = fit,
  response = "political_authority",
  predictor = "religious_authority"
  )

## End(Not run)

Calculate equilibrium trait values (theta) for a fitted coevfit object

Description

Calculate equilibrium trait values θ\theta for one or more traits given a set of intervention values for other traits from a fitted coevfit object.

Usage

coev_calculate_theta(object, intervention_values = NULL)

Arguments

object

An object of class coevfit
intervention_values

Either NULL (the default) or a named list of variables and associated intervention values for calculating equilibrium trait values. If NULL, calculates the equilibrium states when all parameters are free to vary. Otherwise, all coevolving variables must be declared separately in the named list without repetition. If the intervention value for a particular variable is set to NA, this variable is treated as a free variable. If the intervention value for a particular variable is specified, the variable is held constant at this trait value in the calculation. At least one variable must be declared as a free variable and at least one variable must be held constant (e.g., list(var1 = NA, var2 = 0)).

Details

The equilibrium trait values for freely evolving traits η\mathbf{\eta} are calculated using the following formula:

θ=−A−1b\mathbf{\theta} = -\mathbf{A}^{-1}\mathbf{b}

If we hold some variables constant at some value(s) (denoted ηh\eta_{h}) and let others evolve freely (denoted ηf\eta_{f}), we can partition the parameters as follows:

  • Aff\mathbf{A}_{ff}: selection coefficients to/from the free variables

  • Afh\mathbf{A}_{fh}: selection coefficients to the free variables from the held variables

  • bf\mathbf{b}_{f}: continuous time intercepts for the free variables

We can then calculate the equilibrium values for the free variables:

θf=−Aff−1(Afhηh+bf)\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} \mathbf{\eta}_h + \mathbf{b}_f \right)

With the overall equilibrium vector being a mix of the free equilibria and the held values:

θ∣ηh=[θfηh]\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}

Value

Posterior samples in matrix format

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_calculate_delta_theta, coev_plot_delta_theta

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# calculate theta with no interventions
coev_calculate_theta(
  object = fit
 )

# calculate theta given intervention values
coev_calculate_theta(
  object = fit,
  intervention_values = list(
    political_authority = NA,
    religious_authority = 0
    )
  )

## End(Not run)

Fit Bayesian dynamic coevolutionary model in Stan

Description

Fit Bayesian dynamic coevolutionary models in Stan for full Bayesian inference. The model allows users to assess causal directionality (X -> Y vs. Y -> X) and contingencies (X, then Y) in evolution. Several data types are supported, including binary, ordinal, count, continuous, and positive real variables. The model can additionally account for missing data, repeated observations, and controls for spatial proximity. Model fit can be assessed and compared with posterior predictive checks and cross-validation.

Usage

coev_fit(
  data,
  variables,
  id,
  tree,
  effects_mat = NULL,
  complete_cases = FALSE,
  lon_lat = NULL,
  dist_k = NA,
  dist_cov = "exp_quad",
  measurement_error = NULL,
  prior = NULL,
  scale = TRUE,
  estimate_correlated_drift = TRUE,
  estimate_residual = TRUE,
  log_lik = FALSE,
  prior_only = FALSE,
  adapt_delta = 0.95,
  nuts_sampler = "stan",
  backend = deprecated(),
  dist_mat = deprecated(),
  ...
)

Arguments

data

An object of class data.frame (or one that can be coerced to that class) containing data of all variables used in the model.
variables

A named list identifying variables that should coevolve in the model and their associated response distributions as character strings (e.g. list(var1 = "bernoulli_logit", var2 = "ordered_logistic")). Variable names must refer to valid column names in data. Currently, the only supported response distributions are bernoulli_logit, ordered_logistic, poisson_softplus, negative_binomial_softplus, normal, and gamma_log. Bernoulli variables must be 0/1 integers, ordered variables must be ordered factors, Poisson and negative binomial variables must be positive integers, normal variables must be continuous numeric, and gamma variables must be positive numeric.
id

A character of length one identifying the variable in the data that links rows to tips on the phylogeny (strictly case-sensitive). Must refer to a valid column name in the data. The id column must exactly match the tip labels in the phylogeny.
tree

A phylogenetic tree object of class phylo or multiPhylo. The tree(s) must be rooted and must include positive non-zero branch lengths. All trees in multiPhylo objects must have the same number of internal nodes and branches.
effects_mat

(optional) A boolean matrix with row and column names exactly matching the variables declared for the model. If not specified, all cross-lagged effects will be estimated in the model. If specified, the model will only estimate cross-lagged effects where cells in the matrix are TRUE and will ignore cross-lagged effects where cells in the matrix are FALSE. In the matrix, columns represent predictor variables and rows represent outcome variables. All autoregressive effects (e.g., X -> X) must be TRUE in the matrix.
complete_cases

(optional) Logical. If FALSE (default), all missing values are imputed by the model. If TRUE, taxa with missing data are excluded.
lon_lat

(optional) A data.frame containing longitude and latitude values for all tips on the phylogeny. The data frame must contain the following columns: a column labelled "id" that exactly matches the tip labels on the phylogeny; a column labelled "longitude" that declares longitude values in decimal degrees; and a column labelled "latitude" that declares latitude values in decimal degrees. If specified, the model will additionally control for spatial location by including a separate Gaussian Process over locations for every coevolving variable in the model.
dist_k

An integer of length one specifying the number of basis functions for computing Hilbert-space approximate Gaussian Processes over locations. If NA (the default), exact Gaussian Processes are computed.
dist_cov

A string of length one specifying the covariance kernel used for Gaussian Processes over locations (strictly case-sensitive). Currently supported are "exp_quad" (exponentiated-quadratic kernel; default), "exponential" (exponential kernel), and "matern32" (Matern 3/2 kernel).
measurement_error

(optional) A named list of coevolving variables and their associated columns in the dataset containing standard errors. Only valid for normally-distributed variables. For example, if we declare variables = list(x = "normal", y = "normal"), then we could set measurement_error = list(x = "x_std_err") to tell the function to include measurement error on x using standard errors from the x_std_err column of the dataset.
prior

(optional) A named list of priors for the model. If not specified, the model uses default priors (see help(coev_fit)). Alternatively, the user can specify a named list of priors. The list must contain non-duplicated entries for any of the following parameters: the autoregressive effects (A_diag), the cross effects (A_offdiag), the Cholesky factor for the drift matrix (L_R), the drift std. dev. parameters (Q_sigma), the continuous time intercepts (b), the ancestral states for the traits (eta_anc), the cutpoints for ordinal variables (c), the overdispersion parameters for negative binomial variables (phi), the shape parameters for gamma variables (shape), the sigma parameters for Gaussian Processes over locations (sigma_dist), the rho parameters for Gaussian Processes over locations (rho_dist), the residual standard deviations when there are repeated observations (sigma_residual), and the Cholesky factor for the residual correlations when there are repeated observations (L_residual). These must be entered with valid prior strings, e.g. list(A_offdiag = "normal(0, 2)"). Invalid prior strings will throw an error when the function internally checks the syntax of resulting Stan code.
scale

Logical. If TRUE (default), variables following the normal and gamma_log response distributions are scaled before fitting the model. Continuous variables following the normal distribution are standardised (e.g., mean centered and divided by their standard deviation) and positive real variables following the gamma_log distribution are divided by the mean value without centering. This approach is recommended when using default priors to improve efficiency and ensure accurate inferences. If FALSE, variables are left unscaled for model fitting. In this case, users should take care to set sensible priors on variables.
estimate_correlated_drift

Logical. If TRUE (default), the model estimates the off-diagonals for the QQ drift matrix (i.e., correlated drift). If FALSE, the off-diagonals for the QQ drift matrix are set to zero.
estimate_residual

Logical. If TRUE (default), the model estimates residual standard deviations and residual correlations when there are repeated observations for taxa. If FALSE, residual standard deviations and residual correlations are not estimated. The latter may be preferable in cases where repeated observations are sparse (e.g., only some taxa have only few repeated observations). This argument only applies when repeated observations are present in the data.
log_lik

Logical. Set to FALSE by default. If TRUE, the model returns the pointwise log likelihood, which can be used to calculate WAIC and LOO.
prior_only

Logical. If FALSE (default), the model is fitted to the data and returns a posterior distribution. If TRUE, the model samples from the prior only, ignoring the likelihood.
adapt_delta

Target acceptance probability for the NUTS sampler (default 0.95). Passed to cmdstanr::sample() as adapt_delta, or to nutpie as target_accept.
nuts_sampler

Which NUTS implementation to use. "stan" (default) fits the Stan model with CmdStan via cmdstanr. "nutpie" fits a pure JAX model using nutpie's Rust NUTS sampler with JAX gradients via jax.value_and_grad.
backend

Deprecated. Use the nuts_sampler argument instead.
dist_mat

Deprecated. Use the lon_lat argument instead.
...

Additional arguments. For nuts_sampler = "stan": passed to cmdstanr::sample() (and compilation via extra_stanc_args, extra_compile_args). For nuts_sampler = "nutpie": sampling arguments such as chains, iter_sampling, iter_warmup, seed.

Details

The coev_fit function generates the Stan code for the model, generates the data list, and then compiles and fits the model using the cmdstanr package.

Fit a Bayesian dynamic coevolutionary model in Stan. A general overview is provided in the vignette vignette("coevolve")

Available response distributions

Variables that should coevolve in the model are declared in the variables argument, using a named list with associated response distributions. For example: list(x = "bernoulli_logit", y = "ordered_logistic"). Currently, the only supported response distributions are bernoulli_logit, ordered_logistic, poisson_softplus, negative_binomial_softplus, normal, and gamma_log.

Default prior distributions

If priors are not explicitly declared by the user, default priors are used. Default priors were chosen to be weakly regularising, which improves model fitting and conservatism in parameter estimates. Default priors for Stan parameters are as follows:

  • A_diag (autoregressive effects) = std_normal()

  • A_offdiag (cross effects) = std_normal()

  • L_R (Cholesky factor for drift matrix) = lkj_corr_cholesky(4)

  • Q_sigma (drift std. dev. parameters) = std_normal()

  • b (continuous time intercepts) = std_normal()

  • eta_anc (trait ancestral states) = std_normal()

  • c (ordinal cutpoints) = normal(0, 2)

  • shape (shape parameters for gamma distributions) = gamma(0.01, 0.01)

  • sigma_dist (sigma for Gaussian process over locations) = exponential(1)

  • rho_dist (rho for Gaussian process over locations) = exponential(5)

  • sigma_residual (residual standard deviations) = exponential(1)

  • L_residual (Cholesky factor for residual correlations) = lkj_corr_cholesky(4)

The default prior for phi (the overdispersion parameter for the negative-binomial distribution) is scaled automatically based on the variance of the data.

We recommend that users assess the suitability of these default priors by fitting the model with prior_only = TRUE and then plotting prior predictive checks for all variables using the coev_plot_predictive_check() function.

Handling missing data

In order to retain the most information, the coev_fit() function automatically imputes all missing values. To turn off this behaviour and exclude taxa with missing data, set complete_cases = FALSE.

Dealing with repeated observations

If taxa appear in the dataset multiple times (i.e., there are repeated observations), the model will automatically estimate residual standard deviations and residual correlations that capture the within-taxa variation that is not due to the coevolutionary process. To turn off this behaviour, set estimate_residual = FALSE.

We do not recommend modelling residual standard deviations in cases where repeated observations are present for some traits but not others, as there will be no within-species variance for traits without repeated observations.

Note that models with repeated observations for mixtures of continuous and non-continuous traits are currently considered experimental and should be treated with caution.

Incorporating measurement error

If any normally-distributed coevolving variables are measured with error, the user can pass standard errors for one or more of these variables into the model with the measurement_error argument.

Controlling for spatial location

If users declare longitude and latitude values with the lon_lat argument, the model will include several Gaussian processes (one per coevolving variable) over spatial locations.

Value

Fitted model of class coevfit

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_make_stancode, coev_make_standata

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# print model summary
summary(fit)

## End(Not run)

Make Stan code for dynamic coevolutionary model

Description

Make the Stan code for the Bayesian dynamic coevolutionary model. Stan code is generated, checked for syntactical errors, and then returned as a character string.

Usage

coev_make_stancode(
  data,
  variables,
  id,
  tree,
  effects_mat = NULL,
  complete_cases = FALSE,
  lon_lat = NULL,
  dist_k = NA,
  dist_cov = "exp_quad",
  measurement_error = NULL,
  prior = NULL,
  scale = TRUE,
  estimate_correlated_drift = TRUE,
  estimate_residual = TRUE,
  log_lik = FALSE,
  prior_only = FALSE,
  dist_mat = deprecated()
)

Arguments

data

An object of class data.frame (or one that can be coerced to that class) containing data of all variables used in the model.
variables

A named list identifying variables that should coevolve in the model and their associated response distributions as character strings (e.g. list(var1 = "bernoulli_logit", var2 = "ordered_logistic")). Variable names must refer to valid column names in data. Currently, the only supported response distributions are bernoulli_logit, ordered_logistic, poisson_softplus, negative_binomial_softplus, normal, and gamma_log. Bernoulli variables must be 0/1 integers, ordered variables must be ordered factors, Poisson and negative binomial variables must be positive integers, normal variables must be continuous numeric, and gamma variables must be positive numeric.
id

A character of length one identifying the variable in the data that links rows to tips on the phylogeny (strictly case-sensitive). Must refer to a valid column name in the data. The id column must exactly match the tip labels in the phylogeny.
tree

A phylogenetic tree object of class phylo or multiPhylo. The tree(s) must be rooted and must include positive non-zero branch lengths. All trees in multiPhylo objects must have the same number of internal nodes and branches.
effects_mat

(optional) A boolean matrix with row and column names exactly matching the variables declared for the model. If not specified, all cross-lagged effects will be estimated in the model. If specified, the model will only estimate cross-lagged effects where cells in the matrix are TRUE and will ignore cross-lagged effects where cells in the matrix are FALSE. In the matrix, columns represent predictor variables and rows represent outcome variables. All autoregressive effects (e.g., X -> X) must be TRUE in the matrix.
complete_cases

(optional) Logical. If FALSE (default), all missing values are imputed by the model. If TRUE, taxa with missing data are excluded.
lon_lat

(optional) A data.frame containing longitude and latitude values for all tips on the phylogeny. The data frame must contain the following columns: a column labelled "id" that exactly matches the tip labels on the phylogeny; a column labelled "longitude" that declares longitude values in decimal degrees; and a column labelled "latitude" that declares latitude values in decimal degrees. If specified, the model will additionally control for spatial location by including a separate Gaussian Process over locations for every coevolving variable in the model.
dist_k

An integer of length one specifying the number of basis functions for computing Hilbert-space approximate Gaussian Processes over locations. If NA (the default), exact Gaussian Processes are computed.
dist_cov

A string of length one specifying the covariance kernel used for Gaussian Processes over locations (strictly case-sensitive). Currently supported are "exp_quad" (exponentiated-quadratic kernel; default), "exponential" (exponential kernel), and "matern32" (Matern 3/2 kernel).
measurement_error

(optional) A named list of coevolving variables and their associated columns in the dataset containing standard errors. Only valid for normally-distributed variables. For example, if we declare variables = list(x = "normal", y = "normal"), then we could set measurement_error = list(x = "x_std_err") to tell the function to include measurement error on x using standard errors from the x_std_err column of the dataset.
prior

(optional) A named list of priors for the model. If not specified, the model uses default priors (see help(coev_fit)). Alternatively, the user can specify a named list of priors. The list must contain non-duplicated entries for any of the following parameters: the autoregressive effects (A_diag), the cross effects (A_offdiag), the Cholesky factor for the drift matrix (L_R), the drift std. dev. parameters (Q_sigma), the continuous time intercepts (b), the ancestral states for the traits (eta_anc), the cutpoints for ordinal variables (c), the overdispersion parameters for negative binomial variables (phi), the shape parameters for gamma variables (shape), the sigma parameters for Gaussian Processes over locations (sigma_dist), the rho parameters for Gaussian Processes over locations (rho_dist), the residual standard deviations when there are repeated observations (sigma_residual), and the Cholesky factor for the residual correlations when there are repeated observations (L_residual). These must be entered with valid prior strings, e.g. list(A_offdiag = "normal(0, 2)"). Invalid prior strings will throw an error when the function internally checks the syntax of resulting Stan code.
scale

Logical. If TRUE (default), variables following the normal and gamma_log response distributions are scaled before fitting the model. Continuous variables following the normal distribution are standardised (e.g., mean centered and divided by their standard deviation) and positive real variables following the gamma_log distribution are divided by the mean value without centering. This approach is recommended when using default priors to improve efficiency and ensure accurate inferences. If FALSE, variables are left unscaled for model fitting. In this case, users should take care to set sensible priors on variables.
estimate_correlated_drift

Logical. If TRUE (default), the model estimates the off-diagonals for the

QQ

drift matrix (i.e., correlated drift). If FALSE, the off-diagonals for the

QQ

drift matrix are set to zero.
estimate_residual

Logical. If TRUE (default), the model estimates residual standard deviations and residual correlations when there are repeated observations for taxa. If FALSE, residual standard deviations and residual correlations are not estimated. The latter may be preferable in cases where repeated observations are sparse (e.g., only some taxa have only few repeated observations). This argument only applies when repeated observations are present in the data.
log_lik

Logical. Set to FALSE by default. If TRUE, the model returns the pointwise log likelihood, which can be used to calculate WAIC and LOO.
prior_only

Logical. If FALSE (default), the model is fitted to the data and returns a posterior distribution. If TRUE, the model samples from the prior only, ignoring the likelihood.
dist_mat

[Deprecated] Use the lon_lat argument instead.

Details

For further details, see help(coev_fit)

Value

A character string containing the Stan code to fit the dynamic coevolutionary model

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_make_standata, coev_fit

Examples

# make stan code
stan_code <- coev_make_stancode(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny
)

# include effects matrix
effects_mat <-
  matrix(
    c(TRUE, TRUE,
      FALSE, TRUE),
    nrow = 2,
    dimnames = list(
      c("political_authority", "religious_authority"),
      c("political_authority", "religious_authority")
    )
  )
stan_code <- coev_make_stancode(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  effects_mat = effects_mat
)

# include longitude/latitude values
stan_code <- coev_make_stancode(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  lon_lat = authority$coordinates
)

# include measurement error
d <- authority$data
d$x <- rnorm(nrow(d))
d$y <- rnorm(nrow(d))
d$x_std_err <- rexp(nrow(d))
d$y_std_err <- rexp(nrow(d))
stan_code <- coev_make_stancode(
  data = d,
  variables = list(
    x = "normal",
    y = "normal"
  ),
  id = "language",
  tree = authority$phylogeny,
  measurement_error = list(
    x = "x_std_err",
    y = "y_std_err"
  )
)

# set manual priors
stan_code <- coev_make_stancode(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  prior = list(A_offdiag = "normal(0, 2)")
)

Make Stan data for dynamic coevolutionary model

Description

Make the data list for the Stan model. The function takes a dataset, phylogeny, and a set of variables and performs the necessary computations (e.g., cutting up the tree into segments, computing branch lengths, determining parent and child nodes) for the Stan model.

Usage

coev_make_standata(
  data,
  variables,
  id,
  tree,
  effects_mat = NULL,
  complete_cases = FALSE,
  lon_lat = NULL,
  dist_k = NA,
  dist_cov = "exp_quad",
  measurement_error = NULL,
  prior = NULL,
  scale = TRUE,
  estimate_correlated_drift = TRUE,
  estimate_residual = TRUE,
  log_lik = FALSE,
  prior_only = FALSE,
  dist_mat = deprecated()
)

Arguments

data

An object of class data.frame (or one that can be coerced to that class) containing data of all variables used in the model.
variables

A named list identifying variables that should coevolve in the model and their associated response distributions as character strings (e.g. list(var1 = "bernoulli_logit", var2 = "ordered_logistic")). Variable names must refer to valid column names in data. Currently, the only supported response distributions are bernoulli_logit, ordered_logistic, poisson_softplus, negative_binomial_softplus, normal, and gamma_log. Bernoulli variables must be 0/1 integers, ordered variables must be ordered factors, Poisson and negative binomial variables must be positive integers, normal variables must be continuous numeric, and gamma variables must be positive numeric.
id

A character of length one identifying the variable in the data that links rows to tips on the phylogeny (strictly case-sensitive). Must refer to a valid column name in the data. The id column must exactly match the tip labels in the phylogeny.
tree

A phylogenetic tree object of class phylo or multiPhylo. The tree(s) must be rooted and must include positive non-zero branch lengths. All trees in multiPhylo objects must have the same number of internal nodes and branches.
effects_mat

(optional) A boolean matrix with row and column names exactly matching the variables declared for the model. If not specified, all cross-lagged effects will be estimated in the model. If specified, the model will only estimate cross-lagged effects where cells in the matrix are TRUE and will ignore cross-lagged effects where cells in the matrix are FALSE. In the matrix, columns represent predictor variables and rows represent outcome variables. All autoregressive effects (e.g., X -> X) must be TRUE in the matrix.
complete_cases

(optional) Logical. If FALSE (default), all missing values are imputed by the model. If TRUE, taxa with missing data are excluded.
lon_lat

(optional) A data.frame containing longitude and latitude values for all tips on the phylogeny. The data frame must contain the following columns: a column labelled "id" that exactly matches the tip labels on the phylogeny; a column labelled "longitude" that declares longitude values in decimal degrees; and a column labelled "latitude" that declares latitude values in decimal degrees. If specified, the model will additionally control for spatial location by including a separate Gaussian Process over locations for every coevolving variable in the model.
dist_k

An integer of length one specifying the number of basis functions for computing Hilbert-space approximate Gaussian Processes over locations. If NA (the default), exact Gaussian Processes are computed.
dist_cov

A string of length one specifying the covariance kernel used for Gaussian Processes over locations (case-sensitive). Currently supported are "exp_quad" (exponentiated-quadratic kernel; default), "exponential" (exponential kernel), and "matern32" (Matern 3/2 kernel).
measurement_error

(optional) A named list of coevolving variables and their associated columns in the dataset containing standard errors. Only valid for normally-distributed variables. For example, if we declare variables = list(x = "normal", y = "normal"), then we could set measurement_error = list(x = "x_std_err") to tell the function to include measurement error on x using standard errors from the x_std_err column of the dataset.
prior

(optional) A named list of priors for the model. If not specified, the model uses default priors (see help(coev_fit)). Alternatively, the user can specify a named list of priors. The list must contain non-duplicated entries for any of the following parameters: the autoregressive effects (A_diag), the cross effects (A_offdiag), the Cholesky factor for the drift matrix (L_R), the drift std. dev. parameters (Q_sigma), the continuous time intercepts (b), the ancestral states for the traits (eta_anc), the cutpoints for ordinal variables (c), the overdispersion parameters for negative binomial variables (phi), the shape parameters for gamma variables (shape), the sigma parameters for Gaussian Processes over locations (sigma_dist), the rho parameters for Gaussian Processes over locations (rho_dist), the residual standard deviations when there are repeated observations (sigma_residual), and the Cholesky factor for the residual correlations when there are repeated observations (L_residual). These must be entered with valid prior strings, e.g. list(A_offdiag = "normal(0, 2)").
scale

Logical. If TRUE (default), variables following the normal and gamma_log response distributions are scaled before fitting the model. Continuous variables following the normal distribution are standardised (e.g., mean centered and divided by their standard deviation) and positive real variables following the gamma_log distribution are divided by the mean value without centering. This approach is recommended when using default priors to improve efficiency and ensure accurate inferences. If FALSE, variables are left unscaled for model fitting. In this case, users should take care to set sensible priors on variables.
estimate_correlated_drift

Logical. If TRUE (default), the model estimates the off-diagonals for the

QQ

drift matrix (i.e., correlated drift). If FALSE, the off-diagonals for the

QQ

drift matrix are set to zero.
estimate_residual

Logical. If TRUE (default), the model estimates residual standard deviations and residual correlations when there are repeated observations for taxa. If FALSE, residual standard deviations and residual correlations are not estimated. The latter may be preferable in cases where repeated observations are sparse (e.g., only some taxa have only few repeated observations). This argument only applies when repeated observations are present in the data.
log_lik

Logical. Set to FALSE by default. If TRUE, the model returns the pointwise log likelihood, which can be used to calculate WAIC and LOO.
prior_only

Logical. If FALSE (default), the model is fitted to the data and returns a posterior distribution. If TRUE, the model samples from the prior only, ignoring the likelihood.
dist_mat

[Deprecated] Use the lon_lat argument instead.

Details

For further details, see help(coev_fit)

Value

A list containing the data for fitting the dynamic coevolutionary model in Stan

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_make_stancode, coev_fit

Examples

# make stan data
stan_data <- coev_make_standata(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny
)

# include effects matrix
effects_mat <-
  matrix(
    c(TRUE, TRUE,
      FALSE, TRUE),
    nrow = 2,
    dimnames = list(
      c("political_authority", "religious_authority"),
      c("political_authority", "religious_authority")
    )
  )
stan_data <- coev_make_standata(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  effects_mat = effects_mat
)

# include longitude/latitude values
stan_data <- coev_make_standata(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  lon_lat = authority$coordinates
)

# include measurement error
d <- authority$data
d$x <- rnorm(nrow(d))
d$y <- rnorm(nrow(d))
d$x_std_err <- rexp(nrow(d))
d$y_std_err <- rexp(nrow(d))
stan_data <- coev_make_standata(
  data = d,
  variables = list(
    x = "normal",
    y = "normal"
  ),
  id = "language",
  tree = authority$phylogeny,
  measurement_error = list(
    x = "x_std_err",
    y = "y_std_err"
  )
)

# set manual priors
stan_data <- coev_make_standata(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  prior = list(A_offdiag = "normal(0, 2)")
)

Plot delta theta values from a fitted coevfit object

Description

Plot delta theta values for all trait pairs from a fitted coevfit object. This plot can be used to visually assess contingencies and directionality between different variables in the coevolutionary process.

Usage

coev_plot_delta_theta(
  object,
  variables = NULL,
  prob = 0.66,
  prob_outer = 0.95,
  limits = NULL
)

Arguments

object

An object of class coevfit
variables

If NULL (default), the plot includes all coevolving variables from the model. Otherwise, a character vector of length >= 2 declaring the variables to be included in the plot.
prob

Probability mass to include in the inner interval. Default is 0.66 (66% interval).
prob_outer

Probability mass to include in the outer interval. Default is 0.95 (95% interval).
limits

If NULL (default), limits are scaled automatically to include all posterior samples. Otherwise, a numeric vector of length 2 specifying the lower and upper limits for the x-axis.

Details

This function repeatedly uses the coev_calculate_delta_theta function under the hood to generate a pairs plot of Δθ\Delta\theta for all variables in the model. For more details on the definition and calculation of Δθ\Delta\theta, see help(coev_calculate_delta_theta). Note that often the posterior distribution for Δθ\Delta\theta has long tails, meaning that the distribution for different traits can be difficult to visualise in a single pairs plot. If this plot does not produce satisfactory visualisations, the user should either specify narrower limits for the x-axis or use the coev_calculate_delta_theta function to create plots manually.

Value

A ggplot object

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_calculate_delta_theta, coev_calculate_theta

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# plot delta theta values for all effects
coev_plot_delta_theta(fit)

## End(Not run)

Plot flowfield of expected evolutionary change from a fitted coevfit object

Description

Plot a flowfield of the expected evolutionary dynamics from a fitted coevfit object, with help from the phaseR package.

Usage

coev_plot_flowfield(
  object,
  var1,
  var2,
  nullclines = FALSE,
  limits = c(-2.5, 2.5)
)

Arguments

object

An object of class coevfit
var1

A character string of length one equal to one of the coevolving variables in the model (strictly case-sensitive)
var2

A character string of length one equal to one of the coevolving variables in the model (strictly case-sensitive)
nullclines

Logical (defaults to FALSE); whether to show coloured nullclines to indicate where each variable is at equilibrium, depending on the state of the other
limits

A numeric vector of length 2 (defaults to c(-2.5, 2.5)); specifying the lower limit and the upper limit of the x and y axes.

Details

The flowfield visualises the change in two coevolving traits depending on their current states, with the direction and strength of change depicting with the direction and size of arrows. If nullclines are included, they represent the parameter combinations where each trait is at equilibrium, depending on the state of the other trait.

If three or more traits were included in the model, other traits are held at their median values during the computations. Note that flowfield plots can potentially produce misleading pictures of coevolutionary dynamics when other traits are held constant in models with three or more traits.

If the plot does not look right, the user might try zooming out from the default parameter space by setting wider limits. For some variables (e.g., continuous and count variables), the default limits may not be suitable.

Value

A flowfield plot drawn directly to the device

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_plot_selection_gradient

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# plot flow field
coev_plot_flowfield(
  object = fit,
  var1 = "political_authority",
  var2 = "religious_authority"
  )

## End(Not run)

Plot a predicted coevolutionary time series from a fitted coevfit object

Description

This function plots a predicted coevolutionary time series using the estimated parameters from a fitted coevfit model. By default, the plot uses the posterior ancestral states estimated by the model as the starting values, but users can also set their own starting values for traits. Plots can be generated with or without stochastic drift. For more details on the underlying predictive function, see help(coev_pred_series).

Usage

coev_plot_pred_series(object, prob = 0.95, ...)

Arguments

object

An object of class coevfit
prob

A value between 0 and 1 indicating the desired probability to be covered by the uncertainty intervals. The default is 0.95.
...

Additional arguments passed to coev_pred_series

Value

A ggplot object

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

See Also

coev_pred_series

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# simulated trait co-evolution
coev_plot_pred_series(
  object = fit,
  stochastic = TRUE
  )

# expected trait co-evolution, no drift
coev_plot_pred_series(
  object = fit,
  stochastic = FALSE,
  eta_anc = list(political_authority = -2, religious_authority = 1.5)
  )

## End(Not run)

Plot predictive check from a fitted coevfit object

Description

Plot posterior (or prior) predictive checks for variables from a fitted coevfit object, with help from the bayesplot package.

Usage

coev_plot_predictive_check(
  object,
  variables = NULL,
  ndraws = NULL,
  tree_id = NULL
)

Arguments

object

An object of class coevfit
variables

If NULL (default), the function returns a list of plots for all coevolving variables from the model. Otherwise, a character vector declaring the variables to be included in the list of plots.
ndraws

An integer indicating the number of draws to return. The default and maximum number of draws is the size of the posterior sample.
tree_id

An integer indicating the tree ID to use when making posterior predictions. Set to NULL by default, which will use draws from every tree, integrating phylogenetic uncertainty.

Value

A list of ggplot objects

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# plot predictive checks for all variables
coev_plot_predictive_check(fit)

## End(Not run)

Plot selection gradient heatmap from a fitted coevfit object

Description

Plot a heatmap of the selection gradient for two variables from a fitted coevfit object.

Usage

coev_plot_selection_gradient(
  object,
  var1,
  var2,
  contour = FALSE,
  limits = c(-2.5, 2.5)
)

Arguments

object

An object of class coevfit
var1

A character string of length one equal to one of the coevolving variables in the model (strictly case-sensitive)
var2

A character string of length one equal to one of the coevolving variables in the model (strictly case-sensitive)
contour

Logical (defaults to FALSE); whether to show white contour lines to indicate where selection is stronger than drift
limits

A numeric vector of length 2 (defaults to c(-2.5, 2.5)); specifying the lower limit and the upper limit of the x and y axes.

Details

The selection gradient is operationalised as the ratio of the change in a trait due to deterministic selection Δα\Delta\alpha to the change in a trait due to stochastic drift σ\sigma. Values between -1 and 1 indicate parameter space where the change due to drift is greater than change due to selection on the trait. Conversely, values greater than 1 (or less than -1) indicate parameter space where positive (or negative) selection is stronger than drift.

If three or more traits were included in the model, other traits are held at their median values during the computations. Note that selection gradient plots can potentially produce misleading pictures of coevolutionary dynamics when other traits are held constant in models with three or more traits.

If the plot does not look right, the user might try zooming out from the default parameter space by setting wider limits. For some variables (e.g., continuous and count variables), the default limits may not be suitable.

Value

A ggplot object

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

See Also

coev_plot_flowfield

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# plot selection gradient
coev_plot_selection_gradient(
  object = fit,
  var1 = "political_authority",
  var2 = "religious_authority"
  )

## End(Not run)

Plot estimated trait values from a fitted coevfit object

Description

Produce a pairs plot of the estimated trait values for taxa from a fitted coevfit object. The plot includes scatterplot(s) of median estimated trait values, heatmap(s) of median estimated trait values, and density plots of the marginal distributions for each variable with associated posterior uncertainty.

Usage

coev_plot_trait_values(
  object,
  variables = NULL,
  ndraws = 50,
  tree_id = NULL,
  xlim = NULL,
  ylim = NULL
)

Arguments

object

An object of class coevfit
variables

If NULL (default), the function returns a pairs plot including all coevolving variables from the model. Otherwise, a character vector declaring the variables to be included.
ndraws

An integer indicating the number of draws to return in the density plots on the diagonal. The default and maximum number of draws is the size of the posterior sample.
tree_id

An integer indicating the tree ID to use when making posterior predictions. Set to NULL by default, which will use draws from every tree, integrating phylogenetic uncertainty.
xlim

Limits for the x-axis. If NULL (default), limits are set to the minimum and maximum estimated trait values.
ylim

Limits for the y-axis. If NULL (default), limits are set to the minimum and maximum estimated trait values.

Value

A patchwork of ggplot objects

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

References

Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods reveal that resource-use intensification drives the evolution of "complex" societies. EcoEvoRXiv. doi:10.32942/osf.io/wfp95

Sheehan, O., Watts, J., Gray, R. D., Bulbulia, J., Claessens, S., Ringen, E. J., & Atkinson, Q. D. (2023). Coevolution of religious and political authority in Austronesian societies. Nature Human Behaviour, 7(1), 38-45. 10.1038/s41562-022-01471-y

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# pairs plot of trait values
coev_plot_trait_values(fit)

## End(Not run)

Predict a co-evolutionary time series from a fitted coevfit object

Description

This function produces predicted values of traits over a time series, given the estimated parameters for a coevfit model. This function is used under the hood by the plotting function coev_pred_series.

Usage

coev_pred_series(
  object,
  eta_anc = NULL,
  intervention_values = NULL,
  tmax = 1,
  ntimes = 30,
  ndraws = NULL,
  stochastic = FALSE
)

Arguments

object

An object of class coevfit
eta_anc

If NULL (default), the starting values for the latent states η\eta will be set to the estimated ancestral states at the root from the coevfit model. Otherwise, a named list with length equal to the number of variables specifying the initial η\eta values. All variable names must be included in this list.
intervention_values

Either NULL (the default) or a named list of variables and associated intervention values. If NULL, all traits are free to vary. Otherwise, all coevolving variables must be declared separately in the named list without repetition. If the intervention value for a particular variable is set to NA, this variable is treated as a free variable. If the intervention value for a particular variable is specified, the variable is held constant at this trait value in the calculation. At least one variable must be declared as a free variable and at least one variable must be held constant (e.g., list(var1 = NA, var2 = 0)).
tmax

A positive number indicating the total duration of time for which to predict. Set to 1 by default, corresponding to the entire time depth of the original phylogenetic tree(s).
ntimes

A positive integer indicating the total number of discrete time steps over which to make predictions. Each step corresponds to a time difference of dt = tmax/ntimes. Set to 30 by default.
ndraws

An integer indicating the number of draws to return. The default and maximum number of draws is the size of the posterior sample.
stochastic

Logical (defaults to FALSE); indicator of whether predictions should include only the expected co-evolutionary change due to deterministic selection (FALSE) or also stochastic drift (TRUE). Incompatible with the 'intervention_values' argument.

Value

An [ndraw, ntimes, nvariables] array of predicted η\eta values

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

See Also

coev_plot_pred_series

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# simulate trait co-evolution
sims <- coev_pred_series(
  object = fit,
  stochastic = TRUE
  )

# expected trait co-evolution, no drift
epreds <- coev_pred_series(
  object = fit,
  stochastic = FALSE,
  eta_anc = list(political_authority = -2, religious_authority = 1.5)
  )

# expected trait co-evolution under intervention, no drift
epreds_intervention <- coev_pred_series(
  object = fit,
  stochastic = FALSE,
  intervention_values = list(political_authority = NA,
                             religious_authority = 1.5)
  )

## End(Not run)

Simulate the coevolution of multiple variables in discrete time steps

Description

This function simulates the coevolution of multiple continuous variables in discrete time steps following a simple VAR(1) autoregressive model. Users set the sample size, the variable names, the strength of selection and drift, and the probability of a speciation event in a given time step. The function returns a phylogeny, the results of the simulation run, and a dataset of contemporary trait values.

Usage

coev_simulate_coevolution(
  n,
  variables,
  selection_matrix,
  drift,
  prob_split,
  intercepts = NULL,
  ancestral_states = NULL
)

Arguments

n

Number of data points in the resulting data frame.
variables

A character vector of variable names (e.g., c("x","y"))
selection_matrix

A numeric matrix determining the strength of selection between variables. The matrix must have a number of rows and columns equal to the number of variables and its row and column names must contain all specified variables. Each cell determines the strength of selection from the column variable to the row variable. For example, the cell on the "x" column and the "x" row indicates how much previous values of x influence future values of x (autocorrelation). By contrast, the cell on the "x" column and the "y" row indicates how much previous values of x influence future values of y (cross-lagged effect).
drift

A named numeric vector specifying the strength of drift for different variables. Names must include all specified variables.
prob_split

A numeric probability of a species split in any given timestep.
intercepts

Intercepts for the VAR(1) model. If NULL (default), intercepts are set to zero for all variables. Otherwise, a named numeric vector specifying the intercepts for different variables. Names must include all specified variables.
ancestral_states

Ancestral states for different variables. If NULL (default), ancestral states are set to zero for all variables. Otherwise, a named numeric vector specifying the ancestral states for different variables. Names must include all specified variables.

Details

The model underlying this simulation is a simple VAR(1) autoregressive model, where values of all variables at the previous timestep predict values at the current timestep. In the case of two variables, the model is as follows:

Yt=αy+βy,yYt−1+βy,xXt−1+N(0,ϵy)Y_t = \alpha_{y}+\beta_{y,y}Y_{t-1}+\beta_{y,x}X_{t-1} + \mathcal{N}(0,\epsilon_{y})

Xt=αx+βx,xXt−1+βx,yYt−1+N(0,ϵx)X_t = \alpha_{x}+\beta_{x,x}X_{t-1}+\beta_{x,y}Y_{t-1} + \mathcal{N}(0,\epsilon_{x})

where α\alpha represents the intercepts, β\beta represents the selection matrix, and ϵ\epsilon represents the vector of drift parameters. With some probability pp, a speciation event creates two independent evolutionary branches. This simulation continues until the intended sample size of species has been reached. Note: this VAR(1) model is not the same as the generalized dynamic phylogenetic model that is estimated by the coev_fit() function (see the introductory vignette for mathematical details of generalized dynamic phylogenetic models).

Value

List with dataset at final timestep (data), full simulation log (simulation), and pruned phylogenetic tree (tree)

Author(s)

Scott Claessens [email protected], Erik Ringen [email protected]

Examples

# simulate coevolution of x and y
n <- 100
variables <- c("x", "y")
# x -> y but not vice versa
selection_matrix <- matrix(
  c(
    0.95, 0.00,
    0.80, 0.95
  ),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(variables, variables)
)
drift <- c("x" = 0.05, "y" = 0.05)
prob_split <- 0.05

# run simulation with default parameters
coev_simulate_coevolution(
  n, variables, selection_matrix,
  drift, prob_split
)

# set intercepts manually
coev_simulate_coevolution(
  n, variables, selection_matrix,
  drift, prob_split,
  intercepts = c("x" = 1, "y" = -1)
)

# set ancestral states manually
coev_simulate_coevolution(
  n, variables, selection_matrix,
  drift, prob_split,
  ancestral_states = c("x" = 1, "y" = -1)
)

Extract draws from jax_fit object

Description

Internal S3 method. Extracts draws from a jax_fit object, compatible with cmdstanr's draws() method.

Usage

## S3 method for class 'jax_fit'
draws(x, variables = NULL, ...)

Arguments

x

A jax_fit object.
variables

Character vector of variable names to extract.
...

Additional arguments (ignored).

Value

A draws_array object.

Examples

TRUE

Extract draws from nutpie_fit object

Description

Internal S3 method. Extracts draws from a nutpie_fit object, compatible with cmdstanr's draws() method.

Usage

## S3 method for class 'nutpie_fit'
draws(x, variables = NULL, ...)

Arguments

x

A nutpie_fit object.
variables

Character vector of variable names to extract. If NULL, extracts all variables.
...

Additional arguments (currently ignored).

Value

A draws_array object.

Examples

TRUE

Extract samples (draws) from a fitted coevfit object

Description

Extract samples (draws) from a fitted coevfit object

Usage

extract_samples(object)

Arguments

object

An object of class coevfit

Value

Samples in 'rethinking' style list format

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# extract samples from fitted model
extract_samples(fit)

## End(Not run)

Extract metadata from jax_fit object

Description

Internal S3 method. Extracts metadata from a jax_fit object, compatible with cmdstanr's metadata() method.

Usage

## S3 method for class 'jax_fit'
metadata(x)

Arguments

x

A jax_fit object.

Value

A list containing metadata.

Examples

TRUE

Extract metadata from nutpie_fit object

Description

Internal S3 method. Extracts metadata from a nutpie_fit object, compatible with cmdstanr's metadata() method.

Usage

## S3 method for class 'nutpie_fit'
metadata(x)

Arguments

x

A nutpie_fit object.

Value

A list containing metadata.

Examples

TRUE

Create a matrix of output plots from a coevfit object

Description

A pairs method that is customized for MCMC output.

Usage

## S3 method for class 'coevfit'
pairs(x, parameters = NULL, ...)

Arguments

x

An object of class coevfit
parameters

If NULL (default), the function returns plots for the main parameters from the model (the selection matrix "A", the drift matrix "Q", and the continuous time intercepts "b"). Otherwise, a character vector declaring the parameters to plot (strictly case-sensitive).
...

Further arguments to be passed to mcmc_pairs.

Details

For a detailed description see mcmc_pairs.

Value

A matrix of output plots

Author(s)

Scott Claessens [email protected]

See Also

plot.coevfit

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# print pairs plot
pairs(fit)

## End(Not run)

Trace and density plots for coevfit objects

Description

Trace and density plots for coevfit objects

Usage

## S3 method for class 'coevfit'
plot(
  x,
  parameters = NULL,
  combo = c("dens", "trace"),
  npars = 5,
  plot = TRUE,
  ...
)

Arguments

x

An object of class coevfit.
parameters

If NULL (default), the function returns plots for the main parameters from the model (the selection matrix "A", the drift matrix "Q", and the continuous time intercepts "b"). Otherwise, a character vector declaring the parameters to plot (strictly case-sensitive).
combo

A character vector with at least two elements. Each element of combo corresponds to a column in the resulting plot and should be the name of one of the available MCMC functions (omitting the mcmc_ prefix).
npars

Number of parameters displayed per page (default = 5).
plot

Logical. If TRUE (default), plots are plotted in the active graphic device. If FALSE, plots are not plotted.
...

Not used.

Value

A list of ggplot objects

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# print trace and density plots
plot(fit)

## End(Not run)

Brain weight, body weight, diet, and sociality in primates

Description

This dataset and associated phylogeny describes brain weights, body weights, dietary categories, social systems, mating systems, and group sizes for 143 primate species. These data were originally compiled in DeCasien et al. (2017). The phylogeny is a 50% majority rule consensus tree from the 10k Trees website.

Usage

primates

Format

A list containing a dataset and an associated phylogeny. The phylogeny is a pruned consensus tree for 143 species of primate (excluding humans). The dataset is a data frame with 143 observations and 11 variables:

species

The primate species

clade

The clade of the primate species

brain_weight

Average brain weight, in grams

body_weight

Average body weight, in grams

diet

A factor with four levels indicating whether the species' dietary category is folivore, frugivore, frugivore/folivore, or omnivore

percent_fruit

Percentage of fruit in diet

social_system

A factor with four levels indicating whether the species' social system is solitary, pair-living, harem polygyny, or polygynandry

mating_system

A factor with five levels indicating whether the species' mating system is spatial polygyny, monogamy, polyandry, harem polygyny, or polygynandry

group_size

Average group size

Source

DeCasien, A. R., Williams, S. A., & Higham, J. P. (2017). Primate brain size is predicted by diet but not sociality. Nature Ecology & Evolution, 1(5), 0112.

Arnold, C., Matthews, L. J., & Nunn, C. L. (2010). The 10kTrees website: a new online resource for primate phylogeny. Evolutionary Anthropology, 19(3), 114-118.

https://10ktrees.nunn-lab.org/Primates/index.html

Examples

## Not run: 
# fit model to primates data
m <-
  coev_fit(
    data = primates$data,
    variables = list(
      brain_weight = "lognormal",
      body_weight = "lognormal"
    ),
    id = "species",
    tree = primates$phylogeny,
    # arguments to cmdstanr::sample()
    parallel_chains = 4,
    seed = 1
  )
# print model summary
summary(m)
# plot delta theta
coev_plot_delta_theta(m)

## End(Not run)

Print a summary of a fitted model represented by a coevfit object

Description

Print a summary of a fitted model represented by a coevfit object

Usage

## S3 method for class 'coevfit'
print(x, digits = 2, ...)

Arguments

x

An object of class coevfit
digits

The number of significant digits for printing out the summary; defaults to 2
...

Additional arguments that would be passed to method summary of coevfit.

Value

A printed summary of a coevfit model object

See Also

summary.coevfit

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# print summary of model
print(fit)

## End(Not run)

Example dataset with repeated observations

Description

This example dataset and associated phylogeny are used as an example of repeated observations in phylogenetic modelling. The data are adapted from de Villemeruil & Nakagawa (2014) to include only the first twenty species.

Usage

repeated

Format

A list containing a dataset and an associated phylogeny. The dataset is a data frame with 100 observations and 3 variables:

species

The name of the species

x

An example continuous variable

y

An example continuous variable

Source

de Villemeruil P. & Nakagawa, S. (2014). General quantitative genetic methods for comparative biology. In L. Garamszegi (Ed.), Modern phylogenetic comparative methods and their application in evolutionary biology: concepts and practice (pp. 287-303). Springer, New York.

Examples

## Not run: 
# fit model to repeated data
m <-
  coev_fit(
    data = repeated$data,
    variables = list(
      x = "normal",
      y = "normal"
    ),
    id = "species",
    tree = repeated$phylogeny,
    # arguments to cmdstanr::sample()
    parallel_chains = 4,
    seed = 1
  )
# print model summary
summary(m)
# plot delta theta
coev_plot_delta_theta(m)

## End(Not run)

Save a fitted coevfit object to a file

Description

A wrapper around saveRDS that saves a fitted coevfit model object while ensuring that all posterior draws and diagnostics are correctly saved by cmdstanr.

Usage

save_coevfit(object, file, ...)

Arguments

object

An object of class coevfit
file

A string declaring the path where the file should be saved
...

Other arguments to pass to saveRDS besides object and file

Value

An .RDS file containing the fitted coevfit model object

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# save model to .rds file
save_coevfit(fit, file = "model.rds")

## End(Not run)

Expose the Stan code of a fitted model represented by a coevfit object

Description

Expose the Stan code of a fitted model represented by a coevfit object

Usage

stancode(object)

Arguments

object

An object of class coevfit.

Value

Printed Stan code

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# expose stan code
stancode(fit)

## End(Not run)

Expose the Stan data list for a fitted model represented by a coevfit object

Description

Expose the Stan data list for a fitted model represented by a coevfit object

Usage

standata(object)

Arguments

object

An object of class coevfit.

Value

Named list of data for Stan

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# expose stan data
standata(fit)

## End(Not run)

Create a summary of a fitted model represented by a coevfit object

Description

Create a summary of a fitted model represented by a coevfit object

Usage

## S3 method for class 'coevfit'
summary(object, prob = 0.95, robust = FALSE, ...)

Arguments

object

An object of class coevfit.
prob

A value between 0 and 1 indicating the desired probability to be covered by the uncertainty intervals. The default is 0.95.
robust

If FALSE (the default) the mean is used as the measure of central tendency and the standard deviation as the measure of variability. If TRUE, the median and the median absolute deviation are applied instead.
...

Other potential arguments

Details

The convergence diagnostics rhat, ess_bulk, and ess_tail are described in detail in Vehtari et al. (2020).

Value

List of class coevsummary

References

Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner (2020). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian Analysis. 1–28. dpi:10.1214/20-BA1221

Examples

## Not run: 
# fit dynamic coevolutionary model
fit <- coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
  )

# print summary of model
summary(fit)

## End(Not run)

Summary statistics for jax_fit object

Description

Internal S3 method. Computes summary statistics from a jax_fit object, compatible with cmdstanr's summary() method. This method mimics cmdstanr's interface where summary() is called with named arguments where names are column names and values are function names or formulas.

Usage

## S3 method for class 'jax_fit'
summary(object, variables = NULL, ...)

Arguments

object

A jax_fit object.
variables

Character vector of variable names to summarize.
...

Named summary spec arguments (cmdstanr style).

Value

A data.frame with summary statistics.

Examples

TRUE

Compute summary statistics for nutpie_fit object

Description

Internal S3 method. Computes summary statistics from a nutpie_fit object, compatible with cmdstanr's summary() method. This method mimics cmdstanr's interface where summary() is called with named arguments where names are column names and values are function names or formulas.

Usage

## S3 method for class 'nutpie_fit'
summary(object, variables = NULL, ...)

Arguments

object

A nutpie_fit object.
variables

Character vector of variable names to summarize. If NULL, summarizes all variables.
...

Named arguments where names are output column names and values are function names (strings) or formulas for computing statistics.

Value

A data.frame with summary statistics, compatible with cmdstanr format.

Examples

TRUE