Climate change vulnerability assessment with ccva

A climate change vulnerability assessment (CCVA) asks which species are most at risk from climate change. Almost every CCVA rests on the same framework:

\[\text{Vulnerability} = f(\text{Exposure}, \text{Sensitivity}, \text{Adaptive capacity}).\]

In practice \(f\) is usually a sum of expert-scored traits with arbitrary weights, never tested for internal consistency, never validated against an observable outcome, and reported without uncertainty. ccva keeps the framework but puts it on a statistical footing:

Every step has a plot. All data in this vignette are simulated and have no biological meaning.

1. Data

ccva_sim holds 300 simulated species generated from a known latent structure: correlated latent exposure, sensitivity and adaptive capacity drive a latent vulnerability \(V^*\). Each component is observed only through its indicators, and vulnerability only through three proxies – an IUCN-like status, a population trend and range change. Because the truth is known, we can check which method recovers it.

as_ccva() is the data contract. It tags each column to a component and records its direction: whether a higher value means more of the component (or, for an outcome, more vulnerable). A declining population trend and a range loss are bad news, so both get direction -1:

d <- as_ccva(ccva_sim,
  exposure    = c("E1_risk", "E2_risk", "E3_risk"),
  sensitivity = grep("^s[0-9]+_", names(ccva_sim), value = TRUE),
  adaptive    = grep("^a[0-9]+_", names(ccva_sim), value = TRUE),
  outcomes    = c("iucn", "pop_trend", "range_change"),
  direction   = c(pop_trend = -1, range_change = -1),
  id = "species")
d$truth <- ccva_sim[c("E_star", "S_star", "A_star", "V_star")]
d
#> <ccva_data> 300 units
#>   Exposure            3 indicators: E1_risk, E2_risk, E3_risk
#>   Sensitivity        11 indicators: s1_thermal, s2_range, s3_climzones, s4_habitat, ... +7
#>   Adaptive capacity   4 indicators: a1_genetic, a2_bodysize, a3_lowfrag, a4_microhab
#>   Outcomes            3 proxys:     iucn, pop_trend (-), range_change (-)
#>   (-) marks indicators where a higher value means less of the component.

(ccva_sim_data() does exactly this in one call.) With real data, the exposure indicators usually come from a species distribution model. ccva_sdm_exposure() turns the cell counts of thresholded current and future habitat maps into Ofori’s three exposure indices, already recoded so that higher means more exposed:

ccva_sdm_exposure(n_current = c(120, 80), n_future = c(60, 95),
                  n_overlap = c(50, 70), n_protected = c(30, 10))
#>   E1_change E2_overlap E3_protected E1_risk   E2_risk E3_risk
#> 1   -0.5000  0.4166667        0.250  0.5000 0.5833333   0.750
#> 2    0.1875  0.8750000        0.125 -0.1875 0.1250000   0.875

ccva_plot(d, "indicators")

Before combining anything: diagnostics

Summing indicators assumes they measure one construct; adding components assumes they are not redundant. ccva_check() tests both.

chk <- ccva_check(d)
chk
#> <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  good
#>   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

ccva_plot(chk, "alpha")

2. The conventional baseline: Ofori et al. (2017)

Ofori and colleagues scored every indicator 1–3, summed the scores within each component, and combined the components. component_method = "ordinal" reproduces the scoring; the weighted method uses the corrected normalisation \(V = (E + S - A + 1)/3\), which maps the full range of \(E + S - A\) onto \([0, 1]\) without clipping.

ofori <- ccva_combine(d, "weighted", component_method = "ordinal")
ofori
#> <ccva_fit> weighted index | 300 units | constraints: 4 rule(s)
#>   components: exposure (ordinal), sensitivity (ordinal), adaptive (ordinal) 
#>   classes:    Low 19 | Moderate 217 | High 64

ccva_plot(ofori, "structure")

ccva_plot(ofori, "components")

The weighted index makes it easy to see what drives each species’ score:

ccva_plot(ofori, "contributions", n = 8)

3. Vulnerability as a latent variable

Vulnerability is never observed directly; we observe its consequences. The latent model treats it that way: \(V^*\) is measured by the three outcome proxies and caused by latent exposure, sensitivity and adaptive capacity, each measured by its own indicators.

The three exposure indices share a denominator (current habitat area), so their errors are correlated. We declare that in shared_error; it is modelled by an orthogonal method factor.

latent <- ccva_combine(d, "latent",
  shared_error = list(c("E1_risk", "E2_risk", "E3_risk")))
latent
#> <ccva_fit> latent SEM | 300 units | constraints: 4 rule(s)
#>   components: exposure (mean), sensitivity (mean), adaptive (mean) 
#>   classes:    Low 39 | Moderate 200 | High 61 
#>   fit: CFI 1.000, TLI 1.000, RMSEA 0.004, SRMR 0.037
#>   effects on V*: exposure +0.66, sensitivity +0.41, adaptive -0.25 (standardised)

