A statistical CCVA workflow, step by step

This vignette walks through a statistically grounded CCVA workflow, one recommendation at a time: match the method to the data, control trait measurements for methodology, test how well the data measure each component, make causal hypotheses explicit, impose ecologically justified shapes, turn rankings into calibrated risk, and aggregate to sites for spatial prioritisation. It applies to any taxon. All data are simulated.

1. Match the method to the data

A sensible default is shape-constrained parametric models for small samples, regularised regression for regional compilations, and flexible non-parametric methods for global databases. ccva_recommend() reads the number of species and of predictors per component and gives concrete steps:

ccva_recommend(ccva_simulate(n = 40, seed = 1))
#> <ccva_recommendation> small sample (< 50 units) 
#>   predictors: exposure 3 (p/n 0.07), sensitivity 11 (p/n 0.28), adaptive 4 (p/n 0.10) 
#>   outcomes:   iucn, pop_trend, range_change 
#> 
#>   1. [diagnose] Check internal consistency of each component (Cronbach's alpha, mean inter-item r); declare complementary-facet components as formative
#>       ccva_check(x) 
#>   2. [vulnerability] Small sample: shape-constrained parametric model with ecological sign constraints; weighted index as baseline
#>       ccva_combine(x, "shape") 
#>   3. [dependencies] Model dependencies among components you hypothesise (e.g. adaptive capacity altering exposure via range shifts, or sensitivity via plasticity)
#>       constraints = ccva_constraints(constrain_link("adaptive", "exposure")) 
#>   4. [uncertainty] Propagate uncertainty to scores and ranks; report intervals, not a bare ranking
#>       ccva_uncertainty(fit, "bootstrap") 
#>   5. [calibration] Turn scores into calibrated risk (probability of an observed outcome) with intervals
#>       ccva_calibrate(fit, uncertainty = u) 
#>   6. [caution] With fewer than 50 units expect wide intervals and few robust rankings; constraints are what keep latent fits stable

A regional compilation where genomic data make adaptive capacity high-dimensional gets a different route:

regional <- ccva_simulate(n = 120, n_genomic = 250, seed = 2)
ccva_recommend(regional)
#> <ccva_recommendation> regional compilation (50-149 units) 
#>   predictors: exposure 3 (p/n 0.03), sensitivity 11 (p/n 0.09), adaptive 254 (p/n 2.12) 
#>   outcomes:   iucn, pop_trend, range_change 
#> 
#>   1. [diagnose] Check internal consistency of each component (Cronbach's alpha, mean inter-item r); declare complementary-facet components as formative
#>       ccva_check(x) 
#>   2. [components] Predictors outnumber or approach units for adaptive: use regularised (penalised) scoring
#>       component_method = list(adaptive = "penalised") 
#>   3. [vulnerability] Regional compilation with many predictors: regularised two-stage model
#>       ccva_combine(x, "penalised") 
#>   4. [dependencies] Model dependencies among components you hypothesise (e.g. adaptive capacity altering exposure via range shifts, or sensitivity via plasticity)
#>       constraints = ccva_constraints(constrain_link("adaptive", "exposure")) 
#>   5. [uncertainty] Propagate uncertainty to scores and ranks; report intervals, not a bare ranking
#>       ccva_uncertainty(fit, "bootstrap") 
#>   6. [calibration] Turn scores into calibrated risk (probability of an observed outcome) with intervals
#>       ccva_calibrate(fit, uncertainty = u)

2. Control sensitivity traits for methodology

Heat tolerance depends on how it was measured: static or dynamic assays, ramping rate, starting temperature. Before comparing species, the methodological effect should be removed. ccva_adjust_method() estimates it within species and adjusts every measurement to a common reference protocol:

assays <- data.frame(species = rep(sprintf("sp_%02d", 1:25), each = 3),
                     ramp = rep(c(0.1, 0.25, 0.5), 25),
                     assay = rep(c("dynamic", "dynamic", "static"), 25))
true_ctmax <- rnorm(25, 44, 2)
assays$ctmax <- rep(true_ctmax, each = 3) + 3 * assays$ramp -
  1.2 * (assays$assay == "static") + rnorm(75, 0, 0.3)

adj <- ccva_adjust_method(assays, "ctmax", c("ramp", "assay"),
                          unit = "species")
