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.
Introduction to recruitlag: a worked
analysis, step by step
This vignette is for a reader new to the package. There is no
mathematics in it. Every step is shown, and after each step you should
look at what R printed and check that it matches. If you want to know
why the steps work, read
vignette("mathematics", package = "recruitlag")
afterwards.
The question the package answers is this. A plant population
reproduces (flowers, fruits, seeds) and later new seedlings appear
(recruits). If the seedlings are simply the seeds from some earlier
time, then a record of reproduction and a record of seedlings should
look like each other, shifted in time. The package checks whether that
is possible.
Step 1. Load the packages
library(recruitlag)
library(dplyr)
library(ggplot2)
Loading a package prints nothing. dplyr is for picking
rows and columns out of tables, ggplot2 for drawing.
Step 2. Look at a data set
The package includes a made-up data set called regimes.
It is one table with six census records stacked in it, one row per
record, unit (think “host tree”) and census period.
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
The first line says it has 2,880 rows and 8 columns. We take one
record out of it, the one in which reproduction is steady through the
year but the seedlings arrive in bursts.
r <- regimes |> filter(record == "aseasonal_gated")
r
#> # A tibble: 480 × 8
#> record label seasonality mechanism unit period reproduction recruits
#> <chr> <chr> <dbl> <chr> <chr> <int> <int> <dbl>
#> 1 aseasonal_gat… asea… 0 gated u01 1 2 NA
#> 2 aseasonal_gat… asea… 0 gated u01 2 5 NA
#> 3 aseasonal_gat… asea… 0 gated u01 3 4 NA
#> 4 aseasonal_gat… asea… 0 gated u01 4 5 NA
#> 5 aseasonal_gat… asea… 0 gated u01 5 5 NA
#> 6 aseasonal_gat… asea… 0 gated u01 6 5 0
#> 7 aseasonal_gat… asea… 0 gated u01 7 7 0
#> 8 aseasonal_gat… asea… 0 gated u01 8 6 0
#> 9 aseasonal_gat… asea… 0 gated u01 9 6 0
#> 10 aseasonal_gat… asea… 0 gated u01 10 9 0
#> # ℹ 470 more rows
filter() keeps only the rows of that record: 20 units
times 24 periods is 480 rows. Each row has the reproduction counted on
one unit at one period and the recruits (new seedlings) seen there.
Recruits are NA at periods 1 to 5, which means “not
scored”: the test only uses periods that have a full lag window behind
them.
Step 3. Add up over units
The first test works on the whole population, so add the units
together.
totals <- r |>
group_by(period) |>
summarise(reproduction = sum(reproduction),
recruits = if (all(is.na(recruits))) NA else sum(recruits))
totals
#> # A tibble: 24 × 3
#> period reproduction recruits
#> <int> <int> <dbl>
#> 1 1 332 NA
#> 2 2 330 NA
#> 3 3 323 NA
#> 4 4 329 NA
#> 5 5 313 NA
#> 6 6 315 3
#> 7 7 327 4
#> 8 8 315 7
#> 9 9 345 2
#> 10 10 300 4
#> # ℹ 14 more rows
Look down the reproduction column: the numbers wander
around a few hundred and never do anything dramatic. Now look at
recruits: mostly small numbers, then three huge ones. That
difference is what the test is about. Draw it:
totals |>
tidyr::pivot_longer(c(reproduction, recruits), names_to = "series", values_to = "count") |>
filter(!is.na(count)) |>
ggplot(aes(period, count)) +
geom_col(fill = "grey40", na.rm = TRUE) +
facet_wrap(~ series, ncol = 1, scales = "free_y") +
theme_minimal()
Step 4. Measure how “bursty” each series is
The package measures burstiness with the Gini coefficient. 0 means
every census had the same count; the maximum, which is \((n - 1)/n\) on \(n\) censuses, means everything happened at
one census.
rain_gini(totals$reproduction)
#> [1] 0.02493477
rain_gini(totals$recruits)
#> [1] 0.7266465
(rain_gini() ignores the NAs.) Reproduction
scores very low: it is spread evenly. Recruits score high: they are
concentrated in a few censuses. Note those two numbers.
Step 5. Ask the key question
Could a delay turn the smooth series into the bursty one? The answer
is found with lag_ceiling(). Give it the table of totals
and tell it how far back in time a delay is allowed to reach
(K = 4 means five census intervals). It finds the
period, reproduction and recruits
columns by name.
ce <- lag_ceiling(totals, K = 4, kernels = lag_kernels(4))
ce
#> Ceiling on the concentration of expected recruitment
#> reproductive record: 24 periods; recruits scored at 19 periods; K = 4, lag0 = 1
#> reproduction (theorem bound) Gini 0.025 CV 0.046
#> searched ceiling Gini 0.027 CV 0.050 (pure delay bins 0-0, 20 kernels)
#> observed recruitment Gini 0.727 CV 1.960
#> exceedance (obs / ceiling) Gini 27.21 CV 39.30
#> The exceedance bounds the expected series, not the counts; see ceiling_calibration() and host_lag_test().
Read the printout line by line. The line “reproduction (theorem
bound)” repeats the Gini of reproduction. That number is the
ceiling: no delay, of any shape, can make expected recruitment
burstier than that. The line “searched ceiling” is the largest Gini any
of the delays actually tried produced on the scored censuses; it is
close to the bound. The line “observed recruitment” is the Gini of
recruits. The line “exceedance” divides the observed Gini by the
searched ceiling: observed recruitment is many times burstier than any
delay could make it.
Step 6. Draw it
rain_strip(`aseasonal, gated` = ce)
The strip is a ruler from 0 to 1. The blue tick is the ceiling. The
orange dot is what was observed. The dot sits far to the right of the
tick, in the region no delay can reach.
Step 7. See what a delay would actually predict
Here are four possible delays. The first one is what most population
models assume: seedlings appear one census after reproduction. The
others spread the delay in different ways.
delays <- list(
`one interval only` = c(1, 0, 0, 0, 0),
`early` = c(.7, .2, .07, .03, 0),
`evenly spread` = rep(.2, 5),
`late` = c(0, .03, .07, .2, .7))
rain_profiles(delays)
Now apply each of them to the real reproduction series and compare
with the recruits that were counted:
rain_expected(ce, delays)
The bars are the seedlings counted. The lines are what each delay
predicts. All four lines are nearly flat and nearly on top of one
another, and none of them goes anywhere near the three tall bars. This
is the picture that the numbers in Step 5 summarise.
Step 8. The careful version of the test
Step 5 compared the average that a delay predicts with the
counts that were observed. Counts are noisy, and a noisy count
series can look bursty even when its average is flat. The careful test
takes that into account and also uses the fact that recruits were
counted separately on 20 units. It works from the unit-by-period table
r directly. It takes a few seconds because it simulates
many count series.
set.seed(1)
ht <- host_lag_test(r, K = 4, nsim = 2000, kernels = lag_kernels(4))
ht
#> Host-level test of the lag hypothesis: 20 units, 19 censuses, 20 profiles, 2000 simulations each
#> observed census totals: 0 silent, Gini 0.727, CV 1.960, largest 123
#> flat reference (phi = 0.23): P(silent) 1.0000, P(Gini) 0.0000, P(CV) 0.0000
#> most favourable profile for the lag hypothesis, by statistic:
#> moment silent P = 1.0000 (pure delay bins 0-0)
#> moment gini P = 0.0000 (pure delay bins 0-0)
#> moment cv P = 0.0000 (pure delay bins 0-0)
#> moment max P = 0.0030 (pure delay bins 1-1)
The last lines are the verdict. For each way of measuring burstiness
the test tried every delay it knows and reports the best case for
the delay hypothesis: the largest probability, over all delays, of
getting a record as bursty as the one observed. If even the best case is
a very small number, no delay can explain the record. The full table,
one row per delay, is in ht$table:
ht$table |> select(label, phi, p_silent, p_gini, p_cv) |> arrange(desc(p_gini)) |> head()
#> # A tibble: 6 × 5
#> label phi p_silent p_gini p_cv
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 flat reference 0.234 1 0 0
#> 2 pure delay bins 0-0 0.242 1 0 0
#> 3 pure delay bins 1-1 0.247 1 0 0
#> 4 pure delay bins 2-2 0.238 1 0 0
#> 5 pure delay bins 3-3 0.238 1 0 0
#> 6 pure delay bins 4-4 0.243 1 0 0
Step 9. A check with a record where a delay is the truth
The same data set has a record in which recruits really are a delayed
consequence of reproduction. Repeat Steps 2, 3, 5 and 8 on it.
r2 <- regimes |> filter(record == "aseasonal_delayed")
ce2 <- lag_ceiling(r2, K = 4, kernels = lag_kernels(4)) # a unit-by-period table is summed for you
ce2
#> Ceiling on the concentration of expected recruitment
#> reproductive record: 24 periods; recruits scored at 19 periods; K = 4, lag0 = 1
#> reproduction (theorem bound) Gini 0.027 CV 0.048
#> searched ceiling Gini 0.027 CV 0.050 (pure delay bins 3-3, 20 kernels)
#> observed recruitment Gini 0.156 CV 0.282
#> exceedance (obs / ceiling) Gini 5.75 CV 5.69
#> The exceedance bounds the expected series, not the counts; see ceiling_calibration() and host_lag_test().
set.seed(1)
ht2 <- host_lag_test(r2, K = 4, nsim = 2000, kernels = lag_kernels(4))
ht2
#> Host-level test of the lag hypothesis: 20 units, 19 censuses, 20 profiles, 2000 simulations each
#> observed census totals: 0 silent, Gini 0.156, CV 0.282, largest 38
#> flat reference (phi = 1.95): P(silent) 1.0000, P(Gini) 0.3250, P(CV) 0.3705
#> most favourable profile for the lag hypothesis, by statistic:
#> moment silent P = 1.0000 (pure delay bins 0-0)
#> moment gini P = 0.4315 (pure delay bins 0-0)
#> moment cv P = 0.4820 (pure delay bins 0-0)
#> moment max P = 0.7080 (pure delay bins 0-0)
Two things to notice. First, in the ceiling printout the observed
Gini is again above the ceiling: counts are noisier than averages, so
this on its own proves nothing (that is exactly why Step 8 exists).
Second, in the careful test the best-case probabilities are now ordinary
numbers, not tiny ones: the careful test does not reject a delay when a
delay is the truth.
Step 10. The other way round: a seasonal plant
Everything so far used a plant that reproduces steadily all year.
Many plants do not: they flower in one season. Repeat Steps 2 to 5 on
such a record and watch what changes.
r3 <- regimes |> filter(record == "seasonal_delayed")
totals3 <- r3 |> group_by(period) |>
summarise(reproduction = sum(reproduction),
recruits = if (all(is.na(recruits))) NA else sum(recruits))
totals3$reproduction
#> [1] 42 308 933 1261 878 342 64 12 1 0 0 9 57 299 919
#> [16] 1308 895 319 52 7 0 0 0 2
Large numbers for a few periods, then almost nothing, then large
again a year later. That is seasonal reproduction. Measure it:
rain_gini(totals3$reproduction)
#> [1] 0.6750995
rain_gini(totals3$recruits)
#> [1] 0.5343037
Reproduction now scores high on the burstiness scale, because it is
concentrated in a season, and recruits score a little lower
than reproduction. Now the ceiling:
ce3 <- lag_ceiling(totals3, K = 4, kernels = lag_kernels(4))
ce3
#> Ceiling on the concentration of expected recruitment
#> reproductive record: 24 periods; recruits scored at 19 periods; K = 4, lag0 = 1
#> reproduction (theorem bound) Gini 0.675 CV 1.381
#> searched ceiling Gini 0.707 CV 1.513 (pure delay bins 0-0, 20 kernels)
#> observed recruitment Gini 0.534 CV 1.002
#> exceedance (obs / ceiling) Gini 0.76 CV 0.66
#> The exceedance bounds the expected series, not the counts; see ceiling_calibration() and host_lag_test().
rain_strip(`aseasonal, gated` = ce, `seasonal, true delay` = ce3)
In the strip for the seasonal plant the orange dot sits to the
left of the blue tick, inside the region a delay can reach. So
a delay could have produced this record, and in this case one did: this
record was built with a true delay.
This is the important point. A dot to the left of the tick does not
prove there was a delay. It only says the test cannot rule one out. Try
the seasonal record with the gate instead:
r4 <- regimes |> filter(record == "seasonal_gated")
ce4 <- lag_ceiling(r4, K = 4, kernels = lag_kernels(4))
ce4
#> Ceiling on the concentration of expected recruitment
#> reproductive record: 24 periods; recruits scored at 19 periods; K = 4, lag0 = 1
#> reproduction (theorem bound) Gini 0.668 CV 1.357
#> searched ceiling Gini 0.702 CV 1.491 (pure delay bins 0-0, 20 kernels)
#> observed recruitment Gini 0.700 CV 1.793
#> exceedance (obs / ceiling) Gini 1.00 CV 1.20
#> The exceedance bounds the expected series, not the counts; see ceiling_calibration() and host_lag_test().
rain_strip(`seasonal, true delay` = ce3, `seasonal, gated` = ce4)
The gated record lands almost on the tick. With a seasonal plant, the
true delay and the gate look nearly the same on the ruler, because the
ceiling is so high that almost anything fits under it. The careful test
of Step 8 does not help much here either:
set.seed(1)
host_lag_test(r4, K = 4, nsim = 2000, kernels = lag_kernels(4))
#> Host-level test of the lag hypothesis: 20 units, 19 censuses, 20 profiles, 2000 simulations each
#> observed census totals: 0 silent, Gini 0.700, CV 1.793, largest 90
#> flat reference (phi = 0.25): P(silent) 1.0000, P(Gini) 0.0000, P(CV) 0.0000
#> most favourable profile for the lag hypothesis, by statistic:
#> moment silent P = 1.0000 (pure delay bins 0-0)
#> moment gini P = 0.8305 (pure delay bins 0-0)
#> moment cv P = 0.1185 (pure delay bins 0-0)
#> moment max P = 0.4405 (pure delay bins 0-0)
The best-case probabilities are ordinary numbers, even though we know
this record was gated. That is not a failure of the test; it is the
honest answer. When reproduction is seasonal, a burst of seedlings can
always be explained as “the seeds from last season”, and no census
record can say otherwise.
Step 11. What the four cases teach
You have now seen all four combinations. Keep this table in mind when
you read your own results.
| steady all year |
bursts (gated) |
far right of the tick |
tiny probabilities: no delay can explain it |
| steady all year |
true delay |
right of the tick, but only because counts are noisy |
ordinary probabilities: a delay is possible |
| seasonal |
true delay |
left of the tick |
ordinary: a delay is possible |
| seasonal |
bursts (gated) |
near the tick |
ordinary: cannot tell |
The test can only say “no delay can explain this” when reproduction
is steady. When reproduction is seasonal, it cannot say much either way.
That is the whole reason the Lepanthes orchid in the
lepanthes_* data sets is a good case: it flowers every
month of the year.
Step 12. Your own data
Put your census record in a table with one row per unit (host tree,
plot, population) and period, and these columns: unit,
period (an integer starting at 1),
reproduction and recruits (with
NA at periods where recruits were not scored). If your
columns have other names, tell the functions with unit =,
period =, reproduction = and
recruits =, as the Lepanthes vignette does. Then
repeat Steps 3 to 8 with your table in place of r, and read
your result against the table in Step 11. If your first lagged bin
should be the census itself rather than the previous census, add
lag0 = 0 to lag_ceiling() and
host_lag_test().
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.