Package {phylopred}

Type:	Package
Title:	Predicting Species Interactions from Phylogenetic Distance
Version:	0.1.0
Description:	Tools to model and predict host-parasite (and other pairwise species) interactions from phylogenetic distance under presence-only data. Provides pairwise host-sharing logistic models with cluster-aware bootstrap and slope comparison between clades, presence-only performance metrics (a threshold efficiency metric and a phylogenetic partial receiver operating characteristic (ROC) adapted from Peterson, Papes and Soberon (2008) <doi:10.1016/j.ecolmodel.2007.11.008>), positive-unlabeled host prediction, and leave-hosts-out validation. Methods build on Robles-Fernandez and Lira-Noriega (2017) <doi:10.3389/fams.2017.00017>.
License:	MIT + file LICENSE
URL:	https://github.com/alrobles/phylopred, https://alrobles.github.io/phylopred/
BugReports:	https://github.com/alrobles/phylopred/issues
Encoding:	UTF-8
LazyData:	true
Imports:	graphics, stats, utils, xplus
Suggests:	ape, cofid, knitr, rmarkdown, testthat (≥ 3.0.0)
Depends:	R (≥ 4.1.0)
Config/testthat/edition:	3
Config/roxygen2/version:	8.0.0
NeedsCompilation:	no
Packaged:	2026-09-30 15:18:57 UTC; alrobles
Author:	Angel Robles [aut, cre]
Maintainer:	Angel Robles <a.l.robles.fernandez@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-14 12:12:07 UTC

Align incidence matrix and phylogenetic distance matrix

Description

Validates both inputs and reorders the columns of incidence to match the column order of phydist.

Usage

align_phylopred_inputs(incidence, phydist)

Arguments

Value

A named list with:

incidence

The incidence matrix with columns reordered to match phydist.

phydist

The (unchanged) phylogenetic distance matrix.

Examples

incidence <- prepare_incidence_matrix(beetleTreeInteractions)
aligned <- align_phylopred_inputs(incidence, phy_dist)

Beetle tree interaction table.

Description

A dataset of bark beetle species and their host tree genera. Original information is available in https://www.barkbeetles.info/

Usage

beetleTreeInteractions

Format

A data frame with 565 rows and 2 variables:

beetle

Species of bark beetles

genus

Genus of host tree of bark beetles

Details

We use this interaction tables to construct an incidence matrix. This is a table representation of a bipartite graph

Source

https://www.barkbeetles.info/

Bootstrap a phylopred logistic regression model

Description

Runs the phylopred logistic regression n times using random focal host sampling, returning a distribution of coefficient estimates.

Usage

bootstrap_phylopred_model(incidence, phydist, n = 1000, seed = NULL)

Arguments

Value

An S3 object of class "phylopred_bootstrap", a list with:

call

The matched call.

n

Number of bootstrap iterations.

seed

The seed used.

coefficients

A matrix of n rows x 12 columns (coefficient info from each iteration).

summary

A data.frame with the column-wise means (same 12 columns).

metadata

List with n_parasites and n_hosts.

See Also

fit_phylopred_model, predict_phylopred_probability

Examples

incidence <- prepare_incidence_matrix(beetleTreeInteractions)
aligned   <- align_phylopred_inputs(incidence, phy_dist)
bootstrap_phylopred_model(aligned$incidence, aligned$phydist, n = 5, seed = 42)

Cluster bootstrap for pairwise host-sharing models

Description

Nonparametric bootstrap that resamples whole parasites (clusters) with replacement and refits the logistic regression on each replicate. This respects the non-independence of rows that share a parasite in a prepare_pair_data data set, unlike a naive row-level bootstrap. When a group column is present, resampling is stratified within groups and the model includes a group-by-distance interaction.

Usage

cluster_bootstrap_phylopred(
  pairs,
  n_boot = 1000,
  seed = NULL,
  conf_level = 0.95
)

Arguments

Value

An object of class "phylopred_cluster_bootstrap", a list with:

draws

Matrix of bootstrap coefficient draws (replicates x terms).

estimates

Coefficients of the model fitted to the full data.

confidence_intervals

Percentile confidence intervals per term.

n_boot, seed, conf_level

Bootstrap settings.

See Also

prepare_pair_data, compare_phylopred_slopes

