Getting started with SporeLag

SporeLag turns a daily exposure series (pollen or spore counts, ozone, particulate matter) into lagged and moving-average exposure variables for epidemiological models. This vignette walks through the full pipeline on the bundled pollen_demo data and, along the way, shows the ways the pipeline can go wrong and how SporeLag stops them from going wrong silently.

The data: gaps are not missing values

pollen_demo is a synthetic dataset of daily pollen counts at two monitoring sites over one spring season. It is not real surveillance data.

str(pollen_demo)
#> 'data.frame':    266 obs. of  3 variables:
#>  $ site : chr  "North" "North" "North" "North" ...
#>  $ date : Date, format: "2024-02-15" "2024-02-16" ...
#>  $ count: int  NA 98 134 51 41 210 48 93 251 338 ...

The season runs from 2024-02-15 to 2024-06-30, which is 137 days. Neither site has 137 rows:

table(pollen_demo$site)
#> 
#> North South 
#>   132   134

The data contain two different defects, and it is worth being precise about the difference:

a missing value is a row whose count is NA, for example a sample that failed on a day the station was running;
a gap is a day with no row at all, for example a week when the sampler was offline.

is.na() sees the first kind and is blind to the second:

tapply(is.na(pollen_demo$count), pollen_demo$site, sum)
#> North South 
#>     8     8

Here is the North site around its gap. The rows jump from 17 March straight to 23 March:

north <- pollen_demo[pollen_demo$site == "North", ]
north[north$date >= as.Date("2024-03-15") & north$date <= as.Date("2024-03-25"), ]
#>     site       date count
#> 30 North 2024-03-15   224
#> 31 North 2024-03-16   388
#> 32 North 2024-03-17   759
#> 33 North 2024-03-23    NA
#> 34 North 2024-03-24   274
#> 35 North 2024-03-25   510

Lags on a gapped series fail loudly

A lag computed by row position on this series would be wrong. “Yesterday” for 23 March would be 17 March, six days earlier. The resulting column would look plausible, run cleanly through a regression, and give a misaligned lag-response estimate that nothing downstream would flag.

For that reason the two temporal functions, apply_lag() and build_moving_average(), check that each group is on a complete daily grid and refuse to run if it is not:

apply_lag(pollen_demo, value = "count", lags = 1, date = "date", by = "site")
#> Error in `apply_lag()`:
#> ! Date sequence is not a complete daily grid.
#> ✖ Found 8 missing days across 2 groups.
#> ℹ Lagging a gapped series shifts by position, not by time, which silently
#>   misaligns exposure and outcome.
#> → Call `complete_daily_grid()` first, then re-run.

The error has class sporelag_error_gaps, so code that wraps SporeLag can catch it specifically rather than matching on the message text:

tryCatch(
  apply_lag(pollen_demo, value = "count", lags = 1, date = "date", by = "site"),
  sporelag_error_gaps = function(e) "gapped grid detected"
)
#> [1] "gapped grid detected"

There is deliberately no fill_gaps = TRUE argument. Closing gaps is a separate, visible step.

Step 1: complete the daily grid

complete_daily_grid() inserts a row for every absent calendar day within each group. Inserted rows carry NA in every column other than the date and the grouping columns; nothing is guessed.

grid <- complete_daily_grid(pollen_demo, date = "date", by = "site")
table(grid$site)
#> 
#> North South 
#>   137   137
sum(is.na(grid$count))
#> [1] 24

Each site now has all 137 days. The gaps have become ordinary missing values, which is a problem we can reason about explicitly.

Step 2: calendar features

ISO weeks

assign_iso_week() appends the ISO 8601 week and the ISO year. Both are needed, because the ISO year of a date near New Year is often not its calendar year:

new_year <- data.frame(date = as.Date(c("2024-12-29", "2024-12-30", "2025-01-01")))
assign_iso_week(new_year, date = "date")
#>         date iso_week iso_year
#> 1 2024-12-29       52     2024
#> 2 2024-12-30        1     2025
#> 3 2025-01-01        1     2025

30 December 2024 belongs to week 1 of ISO year 2025. Grouping a multi-year series by iso_week alone would pool unrelated days from different years, so always group by both columns. The weeks are computed in base R, giving identical results on every operating system.

Seasons

assign_season() appends a season label as an ordered factor. The default is the meteorological calendar:

grid <- grid |>
  assign_iso_week(date = "date") |>
  assign_season(date = "date")

table(grid$season, grid$site)
#>         
#>          North South
#>   Winter    15    15
#>   Spring    92    92
#>   Summer    30    30
#>   Autumn     0     0

A calendar season is not a pollen season. Pollen seasons are taxon- and region-specific, so for exposure work you will usually want to supply your own season start dates with definition = "custom". The dates below are illustrative, not a recommendation:

custom <- assign_season(
  grid[, c("site", "date")],
  date = "date",
  definition = "custom",
  breaks = c(Dormant = "11-01", Tree = "02-15", Grass = "05-01", Weed = "08-01")
)
table(custom$season, custom$site)
#>          
#>           North South
#>   Tree       76    76
#>   Grass      61    61
#>   Weed        0     0
#>   Dormant     0     0

Step 3: impute missing values, and keep track of what you imputed

