28. Hypothesis tests and paired data

Language: R (Quarto) - Python version

Two surfaces that had no R or Python binding until this release: the twelve Numerics HypothesisTests statics, reached directly through hypothesis_test() or, for a flood-frequency record, through the DataFrame facades analysis_data_hypothesis_test()/analysis_data_statistics(); and the RMC.BestFit “Paired Data” subsystem – an x-y curve container with interpolation, area, three curve-simplification algorithms, and an uncertain twin whose y-coordinate is a whole distribution per point. This page works through both, on a real flood record and a real reservoir curve.

What you’ll learn

  • What four hypothesis tests answer about a flood-frequency record – normality, independence, and homogeneity – and how to read a p-value from each.
  • The DataFrame path through analysis_data(): the same tests over a censored-data-aware container, the use_log10 switch, and the two-sample split_index split.
  • correlation()’s matrix form over more than two series at once.
  • A reservoir stage-storage curve interpolated in linear and log space, and its trapezoidal area.
  • How the three curve-simplification algorithms differ on the same curve, including one that drops the curve’s last point without warning.
  • Sampling an uncertain curve at its mean versus its median.

Setup

library(corehydror)

A flood-frequency record

Sixty-nine years of annual peak discharge (m3/s) for the Harricana River, Quebec (Bobee & Ashkar 1991, Multivariate Flood Frequency Analysis, Table 1.2):

harricana <- c(
  122, 244, 214, 173, 229, 156, 212, 263, 146, 183, 161, 205, 135, 331, 225,
  174, 98.8, 149, 238, 262, 132, 235, 216, 240, 230, 192, 195, 172, 173, 172,
  153, 142, 317, 161, 201, 204, 194, 164, 183, 161, 167, 179, 185, 117, 192,
  337, 125, 166, 99.1, 202, 230, 158, 262, 154, 164, 182, 164, 183, 171, 250,
  184, 205, 237, 177, 239, 187, 180, 173, 174
)
length(harricana)
[1] 69

Four questions a flood-frequency analysis usually asks of a record before fitting anything to it: is it normally distributed, is it independent from year to year, and is it homogeneous (no jump, no trend)? hypothesis_test() answers each with a p-value.

c(
  jarque_bera    = hypothesis_test(harricana, method = "jarque_bera")[["p_value"]],
  wald_wolfowitz = hypothesis_test(harricana, method = "wald_wolfowitz")[["p_value"]],
  ljung_box      = hypothesis_test(harricana, method = "ljung_box")[["p_value"]],
  mann_kendall   = hypothesis_test(harricana, method = "mann_kendall")[["p_value"]]
)
   jarque_bera wald_wolfowitz      ljung_box   mann_kendall 
    0.00105658     0.24359142     0.82041073     0.77572941

Jarque-Bera tests normality, and its p-value here, 0.001, rejects it firmly – unsurprising for a right-skewed peak-flow record (skewness 0.86 below), which is exactly why flood-frequency work fits a skewed family instead of a Normal one. Wald-Wolfowitz (a runs test) and Ljung-Box (serial correlation up to a default lag) both test independence between consecutive years; neither rejects it. Mann-Kendall tests for a monotonic trend – homogeneity in the “no trend” sense – and its own large p-value agrees: no trend detected over these 69 years.

Mann-Whitney tests the other kind of homogeneity, a jump: are the first and second halves of the record drawn from the same distribution? It is the one two-sample test above, so it needs both samples split out by hand when called directly (analysis_data_hypothesis_test() below does the split for you, from a split_index):

first_half <- harricana[1:34]
second_half <- harricana[35:69]
hypothesis_test(first_half, second_half, method = "mann_whitney")
 p_value 
0.687614

The DataFrame path

analysis_data() wraps a record the way an actual flood-frequency analysis would – with room for censored observations the plain vector above has none of. Its own hypothesis-test and summary-statistics facades read the same twelve HypothesisTests statics, plus a use_log10 switch and a split-index-based two-sample split:

d <- analysis_data(harricana)
d
<corehydro_data> 69 exact
analysis_data_hypothesis_test(d, "mann_kendall")
mann_kendall 
   0.7757294 
analysis_data_hypothesis_test(d, "mann_whitney", split_index = 34)
mann_whitney 
    0.687614

The split_index split is on the record’s INDEX (0-based, so split_index = 34 puts the observations with index below 34 – the first 34 – in the first sample and the rest in the second) – matching the manual harricana[1:34]/harricana[35:69] split above exactly.

use_log10 = TRUE reruns the same test against log-space values. For a rank-based test like Mann-Kendall this changes nothing – log10 is monotonic, so it cannot change which of two years’ peaks is larger, and Mann-Kendall depends only on those pairwise orderings:

c(
  real_space = analysis_data_hypothesis_test(d, "mann_kendall"),
  log_space  = analysis_data_hypothesis_test(d, "mann_kendall", use_log10 = TRUE)
)
real_space.mann_kendall  log_space.mann_kendall 
              0.7757294               0.7757294

Jarque-Bera is not rank-based, so use_log10 changes its answer outright – this is the log-normal version of the normality question asked above, and it comes back with a very different verdict:

c(
  real_space = analysis_data_hypothesis_test(d, "jarque_bera"),
  log_space  = analysis_data_hypothesis_test(d, "jarque_bera", use_log10 = TRUE)
)
real_space.jarque_bera  log_space.jarque_bera 
            0.00105658             0.57562570

Log10-transformed, the record no longer rejects normality – the log-Normal (or log-Pearson) family a flood-frequency analysis usually reaches for is a reasonable fit exactly because this record’s SHAPE, not just its scale, is closer to Normal after the transform.

analysis_data_statistics() reports the twenty summary-statistics keys DataFrame computes – moments, percentiles, both in real space and log space:

s <- analysis_data_statistics(d)
s$value[c("Record Length", "Mean", "Std Dev", "Skewness", "Kurtosis", "50%")]
Record Length          Mean       Std Dev      Skewness      Kurtosis 
   69.0000000   191.3173913    47.9616114     0.8605451     4.3434868 
          50% 
  183.0000000

A correlation matrix

correlation() takes a matrix (or a data frame) instead of two vectors and returns every pairwise coefficient at once. Ten years of peak flow at three tributary gauges:

gauge_a <- c(122, 244, 214, 173, 229, 156, 212, 263, 146, 183)
gauge_b <- c(98, 210, 187, 145, 201, 130, 178, 240, 120, 155)
gauge_c <- c(310, 180, 405, 220, 350, 190, 300, 410, 175, 260)

m <- correlation(cbind(a = gauge_a, b = gauge_b, c = gauge_c))
m
          a         b         c
a 1.0000000 0.9969634 0.4653750
b 0.9969634 1.0000000 0.5043234
c 0.4653750 0.5043234 1.0000000

Gauges A and B track closely (a shared upstream driver, high peak years at one are high peak years at the other); gauge C answers to something else – its correlation with either is much weaker.

A reservoir stage-storage curve

The rest of this page is the Paired Data subsystem: OrderedPairedData, an x-y curve that keeps itself sorted and validated, with interpolation, trapezoidal area, and three ways to reduce its point count. The curve below is upstream’s own 15-point reservoir fixture – storage in acre-feet against pool elevation in feet:

storage <- c(230408, 288010, 345611, 403213, 460815, 518417, 576019, 633612,
             691223, 748825, 806427, 864029, 921631, 1036834, 1152038)
elevation <- c(1519.7, 1520.5, 1520.9, 1521.7, 1523.5, 1525.9, 1528.4, 1530.9,
               1533.2, 1534.7, 1535.9, 1538, 1541.3, 1547.7, 1552.7)

curve_interpolate() reads the elevation at a given storage, linearly by default:

target <- 500000
lin <- curve_interpolate(storage, elevation, xout = target)
log_space <- curve_interpolate(storage, elevation, xout = target, x_transform = "logarithmic")
c(linear = lin, log_x = log_space)
  linear    log_x 
1525.133 1525.163

The two differ by a small fraction of a foot here – the curve is nearly linear over this stretch, so which space you interpolate in barely matters; it would matter much more over a stretch where storage grows convexly with elevation.

curve_area() is the trapezoidal-rule area between the curve and the x-axis:

curve_area(storage, elevation)
[1] 1413175623

Three ways to simplify the same curve

curve_simplify() reduces a curve’s point count three ways – Ramer-Douglas-Peucker, Visvalingam- Whyatt, and Lang – each trading a tolerance (or a target count) for fewer points:

rdp <- curve_simplify(storage, elevation, method = "rdp", tolerance = 0.5)
vis <- curve_simplify(storage, elevation, method = "visvalingam", num_to_keep = 7)
lang <- curve_simplify(storage, elevation, method = "lang", tolerance = 0.5, look_ahead = 2)

data.frame(
  method = c("rdp", "visvalingam", "lang"),
  retained_points = c(nrow(rdp), nrow(vis), nrow(lang)),
  last_x_matches = c(
    tail(rdp$x, 1) == tail(storage, 1),
    tail(vis$x, 1) == tail(storage, 1),
    tail(lang$x, 1) == tail(storage, 1)
  )
)
       method retained_points last_x_matches
1         rdp               7           TRUE
2 visvalingam               7           TRUE
3        lang               8          FALSE