Examples

# Three parasites keep the example fast
top <- names(sort(table(beetleTreeInteractions[[1]]),
                  decreasing = TRUE))[1:3]
sub <- beetleTreeInteractions[
  as.character(beetleTreeInteractions[[1]]) %in% top, ]
pairs <- prepare_pair_data(sub, phy_dist)
cluster_bootstrap_phylopred(pairs, n_boot = 5, seed = 42)

Compare host-sharing slopes between two clades

Description

Fits a single logistic regression with a group-by-phylogenetic-distance interaction to a pairwise data set built with prepare_pair_data and tests whether the slopes of the two groups differ. A more negative slope means the probability of sharing a parasite decays faster with phylogenetic distance, i.e. the clade is more of a phylogenetic specialist. Uncertainty for the slope difference comes from a parasite-level cluster bootstrap.

Usage

compare_phylopred_slopes(pairs, n_boot = 1000, seed = NULL, conf_level = 0.95)

Arguments

Value

An object of class "phylopred_slope_comparison", a list with:

groups

The two group labels (reference first).

slopes

Named vector with the slope of each group.

slope_difference

Slope of the second group minus the reference.

interaction_p_value

Wald p-value of the interaction term.

bootstrap

The cluster_bootstrap_phylopred object.

slope_cis

Bootstrap percentile CIs for each group's slope and for the difference.

fitted_model

The underlying glm object.

See Also

prepare_pair_data, cluster_bootstrap_phylopred

Examples

# Compare Xyleborus beetles against other genera (subset kept small
# so the bootstrap example runs quickly)
parasites <- unique(as.character(beetleTreeInteractions[[1]]))
keep <- c(grep("^Xyleborus", parasites, value = TRUE)[1:3],
          parasites[!grepl("^Xyleborus", parasites)][1:3])
sub <- beetleTreeInteractions[
  as.character(beetleTreeInteractions[[1]]) %in% keep, ]
pairs <- prepare_pair_data(
  sub, phy_dist,
  group = function(x) ifelse(grepl("^Xyleborus", x), "Xyleborus", "other")
)
compare_phylopred_slopes(pairs, n_boot = 5, seed = 42)

Evaluate the predictive capacity of a pairwise host-sharing model

Description

Computes goodness-of-fit and discrimination metrics for a logistic host-sharing model fitted to a prepare_pair_data data set: McFadden and Tjur pseudo-R-squared, the area under the ROC curve (AUC), and optionally a grouped cross-validated AUC where whole parasites are held out (parasite-grouped k-fold), which measures the ability to predict the host range of unseen parasites.

Usage

evaluate_phylopred_model(pairs, cv = FALSE, k = 5, seed = NULL)

Arguments

Value

An object of class "phylopred_evaluation", a list with:

mcfadden_r2

McFadden pseudo-R-squared, 1 - logLik(model)/logLik(null).

tjur_r2

Tjur coefficient of discrimination (mean fitted probability for 1s minus mean for 0s).

auc

In-sample area under the ROC curve.

cv_auc

Parasite-grouped k-fold cross-validated AUC (NA if cv = FALSE).

fitted_model

The underlying glm object.

n_pairs, n_parasites

Data set size.

See Also

prepare_pair_data, cluster_bootstrap_phylopred

Examples

# Three parasites keep the example fast
top <- names(sort(table(beetleTreeInteractions[[1]]),
                  decreasing = TRUE))[1:3]
sub <- beetleTreeInteractions[
  as.character(beetleTreeInteractions[[1]]) %in% top, ]
pairs <- prepare_pair_data(sub, phy_dist)
evaluate_phylopred_model(pairs)
evaluate_phylopred_model(pairs, cv = TRUE, k = 3, seed = 42)

Extract coefficients from a phylopred model or bootstrap object

Description

Returns a named list of model coefficients and their associated statistics from a phylopred_model or phylopred_bootstrap object.

Usage

extract_phylopred_coefficients(model)

Arguments

Value

A named list with elements:

intercept

Intercept coefficient (mean for bootstrap).

slope

Slope coefficient (mean for bootstrap).

intercept_se

Standard error of intercept.

slope_se

Standard error of slope.

intercept_ci

Two-element vector with 2.5% and 97.5% CI for intercept.

slope_ci

