Probabilistic establishment from invasion fitness — compute_establishment

Invasion fitness $\lambda_{is}$ integrates trait space geometry (distances, overlaps, convex hull, cloud centroid) with abiotic suitability (alignment of invader traits to local environment), niche crowding (overlap with resident trait space weighted by composition), and resident competition (site saturation). compute_establishment_probability() maps $\lambda_{is}$ to probabilities of establishment using a unified interface:

Probit: $P = \Phi(\lambda/\sigma)$ with $\sigma$ as a scalar, the residual SD from a fitted auxiliary GLMM, or a cell-wise predictive SD.
Logistic: $P = \mathrm{logit}^{-1}(\lambda / \tau)$ with scale $\tau$.
Hard rule: $P = \mathbb{I}\{\lambda>0\}$ for a binary map.

If lambda_is is not supplied, the function will build it from standardized predictors via $\lambda_{is} = \gamma \, r^{(z)}_{is} - \alpha \, C^{(z)}_{is} - \beta \, S^{(z)}_{is} + \kappa$.

Usage

compute_establishment_probability(
  lambda_is = NULL,
  r_is_z = NULL,
  C_is_z = NULL,
  S_is_z = NULL,
  gamma = 1,
  alpha = NULL,
  beta = NULL,
  kappa = 0,
  method = c("probit", "logit", "hard"),
  sigma = NULL,
  tau = 1,
  fit = NULL,
  predictive = FALSE,
  sigma_mat = NULL,
  use_vcov = FALSE,
  Q_inv = NULL,
  site_df = NULL,
  return_long = TRUE,
  make_plots = TRUE,
  option_label = NULL
)

Arguments

lambda_is: Optional matrix $S\times I$ of invasion fitness (rows = sites, cols = invaders). If NULL, it is computed from components (see below).
r_is_z, C_is_z, S_is_z: Optional matrices $S\times I$ of standardized abiotic suitability, niche crowding, and saturation. Used only if lambda_is=NULL.
gamma: Optional vector (length $I$) or matrix $S\times I$ for the abiotic slope $\gamma$; recycled to $S\times I$ as needed.
alpha: Optional vector (length $I$) or matrix $S\times I$ of crowding penalties $\alpha \ge 0$ by convention); recycled if needed.
beta: Optional vector (length $I$) of saturation penalties $\beta$. If you allow facilitation, pass signed values here.
kappa: Optional scalar calibration offset added to $\lambda$; default 0.
method: Character, one of c("probit","logit","hard").
sigma: Numeric scalar SD for probit (ignored for other methods). If not given and fit is supplied, sigma defaults to sigma(fit).
tau: Numeric scalar scale for logistic (denominator inside logit); default 1.
fit: Optional fitted model (e.g., glmmTMB), used to obtain sigma(fit) when method="probit" and sigma is missing.
predictive: Logical; for method="probit", if TRUE the function uses a predictive SD matrix instead of a scalar. Provide it via sigma_mat, or set use_vcov=TRUE to compute it from fit (see below).
sigma_mat: Optional matrix $S\times I$ of SDs (e.g., predictive SD per cell).
use_vcov: Logical; if TRUE and method="probit" with predictive=TRUE, the function will attempt to compute mean-SE from vcov(fit) and add the residual SD to yield a predictive SD. Requires fit, Q_inv, and the standardized matrices.
Q_inv: Optional data frame with invader trait-plane scores (tr1,tr2), rownames = invader IDs. Required only if use_vcov=TRUE.
site_df: Optional site table with columns site,x,y for mapping; rownames of the matrices must match site_df$site.
return_long: Logical; if TRUE, include a tidy tibble in the output.
make_plots: Logical; if TRUE, return ggplot objects (site mean map, invader ranking, and site×invader heatmap).
option_label: Optional label attached to the tidy output and plot titles.

Value

A list with:

p_is: matrix $S\times I$ of probabilities (or 0/1 for hard);
lambda_is: the $S\times I$ fitness matrix used;
sigma_used: scalar or matrix actually used for probit (else NULL);
method, option_label;
prob_long: tidy tibble (if return_long=TRUE);
plots: list of ggplots (if make_plots=TRUE): site_mean, invader_mean, heatmap.

Details

Convert invasion fitness to probabilistic establishment (probit/logit/hard)

Probit: $P_{is}=\Phi(\lambda_{is}/\sigma)$. Use a scalar sigma (fast), the residual SD from fit, or set predictive=TRUE to inject a cell-wise SD: either supply sigma_mat directly or compute it with use_vcov=TRUE, which uses the fixed-effect covariance from vcov(fit) plus the model residual SD.

Logistic: $P_{is}=\mathrm{logit}^{-1}(\lambda_{is}/\tau)$. The scale tau controls how sharply probabilities switch around $\lambda=0$.

Hard rule: $P_{is} = \mathbb{I}\{\lambda_{is}>0\}$ yields a binary map useful for thresholds and counts (invasibility = \#invaders with $P=1$ at a site).

Examples

## Minimal example (toy shapes)
set.seed(1)
S = 6; I = 4
sites   = paste0("s", 1:S)
inv     = paste0("i", 1:I)
r_is_z  = matrix(rnorm(S*I), S, I, dimnames=list(sites, inv))
C_is_z  = matrix(rnorm(S*I), S, I, dimnames=dimnames(r_is_z))
S_is_z  = matrix(rep(scale(rnorm(S)), each=I), S, I, dimnames=dimnames(r_is_z))
gamma   = setNames(runif(I, 0.5, 1.2), inv)
alpha   = setNames(runif(I, 0.2, 1.0), inv)
beta    = setNames(runif(I, 0.1, 0.6), inv)

# Build lambda internally, then get logistic probabilities
out_logit = compute_establishment_probability(
  r_is_z=r_is_z, C_is_z=C_is_z, S_is_z=S_is_z,
  gamma=gamma, alpha=alpha, beta=beta,
  method="logit", tau=1, return_long=FALSE, make_plots=FALSE
)
str(out_logit$p_is)
#>  num [1:6, 1:4] 0.274 0.57 0.302 0.93 0.509 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:6] "s1" "s2" "s3" "s4" ...
#>   ..$ : chr [1:4] "i1" "i2" "i3" "i4"

# Probit with a scalar sigma
out_probit = compute_establishment_probability(
  r_is_z=r_is_z, C_is_z=C_is_z, S_is_z=S_is_z,
  gamma=gamma, alpha=alpha, beta=beta,
  method="probit", sigma=1
)
# View site-mean probability map (requires ggplot2)
if (requireNamespace("ggplot2", quietly=TRUE)) print(out_probit$plots$site_mean)
#> NULL

# Hard rule (λ>0)
out_hard = compute_establishment_probability(
  r_is_z=r_is_z, C_is_z=C_is_z, S_is_z=S_is_z,
  gamma=gamma, alpha=alpha, beta=beta,
  method="hard", return_long=TRUE, make_plots=FALSE
)
table(out_hard$prob_long$val)  # 0/1 counts
#> 
#>  0  1 
#> 10 14