adj
#> <ccva_adjustment> ctmax adjusted to reference protocol (ramp = 0.25, assay = dynamic)
#>   ramp                     +3.518 (se 0.485, p 3.08e-09)
#>   assaystatic              -1.402 (se 0.170, p 9.27e-11)
#>   25 units; use $by_unit as an indicator in as_ccva()
head(adj$by_unit)
#>    unit ctmax_adj    ctmax n
#> 1 sp_01  44.90867 44.55870 3
#> 2 sp_02  44.62003 44.27005 3
#> 3 sp_03  45.77234 45.42236 3
#> 4 sp_04  42.89313 42.54316 3
#> 5 sp_05  45.29218 44.94221 3
#> 6 sp_06  48.17168 47.82170 3

adj$by_unit$ctmax_adj is then a sensitivity indicator for as_ccva(). Warming tolerance and thermal safety margins against measured or modelled microclimates come from ccva_thermal_margin(); sublethal performance traits, biotic-interaction metrics and effective population size enter as indicators in the same way.

3. Test how well the data measure E, S and A

Cronbach’s alpha checks whether a component’s indicators covary like measures of one construct. It is only expected to be high when they are interchangeable measures (reflective). When they are complementary facets – dispersal ability and genetic diversity both contribute to adaptive capacity, through distinct mechanisms – a low alpha is expected and is not a problem. Declare such components formative:

d <- ccva_sim_data()
d$formative <- "adaptive"
ccva_check(d)
#> <ccva_check> internal consistency (Cronbach's alpha)
#>   exposure       3 indicators  alpha 0.868  (std 0.868)  mean r 0.69  good
#>   sensitivity   11 indicators  alpha 0.871  (std 0.871)  mean r 0.38  good
#>   adaptive       4 indicators  alpha 0.801  (std 0.802)  mean r 0.50  formative: alpha not required
#>   guide (reflective components): >= 0.8 good, >= 0.7 acceptable, < 0.6 questionable;
#>   alpha rises with the number of indicators, so compare mean r too
#> 
#> cross-component Spearman correlations
#>             exposure sensitivity adaptive
#> exposure       1.000       0.269   -0.157
#> sensitivity    0.269       1.000   -0.172
#> adaptive      -0.157      -0.172    1.000
#> 
#> no indicator pairs with |r| > 0.90

Alpha rises with the number of indicators, so the mean inter-item correlation is reported alongside it. In the latent model a formative component enters as a composite of its indicators rather than as a reflective factor. With formative_weights = "estimated" every indicator of a formative component gets its own sign-constrained weight on latent vulnerability (a MIMIC model), so the data decide which facets matter:

fw <- ccva_combine(d, "latent", formative_weights = "estimated",
                   shared_error = list(c("E1_risk", "E2_risk", "E3_risk")))
fw$details$indicator_weights
#>   component   indicator   estimate        se         std    pvalue
#> 1  adaptive  a1_genetic -0.5166512 0.4218689 -0.11451676 0.2206987
#> 2  adaptive a2_bodysize -0.5953934 0.4551231 -0.13197012 0.1908045
#> 3  adaptive  a3_lowfrag -0.2918731 0.3076811 -0.06469426 0.3428128
#> 4  adaptive a4_microhab -0.1013481 0.2828641 -0.02246402 0.7201243

Incomplete trait data

Real trait databases have gaps, and the species without a Red List category are often those that most need an assessment. Build the data with na_action = "keep": the latent model then uses full-information maximum likelihood, and units with a missing outcome are scored from their indicators.

dm <- ccva_sim
dm$s1_thermal[sample(300, 40)] <- NA
dm$iucn[1:30] <- NA                       # e.g. Data Deficient species
xm <- as_ccva(dm, exposure = c("E1_risk", "E2_risk", "E3_risk"),
  sensitivity = grep("^s[0-9]+_", names(dm), value = TRUE),
  adaptive = grep("^a[0-9]+_", names(dm), value = TRUE),
  outcomes = c("iucn", "pop_trend", "range_change"),
  direction = c(pop_trend = -1, range_change = -1), id = "species",
  na_action = "keep")
fm <- ccva_combine(xm, "latent",
                   shared_error = list(c("E1_risk", "E2_risk", "E3_risk")))
cor(fm$score[1:30], ccva_sim$V_star[1:30], method = "spearman")
#> [1] 0.8117909

Rule-based assessments