Two-element vector with 2.5% and 97.5% CI for slope.

Examples

inc <- prepare_incidence_matrix(beetleTreeInteractions)
al  <- align_phylopred_inputs(inc, phy_dist)
model <- fit_phylopred_model(al$incidence, al$phydist, seed = 42)
extract_phylopred_coefficients(model)

Phylogenetic tree of marine fish hosts of Caligus and Lepeophtheirus

Description

A phylogenetic tree (class phylo) of 625 marine fish species that are recorded hosts of sea lice of the genera Caligus and Lepeophtheirus in the cofid copepod-fish interaction database. The tree is pruned from the Fish Tree of Life actinopterygian phylogeny (Rabosky and colleagues 2018, actinopt_12k_raxml). Tip labels are binomial species names with spaces.

Usage

fish_tree

Format

A phylo object with 625 tips.

Source

https://fishtreeoflife.org/. The underlying data are released under the CC0 1.0 public domain dedication (Dryad, doi:10.5061/dryad.fc71cp4). Please cite Rabosky, D. L. and colleagues (2018) An inverse latitudinal gradient in speciation rate for marine fishes. Nature 559, 392-395. doi:10.1038/s41586-018-0273-1

See Also

The cofid package for the interaction records.

Fit a phylopred logistic regression model

Description

Fits a logistic regression model relating host susceptibility to phylogenetic distance.

Usage

fit_phylopred_model(incidence, phydist, seed = NULL)

Arguments

Value

An S3 object of class "phylopred_model", a list with:

call

The matched call.

coefficients

Named numeric vector c(intercept, slope).

standard_errors

Named numeric vector c(intercept, slope).

confidence_intervals

Matrix with rows c("intercept","slope") and cols c("2.5%","97.5%").

fitted_model

The underlying glm object.

training_data

The data.frame used for model fitting.

convergence

Logical; whether the GLM converged.

metadata

List with n_parasites, n_hosts, and seed.

See Also

bootstrap_phylopred_model, predict_phylopred_probability

Examples

incidence <- prepare_incidence_matrix(beetleTreeInteractions)
aligned   <- align_phylopred_inputs(incidence, phy_dist)
fit_phylopred_model(aligned$incidence, aligned$phydist, seed = 42)

Positive-unlabeled host prediction via xplus (PLUS algorithm)

Description

Predicts the hosts of a parasite from the phylogeny using positive-unlabeled (PU) learning, as a parallel alternative to the pairwise logistic model. Known hosts are treated as positives and every other species in the phylogeny as unlabeled (never as a true absence). Each species is described by phylogenetic features and the PLUS algorithm (xplus::xplus) iteratively re-labels the unlabeled species while fitting a penalized logistic model.

Usage

fit_phylopred_pu(
  interactions,
  phydist,
  parasite,
  n_eigen = 10,
  seed = NULL,
  ...
)

Arguments

Details

The feature matrix contains, for every species in phydist:

• the minimum and mean phylogenetic distance to the parasite's known hosts (excluding the species itself), and

• the first n_eigen principal coordinates (PCoA) of the phylogenetic distance matrix, describing the global position of the species in the phylogeny (the analogue of environmental layers in an SDM).

Value

An object of class "phylopred_pu": a list with

prediction

data.frame with species, pred (predicted host suitability) and is_host.

fit

The underlying xplus object.

parasite

The parasite name.

The prediction columns can be passed directly to host_threshold_metric and phylo_partial_roc.

See Also

host_threshold_metric, phylo_partial_roc, prepare_pair_data

Examples

pu <- fit_phylopred_pu(beetleTreeInteractions, phy_dist,
                     parasite = beetleTreeInteractions[[1]][1],
                     max_iter = 20, seed = 1)
host_threshold_metric(pu$prediction$pred, pu$prediction$is_host)

Host specificity index (S_TD) for every parasite in an interaction table

Description

Applies poulin_std to every parasite in an interaction table, following the same expansion logic as prepare_pair_data. Useful for contrasting the Poulin and Mouillot (2003) specificity index against phylopred's model-based host-range metric b* (host_threshold_metric) across many parasites.

Usage

host_specificity_index(interactions, phydist, group = NULL)

Arguments

Value

