Checking a phenology analysis

Three checks for a phenology analysis in R: effort sensitivity by thinning, survey-window censoring, and whether first, mean and median dates agree on the trend.

The first three tutorials in this series argued that phenological summaries are estimators with sampling properties, and that the first date is the shakiest of them. This one turns that argument into a short checklist. Before you trust a reported onset, sensitivity or trend, three quick diagnostics tell you whether the number reflects the biology or the sampling: does it move when you thin the data, is the season censored by the survey window, and do different summaries agree.

Check 1: does the estimate survive thinning?

Take a well-monitored season and throw away records at random. A summary that tracks the true phenology should barely move; a summary that tracks effort will shift. This is the single most useful check, because it turns the effort bias into something you can measure on data you already have.

set.seed(505)
mu <- 128; sdev <- 12
full <- round(rnorm(300, mu, sdev))              # a well-monitored season
metrics <- function(x) c(first = min(x), mean = round(mean(x), 1), median = median(x))
m_full <- metrics(full)
m_50   <- metrics(sample(full, 150))             # keep half the records
m_20   <- metrics(sample(full, 60))              # keep a fifth
rbind(full = m_full, half = m_50, fifth = m_20)

      first  mean median
full     86 127.5    127
half     94 127.8    127
fifth   104 124.9    125

Halving the records pushes the first date +8 days later while the mean moves only +0.3 days. Keep only a fifth and the first date is +18 days later still. The interpretation is direct: if your onset estimate lurches when you drop data, it was measuring your effort, not the plant. The mean and median pass this check; the first date fails it.

tdf <- data.frame(
  frac = rep(c(100, 50, 20), 3),
  day  = c(m_full["first"], m_50["first"], m_20["first"],
           m_full["mean"],  m_50["mean"],  m_20["mean"],
           m_full["median"],m_50["median"],m_20["median"]),
  metric = rep(c("first date", "mean date", "median date"), each = 3))
ggplot(tdf, aes(frac, day, colour = metric)) +
  geom_line(linewidth = 0.9) + geom_point(size = 2.4) +
  scale_x_reverse(breaks = c(100, 50, 20)) +
  scale_colour_manual(values = c("first date" = terra, "mean date" = forest,
                                 "median date" = gold)) +
  labs(x = "percent of records retained", y = "day of year", colour = NULL)

Figure 1: Each summary under three levels of retained effort. The first date drifts steadily later as records are removed, while the mean and median hold roughly steady.

Check 2: is the season censored by the survey window?

Surveys start and end on fixed dates. If flowering begins before the survey opens, the early plants are simply invisible, and the first observed date is pinned to the opening day rather than measured. This is left censoring, and it is easy to miss because the observed first date looks like a real number.

set.seed(606)
win_lo <- 100; win_hi <- 165                     # the survey runs day 100 to 165
cold  <- round(rnorm(200, 116, 12))              # an early (cold-adapted) year
true_first <- min(cold)
seen  <- cold[cold >= win_lo & cold <= win_hi]   # only what falls inside the window
obs_first <- min(seen)
below <- mean(cold < win_lo)                      # fraction that flowered before the window
c(true_first = true_first, observed_first = obs_first,
  frac_before_window = round(below, 3))

        true_first     observed_first frac_before_window 
             84.00             100.00               0.11 

The true first flowering is day 84, but the survey does not open until day 100, so the observed first date is day 100, exactly the opening day. Fully 11 percent of the population had already flowered before anyone was watching. The tell is that the observed first date equals the survey start (TRUE): when that happens you cannot separate a real onset on the opening day from an onset hidden before it. The diagnostic is to check whether records pile up against the window edge, and to be honest that the onset is a lower bound, not a measurement.

cdf <- data.frame(day = cold, visible = cold >= win_lo)
ggplot(cdf, aes(day, fill = visible)) +
  annotate("rect", xmin = win_lo, xmax = win_hi, ymin = 0, ymax = Inf,
           fill = "#eef1e9") +
  geom_histogram(binwidth = 2, colour = "white") +
  geom_vline(xintercept = win_lo, colour = forest, linewidth = 0.8) +
  annotate("text", x = win_lo, y = 20, label = "survey opens", angle = 90,
           vjust = -0.3, hjust = 1, size = 3.3, colour = forest) +
  scale_fill_manual(values = c("FALSE" = terra, "TRUE" = "#b9c6b0"),
                    labels = c("FALSE" = "flowered unseen", "TRUE" = "recorded")) +
  labs(x = "day of year", y = "plants flowering", fill = NULL) +
  coord_cartesian(clip = "off") + theme(plot.margin = margin(12, 12, 6, 6))

Figure 2: A season whose onset falls before the survey opens. The shaded band is the survey window; the plants to its left flowered unseen, so the observed first date is pinned to the opening day rather than measured.

Check 3: do the summaries agree on the trend?

Run the trend three ways, with the first date, the mean and the median. If they agree, the reported shift does not depend on which summary you chose. If they diverge, either the distribution is changing shape or effort is confounding the result, and the single number you were about to report is not safe.

set.seed(707)
Y <- 15; sdev <- 12
temp <- rnorm(Y, 0, 1)
mean_true <- 130 - 3 * temp                       # true 3 days per degree
n_year <- round(exp(3.2 + 0.8 * temp))            # warm years better monitored
first <- mn <- med <- numeric(Y)
for (y in 1:Y) {
  x <- round(rnorm(n_year[y], mean_true[y], sdev))
  first[y] <- min(x); mn[y] <- mean(x); med[y] <- median(x)
}
slope <- function(z) coef(lm(z ~ temp))[2]
c(first = round(slope(first), 2), mean = round(slope(mn), 2),
  median = round(slope(med), 2))

 first.temp   mean.temp median.temp 
      -9.13       -3.06       -3.08 

The mean and median agree closely, at -3.06 and -3.08 days per degree, both near the true value of -3. The first date reports -9.13, far steeper. That disagreement is the warning: here it flags the effort confound, but the same check catches a season whose shape is changing over time. CaraDonna and colleagues showed with a long flowering record that first, peak and last dates often shift independently of one another, so no single statistic captures the change, which is exactly why cross-checking summaries is worth the few extra lines.

Bringing it together

None of these checks needs a new model or a special package. Thin the data and watch whether the estimate holds; look at the season edges for a censoring pile-up; compute the trend with more than one summary and see whether they agree. A phenological shift that survives all three is worth reporting. One that fails any of them is a statement about the sampling, dressed up as a statement about the season. Every summary is a decision, and these three checks make the decision defensible; the checking a circular analysis tutorial applies the same spirit to directional and time-of-day data.

References

Moussus J-P, Julliard R, Jiguet F 2010. Methods in Ecology and Evolution 1(2):140-150 (10.1111/j.2041-210X.2010.00020.x)

CaraDonna PJ, Iler AM, Inouye DW 2014. Proceedings of the National Academy of Sciences 111(13):4916-4921 (10.1073/pnas.1323073111)

Miller-Rushing AJ, Inouye DW, Primack RB 2008. Journal of Ecology 96(6):1289-1296 (10.1111/j.1365-2745.2008.01436.x)

Related tutorials

• First flowering dates and sampling effort
• Phenological trends and temperature
• Phenology in R: day of year and event timing
• Checking a circular analysis