Invented records: where the test has power and where it does not

library(recruitlag)
library(dplyr)
library(tidyr)
library(ggplot2)

The regimes tibble holds six constructed census records, one row per record, unit and period: three reproductive regimes (seasonal, intermediate, aseasonal) crossed with two recruitment mechanisms. Under delayed, recruits on each of 20 units are a true delayed consequence of that unit’s own reproduction, through a fixed five-bin lag profile with negative-binomial counts. Under gated, recruits arrive only at three of 24 periods, when an external gate opens, at a rate proportional to each unit’s average reproduction, so the record is episodic and does not track reproduction. Nothing here is data; the generating rules are in regimes_truth and in the script data-raw/regimes.R of the source repository (https://github.com/RaymondLTremblay/recruitlag).

regimes
#> # A tibble: 2,880 × 8
#>    record         label seasonality mechanism unit  period reproduction recruits
#>    <chr>          <chr>       <dbl> <chr>     <chr>  <int>        <int>    <dbl>
#>  1 seasonal_dela… seas…           1 delayed   u01        1            4       NA
#>  2 seasonal_dela… seas…           1 delayed   u01        2            6       NA
#>  3 seasonal_dela… seas…           1 delayed   u01        3           14       NA
#>  4 seasonal_dela… seas…           1 delayed   u01        4           31       NA
#>  5 seasonal_dela… seas…           1 delayed   u01        5           17       NA
#>  6 seasonal_dela… seas…           1 delayed   u01        6            2        0
#>  7 seasonal_dela… seas…           1 delayed   u01        7            0        3
#>  8 seasonal_dela… seas…           1 delayed   u01        8            0        3
#>  9 seasonal_dela… seas…           1 delayed   u01        9            0        1
#> 10 seasonal_dela… seas…           1 delayed   u01       10            0        0
#> # ℹ 2,870 more rows
regimes_truth
#> # A tibble: 29 × 4
#>    mechanism index what                                        value
#>    <chr>     <int> <chr>                                       <dbl>
#>  1 delayed       1 lag bin (periods before the recruit census)  0.1 
#>  2 delayed       2 lag bin (periods before the recruit census)  0.25
#>  3 delayed       3 lag bin (periods before the recruit census)  0.35
#>  4 delayed       4 lag bin (periods before the recruit census)  0.2 
#>  5 delayed       5 lag bin (periods before the recruit census)  0.1 
#>  6 gated         1 period                                       0.03
#>  7 gated         2 period                                       0.03
#>  8 gated         3 period                                       0.03
#>  9 gated         4 period                                       0.03
#> 10 gated         5 period                                       0.03
#> # ℹ 19 more rows

The point of the set is to show three things: the ceiling is high where reproduction is seasonal and low where it is aseasonal; the pooled comparison of observed concentration with the ceiling cannot decide the question on its own; and the host-level test separates a true delay from a gate where the ceiling is low, and loses that power as reproduction becomes seasonal.

The six records

System totals by period, then reproduction and recruits side by side, each on its own scale from zero (recruits are an order of magnitude fewer and would vanish on a shared axis).

totals <- regimes |>
  group_by(record, label, period) |>
  summarise(reproduction = sum(reproduction),
            recruits = if (all(is.na(recruits))) NA_real_ else sum(recruits),
            .groups = "drop")
totals
#> # A tibble: 144 × 5
#>    record            label                          period reproduction recruits
#>    <chr>             <chr>                           <int>        <int>    <dbl>
#>  1 aseasonal_delayed aseasonal reproduction, true …      1          334       NA
#>  2 aseasonal_delayed aseasonal reproduction, true …      2          327       NA
#>  3 aseasonal_delayed aseasonal reproduction, true …      3          316       NA
#>  4 aseasonal_delayed aseasonal reproduction, true …      4          314       NA
#>  5 aseasonal_delayed aseasonal reproduction, true …      5          321       NA
#>  6 aseasonal_delayed aseasonal reproduction, true …      6          294       28
#>  7 aseasonal_delayed aseasonal reproduction, true …      7          299       26
#>  8 aseasonal_delayed aseasonal reproduction, true …      8          308       12
#>  9 aseasonal_delayed aseasonal reproduction, true …      9          291       30
#> 10 aseasonal_delayed aseasonal reproduction, true …     10          299       16
#> # ℹ 134 more rows
totals |>
  pivot_longer(c(reproduction, recruits), names_to = "series", values_to = "count") |>
  filter(!is.na(count)) |>                      # unscored periods are not drawn
  mutate(series = factor(series, c("reproduction", "recruits")),
         label = factor(label, unique(regimes$label))) |>
  mutate(panel = factor(paste(label, series, sep = ": "),
                        paste(rep(unique(regimes$label), each = 2), c("reproduction", "recruits"), sep = ": "))) |>
  ggplot(aes(period, count, fill = series)) +
  geom_col(width = 0.7, na.rm = TRUE) +
  facet_wrap(~ panel, ncol = 2, scales = "free_y") +
  scale_fill_manual(values = c(reproduction = "#0072B2", recruits = "#D55E00"), guide = "none") +
  theme_minimal(base_size = 8.5) +
  theme(strip.text = element_text(hjust = 0))

Going down the left column, reproduction goes from strongly seasonal to flat. In the right column, the delayed records track reproduction as a smoothed, shifted copy, and the gated records arrive in bursts at the same three periods whatever the reproductive regime.

The ceiling

For each record, the ceiling is the concentration of reproduction itself, and lag_ceiling() confirms by search that no lag profile exceeds it on the scored periods. A unit-by-period table is summed to the system series automatically, so each record can be passed straight in.

k4 <- lag_kernels(4)
ceilings <- regimes |>
  group_split(record) |>
  setNames(sort(unique(regimes$record))) |>
  lapply(lag_ceiling, K = 4, kernels = k4)

tibble(record = names(ceilings),
       ceiling = sapply(ceilings, function(x) x$theorem["gini"]),
       observed = sapply(ceilings, function(x) x$obs["gini"]),
       exceedance = sapply(ceilings, function(x) x$exceedance["gini"])) |>
  mutate(across(where(is.numeric), \(x) round(x, 2)))
#> # A tibble: 6 × 4
#>   record               ceiling observed exceedance
#>   <chr>                  <dbl>    <dbl>      <dbl>
#> 1 aseasonal_delayed       0.03     0.16       5.75
#> 2 aseasonal_gated         0.02     0.73      27.2 
#> 3 intermediate_delayed    0.35     0.33       0.93
#> 4 intermediate_gated      0.34     0.72       2.07
#> 5 seasonal_delayed        0.68     0.53       0.76
#> 6 seasonal_gated          0.67     0.7        1
rain_strip(ceilings)

Two rows of the table carry the lesson. Under seasonal reproduction the ceiling is about 0.68, and both records, the true delay and the gate, sit near it: the gated record exceeds the ceiling only slightly, and a reader could not tell the two mechanisms apart from the pooled series. Under aseasonal reproduction the ceiling is near zero, so the gated record exceeds it by an enormous factor, but so does the true delay, by a large one. That second fact is the reason the pooled comparison cannot be the test: counts are noisier than their expectation, and when the expectation is nearly constant, any count series at all looks concentrated by comparison.

Why the pooled comparison has no power

ceiling_calibration() simulates negative-binomial counts about two means that both obey the ceiling, a flat mean and the best-kernel mean, and asks how often the simulated series is at least as concentrated as the observed one.

set.seed(1)
cal <- lapply(ceilings, ceiling_calibration, nsim = 3000)
bind_rows(cal, .id = "record") |> select(record, null, phi, p_gini, p_cv)
#> Observed recruitment: Gini 0.156, CV 0.282, 0 zero periods of 19
#> Simulated counts about a mean that obeys the ceiling:
#> # A tibble: 12 × 5
#>    record               null                   phi p_gini  p_cv
#>    <chr>                <chr>                <dbl>  <dbl> <dbl>
#>  1 aseasonal_delayed    N0 flat mean        24.6    0.395 0.453
#>  2 aseasonal_delayed    N1 best-kernel mean 24.6    0.431 0.491
#>  3 aseasonal_gated      N0 flat mean         0.264  0.482 0.259
#>  4 aseasonal_gated      N1 best-kernel mean  0.264  0.468 0.251
#>  5 intermediate_delayed N0 flat mean         2.93   0.371 0.426
#>  6 intermediate_delayed N1 best-kernel mean  2.93   0.979 0.985
#>  7 intermediate_gated   N0 flat mean         0.274  0.499 0.268
#>  8 intermediate_gated   N1 best-kernel mean  0.274  0.643 0.419
#>  9 seasonal_delayed     N0 flat mean         1.03   0.185 0.352
#> 10 seasonal_delayed     N1 best-kernel mean  1.03   1     1    
#> 11 seasonal_gated       N0 flat mean         0.317  0.467 0.293
#> 12 seasonal_gated       N1 best-kernel mean  0.317  0.995 0.937
#> p_* is the probability that a simulated series is at least as concentrated as the observed one.

Every probability is ordinary. The calibration says that a mean obeying the ceiling, once sampled as clumped counts, routinely produces series as concentrated as any of the six, gated or not. Pooled over units, the test cannot refute anything.

The host-level test

The null stated at the level at which the lag hypothesis is asserted: expected recruitment on each unit is that unit’s own lagged reproduction, scaled to its own total, with clumping estimated by the moment identity under that heterogeneous mean. Simulated counts are summed over units and compared with the observed census totals on four statistics.

set.seed(2)
tests <- regimes |>
  group_split(record) |>
  setNames(sort(unique(regimes$record))) |>
  lapply(host_lag_test, K = 4, nsim = 2000, kernels = k4)

verdict <- bind_rows(lapply(tests, `[[`, "verdict"), .id = "record") |>
  select(record, statistic, max_p) |>
  pivot_wider(names_from = statistic, values_from = max_p, names_prefix = "P_")
verdict
#> # A tibble: 6 × 5
#>   record               P_silent P_gini   P_cv  P_max
#>   <chr>                   <dbl>  <dbl>  <dbl>  <dbl>
#> 1 aseasonal_delayed       1      0.426 0.484  0.708 
#> 2 aseasonal_gated         1      0     0      0.0035
#> 3 intermediate_delayed    1      0.952 0.984  0.978 
#> 4 intermediate_gated      0.118  0     0.0015 0.109 
#> 5 seasonal_delayed        1      1     1      0.995 
#> 6 seasonal_gated          1      0.838 0.116  0.457

Read the rows in pairs. Under seasonal reproduction the true delay is not refuted (P near 1), but neither is the gate: a lagged seasonal input itself predicts a concentrated census series, and counts about that expectation are as concentrated as the gated record, so the host-level test has little power where the ceiling is high, exactly as the pooled comparison had none. Under intermediate seasonality the true delay is not refuted and the gate is refuted on the Gini and CV of the census totals, because its bursts do not fall where the lagged peaks predict them. Under aseasonal reproduction the true delay is again not refuted, with P in the middle of its range, and the gate is refuted decisively. The silent-census statistic is uninformative in these records because the gate keeps a small background rate open at every period; in the Lepanthes record it is the most interpretable statistic because the gate there closed completely.

plot(tests$aseasonal_gated, statistic = "gini")

Is the host-level test calibrated?

When the true profile is among the kernels, the probability it receives should be uniformly distributed over replicate draws of the data. The following draws new data sets from the aseasonal generating rule, evaluates the true profile on each with a small simulation budget, and shows the distribution of P. It is kept small so that the vignette builds quickly; increase n_rep and nsim to tighten it.

set.seed(11)
H <- 20; Tn <- 24; K <- 4; t_R <- 6:24
w_true <- regimes_truth |> filter(mechanism == "delayed") |> pull(value)
base <- round(exp(rnorm(H, log(12), 0.6)))
n_rep <- 12
p_true <- replicate(n_rep, {
  d <- tibble(unit = rep(sprintf("u%02d", 1:H), each = Tn), period = rep(1:Tn, H),
              reproduction = as.vector(t(sapply(base, function(b) rpois(Tn, b)))))
  d <- d |> group_by(unit) |>
    mutate(mu = c(rep(NA, 5), 0.08 * convolve_lag(reproduction, w_true, t_R, 1)),
           recruits = ifelse(is.na(mu), NA, rnbinom(n(), size = 2, mu = ifelse(is.na(mu), 0, mu)))) |>
    ungroup()
  ht <- host_lag_test(d, K = K, nsim = 300, flat = FALSE,
                      kernels = lag_kernels(K, families = "delay", extra = list(truth = w_true)))
  ht$table$p_gini[ht$table$label == "truth"]
})
round(sort(p_true), 2)
#>  [1] 0.01 0.15 0.17 0.23 0.29 0.37 0.44 0.52 0.63 0.63 0.69 0.91

The values spread across the unit interval rather than piling up near zero, which is what a calibrated test does when the null is true. A single draw can of course land in the tail; the seed of the generating script was chosen so that the delayed records in regimes are typical draws rather than tail ones.

What to take from the set

The ceiling is a statement about expectations, and where reproduction is seasonal it is too high to bite. Where reproduction is aseasonal the ceiling is low, but the pooled comparison still cannot decide the question, because counts about a flat mean look concentrated. The decision comes from replicate units: a true delay predicts, for each unit, a series shaped by that unit’s own reproduction, and a gate does not.