A data.frame of class "phylopred_specificity" with columns parasite, n_hosts (known hosts present in phydist), std (S_{TD}), var_std (Var(S_{TD})), and group (only if group is supplied).

References

Poulin, R. and Mouillot, D. (2003) Parasite specialization from a phylogenetic perspective: a new index of host specificity. Parasitology 126, 473-480. doi:10.1017/S0031182003002993

See Also

poulin_std, prepare_pair_data, host_threshold_metric

Examples

std_by_parasite <- host_specificity_index(beetleTreeInteractions, phy_dist)
head(std_by_parasite)

Threshold metric f_gamma for predicted host suitability

Description

Computes the threshold-dependent performance metric

f_\gamma(t) = a(t) / (\gamma b(t) + c(t))

where, for a cutoff t on the predicted host-suitability scores:

a(t)

proportion of known hosts with prediction \ge t (sensitivity / recall);

b(t)

proportion of ALL species in the phylogeny with prediction \ge t (predicted host breadth, the phylogenetic analogue of predicted area in SDM threshold metrics);

c(t) = 1 - a(t)

proportion of known hosts lost (omission).

The metric is designed for presence-only / positive-unlabeled settings: it never uses absences, and gamma tunes the trade-off between retaining known hosts and predicting a compact host set. The maximum of f_\gamma(t) gives an optimal cutoff t^*.

Usage

host_threshold_metric(pred, is_host, gamma = 1)

Arguments

Value

An object of class "phylopred_threshold_metric": a data.frame with columns threshold, a, b, c, and f_gamma, evaluated at every distinct prediction value. The optimal cutoff is stored in attributes "optimal_threshold" and "optimal_f", and gamma in attribute "gamma".

See Also

phylo_partial_roc, evaluate_phylopred_model

Examples

pairs <- prepare_pair_data(beetleTreeInteractions, phy_dist)
fit <- stats::glm(suscept ~ phydist, data = pairs, family = stats::binomial())
# Suitability of every host species for one focal parasite
one <- pairs[pairs$parasite == pairs$parasite[1], ]
one <- one[!duplicated(one$target), ]
pred <- stats::predict(fit, newdata = one, type = "response")
m <- host_threshold_metric(pred, one$suscept, gamma = 1)
plot(m)

Phylogenetic distance matrix among host data on phylopred

Description

Data of phylogenetic distance matrix for bark beetles host species

Usage

phy_dist

Format

A matrix with 235 rows and 235 columns. Contain Phylogenetic distance matrix between bark beetle host genera

Source

https://www.barkbeetles.info/

Phylogenetic partial ROC (Peterson, Papes & Soberon 2008 analogue)

Description

Partial ROC test for presence-only host-suitability predictions, following Peterson, Papes & Soberon (2008, Ecological Modelling 213, 63-72). The ROC space is redefined for presence-only data: the y axis is the sensitivity a(t) (proportion of known hosts predicted at cutoff t) and the x axis is the proportion of ALL species in the phylogeny predicted as hosts, b(t) (the analogue of proportional predicted area). The curve is restricted to the low-omission region a(t) \ge 1 - E for a user-chosen omission tolerance error_rate E, and performance is summarized as the ratio of the partial AUC of the observed curve to the partial AUC of the null (random-ranking) expectation, the diagonal a = b. A ratio of 1 equals random performance; ratios above 1 indicate predictive signal. Significance is assessed by bootstrap resampling of the known hosts.

Usage

phylo_partial_roc(
  pred,
  is_host,
  error_rate = 0.05,
  n_boot = 500,
  boot_prop = 0.5,
  seed = NULL
)

Arguments

Value

An object of class "phylopred_partial_roc", a list with:

auc_ratio

Observed partial AUC ratio (observed / null).

p_value

Bootstrap proportion of replicates with ratio \le 1.

boot_ratios

Vector of bootstrap AUC ratios.

curve

data.frame with threshold, a (sensitivity) and b (predicted breadth).

error_rate, n_boot, boot_prop

The configuration used.

See Also

host_threshold_metric

Examples

pairs <- prepare_pair_data(beetleTreeInteractions, phy_dist)
fit <- stats::glm(suscept ~ phydist, data = pairs, family = stats::binomial())
one <- pairs[pairs$parasite == pairs$parasite[1], ]
one <- one[!duplicated(one$target), ]
pred <- stats::predict(fit, newdata = one, type = "response")
phylo_partial_roc(pred, one$suscept, n_boot = 50, seed = 42)