ccva_plot(latent, "structure")

The measurement model shows how strongly each indicator reflects its component – a data-driven replacement for expert weights:

ccva_plot(latent, "effects")

Which method recovers the truth?

cmp <- ccva_compare(ofori = ofori, latent = latent)
cmp
#> <ccva_comparison> Spearman correlation of scores with each target
#>  method  iucn pop_trend range_change V* (truth) class_agreement
#>   ofori 0.580     0.584        0.451      0.812            1.00
#>  latent 0.745     0.748        0.481      0.888            0.85
#> 
#> rank correlation between methods
#>        ofori latent
#> ofori  1.000  0.897
#> latent 0.897  1.000

ccva_plot(cmp, "truth")

The latent model recovers the true vulnerability better than the Ofori index, because it uses the outcome proxies to infer \(V^*\) instead of computing it from discretised scores. Rankings change between methods:

ccva_plot(cmp, "ranks")

4. Constraints: imposing ecological knowledge

Free estimation ignores what ecologists already know. In ccva that knowledge is imposed explicitly, at two levels:

predictor level – the sign of each indicator’s loading, bounds, equality between indicators, or a fixed value (0 drops it);
component level – exposure and sensitivity raise vulnerability, adaptive capacity lowers it; bounds on the effects; or adaptive capacity acting only as a moderator.

The default constraints encode the standard ecological directions:

ccva_default_constraints()
#> <ccva_constraints> 4 rule(s)
#>   [component] exposure: positive effect
#>   [component] sensitivity: positive effect
#>   [component] adaptive: negative effect
#>   [predictor] all indicators & outcomes: sign follows declared direction

Rules are added on top of the defaults:

cons <- ccva_constraints(
  constrain_equal(c("a1_genetic", "a2_bodysize")),
  constrain_fixed("s7_repromode", 0),
  constrain_bounds("exposure", lower = 0, upper = 2))
cons
#> <ccva_constraints> 7 rule(s)
#>   [component] exposure: positive effect
#>   [component] sensitivity: positive effect
#>   [component] adaptive: negative effect
#>   [predictor] all indicators & outcomes: sign follows declared direction
#>   [predictor] equal: a1_genetic = a2_bodysize
#>   [predictor] s7_repromode fixed at 0 (dropped)
#>   [component] exposure in [0, 2]
fit_c <- ccva_combine(d, "latent", constraints = cons,
  shared_error = list(c("E1_risk", "E2_risk", "E3_risk")))
fit_c$details$effects
#>     component   estimate        se        std       pvalue
#> 1    exposure  2.0000000 0.0000000  0.6265393           NA
#> 2 sensitivity  1.3047784 0.2307204  0.4087475 1.556320e-08
#> 3    adaptive -0.7932589 0.1926316 -0.2485039 3.821578e-05

When do constraints matter?

With 300 species the free estimates already have the ecological signs, so constraints change little. CCVAs, however, are usually run on a few dozen species. Here we fit the latent model with and without constraints to 20 simulated datasets of 50 species each:

small_n <- do.call(rbind, lapply(1:20, function(i) {
  s <- ccva_simulate(n = 50, seed = 100 + i)
  run <- function(cons) {
    f <- tryCatch(suppressWarnings(ccva_combine(s, "latent",
           constraints = cons,
           shared_error = list(c("E1_risk", "E2_risk", "E3_risk")))),
           error = function(e) NULL)
    if (is.null(f)) return(c(converged = 0, rho = NA))
    c(converged = f$details$converged,
      rho = cor(f$score, s$truth$V_star, method = "spearman"))
  }
  rbind(data.frame(model = "constrained", t(run(NULL))),
        data.frame(model = "free", t(run(ccva_unconstrained()))))
}))
aggregate(cbind(converged, rho) ~ model, small_n, mean, na.action = na.pass)
#>         model converged       rho
#> 1 constrained      0.95 0.9118463
#> 2        free      0.95 0.9118559

Imposing the ecological directions is what lets the latent model converge reliably at realistic sample sizes. Sign constraints on loadings have a second benefit: the sign of a latent factor is otherwise arbitrary, and an unconstrained fit can return vulnerability upside down. With every loading pinned to its ecological direction that cannot happen.

Adaptive capacity as a moderator

Biologically, adaptive capacity does nothing on its own: it modulates the impact of exposure and sensitivity. constrain_moderator() removes its main effect. In the shape-constrained model this fits monotone interaction surfaces \(f_{EA}(E, A) + f_{SA}(S, A)\); in the weighted index it gives \(V \propto E(1 - \gamma_E A) + S(1 - \gamma_S A)\).