impute_weekly_mean() fills each missing day with the mean of the observed days in the same ISO week, within site. It never overwrites the input column. Instead, it appends the completed series and a logical flag:

grid <- impute_weekly_mean(grid, value = "count", by = "site")

grid[grid$site == "North" &
       grid$date >= as.Date("2024-03-16") & grid$date <= as.Date("2024-03-25"),
     c("date", "iso_week", "count", "count_imputed", "count_imputed_flag")]
#>          date iso_week count count_imputed count_imputed_flag
#> 31 2024-03-16       11   388           388              FALSE
#> 32 2024-03-17       11   759           759              FALSE
#> 33 2024-03-18       12    NA           274               TRUE
#> 34 2024-03-19       12    NA           274               TRUE
#> 35 2024-03-20       12    NA           274               TRUE
#> 36 2024-03-21       12    NA           274               TRUE
#> 37 2024-03-22       12    NA           274               TRUE
#> 38 2024-03-23       12    NA           274               TRUE
#> 39 2024-03-24       12   274           274              FALSE
#> 40 2024-03-25       13   510           510              FALSE

This output exposes a real weakness of the default. North’s offline period fills almost all of ISO week 12, so six days have been filled from a single observation (24 March). The min_obs argument sets how many observed days a week needs before its mean is used:

strict <- impute_weekly_mean(
  grid[, c("site", "date", "count", "iso_week", "iso_year")],
  value = "count", by = "site", min_obs = 4
)
sum(strict$count_imputed_flag)
#> [1] 12
sum(is.na(strict$count_imputed))
#> [1] 12

With min_obs = 4, weeks with fewer than four observed days are left NA rather than extrapolated.

Whatever threshold you choose, imputation replaces day-to-day variation with a constant. That shrinks the variance of the exposure and tends to bias effect estimates toward the null. Use the flag column to report the proportion imputed and to run a complete-case sensitivity analysis:

tapply(grid$count_imputed_flag, grid$site, mean)
#>      North      South 
#> 0.09489051 0.08029197

Step 4: moving averages

build_moving_average() appends one column per window. By default the window is trailing and includes the current day: window = 3 averages today and the two previous days.

grid <- build_moving_average(
  grid, value = "count_imputed", window = c(3, 7), date = "date", by = "site"
)

Three defaults are worth knowing, because each changes what the variable means:

Edge windows are NA. The first window - 1 days of each site have no complete window, so a 7-day mean is NA for the first six days rather than being computed from fewer days.
min_obs = window. Any NA inside a window makes the mean NA. You can lower it to tolerate missing days, but that has the same variance-shrinking effect as imputation.
align = "center" uses future exposure. It is useful for smoothing a plot, but it should not be used to build an exposure for an outcome measured on the current day.

To see what min_obs does, apply a 7-day window to the raw (non-imputed) counts:

raw_ma <- build_moving_average(grid[, c("site", "date", "count")],
                               value = "count", window = 7,
                               date = "date", by = "site")
sum(is.na(raw_ma$count_ma7))
#> [1] 106

raw_ma5 <- build_moving_average(grid[, c("site", "date", "count")],
                                value = "count", window = 7, min_obs = 5,
                                date = "date", by = "site")
sum(is.na(raw_ma5$count_ma7))
#> [1] 30

Step 5: lags

apply_lag() appends one column per lag. Lag n places the value from n days earlier alongside the current day. lags = 0 gives a copy of the same-day value, so you can build a uniform lag0, lag1, ... model matrix. Negative lags are an error, since a negative lag would pair an outcome with exposure measured after it.

grid <- apply_lag(
  grid, value = "count_imputed", lags = 0:3, date = "date", by = "site"
)

Lags are computed strictly within each site. The first rows of the South series are NA, not the last values of the North series:

south <- grid[grid$site == "South", ]
head(south[, c("site", "date", "count_imputed",
               "count_imputed_lag1", "count_imputed_lag3")], 4)
#>      site       date count_imputed count_imputed_lag1 count_imputed_lag3
#> 138 South 2024-02-15      53.00000                 NA                 NA
#> 139 South 2024-02-16      55.00000                 53                 NA
#> 140 South 2024-02-17      28.00000                 55                 NA
#> 141 South 2024-02-18      45.33333                 28                 53

The whole pipeline

Each function takes a data frame and returns the same data frame with new columns appended, so the steps chain naturally:

model_ready <- pollen_demo |>
  complete_daily_grid(date = "date", by = "site") |>
  assign_iso_week(date = "date") |>
  assign_season(date = "date") |>
  impute_weekly_mean(value = "count", by = "site") |>
  build_moving_average(value = "count_imputed", window = c(3, 7),
                       date = "date", by = "site") |>
  apply_lag(value = "count_imputed", lags = 0:3,
            date = "date", by = "site")

dim(model_ready)
#> [1] 274  14
names(model_ready)
#>  [1] "site"               "date"               "count"             
#>  [4] "iso_week"           "iso_year"           "season"            
#>  [7] "count_imputed"      "count_imputed_flag" "count_imputed_ma3" 
#> [10] "count_imputed_ma7"  "count_imputed_lag0" "count_imputed_lag1"
#> [13] "count_imputed_lag2" "count_imputed_lag3"

The input columns site, date, and count are unchanged; everything else was appended. If you leave out complete_daily_grid(), the last two steps raise an error instead of returning misaligned lags.