Joint confidence ellipse for intercept and slope

Description

Computes the boundary of the joint confidence region ("ellipse") for a pair of coefficients, either from a fitted glm (using the asymptotic covariance matrix) or from a cluster_bootstrap_phylopred object (using the empirical covariance of the bootstrap draws, which respects parasite-level clustering).

Usage

phylopred_confidence_ellipse(
  object,
  terms = c("(Intercept)", "phydist"),
  level = 0.95,
  n_points = 200
)

Arguments

Value

A data.frame with n_points rows and two columns named after terms, tracing the ellipse boundary. The center (point estimates) is returned as the "center" attribute.

See Also

cluster_bootstrap_phylopred

Examples

# Three parasites keep the example fast
top <- names(sort(table(beetleTreeInteractions[[1]]),
                  decreasing = TRUE))[1:3]
sub <- beetleTreeInteractions[
  as.character(beetleTreeInteractions[[1]]) %in% top, ]
pairs <- prepare_pair_data(sub, phy_dist)

# From a glm: asymptotic ellipse
fit <- glm(suscept ~ phydist, data = pairs, family = binomial())
ell <- phylopred_confidence_ellipse(fit)
plot(ell, type = "l")
points(t(attr(ell, "center")), pch = 19)

# From a cluster bootstrap
boot <- cluster_bootstrap_phylopred(pairs, n_boot = 5, seed = 42)
ell_b <- phylopred_confidence_ellipse(boot)
lines(ell_b, lty = 2)

Poulin and Mouillot taxonomic/phylogenetic distinctness index (S_TD)

Description

Computes the average taxonomic (or phylogenetic) distinctness S_{TD} of Poulin and Mouillot (2003) for a single parasite's known host set:

S_{TD} = \frac{2 \sum_{i<j} \omega_{ij}}{s(s-1)}

where \omega_{ij} is the distance between hosts i and j taken from phydist, and s is the number of known hosts. When phydist is a taxonomic step matrix (see taxonomic_step_matrix), S_{TD} is the original 1-5 taxonomic index; when phydist is a matrix of phylogenetic (patristic) distances, this is the continuous phylogenetic generalization noted as possible by Poulin and Mouillot (2003, after Clarke and Warwick, 2001), expressed directly in branch-length units. The corresponding sampling variance Var(S_{TD}) = \sum_{i \ne j} (\omega_{ij} - S_{TD})^2 / [s(s-1)] is also returned.

Usage

poulin_std(phydist, hosts)

Arguments

Value

An object of class "phylopred_std": a list with elements n_hosts, std (S_{TD}), and var_std (Var(S_{TD}), NA when n_hosts < 3). S_{TD} is undefined (NA, with a warning) when n_hosts < 2, as noted by Poulin and Mouillot (2003) for single-host parasites.

References

Poulin, R. and Mouillot, D. (2003) Parasite specialization from a phylogenetic perspective: a new index of host specificity. Parasitology 126, 473-480. doi:10.1017/S0031182003002993

Clarke, K.R. and Warwick, R.M. (2001) A further biodiversity index applicable to species lists: variation in taxonomic distinctness. Marine Ecology Progress Series 216, 265-278.

See Also

host_specificity_index, taxonomic_step_matrix, host_threshold_metric

Examples

hosts <- unique(beetleTreeInteractions$genus[
  beetleTreeInteractions$beetle == beetleTreeInteractions$beetle[1]])
poulin_std(phy_dist, hosts)

Predict host susceptibility probabilities from logistic regression coefficients

Description

Applies a logit function using intercept and slope coefficients to a vector of phylogenetic distances. Accepts a numeric vector, a phylopred_model, or a phylopred_bootstrap object.

Usage

predict_phylopred_probability(coef, dist)

Arguments

Value

A named numeric vector of predicted probabilities of the same length as dist.

See Also

fit_phylopred_model

Examples

inc <- prepare_incidence_matrix(beetleTreeInteractions)
al  <- align_phylopred_inputs(inc, phy_dist)
model <- fit_phylopred_model(al$incidence, al$phydist, seed = 42)
dist_vec <- seq(0, max(al$phydist), length.out = 50)
predict_phylopred_probability(model, dist_vec)