mod <- ccva_combine(d, "shape",
                    constraints = ccva_constraints(constrain_moderator()))
ccva_plot(mod, "structure")

Without the moderator constraint, the shape-constrained model learns the shape of each component’s effect from data while the direction is fixed:

shape <- ccva_combine(d, "shape")
ccva_plot(shape, "effects")

5. High-dimensional predictors

Some components arrive with more predictors than species – genomic markers for adaptive capacity are the obvious case. The penalised two-stage method estimates latent vulnerability from the outcomes, then regresses it on every predictor with a LASSO whose coefficient signs follow the constraints.

g <- ccva_simulate(n = 150, n_genomic = 300, n_causal = 10, seed = 3)
length(unlist(g$indicators))   # predictors
#> [1] 318
pen <- ccva_combine(g, "penalised")
pen
#> <ccva_fit> penalised two-stage | 150 units | constraints: 4 rule(s)
#>   components: exposure (mean), sensitivity (mean), adaptive (mean) 
#>   classes:    Low 19 | Moderate 90 | High 41 
#>   selected predictors: 23 of 318

ccva_plot(pen, "effects")

6. How certain is the ranking?

A sorted list of point estimates implies a precision the data rarely support. ccva_uncertainty() refits the whole pipeline on bootstrap resamples; for the latent model each species’ factor-score error is added as well.

unc <- ccva_uncertainty(latent, "bootstrap", R = 40, seed = 1,
                        progress = FALSE)
unc
#> <ccva_uncertainty> bootstrap (sampling), 40 replicates, 95% intervals
#>   mean interval width: 0.302 (score range 1.000)
#>   pairs with overlapping intervals: 71.1%
#>   robust pairwise rankings (P > 0.95): 51.8%
#>   most vulnerable:
#>     Species_241      1.000 [0.884, 1.209]  rank 1 (1-22)
#>     Species_280      0.989 [0.865, 1.144]  rank 2 (1-16)
#>     Species_013      0.953 [0.839, 1.151]  rank 3 (1-19)
#>     Species_287      0.935 [0.776, 1.047]  rank 4 (1-22)
#>     Species_022      0.921 [0.796, 1.043]  rank 5 (1-27)

ccva_plot(unc, "intervals", n = 25)

For the weighted index, a second question is how much the ranking depends on the arbitrary choice of weights:

ccva_uncertainty(ofori, "weights", R = 1000, seed = 1)
#> <ccva_uncertainty> random weights (weighting), 1000 replicates, 95% intervals
#>   mean interval width: 0.365 (score range 0.646)
#>   pairs with overlapping intervals: 87.1%
#>   robust pairwise rankings (P > 0.95): 48.1%
#>   most vulnerable:
#>     Species_013      0.889 [0.712, 0.996]  rank 1 (1-1)
#>     Species_250      0.889 [0.712, 0.996]  rank 1 (1-1)
#>     Species_280      0.879 [0.710, 0.985]  rank 3 (3-3)
#>     Species_112      0.859 [0.704, 0.973]  rank 4 (4-12)
#>     Species_171      0.859 [0.704, 0.973]  rank 4 (4-12)

7. Classification and outputs

Cut-offs at 1/3 and 2/3 assume uniformly distributed scores. Alternatives calibrate them to an observed outcome or to natural clusters:

sapply(c("uniform", "quantile", "mixture", "roc"), function(m)
  table(ccva_classify(latent$score, m, data = d)))
#>          uniform quantile mixture roc
#> Low           39      100      65  99
#> Moderate     200      100     199  58
#> High          61      100      36 143

ccva_plot(latent, "ranking", n = 8)

head(ccva_rank(latent, unc))
#>            unit     score rank class  exposure sensitivity  adaptive     lower
#> 241 Species_241 1.0000000    1  High 0.8919842   0.6868433 0.3320612 0.8835635
#> 280 Species_280 0.9890865    2  High 0.9038663   0.7170390 0.2866167 0.8646315
#> 13  Species_013 0.9532944    3  High 0.7927621   0.7350135 0.2202469 0.8391743
#> 287 Species_287 0.9351071    4  High 0.8846028   0.5437302 0.2262880 0.7762848
#> 22  Species_022 0.9214363    5  High 0.6583220   0.6789288 0.2349839 0.7960000
#> 159 Species_159 0.9196621    6  High 0.8824058   0.6684904 0.5109389 0.7789079
#>        upper rank_low rank_high
#> 241 1.208669    1.000    22.250
#> 280 1.143911    1.000    16.000
#> 13  1.150911    1.000    19.025
#> 287 1.046994    1.000    22.400
#> 22  1.042746    1.000    27.025
#> 159 1.110225    1.975    27.075

ccva_report(latent, "report.html", uncertainty = unc, comparison = cmp) writes all of this into one reproducible HTML report.