Ramer-Douglas-Peucker and Visvalingam-Whyatt both always keep the curve’s first and last point. Lang does not, and this is not merely a theoretical concern raised elsewhere on this site – at this tolerance and look-ahead it drops the reservoir curve’s own last point right here, the same defect docs/upstream-csharp-issues.md documents against the sin curve upstream’s own Test_LangSimplify uses. last_x_matches is FALSE for Lang and TRUE for the other two.

Sampling an uncertain curve

UncertainOrderedPairedData is the same idea with a distribution at every y instead of a number – a stage-storage curve, say, whose elevation at each storage level is only known within a range. uncertain_curve_sample() collapses it to a plain curve at a chosen quantile, or at the mean:

storage_pts <- storage[1:5]
elev_dists <- lapply(elevation[1:5], function(e) distribution("Triangular", c(e - 0.6, e, e + 1.2)))

median_curve <- uncertain_curve_sample(storage_pts, elev_dists, probability = 0.5)
mean_curve <- uncertain_curve_sample(storage_pts, elev_dists)

data.frame(storage = storage_pts, median = median_curve$y, mean = mean_curve$y)
  storage   median   mean
1  230408 1519.861 1519.9
2  288010 1520.661 1520.7
3  345611 1521.061 1521.1
4  403213 1521.861 1521.9
5  460815 1523.661 1523.7

The Triangular distributions here are not symmetric (each is skewed further above its mode than below), so the mean and median curves differ by a consistent small offset – exactly what you would expect from a right-skewed uncertainty band at every point, and exactly what collapsing to a single number by mean instead of by median would hide.

Key takeaways

  1. hypothesis_test() reaches all twelve HypothesisTests statics directly; analysis_data_*() reaches nine of them through a DataFrame, with the censored-data and use_log10/split_index machinery a raw vector does not have.
  2. correlation() without a second argument treats its input as a table and returns the whole pairwise matrix in one call.
  3. curve_interpolate()/curve_area() work in linear or log-transformed space on the same curve.
  4. Ramer-Douglas-Peucker and Visvalingam-Whyatt always keep a curve’s endpoints; Lang does not – an upstream defect this port reproduces on purpose rather than silently fixing.
  5. uncertain_curve_sample() collapses a distribution-valued curve to a plain one at any quantile, including the mean and the median, which need not agree.

Reproduction check

All values below are deterministic given the fixed inputs above; checked at 1e-15 relative tolerance except the Mann-Kendall p-value, which is additionally pinned against the upstream C# test literal (ExactDataHypothesisTests.Test_MannKendall, fixtures/data/data_frame_facades.json) at that test’s own tolerance. Mann-Whitney has no such upstream literal pinned in this repo – it is checked here only for internal bit-for-bit reproducibility (the analysis_data_*() path against the direct hypothesis_test() call on the manually-split halves), the same standard every other non-oracle value on this page is held to.

near <- \(x, literal, tol = 1e-15) abs(x / literal - 1) < tol

mk <- analysis_data_hypothesis_test(d, "mann_kendall")
mw <- analysis_data_hypothesis_test(d, "mann_whitney", split_index = 34)
jb_log <- analysis_data_hypothesis_test(d, "jarque_bera", use_log10 = TRUE)

stopifnot(
  # Pinned against the upstream C# oracle.
  abs(unname(mk) - 0.7757) < 1e-3,

  # This page's own deterministic numbers, checked bit-for-bit against themselves.
  near(unname(mk), unname(hypothesis_test(harricana, method = "mann_kendall")[["p_value"]])),
  near(unname(mw), unname(hypothesis_test(first_half, second_half, method = "mann_whitney")[["p_value"]])),

  # The correlation matrix, pinned to this page's own executed values (not merely checked for
  # symmetry, which is guaranteed by construction and proves nothing about the actual numbers).
  near(m["a", "b"], 0.99696337794688294),
  near(m["a", "c"], 0.46537498754641160),
  near(m["b", "c"], 0.50432343340073482),

  near(curve_area(storage, elevation), 1413175623, tol = 1e-9),

  # The two interpolations, pinned to this page's own executed values (not merely recomputed and
  # compared against themselves, which is vacuous).
  near(lin, 1525.1326516440402),
  near(log_space, 1525.1629478762325),

  nrow(rdp) == 7L,
  nrow(vis) == 7L,
  nrow(lang) == 8L,
  isTRUE(tail(rdp$x, 1) == tail(storage, 1)),
  isTRUE(tail(vis$x, 1) == tail(storage, 1)),
  isFALSE(tail(lang$x, 1) == tail(storage, 1)),   # the LangSimplify drop, reproduced on purpose
  jb_log > 0.5,                                    # log-space does not reject normality
  near(mean_curve$y[1], 1519.8999999999999),
  near(median_curve$y[1], 1519.8607695154587),
  near(mean_curve$y[1] - median_curve$y[1], 0.039230484541121768)
)
cat("All reproduction checks passed.\n")
All reproduction checks passed.