Prepare an incidence matrix from an interaction table

Description

Builds an incidence (presence/absence) matrix from a two-column data.frame of species interactions.

Usage

prepare_incidence_matrix(db)

Arguments

Value

An incidence matrix (rows = first column unique values, cols = second column unique values) filled with 0/1.

Examples

prepare_incidence_matrix(beetleTreeInteractions)

Build a pairwise host-sharing data set from an interaction table

Description

Expands an interaction table into all (focal host, target host) pairs per parasite. For every parasite, each of its known hosts is used as a focal host, and every other host in phydist is a target host whose response is 1 if it is also a known host of that parasite and 0 otherwise.

Usage

prepare_pair_data(interactions, phydist, group = NULL)

Arguments

Details

Unlike the legacy workflow (which samples a single focal host per parasite), this uses every known host as a focal host, so no information is discarded and results are deterministic. Rows belonging to the same parasite are not independent; account for this with cluster_bootstrap_phylopred when computing uncertainty.

Value

A data.frame of class "phylopred_pairs" with columns:

parasite

Parasite species (cluster identifier).

focal

Focal (known) host.

target

Target host.

phydist

Distance between focal and target host.

suscept

1 if the target host is a known host of the parasite.

group

Grouping factor (only if group is supplied).

See Also

fit_phylopred_model, cluster_bootstrap_phylopred, compare_phylopred_slopes

Examples

pairs <- prepare_pair_data(beetleTreeInteractions, phy_dist)
head(pairs)

Build a taxonomic step-distance matrix (Poulin & Mouillot scale)

Description

Constructs a symmetric taxonomic distance matrix among species from a taxonomic hierarchy table, following the step scale of Poulin and Mouillot (2003): 1 if two species share a genus, 2 if they share a family but not a genus, 3 if they share an order but not a family, 4 if they share a class but not an order, and 5 if they share none of the above (different class / phylum). Any subset of genus, family, order, class may be supplied; missing ranks are skipped and the maximum step for that pair falls back to length(ranks) + 1.

Usage

taxonomic_step_matrix(tax_table)

Arguments

Value

A square numeric matrix with row/column names equal to tax_table$species, diagonal 0, suitable as the phydist argument of poulin_std or host_specificity_index.

References

Poulin, R. and Mouillot, D. (2003) Parasite specialization from a phylogenetic perspective: a new index of host specificity. Parasitology 126, 473-480. doi:10.1017/S0031182003002993

See Also

poulin_std, host_specificity_index

Examples

tax <- data.frame(
  species = c("Caligus elongatus", "Caligus curtus",
              "Lepeophtheirus salmonis", "Ergasilus sieboldi"),
  genus  = c("Caligus", "Caligus", "Lepeophtheirus", "Ergasilus"),
  family = c("Caligidae", "Caligidae", "Caligidae", "Ergasilidae"),
  order  = c("Siphonostomatoida", "Siphonostomatoida",
              "Siphonostomatoida", "Cyclopoida")
)
taxonomic_step_matrix(tax)

Leave-hosts-out validation of host predictions

Description

Presence-only predictive validation: a fraction of the parasite's known hosts is hidden, the model is refitted on the remaining hosts, and the hidden hosts are scored against the unlabeled species. Because true absences do not exist, performance is measured as the AUC of the hidden (held-out) hosts versus the unlabeled species — the probability that a truly interacting species outranks a random unlabeled one. This measures actual predictive ability, not in-sample fit.

Usage

validate_host_prediction(
  interactions,
  phydist,
  parasite,
  method = c("glm", "pu"),
  prop_holdout = 0.25,
  n_rep = 10,
  seed = NULL,
  ...
)

Arguments

Value

An object of class "phylopred_holdout": a list with auc (vector of holdout AUCs, one per replicate), mean_auc, method, parasite, prop_holdout, and n_rep.

See Also

fit_phylopred_pu, host_threshold_metric, phylo_partial_roc

Examples

parasite <- names(sort(table(beetleTreeInteractions[[1]]),
                       decreasing = TRUE))[1]
validate_host_prediction(beetleTreeInteractions, phy_dist, parasite,
                         method = "glm", n_rep = 3, seed = 1)