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.

library(ccva)

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.

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). With real data, calibrate against an outcome the model did not use.

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.578 0.491 0.615      0.214    1.000     5.000
#>  site_06        29              0.567 0.470 0.613      0.241    1.000     7.000
#>  site_13        28              0.555 0.450 0.611      0.179    1.000    12.025
#>  site_02        28              0.549 0.451 0.608      0.214    2.000    10.025
#>  site_09        28              0.547 0.442 0.602      0.214    3.000    13.000
#>  site_07        29              0.547 0.458 0.591      0.172    3.000    14.025
#>  site_15        30              0.546 0.454 0.587      0.267    2.975    12.000
#>  site_08        28              0.546 0.443 0.594      0.143    1.975    12.025
#>  site_11        29              0.542 0.437 0.589      0.103    2.975    13.000
#>  site_05        28              0.537 0.444 0.583      0.250    3.000    13.000
#>   ... 5 more sites
ccva_plot(sites)

What is not (yet) in the package