7.5. Ordinal Fragility Functions

postprocessor.calculate_ordinal_fragility(imls, edps, damage_thresholds, intensities=np.round(np.geomspace(0.05, 10.0, 50), 3), dispersion_type='constant')[source]

Fits a hierarchical (cumulative link) probit model to estimate fragility curves for ordered damage states.

This function implements the hierarchical fragility modelling framework for ordinal damage levels, in which the probability of equalling or exceeding damage state \(i\) given intensity measure level IM is expressed as a lognormal CDF:

\[P(\mathrm{DS} \geq ds_i \mid \mathrm{IM}) = \Phi\!\left(\frac{\ln(\mathrm{IM}) - \ln(\theta_i)}{\beta_i}\right), \quad i = 1, \ldots, k\]

where \(\theta_i\) is the median IM capacity for damage state \(i\) and \(\beta_i\) is the associated dispersion. The hierarchical property—i.e., that curves do not cross at the median—is maintained through the constraint \(\theta_1 \leq \theta_2 \leq \cdots \leq \theta_k\).

Two dispersion models are supported via dispersion_type:

Constant dispersion (dispersion_type='constant', default)

All damage states share a single slope \(\beta\), i.e. \(\beta_i = \beta\) for all \(i\). This is the “fully ordered” special case. The model is parameterised by \(k\) ordered intercepts \(\alpha_i\) and one common slope:

\[P(\mathrm{DS} \geq ds_i \mid \mathrm{IM}) = \Phi(\alpha_i - \beta\,\ln \mathrm{IM}), \quad \alpha_1 \leq \cdots \leq \alpha_k\]

The shared slope is enforced by fitting a single statsmodels.OrderedModel (probit distribution) to all damage state observations simultaneously. This prevents fragility curve crossings everywhere along the IM axis.

Variable dispersion (dispersion_type='variable')

Each damage state \(i\) has its own dispersion \(\beta_i\), making the model more realistic when damage states differ in their sensitivity to ground-motion intensity. \(k\) independent binary probit regressions are fitted, one per damage state, yielding damage-state-specific medians \(\theta_i\) and dispersions \(\beta_i\). To restore the hierarchical ordering constraint, isotonic regression (Pool-Adjacent Violators Algorithm) is applied to \(\ln(\theta_i)\) after fitting, ensuring \(\theta_1 \leq \cdots \leq \theta_k\). Curves may cross in the tails (beyond the median), but the ordering of 50th-percentile capacities is guaranteed.

Parameters:

  • imls (array-like) – Intensity measure levels corresponding to the observed EDPs.

  • edps (array-like) – Engineering Demand Parameters (EDPs) representing structural responses.

  • damage_thresholds (array-like) – Damage state thresholds for classifying exceedance levels.

  • intensities (array-like, optional) – Intensity measure levels for which fragility curves are evaluated (default: np.geomspace(0.05, 10.0, 50)).

  • dispersion_type ({'constant', 'variable'}, optional) – Controls whether a single dispersion parameter is shared across all damage states (‘constant’, default) or each damage state receives its own dispersion estimated independently (‘variable’).

Returns:

poes – A 2D array of exceedance probabilities (CDF values) for each intensity level.

  • dispersion_type='constant': shape (len(intensities), n_categories) where n_categories is the number of distinct damage states observed in the data (typically len(damage_thresholds) + 1 when the no-damage state DS=0 is present).

  • dispersion_type='variable': shape (len(intensities), len(damage_thresholds)).

Return type:

numpy.ndarray

References

1) Jalayer, F., Ebrahimian, H., Trevlopoulos, K., and Bradley, B. (2023). Empirical tsunami fragility modelling for hierarchical damage levels. Natural Hazards and Earth System Sciences, 23(2), 909-931. https://doi.org/10.5194/nhess-23-909-2023

2) Nguyen, M. and Lallemant, D. (2022). Order Matters: The Benefits of Ordinal Fragility Curves for Damage and Loss Estimation. Risk Analysis, 42: 1136-1148. https://doi.org/10.1111/risa.13815

3) Lallemant, D., Kiremidjian, A., and Burton, H. (2015). Statistical procedures for developing earthquake damage fragility curves. Earthquake Engng Struct. Dyn., 44, 1373-1389. doi: 10.1002/eqe.2522.

Theoretical Background

Hierarchical fragility model

The hierarchical (cumulative link) approach treats the probability of equalling or exceeding damage state \(i\) given intensity measure level IM as a lognormal CDF:

\[P(\mathrm{DS} \geq ds_i \mid \mathrm{IM}) = \Phi\!\left(\frac{\ln(\mathrm{IM}) - \ln(\theta_i)}{\beta_i}\right), \quad i = 1, \ldots, k\]

where \(\theta_i\) is the median IM capacity for damage state \(i\) and \(\beta_i\) is its associated dispersion. All damage states are fitted simultaneously (or with an ordering constraint applied post-fit), exploiting the ordinal structure of the damage scale. The hierarchical constraint

\[\theta_1 \leq \theta_2 \leq \cdots \leq \theta_k\]

ensures that the 50th-percentile capacities are ordered, i.e. higher damage states are first exceeded at larger intensity levels (Jalayer et al., 2023; Nguyen & Lallemant, 2022).

Constant dispersion (dispersion_type='constant', default)

When a single dispersion \(\beta_i = \beta\) is assumed for all damage states — the fully ordered special case — the model reduces to the standard ordered-probit parameterisation with \(k\) ordered intercepts \(\alpha_i\) and one common slope:

\[P(\mathrm{DS} \geq ds_i \mid \mathrm{IM}) = \Phi(\alpha_i - \beta\,\ln \mathrm{IM}), \quad \alpha_1 \leq \cdots \leq \alpha_k\]

The shared slope is estimated by fitting a single statsmodels.OrderedModel (probit distribution) to all damage-state observations simultaneously. Because every curve has identical shape and only the threshold differs, crossings are strictly impossible everywhere along the IM axis.

Variable dispersion (dispersion_type='variable')

In general, each damage state may exhibit its own sensitivity to ground-motion intensity, yielding damage-state-specific dispersions \(\beta_i\). This is the more realistic formulation and does not sacrifice the hierarchical property. Here \(k\) independent binary probit regressions are fitted, one per damage state:

\[P(\mathrm{DS} \geq ds_i \mid \mathrm{IM}) = \Phi(b_{0,i} + b_{1,i}\,\ln \mathrm{IM})\]

Median and dispersion for each DS are recovered as \(\theta_i = \exp(-b_{0,i}/b_{1,i})\) and \(\beta_i = 1/b_{1,i}\). After fitting, isotonic regression (Pool-Adjacent Violators Algorithm, PAVA) is applied to \(\ln(\theta_i)\) to restore the hierarchical ordering constraint. Curves may cross in the tails but the ordering of 50th-percentile capacities is guaranteed.

Fragility curves

The probability of being in exactly damage state \(i\) is:

\[P(D = i \mid \mathrm{IM}) = P(D \geq i \mid \mathrm{IM}) - P(D \geq i+1 \mid \mathrm{IM})\]

The exceedance fragility curves \(P(D \geq i \mid \mathrm{IM})\) are obtained directly from the fitted cumulative model.

References

  1. Jalayer, F., Ebrahimian, H., Trevlopoulos, K., and Bradley, B. (2023). Empirical tsunami fragility modelling for hierarchical damage levels. Natural Hazards and Earth System Sciences, 23(2), 909–931. https://doi.org/10.5194/nhess-23-909-2023

  2. Nguyen, M. and Lallemant, D. (2022). Order Matters: The Benefits of Ordinal Fragility Curves for Damage and Loss Estimation. Risk Analysis, 42: 1136–1148. https://doi.org/10.1111/risa.13815

  3. Lallemant, D., Kiremidjian, A., and Burton, H. (2015). Statistical procedures for developing earthquake damage fragility curves. Earthquake Engineering & Structural Dynamics, 44, 1373–1389. https://doi.org/10.1002/eqe.2522

Example

import numpy as np
from openquake.vmtk.postprocessor import postprocessor

pp = postprocessor()
intensities = np.geomspace(0.05, 3.0, 50)
damage_thresholds = [0.005, 0.015, 0.040, 0.080]

# Constant dispersion — fully-ordered special case (default)
# shape: (50, n_categories), where n_categories is the number of
# distinct damage states observed (typically n_DS + 1 when DS=0
# is present in the data)
poes_const = pp.calculate_ordinal_fragility(
    imls=imls,
    edps=edps,
    damage_thresholds=damage_thresholds,
    intensities=intensities,
    dispersion_type='constant',
)

# Variable dispersion — general hierarchical case
# shape: (50, len(damage_thresholds)) = (50, 4)
poes_var = pp.calculate_ordinal_fragility(
    imls=imls,
    edps=edps,
    damage_thresholds=damage_thresholds,
    intensities=intensities,
    dispersion_type='variable',
)