The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
The Lepanthes eltoroensis
analysis
library(recruitlag)
library(dplyr)
This vignette reproduces the analysis of the companion paper on the
epiphytic orchid Lepanthes eltoroensis (Luquillo Mountains,
Puerto Rico), in which reproduction is continuous in time and
recruitment is episodic. The package’s tests check the numbers below
against the archived outputs of the paper’s pipeline; here they are
simply computed.
The two records
lepanthes_monthly
#> # A tibble: 23 × 3
#> period reproductive_adults recruits
#> <int> <int> <int>
#> 1 1 60 NA
#> 2 2 64 0
#> 3 3 84 4
#> 4 4 71 1
#> 5 5 80 0
#> 6 6 81 0
#> 7 7 86 0
#> 8 8 78 0
#> 9 9 61 0
#> 10 10 54 0
#> # ℹ 13 more rows
lepanthes_census |> select(period, census, adults, inflorescences, recruits)
#> # A tibble: 12 × 5
#> period census adults inflorescences recruits
#> <int> <chr> <int> <int> <dbl>
#> 1 1 Mar 1999 264 509 NA
#> 2 2 Sep 1999 290 218 NA
#> 3 3 Mar 2000 311 216 5
#> 4 4 Sep 2000 324 237 3
#> 5 5 Mar 2001 333 225 62
#> 6 6 Sep 2001 303 142 30
#> 7 7 Mar 2002 289 185 0
#> 8 8 Sep 2002 288 171 0
#> 9 9 Mar 2003 274 128 14
#> 10 10 Sep 2003 243 82 32
#> 11 11 Mar 2004 217 118 3
#> 12 12 Sep 2004 179 111 0
lepanthes_hosts
#> # A tibble: 276 × 5
#> host period adults inflorescences recruits
#> <chr> <int> <int> <int> <dbl>
#> 1 209 1 6 10 NA
#> 2 209 2 9 11 NA
#> 3 209 3 9 4 0
#> 4 209 4 9 5 0
#> 5 209 5 8 6 0
#> 6 209 6 7 5 1
#> 7 209 7 7 6 0
#> 8 209 8 7 3 0
#> 9 209 9 8 1 0
#> 10 209 10 10 3 1
#> # ℹ 266 more rows
The monthly survey (1994 to 1996) counts inflorescence-bearing adults
and new seedlings at 23 monthly censuses of four populations; recruits
are defined by first appearance, so they are scored from census 2. The
six-monthly survey (1999 to 2004) counts every plant on 23 host trees at
12 censuses; recruits are scored at censuses 3 to 12.
lepanthes_census holds its system totals and
lepanthes_hosts one row per host tree and census, which is
the input of the host-level test.
The lag conventions differ, and the package makes them explicit. The
monthly lag model of the paper reads lag k at month t -
k over a six-month horizon (K = 6,
lag0 = 0). The six-monthly model reads lag k at
census t - 1 - k over five six-month bins (K = 4,
lag0 = 1). The bins are labelled 0-6, 6-12, 12-18, 18-24
and 24-30 months, counted from the first lagged census, which is the
convention of the kernel labels throughout.
Concentration and the ceiling
k_m <- lag_kernels(6, bin = 1, unit = "mo")
k_6 <- lag_kernels(4, bin = 6, unit = "mo")
c_m <- lag_ceiling(lepanthes_monthly, reproduction = "reproductive_adults", K = 6, lag0 = 0, kernels = k_m)
c_6 <- lag_ceiling(lepanthes_census, reproduction = "inflorescences", K = 4, lag0 = 1, kernels = k_6)
c_m
#> Ceiling on the concentration of expected recruitment
#> reproductive record: 23 periods; recruits scored at 22 periods; K = 6, lag0 = 0
#> reproduction (theorem bound) Gini 0.089 CV 0.159
#> searched ceiling Gini 0.095 CV 0.173 (pure delay 6-7 mo, 33 kernels)
#> observed recruitment Gini 0.859 CV 2.749
#> exceedance (obs / ceiling) Gini 9.00 CV 15.88
#> The exceedance bounds the expected series, not the counts; see ceiling_calibration() and host_lag_test().
c_6
#> Ceiling on the concentration of expected recruitment
#> reproductive record: 12 periods; recruits scored at 10 periods; K = 4, lag0 = 1
#> reproduction (theorem bound) Gini 0.259 CV 0.570
#> searched ceiling Gini 0.254 CV 0.547 (pure delay 12-18 mo, 20 kernels)
#> observed recruitment Gini 0.649 CV 1.377
#> exceedance (obs / ceiling) Gini 2.56 CV 2.51
#> The exceedance bounds the expected series, not the counts; see ceiling_calibration() and host_lag_test().
At both time scales the ceiling is low because the input is flat, the
searched maximum is attained by a pure delay and lies within 0.01 of the
theorem bound (the small difference arises from backfilling pre-record
lags), and observed recruitment exceeds the ceiling by a factor of about
9 (monthly) and 2.6 (six-monthly) on the Gini.
rain_strip(`monthly 1994` = c_m, `six-monthly 1999-2004` = c_6)
What any lag profile predicts
profiles <- list(
`one interval only (projection-matrix null)` = c(1, 0, 0, 0, 0),
`concentrated early (0 to 6 months)` = c(.7, .2, .07, .03, 0),
`evenly spread (0 to 30 months)` = rep(.2, 5),
`concentrated late (24 to 30 months)` = c(0, .03, .07, .2, .7),
`fitted posterior mean` = lepanthes_profile$weight)
rain_profiles(profiles, bin = 6, unit = "months")
rain_expected(c_6, profiles)
The expected series all but coincide and follow the slow decline of
reproduction; none resembles the recruits counted (62 at census 5, none
at censuses 7, 8 and 12, 32 at census 10). The fitted posterior-mean
profile from the paper’s distributed-lag model is among them and is
indistinguishable from the rest, which is the identification problem the
paper describes: with a flat input the lag cannot be recovered, however
long the series.
How firm the exceedance is
The ceiling is arithmetic and carries no uncertainty. The two
concentration indices it is compared with are estimates from ten
censuses, and so is their ratio, and a reader is entitled to an interval
on them. The package gives three, because they answer differently worded
questions and a user should take the one whose logic they accept.
The first resamples the 23 host trees with replacement, rebuilds both
series from the resampled trees, and recomputes the ceiling and the
exceedance in each replicate. It is a confidence interval, BCa by
default, and is read in the usual way. The share of replicates at or
below 1 is reported as what it is, a one-sided bootstrap p-value, and
not as a probability that the exceedance is at most 1.
set.seed(20260913)
k4 <- lag_kernels(4, bin = 6, unit = "mo")
lag_ceiling_boot(lepanthes_hosts, unit = "host", reproduction = "inflorescences",
K = 4, kernels = k4, R = 2000)
#> Cluster bootstrap of the ceiling: 23 units resampled 2000 times, 90% BCa intervals
#> estimand: realised counts; K = 4, lag0 = 1, 20 kernels
#> observed Gini 0.649 [0.578, 0.694]
#> ceiling Gini 0.254 [0.196, 0.291]
#> exceedance Gini 2.56 [2.11, 3.18]
#> observed CV 1.377 [1.143, 1.598]
#> ceiling CV 0.547 [0.426, 0.637]
#> exceedance CV 2.51 [1.94, 3.27]
#> one-sided bootstrap p-value for exceedance <= 1: Gini < 5e-04, CV < 5e-04
#> A confidence interval, read in the usual way; p_boot is not a posterior probability.
The second is Rubin’s Bayesian bootstrap: Dirichlet weights over the
hosts in place of resampling, and a posterior in place of a sampling
distribution. The numbers are close to the first, which is expected; the
difference is that P(exceedance > 1) is now a posterior
probability under a stated prior and may be written as one.
lag_ceiling_bayesboot(lepanthes_hosts, unit = "host", reproduction = "inflorescences",
K = 4, kernels = k4, draws = 4000)
#> Bayesian bootstrap of the ceiling: 23 units, 4000 posterior draws, 90% equal-tailed credible intervals
#> estimand: realised counts; K = 4, lag0 = 1, 20 kernels
#> observed Gini 0.649 [0.609, 0.731]
#> ceiling Gini 0.254 [0.222, 0.318]
#> exceedance Gini 2.56 [2.06, 3.07]
#> observed CV 1.377 [1.228, 1.809]
#> ceiling CV 0.547 [0.462, 0.669]
#> exceedance CV 2.51 [2.04, 3.46]
#> P(exceedance > 1): Gini > 0.9998, CV > 0.9998
The third fits a negative-binomial model to the counts on every host
at every census, with host and census effects, and puts the posterior of
the expected series through the same ceiling
arithmetic. That is the object the bound applies to, and its exceedance
is lower than the raw one, because part of the concentration of small
counts is sampling variation about a smoother mean. It needs
cmdstanr and CmdStan, so it is not run here.
lag_ceiling_stan(lepanthes_hosts, unit = "host", reproduction = "inflorescences",
K = 4, kernels = k4)
Why the pooled comparison cannot decide
set.seed(20260903)
cal <- ceiling_calibration(c_6)
cal
#> Observed recruitment: Gini 0.649, CV 1.377, 3 zero periods of 10
#> Simulated counts about a mean that obeys the ceiling:
#> # A tibble: 2 × 11
#> null phi gini_med gini_q05 gini_q95 p_gini cv_med cv_q05
#> * <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 N0 flat mean 0.547 0.570 0.408 0.72 0.207 1.19 0.781
#> 2 N1 best-kernel mean 0.547 0.598 0.424 0.756 0.306 1.28 0.818
#> cv_q95 p_cv p_zeros
#> * <dbl> <dbl> <dbl>
#> 1 1.81 0.279 0.203
#> 2 2.05 0.394 0.245
#> p_* is the probability that a simulated series is at least as concentrated as the observed one.
About a mean that obeys the ceiling, negative-binomial counts with
the observed variance-to-mean ratio produce a series as concentrated as
the observed one roughly one time in 3. The exceedance of the ceiling by
the pooled counts is therefore not evidence; the test has to be made
across host trees.
The host-level test
lepanthes_hosts |>
filter(!is.na(recruits)) |>
group_by(period) |>
summarise(recruits = sum(recruits), hosts_recruiting = sum(recruits > 0))
#> # A tibble: 10 × 3
#> period recruits hosts_recruiting
#> <int> <dbl> <int>
#> 1 3 5 1
#> 2 4 3 1
#> 3 5 62 1
#> 4 6 30 1
#> 5 7 0 0
#> 6 8 0 0
#> 7 9 14 1
#> 8 10 32 1
#> 9 11 3 1
#> 10 12 0 0
Three of the ten censuses are silent (no recruit on any of the 23
hosts) and one carries 62. Under the lag hypothesis, expected
recruitment on each host is that host’s own lagged inflorescence record,
scaled to its own total. The kernels are the pure delays, every uniform
window, geometric decay and the fitted profile; the paper adds 800
Dirichlet draws, which take longer and change nothing.
set.seed(20260910)
ht <- host_lag_test(lepanthes_hosts, unit = "host", reproduction = "inflorescences",
K = 4, nsim = 4000,
kernels = lag_kernels(4, bin = 6, unit = "mo",
extra = list(`fitted posterior mean` = lepanthes_profile$weight)),
phi = 0.089)
ht
#> Host-level test of the lag hypothesis: 23 units, 10 censuses, 21 profiles, 4000 simulations each
#> observed census totals: 3 silent, Gini 0.649, CV 1.377, largest 62
#> flat reference (phi = 0.19): P(silent) 0.0000, P(Gini) 0.0003, P(CV) 0.0030
#> most favourable profile for the lag hypothesis, by statistic:
#> moment silent P = 0.0008 (pure delay 12-18 mo)
#> moment gini P = 0.0025 (pure delay 6-12 mo)
#> moment cv P = 0.0320 (pure delay 6-12 mo)
#> moment max P = 0.1403 (pure delay 6-12 mo)
#> fixed 0.089 silent P = 0.0517 (pure delay 6-12 mo)
#> fixed 0.089 gini P = 0.1532 (pure delay 6-12 mo)
#> fixed 0.089 cv P = 0.2650 (pure delay 6-12 mo)
#> fixed 0.089 max P = 0.3320 (pure delay 12-18 mo)
The moment estimate of the clumping parameter under the flat
reference is about 0.19; the paper also evaluates every profile at
0.089, the estimate obtained about a single pooled mean, which
attributes the between-host spread of the means to clumping and is about
twice the clumping the counts support. Even at that value the most
favourable profile does not make three silent censuses ordinary.
plot(ht, statistic = "silent")
plot(ht, statistic = "gini")
The verdict, keyed on labels
ht$verdict
#> # A tibble: 8 × 5
#> phi_source statistic observed max_p kernel
#> <chr> <chr> <dbl> <dbl> <chr>
#> 1 moment silent 3 0.00075 pure delay 12-18 mo
#> 2 moment gini 0.649 0.0025 pure delay 6-12 mo
#> 3 moment cv 1.38 0.032 pure delay 6-12 mo
#> 4 moment max 62 0.140 pure delay 6-12 mo
#> 5 fixed 0.089 silent 3 0.0518 pure delay 6-12 mo
#> 6 fixed 0.089 gini 0.649 0.153 pure delay 6-12 mo
#> 7 fixed 0.089 cv 1.38 0.265 pure delay 6-12 mo
#> 8 fixed 0.089 max 62 0.332 pure delay 12-18 mo
ht$table |>
filter(phi_source == "moment",
label %in% c("flat reference", "pure delay 24-30 mo", "uniform 0-30 mo", "fitted posterior mean")) |>
select(label, phi, gini_expected, p_silent, p_gini, p_cv)
#> # A tibble: 4 × 6
#> label phi gini_expected p_silent p_gini p_cv
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 flat reference 0.190 0 0 0.00025 0.003
#> 2 pure delay 24-30 mo 0.275 0.182 0 0 0.00725
#> 3 uniform 0-30 mo 0.285 0.138 0 0 0.00525
#> 4 fitted posterior mean 0.289 0.137 0 0 0.006
No lag profile, of any shape over the 30-month window, makes the
observed record ordinary. Recruitment in this population is not the
delayed consequence of its own reproduction; something else opens and
closes.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.