Skip to contents

Auxiliary GLMM for trait-varying (and optionally site-varying) slopes

Source: R/fit_auxiliary_residents_glmm.R
fit_auxiliary_residents_glmm.Rd

fit_auxiliary_residents_glmm() assembles a long table of resident (site, species) cells with standardized predictors \(r^{(z)}_{js}\), \(C^{(z)}_{js}\), \(S^{(z)}_{js}\) and 2D trait scores (tr1, tr2), then fits a Gaussian GLMM to estimate how slopes on abiotic suitability, niche crowding, and site saturation vary across the trait plane. You can optionally include site-level random slopes for r_z and C_z. (Do not include (0 + S_z || site), as S_z is site-only and has no within-site variation across residents.)

This version is forgiving: it coerces non-matrix inputs with as.matrix(), repairs missing dimnames where possible, and aligns inputs to comm_res.

Usage

fit_auxiliary_residents_glmm(
  comm_res,
  r_js_z,
  C_js_z,
  S_js_z,
  Q_res,
  use_site_random_slopes = TRUE,
  family = gaussian(),
  control = glmmTMB::glmmTMBControl(optCtrl = list(iter.max = 1e+05, eval.max = 1e+05)),
  na_action = c("drop", "error"),
  verbose = TRUE
)

Arguments

comm_res

Numeric matrix of site × resident abundances; rownames = site IDs, colnames = resident IDs (data.frame/tibble will be coerced with as.matrix()).

r_js_z

Numeric matrix of site × resident standardized abiotic suitability (row-wise z by site); same dimnames as comm_res (coerced with as.matrix() if needed).

C_js_z

Numeric matrix of site × resident standardized crowding (row-wise z by site); same dimnames as comm_res (coerced with as.matrix() if needed).

S_js_z

Numeric matrix of site × resident standardized site saturation (broadcast site-only); same dimnames as comm_res (coerced with as.matrix() if needed).

Q_res

Data frame with resident trait scores; must contain columns tr1 and tr2, and rownames matching colnames(comm_res). If rownames are missing but a column named "species" exists, it is used for rownames.

use_site_random_slopes

Logical; if TRUE, add (0 + r_z || site) and (0 + C_z || site) to allow site-varying slopes. Default TRUE.

family

GLM family for the conditional model. Default gaussian() for log1p(abundance); keep as-is for slope learning.

control

Control list passed to glmmTMB::glmmTMB(). Default increases optimizer limits for stability.

na_action

How to handle rows with any missing values among abundance, r_z, C_z, S_z, tr1, tr2. One of "drop" (default) or "error".

verbose

Logical; if TRUE, prints a compact model summary line.

Value

A named list with:

  • fit — the glmmTMB fitted model (auxiliary residents GLMM).

  • data — the long table used for fitting (one row per site×resident).

  • formula — the exact formula used.

  • args — list of key arguments resolved inside the function.

Details

Fit an auxiliary residents-only GLMM on standardized predictors

The fixed part of the model is $$\log(1+\mathrm{abundance}) \sim (r_z + C_z + S_z)\times(tr1 + tr2),$$ with random intercepts (1 | species) + (1 | site) always included. If use_site_random_slopes = TRUE, the random-slope terms (0 + r_z || site) + (0 + C_z || site) are added. We intentionally omit (0 + S_z || site) because S_z is site-only, therefore constant within a site and not estimable as a random slope.

Why this matters for invasion fitness

Fitting this auxiliary model yields trait-varying slope systems for r_z, C_z, and S_z, which you can map to \(\theta_i\), \(\alpha_i\), and \(\beta_i\) (and optionally site-varying versions via random slopes). These parameters plug directly into the linear invasion-fitness decomposition \(\lambda_{is} = \Gamma_{is}\, r^{(z)}_{is} - \alpha_{is}\, C^{(z)}_{is} - \beta_i\, S^{(z)}_{is}\), linking trait position (trait space, convex hull, centroid) to how abiotic suitability, niche crowding, and resident competition shape establishment.

Examples

if (FALSE) { # \dontrun{
# Toy dimensions
set.seed(1)
S = 8  # sites
J = 6  # residents
sites   = paste0("s", 1:S)
res_ids = paste0("sp", 1:J)

# Site × resident abundance
comm_res = matrix(rpois(S*J, lambda = 2), S, J,
                   dimnames = list(sites, res_ids))

# Standardized predictors (row-wise z by site for r and C; site-only for S)
r_raw = matrix(rnorm(S*J), S, J, dimnames = dimnames(comm_res))
C_raw = matrix(rnorm(S*J), S, J, dimnames = dimnames(comm_res))
S_raw = matrix(rep(scale(rnorm(S))[,1], each = J), S, J,
                dimnames = dimnames(comm_res))
# Simple per-site z for demo
r_js_z = t(scale(t(r_raw))); r_js_z[!is.finite(r_js_z)] = 0
C_js_z = t(scale(t(C_raw))); C_js_z[!is.finite(C_js_z)] = 0
S_js_z = S_raw

# Resident trait-plane scores (PCoA on Gower in real workflow)
Q_res = data.frame(tr1 = rnorm(J), tr2 = rnorm(J))
rownames(Q_res) = res_ids

aux = fit_auxiliary_residents_glmm(
  comm_res   = comm_res,
  r_js_z     = r_js_z,
  C_js_z     = C_js_z,
  S_js_z     = S_js_z,
  Q_res      = Q_res,
  use_site_random_slopes = TRUE
)

summary(aux$fit)
head(aux$data)
aux$formula
} # }