Many published assessments use the rules of Foden et al. (2013): a species is highly vulnerable when at least one trait flags high sensitivity, at least one flags low adaptive capacity and its exposure is in the top quarter. method = "rule" reproduces them and shows how many species each trait flags on its own – often a large share of the classification rests on one trait and one threshold:

r <- ccva_combine(d, "rule")
r
#> <ccva_fit> rule-based (Foden et al. 2013) | 300 units | constraints: 4 rule(s)
#>   components: exposure (mean), sensitivity (mean), adaptive (mean) 
#>   classes:    Low 116 | Moderate 130 | High 54 
#>   categories: highly vulnerable 54 | high latent risk 73 | potential adapter 48 | potential persister 9 | lower risk 116 
#>   flagged:    high exposure 121 | high sensitivity 237 | low adaptive capacity 156 (unknowns treated as optimistic)
head(r$details$drivers[order(-r$details$drivers$n_only), ])
#>    component   indicator     type        cut n_flagged n_only n_unknown
#> 17  adaptive  a3_lowfrag quantile -0.5420646        75     25         0
#> 3   exposure     E3_risk quantile  0.6846644        75     19         0
#> 15  adaptive  a1_genetic quantile -0.6065618        75     19         0
#> 1   exposure     E1_risk quantile  0.6216212        75     18         0
#> 16  adaptive a2_bodysize quantile -0.6052285        75     18         0
#> 2   exposure     E2_risk quantile  0.6247613        75     15         0

4. Make the causal hypotheses explicit

A common causal hypothesis is that adaptive capacity acts on vulnerability partly through the other components: range shifts change the climate a species is exposed to, and plasticity or evolution of heat tolerance changes its sensitivity. Each arrow is a constrain_link(). The simulated data below contain such a dependency (species with more capacity to move face less exposure):

s <- ccva_simulate(n = 400, link_AE = -0.6, seed = 4)
causal <- ccva_combine(s, "latent",
  shared_error = list(c("E1_risk", "E2_risk", "E3_risk")),
  constraints = ccva_constraints(
    constrain_link("adaptive", "exposure", "-"),
    constrain_link("adaptive", "sensitivity", "-")))
causal
#> <ccva_fit> latent SEM | 400 units | constraints: 6 rule(s)
#>   components: exposure (mean), sensitivity (mean), adaptive (mean) 
#>   classes:    Low 59 | Moderate 268 | High 73 
#>   fit: CFI 0.994, TLI 0.994, RMSEA 0.015, SRMR 0.034
#>   effects on V*: exposure +0.49, sensitivity +0.43, adaptive -0.36 (standardised)
#>   link adaptive -> exposure: -0.64 (standardised)
#>   link adaptive -> sensitivity: -0.22 (standardised)
#>   total effect of adaptive on V*: -0.77 = direct -0.36 + indirect -0.41 (via exposure -0.31; via sensitivity -0.09)

ccva_plot(causal, "structure")

The structural equation model estimates every arrow at once and splits the effect of adaptive capacity into a direct part and indirect parts through exposure and sensitivity. Ignoring the dependency would count part of adaptive capacity twice – once in its own score and once inside exposure.

The same links are resolved in two stages by the other models: the part of the target component explained by the source is removed before combining.

ccva_combine(s, "weighted", constraints = ccva_constraints(
  constrain_link("adaptive", "exposure")))
#> <ccva_fit> weighted index | 400 units | constraints: 5 rule(s)
#>   components: exposure (mean), sensitivity (mean), adaptive (mean) 
#>   classes:    Low 27 | Moderate 369 | High 4 
#>   link adaptive -> exposure: two-stage, 23% of exposure explained by adaptive removed
#>   note: two-stage link: exposure adjusted for adaptive (R2 = 0.23)

5. Shapes that ecology justifies

Shape-constrained models keep the data in charge of the details while excluding implausible relationships. Beyond direction, curvature can be imposed: an exposure effect that saturates (concave) and a sensitivity effect that accelerates (convex).

nl <- ccva_simulate(n = 300, nonlinear = TRUE, seed = 7)
curved <- ccva_combine(nl, "shape", constraints = ccva_constraints(
  constrain_shape("exposure", "concave"),
  constrain_shape("sensitivity", "convex")))
ccva_plot(curved, "effects")

Dropping all constraints gives a flexible smooth model, appropriate for large databases: ccva_combine(x, "shape", constraints = ccva_unconstrained()).

6. From rankings to calibrated risk

A ranking says which species are more vulnerable than others; managers also need to know how likely a decline is. ccva_calibrate() turns vulnerability scores into the probability of an observed outcome – here IUCN-like status of Vulnerable or worse – and, given bootstrap draws, carries the uncertainty of the scores into the risk:

fit <- ccva_combine(d, "latent",
  shared_error = list(c("E1_risk", "E2_risk", "E3_risk")))
unc <- ccva_uncertainty(fit, "bootstrap", R = 40, seed = 1,
                        progress = FALSE)
risk <- ccva_calibrate(fit, outcome = "iucn", uncertainty = unc,
                       R = 200, seed = 1)
risk
#> <ccva_calibration> P(iucn >= 2) from latent vulnerability
#>   AUC 0.882 | Brier 0.137 | event rate 0.42 | 200 bootstrap fits (score uncertainty propagated)
#>   note: this outcome was used to fit the model; calibration is in-sample
#>   highest risk (95% intervals):
#>     Species_241      1.00 [0.93, 1.00]
#>     Species_280      1.00 [0.93, 1.00]
#>     Species_013      0.99 [0.92, 1.00]
#>     Species_287      0.99 [0.89, 1.00]
#>     Species_159      0.99 [0.85, 1.00]

ccva_plot(risk, "curve")

ccva_plot(risk, "risk", n = 20)

Here the outcome was also used to fit the latent model, so the calibration is in-sample (the print-out says so). ccva_cv() gives honest scores: each fold is scored by a model fitted without it, with the outcomes of the held-out species hidden, as for a species that has never been assessed. Passing them to ccva_calibrate() makes the calibration out of sample:

cv <- ccva_cv(fit, k = 5, seed = 1, progress = FALSE)
cv
#> <ccva_cv> 5-fold cross-validation of the latent model
#>   hidden for held-out units: iucn, pop_trend, range_change 
#>   Spearman correlation with each hidden outcome:
#>       outcome in_sample cross_validated   n
#>          iucn     0.745           0.609 300
#>     pop_trend     0.748           0.607 300
#>  range_change     0.480           0.435 300
ccva_calibrate(fit, outcome = "iucn", cv = cv, R = 200, seed = 1)
#> <ccva_calibration> P(iucn >= 2) from latent vulnerability
#>   AUC 0.804 | Brier 0.176 | event rate 0.42 | 200 bootstrap fits
#>   scores: cross-validated (outcome hidden for held-out units)
#>   highest risk (95% intervals):
#>     Species_280      1.00 [0.99, 1.00]
#>     Species_013      1.00 [0.99, 1.00]
#>     Species_159      0.99 [0.97, 1.00]
#>     Species_030      0.98 [0.96, 1.00]
#>     Species_222      0.98 [0.96, 1.00]

Related species are not independent replicates. Give a grouping (family, genus, river basin) as cluster to ccva_cv(), so that whole groups are held out together, and to ccva_uncertainty(), so that the bootstrap resamples groups rather than species and the intervals respect the dependence within them.

7. Spatial prioritisation

Conservation is often planned for places and communities rather than single species, especially for cryptic or poorly known taxa. ccva_aggregate() combines species-level results with occurrence records into community-level vulnerability per site, with intervals:

occ <- unique(data.frame(site = rep(sprintf("site_%02d", 1:15), each = 30),
                         species = sample(names(fit$score), 450,
                                          replace = TRUE)))
sites <- ccva_aggregate(fit, occ, uncertainty = unc)
sites
#> <ccva_sites> 15 sites, 430 species occurrences
#>     site n_species mean_vulnerability lower upper share_high rank_low rank_high
#>  site_03        28              0.562 0.472 0.597      0.214    1.000     8.025
#>  site_02        28              0.559 0.462 0.621      0.214    1.000     9.025
#>  site_13        28              0.556 0.454 0.614      0.179    1.000    10.100
#>  site_05        28              0.553 0.462 0.596      0.250    1.000    10.075
#>  site_07        29              0.552 0.464 0.594      0.207    1.000    13.000
#>  site_10        28              0.550 0.458 0.605      0.214    1.000    12.025
#>  site_06        29              0.549 0.443 0.594      0.207    1.000    10.050
#>  site_08        28              0.543 0.443 0.593      0.143    1.000    11.000
#>  site_09        28              0.542 0.435 0.596      0.179    1.975    13.000
#>  site_11        29              0.542 0.437 0.589      0.103    2.000    12.025
#>   ... 5 more sites

ccva_plot(sites)