From 7cc447892f501be6d78e3676feb5a84363ee3fd6 Mon Sep 17 00:00:00 2001 From: "Daniel J. McDonald" Date: Wed, 13 May 2026 14:08:45 -0700 Subject: [PATCH 1/6] upd: tvar delay canada, draft inc. --- vignettes/articles/delay-distributions.Rmd | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/vignettes/articles/delay-distributions.Rmd b/vignettes/articles/delay-distributions.Rmd index 351a04d..8d9aa09 100644 --- a/vignettes/articles/delay-distributions.Rmd +++ b/vignettes/articles/delay-distributions.Rmd @@ -160,7 +160,7 @@ ggplot(cdfs |> pivot_longer(-x), aes(x, value, color = name)) + Finally, we also allow time-varying delay distributions. This is accomplished with a matrix. This requires a bit of work, but is not too challenging. For example, to create the correct matrix using the Baker et al. delays, the -necessary code is the following. +necessary code is the following (note that this is _not_ yet time varying): ```{r backer-matrix, eval=FALSE} # library(Matrix) @@ -396,10 +396,25 @@ delay_dstns_can |> ``` -We use these delays during the period they were most prevalent to -estimate Rt. First, we build the delay matrix. +We use these delays weighted by the probabilities of the circulating variants to +estimate $R_t$. First, we build the delay matrix. ```{r tvar-matrix, message=FALSE} +can_props_voc <- can_props_smoothed |> + filter(Date <= max(cancovid$date)) |> + transmute( + Date = Date, + `Ancestral lineage` = other, + Alpha = Alpha + Beta, # Beta didn't really circulate + Delta = Delta, + Omicron = 1 - (`Ancestral lineage` + Alpha + Delta) + ) + +can_tvar_delay <- can_props_voc |> + pivot_longer(-Date) |> + left_join(delay_dstns_can |> select(name = type, shape, scale), by = "name") |> + + n <- nrow(cancovid) delay_mat <- matrix(0, n, n) for (iter in 1:n) { From f530e61b4c0fa9256bdc45cff1a24e83608b28b9 Mon Sep 17 00:00:00 2001 From: "Daniel J. McDonald" Date: Thu, 14 May 2026 10:15:02 -0700 Subject: [PATCH 2/6] upd: news to reflect tvar delay adjustment --- NEWS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 387764b..be84275 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,7 @@ # rtestim (development version) -# rtestim 1.0.2 +* Correct minor issue in the construction of the time-varying delay example +in `vignettes/articles/delay-distributions.Rmd` # rtestim 1.0.2 From 543c8a3af9dccfd9273459376ead38d0a7b6ff17 Mon Sep 17 00:00:00 2001 From: "Daniel J. McDonald" Date: Thu, 14 May 2026 10:24:29 -0700 Subject: [PATCH 3/6] fix: vignette code corrected --- vignettes/articles/delay-distributions.Rmd | 333 ++++++++++----------- 1 file changed, 166 insertions(+), 167 deletions(-) diff --git a/vignettes/articles/delay-distributions.Rmd b/vignettes/articles/delay-distributions.Rmd index 8d9aa09..71a520c 100644 --- a/vignettes/articles/delay-distributions.Rmd +++ b/vignettes/articles/delay-distributions.Rmd @@ -21,7 +21,6 @@ theme_set(theme_bw()) \newcommand{\argmin}{\mathop{\mathrm{argmin}}} \newcommand{\Rtseq}{\mathop{\{R_t\}_{t=1}^n}} - ```{r, include = FALSE} knitr::opts_chunk$set( dpi = 300, @@ -34,79 +33,56 @@ knitr::opts_chunk$set( ) ``` -This package accommodates 3 different ways of specifying a delay distribution -in the renewal equation, and this vignette illustrates these in turn. The -renewal equation specifies that the expected incidence at time $t$ is a -weighted sum of previous incidence, multiplied by $R_t$: -\begin{equation} +This package accommodates 3 different ways of specifying a delay distribution in the renewal equation, and this vignette illustrates these in turn. The renewal equation specifies that the expected incidence at time $t$ is a weighted sum of previous incidence, multiplied by $R_t$: \begin{equation} \mathbb{E}[I_t \mid I_1,\ldots,I_{t-1}] = R_t \sum_{a=0}^t w_a I_{t-a}, -\end{equation} -calculated by convolving the -preceding $a$ days of new infections with the discretized generation (or serial) -interval distribution $w$ of length $t$. It is the allowable distributions -$\{w\}_{a=0}^n$ that are the focus. +\end{equation} calculated by convolving the preceding $a$ days of new infections with the discretized generation (or serial) interval distribution $w$ of length $t$. It is the allowable distributions $\{w\}_{a=0}^n$ that are the focus. ## Discretization -First, it is important to recognize that the sequence $\{w\}_{a=0}^n$ is usually -a discretization of a probability density function, most frequently gamma or -Weibull. This is because observed incidence happens at discrete time points -like days or weeks, while time is continuous. Using the default (parametric) -delay distributions requires calculating a discrete approximation. +First, it is important to recognize that the sequence $\{w\}_{a=0}^n$ is usually a discretization of a probability density function, most frequently gamma or Weibull. This is because observed incidence happens at discrete time points like days or weeks, while time is continuous. Using the default (parametric) delay distributions requires calculating a discrete approximation. ## Default (parametric) delay distribution -By default, `estimate_rt()` uses a gamma distribution parameterized by the -shape $k$ and scale $\theta$. This density has pdf -\[ -f_W(w) = \frac{1}{\Gamma(k)\theta^k} w^{k-1} e^{-w/\theta} I(w > 0), -\] -where $I$ is the indicator function. The mean of this distribution is $k\theta$, -and the variance is $k\theta^2$. Both $k$ and $\theta$ must be greater than 0. -The figure below shows a few examples densities from this family. +By default, `estimate_rt()` uses a gamma distribution parameterized by the shape $k$ and scale $\theta$. This density has pdf \[ f_W(w) = \frac{1}{\Gamma(k)\theta^k} w\^{k-1} e\^{-w/\theta} I(w \> 0), \] where $I$ is the indicator function. The mean of this distribution is $k\theta$, and the variance is $k\theta^2$. Both $k$ and $\theta$ must be greater than 0. The figure below shows a few examples densities from this family. ```{r gamma-pdf, echo = FALSE} ggplot(data.frame(x = seq(.001, 20, length.out = 100)), aes(x)) + - stat_function(fun = dgamma, args = list(shape = 2.5, scale = 2.5), color = "orange") + - stat_function(fun = dgamma, args = list(shape = 1, scale = 2), color = "cornflowerblue") + - stat_function(fun = dgamma, args = list(shape = 2, scale = 1), color = "forestgreen") + + stat_function( + fun = dgamma, + args = list(shape = 2.5, scale = 2.5), + color = "orange" + ) + + stat_function( + fun = dgamma, + args = list(shape = 1, scale = 2), + color = "cornflowerblue" + ) + + stat_function( + fun = dgamma, + args = list(shape = 2, scale = 1), + color = "forestgreen" + ) + scale_y_continuous(expand = expansion(c(0, 0.05))) ``` -The default is the orange curve shown above. Given an incidence sequence of -length $n$, internally, `delay_calculator()` creates $\{w\}_{a=0}^n$ by -\[ -w_a = F(a + 1) - F(a - 1), -\] -and then renormalizing with `w / sum(w)` to ensure that it sums to 1. Here, $F$ -is the cumulative distribution function for Gamma (though similar works for any -continuous distribution). Note that this formula assumes that the probability -of a 0 delay is _not_ 0. This is related to the phenomenon of "double interval -censoring": a primary infection on day $t$ may have occurred any time between -0:00 and 23:59 and the same for a secondary infection on day $t+1$. So even if -the gap in reported days is only "1 day", those infections may have happened as -little as a few minutes apart, or as much as 48 hours apart. - -When using this default, the delay distribution is necessarily assumed to be the same for all $t$, and we -compute the convolution of `w` and `I` with the Fast Fourier Transform. - -Just to illustrate this behaviour, we show the results for the default setting -on the included `cancovid` data. +The default is the orange curve shown above. Given an incidence sequence of length $n$, internally, `delay_calculator()` creates $\{w\}_{a=0}^n$ by \[ w_a = F(a + 1) - F(a - 1), \] and then renormalizing with `w / sum(w)` to ensure that it sums to 1. Here, $F$ is the cumulative distribution function for Gamma (though similar works for any continuous distribution). Note that this formula assumes that the probability of a 0 delay is *not* 0. This is related to the phenomenon of "double interval censoring": a primary infection on day $t$ may have occurred any time between 0:00 and 23:59 and the same for a secondary infection on day $t+1$. So even if the gap in reported days is only "1 day", those infections may have happened as little as a few minutes apart, or as much as 48 hours apart. + +When using this default, the delay distribution is necessarily assumed to be the same for all $t$, and we compute the convolution of `w` and `I` with the Fast Fourier Transform. + +Just to illustrate this behaviour, we show the results for the default setting on the included `cancovid` data. ```{r can-default} -can_default <- estimate_rt(cancovid$incident_cases, x = cancovid$date, nsol = 20L) +can_default <- estimate_rt( + cancovid$incident_cases, + x = cancovid$date, + nsol = 20L +) plot(can_default) + coord_cartesian(ylim = c(0.5, 2)) ``` ## Constant, non-parametric delay distribution -If we don't believe that the Gamma distribution closely approximates the serial -interval distribution, then we can specify our own distribution manually. For example, -[Backer et al., Table S1](https://doi.org/10.2807/1560-7917.ES.2022.27.6.2200042) -gives observed serial intervals for the Omicron (SGTF) and Delta (non-SGTF) -COVID-19 variants during 2 weeks in 2021 in the Netherlands. For the sake of -illustration, we aggregate these together and use this as our (constant) delay -distribution. +If we don't believe that the Gamma distribution closely approximates the serial interval distribution, then we can specify our own distribution manually. For example, [Backer et al., Table S1](https://doi.org/10.2807/1560-7917.ES.2022.27.6.2200042) gives observed serial intervals for the Omicron (SGTF) and Delta (non-SGTF) COVID-19 variants during 2 weeks in 2021 in the Netherlands. For the sake of illustration, we aggregate these together and use this as our (constant) delay distribution. ```{r cancov-nonpar} # Data from Backer et al. @@ -119,7 +95,7 @@ delay <- delay / sum(delay) ```{r cancov-nonpar-plot, echo=FALSE} ggplot( - data.frame(delay = 1:length(delay), probability = delay), + data.frame(delay = 1:length(delay), probability = delay), aes(delay, probability) ) + geom_point(colour = "cornflowerblue") + @@ -127,21 +103,19 @@ ggplot( scale_y_continuous(expand = expansion(c(0, 0.05))) ``` -This distribution looks something like a gamma, but it has finite support. -We can easily use it instead. +This distribution looks something like a gamma, but it has finite support. We can easily use it instead. ```{r can-nonpar} can_nonpar <- estimate_rt( - cancovid$incident_cases, - x = cancovid$date, + cancovid$incident_cases, + x = cancovid$date, delay_distn = delay, - nsol = 20L) + nsol = 20L +) plot(can_nonpar) + coord_cartesian(ylim = c(0.5, 2)) ``` -The result is much less dramatic than the previous version. This is likely -because the distribution is much more concentrated near short delays. We can see -this by examining the two CDFs. +The result is much less dramatic than the previous version. This is likely because the distribution is much more concentrated near short delays. We can see this by examining the two CDFs. ```{r cdfs, echo = FALSE} cdfs <- data.frame( @@ -151,52 +125,65 @@ cdfs <- data.frame( ) ggplot(cdfs |> pivot_longer(-x), aes(x, value, color = name)) + geom_step() + - scale_y_continuous(expand = expansion(0), limits = c(0, 1)) + + coord_cartesian(ylim = c(0, 1)) + + scale_y_continuous(expand = expansion(0)) + scale_color_manual(values = c("darkorange", "cornflowerblue"), name = "") ``` ## Time-varying delays -Finally, we also allow time-varying delay distributions. This is accomplished -with a matrix. This requires a bit of work, but is not too challenging. For -example, to create the correct matrix using the Baker et al. delays, the -necessary code is the following (note that this is _not_ yet time varying): +Finally, we also allow time-varying delay distributions. This is accomplished with a matrix. This requires a bit of work, but is not too challenging. For example, to create the correct matrix using the Baker et al. delays, the necessary code is the following (note that this is *not* yet time varying): ```{r backer-matrix, eval=FALSE} # library(Matrix) n <- nrow(cancovid) backer_delay <- c(delay, rep(0, n - length(delay) - 1)) delay_mat <- matrix(0, n, n) -for (iter in 1:n) delay_mat[iter, 1:iter] <- rev(backer_delay[1:iter]) +for (iter in 1:n) { + delay_mat[iter, 1:iter] <- rev(backer_delay[1:iter]) +} delay_mat <- drop0(as(delay_mat, "CsparseMatrix")) # make it sparse, not necessary delay_mat <- delay_mat / rowSums(delay_mat) # renormalize ``` -We could then simply pass `delay_mat` to `estimate_rt(..., delay_distn = delay_mat)`. -The result would be the same as that shown above. +We could then simply pass `delay_mat` to `estimate_rt(..., delay_distn = delay_mat)`. The result would be the same as that shown above. ## Variant-specific delays -For a more illustrative version of the previous example, we revisit the Canadian -data, but combine it with some additional information. +For a more illustrative version of the previous example, we revisit the Canadian data, but combine it with some additional information. -We use variant -circulation information from [CoVaRR-Net's Duotang notebook](https://covarr-net.github.io/duotang/duotang.html). -For simplicity, the code for this processing is hidden in the fold. +We use variant circulation information from [CoVaRR-Net's Duotang notebook](https://covarr-net.github.io/duotang/duotang.html). For simplicity, the code for this processing is hidden in the fold.
+ ```{r duotang-processing, eval=FALSE, echo=TRUE} # Run on 19 April 2024 -duotang <- read_tsv("https://github.com/CoVaRR-NET/duotang/raw/main/data_needed/virusseq.metadata.csv.gz") +duotang <- read_tsv( + "https://github.com/CoVaRR-NET/duotang/raw/main/data_needed/virusseq.metadata.csv.gz" +) columnlist <- c( - "fasta_header_name", "province", "host_gender", "host_age_bin", - "sample_collection_date", "sample_collected_by", - "purpose_of_sampling", "purpose_of_sequencing", "lineage", - "raw_lineage", "gisaid_accession", "isolate" + "fasta_header_name", + "province", + "host_gender", + "host_age_bin", + "sample_collection_date", + "sample_collected_by", + "purpose_of_sampling", + "purpose_of_sequencing", + "lineage", + "raw_lineage", + "gisaid_accession", + "isolate" ) unknown.str <- c( - "Undeclared", "Not Provided", "Restricted Access", "Missing", - "Not Applicable", "", "NA", "unknow" + "Undeclared", + "Not Provided", + "Restricted Access", + "Missing", + "Not Applicable", + "", + "NA", + "unknow" ) duotang <- duotang |> rename(province = geo_loc_name_state_province_territory) |> @@ -212,21 +199,35 @@ meta <- meta |> mutate(gisaid_accession = str_replace(gisaid_accession, "EPI_ISL_", "")) |> rename(GID = gisaid_accession) |> rowwise() |> - mutate(raw_lineage = ifelse( - grepl("^X", raw_lineage), - str_replace_all(paste0( - realtorawlineage(substr( - raw_lineage, 1, str_locate(raw_lineage, "\\.") - 1 - )), - ".", - substr(raw_lineage, str_locate(raw_lineage, "\\.") + 1, nchar(raw_lineage)) - ), "[\r\n]", ""), - raw_lineage - )) |> + mutate( + raw_lineage = ifelse( + grepl("^X", raw_lineage), + str_replace_all( + paste0( + realtorawlineage(substr( + raw_lineage, + 1, + str_locate(raw_lineage, "\\.") - 1 + )), + ".", + substr( + raw_lineage, + str_locate(raw_lineage, "\\.") + 1, + nchar(raw_lineage) + ) + ), + "[\r\n]", + "" + ), + raw_lineage + ) + ) |> ungroup() dico <- makepangolindico() # rebuild the lineage dictionary so the correct names gets assigned for XBB descedants not named XBB -VOCVOI <- read_csv("https://raw.githubusercontent.com/CoVaRR-NET/duotang/main/resources/vocvoi.csv") +VOCVOI <- read_csv( + "https://raw.githubusercontent.com/CoVaRR-NET/duotang/main/resources/vocvoi.csv" +) meta$pango_group <- create.pango.group(VOCVOI, meta) meta <- select(meta, province, week, pango_group) |> mutate(week = as.Date(week)) @@ -243,10 +244,10 @@ can_counts <- group_by(meta, week, pango_group) |> counts <- bind_rows(can_counts, counts) saveRDS(counts, "duotang-counts.rds") ``` +
-We smooth the raw data using multinomial logistic regression on a third order orthogonal polynomial to produce the following estimated variant proportions -in Canada. +We smooth the raw data using multinomial logistic regression on a third order orthogonal polynomial to produce the following estimated variant proportions in Canada. ```{r duotang-counts, echo=FALSE, message=FALSE} props <- readRDS("duotang-counts.rds") |> @@ -260,7 +261,7 @@ smooth_it <- function(props_group) { n <- names(z) nn <- gsub(" ", "_", n) names(z) <- nn - form_resp <- paste0("cbind(", paste0(names(z), collapse = ",") ,") ~ ") + form_resp <- paste0("cbind(", paste0(names(z), collapse = ","), ") ~ ") z$time <- as.numeric(props_group$week) form <- as.formula(paste0(form_resp, "poly(time, degree = 3)")) fits <- multinom(form, z, trace = FALSE) @@ -278,38 +279,42 @@ can_props_smoothed |> ggplot(aes(Date, y = value, fill = name)) + geom_area(position = "stack") + ylab("Variants in circulation") + - xlab("") + + xlab("") + theme_bw() + - scale_x_date(name = "", date_breaks = "1 year", date_labels = "%Y", expand = c(0,0)) + - scale_y_continuous(expand = c(0,0)) + + scale_x_date( + name = "", + date_breaks = "1 year", + date_labels = "%Y", + expand = c(0, 0) + ) + + scale_y_continuous(expand = c(0, 0)) + scale_fill_viridis_d(name = "") can_pred_class <- can_props_smoothed |> filter(Date <= max(cancovid$date)) |> - pivot_longer(-Date) |> - group_by(Date) |> - summarise(var = name[which.max(value)], .groups = "drop") |> - mutate(var = case_when( - var == "other" ~ "Ancestral lineage", - var == "Alpha" ~ "Alpha", - var == "Beta" ~ "Beta", - var == "Delta" ~ "Delta", - TRUE ~ "Omicron" - )) - -boot <- data.frame( - Date = seq(min(cancovid$date), min(can_props_smoothed$Date) - 1, by = "day"), - var = "Ancestral lineage" -) -can_pred_class <- bind_rows(boot, can_pred_class) + pivot_longer(-Date) |> + summarise(variant = name[which.max(value)], .by = Date) |> + mutate( + variant = case_when( + variant == "other" ~ "Ancestral lineage", + variant == "Alpha" ~ "Alpha", + variant == "Beta" ~ "Beta", + variant == "Delta" ~ "Delta", + TRUE ~ "Omicron" + ) + ) |> + right_join(cancovid |> select(date), by = join_by(Date == date)) |> + arrange(Date) |> + fill(variant, .direction = "downup") ``` -Using the estimated proportions, we label each date with the dominant variant -at the time (and restrict ourselves to the time period of the included case -data). +Using the estimated proportions, we label each date with the dominant variant at the time (and restrict ourselves to the time period of the included case data). ```{r, echo=FALSE, fig.height=1, dev='png'} -ggplot(can_pred_class, aes(x = Date, fill = fct_relevel(var, "Ancestral lineage"))) + +ggplot( + can_pred_class, + aes(x = Date, fill = fct_relevel(variant, "Ancestral lineage")) +) + geom_ribbon(aes(ymax = 1, ymin = 0)) + scale_y_continuous(expand = expansion(0)) + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) + @@ -318,20 +323,15 @@ ggplot(can_pred_class, aes(x = Date, fill = fct_relevel(var, "Ancestral lineage" ## Creating delay distributions -We use the delay distribution data from the meta analysis in -[Xu, X., Wu, Y., Kummer, A.G. et al.](https://doi.org/10.1186/s12916-023-03070-8). -The raw data for their analysis is available on [GitHub](https://github.com/xxyy0574/COVID-19-transmission-parameters). We use -a simple version of their procedure, taking the median calculations across -studies to find the mean and standard deviation of the delay for each variant -separately. Then we convert these to the $k$ and $\theta$ parameters for the gamma distribution. -This analysis is hidden below the fold. +We use the delay distribution data from the meta analysis in [Xu, X., Wu, Y., Kummer, A.G. et al.](https://doi.org/10.1186/s12916-023-03070-8). The raw data for their analysis is available on [GitHub](https://github.com/xxyy0574/COVID-19-transmission-parameters). We use a simple version of their procedure, taking the median calculations across studies to find the mean and standard deviation of the delay for each variant separately. Then we convert these to the $k$ and $\theta$ parameters for the gamma distribution. This analysis is hidden below the fold.
+ ```{r meta-analysis, eval=FALSE} data_raw <- readxl::read_excel("xu-etal-DATA_RAW.xlsx") |> select(type, para, n = Sample_size, mean, sd, se, median) |> filter(!is.na(type)) |> - mutate(across(-c(type, para), as.numeric)) + mutate(across(-c(type, para), as.numeric)) bonehead_meta <- data_raw |> group_by(type, para) |> @@ -350,16 +350,17 @@ bonehead_meta <- data_raw |> sd = median(sd, na.rm = TRUE), .groups = "drop" ) -## There's only one Beta and only IP. -## We use the corresponding sd for Alpha IP, -## and duplicate Alpha for GT / ST +# There's only one Beta and only IP. +# We use the corresponding sd for Alpha IP, +# and duplicate Alpha for GT / ST Beta_IP <- bonehead_meta |> filter(type == "Beta") -Beta_IP$sd = bonehead_meta |> - filter(type == "Alpha", para == "IP") |> pull(sd) +Beta_IP$sd = bonehead_meta |> + filter(type == "Alpha", para == "IP") |> + pull(sd) Beta <- bind_rows( Beta_IP, - bonehead_meta |> + bonehead_meta |> filter(type == "Alpha", para != "IP") |> mutate(type = "Beta") ) @@ -371,15 +372,14 @@ delay_dstns_byvar <- bonehead_meta |> mutate(shape = mean^2 / sd^2, scale = mean / shape) saveRDS(delay_dstns_byvar, "delay-distns-byvar.rds") ``` +
-We'll use just the Serial Intervals (SI), though the incubation periods (IP), -and generation times (GT) are also available. Below, we visualize the estimated -delays for the 4 variants that were most prevalent in Canada. +We'll use just the Serial Intervals (SI), though the incubation periods (IP), and generation times (GT) are also available. Below, we visualize the estimated delays for the 4 variants that were most prevalent in Canada. ```{r variant-si-plot, echo=FALSE} delay_dstns_can <- readRDS("delay-distns-byvar.rds") |> - filter(para == "SI", type %in% unique(can_pred_class$var)) + filter(para == "SI", type %in% unique(can_pred_class$variant)) delay_dstns_can |> rowwise() |> @@ -390,40 +390,42 @@ delay_dstns_can |> group_by(type) |> mutate(delay = row_number() - 1) |> ggplot(aes(delay, probability, colour = type)) + - geom_line() + + geom_line() + scale_color_brewer(palette = "Dark2", name = "") + scale_y_continuous(expand = expansion(c(0, 0.05))) ``` - -We use these delays weighted by the probabilities of the circulating variants to -estimate $R_t$. First, we build the delay matrix. +We use these delays weighted by the probabilities of the circulating variants to estimate $R_t$. First, we build the delay matrix. ```{r tvar-matrix, message=FALSE} can_props_voc <- can_props_smoothed |> filter(Date <= max(cancovid$date)) |> transmute( - Date = Date, - `Ancestral lineage` = other, - Alpha = Alpha + Beta, # Beta didn't really circulate - Delta = Delta, + Date = Date, + `Ancestral lineage` = other, + Alpha = Alpha + Beta, # Beta didn't really circulate + Delta = Delta, Omicron = 1 - (`Ancestral lineage` + Alpha + Delta) ) -can_tvar_delay <- can_props_voc |> - pivot_longer(-Date) |> - left_join(delay_dstns_can |> select(name = type, shape, scale), by = "name") |> - +# The delays are truncated to have probability 0 after 30 days +can_lagd_delays <- delay_dstns_can |> + rowwise() |> + transmute(name = type, delay = list(discretize_gamma(0:30, shape, scale))) |> + ungroup() -n <- nrow(cancovid) -delay_mat <- matrix(0, n, n) -for (iter in 1:n) { - current_var <- can_pred_class$var[iter] - current_pars <- delay_dstns_can |> filter(type == current_var) - delay <- discretize_gamma(0:(iter - 1), current_pars$shape, current_pars$scale) - delay_mat[iter, 1:iter] <- rev(delay) -} -delay_mat <- drop0(as(delay_mat, "CsparseMatrix")) # make it sparse, not necessary +can_tvar_delay <- can_props_voc |> + pivot_longer(-Date, values_to = "prop") |> + left_join(can_lagd_delays, by = "name") |> + summarise(delay = list(colSums(prop * do.call(rbind, delay))), .by = Date) |> + right_join(cancovid, by = join_by(Date == date)) |> + arrange(Date) |> + fill(delay, .direction = "downup") |> + pull(delay) + +n <- length(can_tvar_delay) +narrow_delay_mat <- do.call(rbind, can_tvar_delay) # n x 31 matrix +delay_mat <- t(bandSparse(n, n, k = 0:30, diagonals = narrow_delay_mat)) delay_mat <- delay_mat / rowSums(delay_mat) # renormalize ``` @@ -431,13 +433,10 @@ Finally, we use this time-varying delay matrix to estimate Rt. ```{r tvar-Rt} can_tvar <- estimate_rt( - cancovid$incident_cases, - x = cancovid$date, + cancovid$incident_cases, + x = cancovid$date, delay_distn = delay_mat, - nsol = 20L) + nsol = 20L +) plot(can_tvar) + coord_cartesian(ylim = c(0.5, 2)) ``` - - - - From 2535382f613e99f7cdcdc462c46568ff8864df44 Mon Sep 17 00:00:00 2001 From: "Daniel J. McDonald" Date: Thu, 14 May 2026 11:28:45 -0700 Subject: [PATCH 4/6] minor textural corrections --- vignettes/articles/delay-distributions.Rmd | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/vignettes/articles/delay-distributions.Rmd b/vignettes/articles/delay-distributions.Rmd index 71a520c..1a7cc96 100644 --- a/vignettes/articles/delay-distributions.Rmd +++ b/vignettes/articles/delay-distributions.Rmd @@ -137,12 +137,18 @@ Finally, we also allow time-varying delay distributions. This is accomplished wi ```{r backer-matrix, eval=FALSE} # library(Matrix) n <- nrow(cancovid) -backer_delay <- c(delay, rep(0, n - length(delay) - 1)) -delay_mat <- matrix(0, n, n) -for (iter in 1:n) { - delay_mat[iter, 1:iter] <- rev(backer_delay[1:iter]) -} -delay_mat <- drop0(as(delay_mat, "CsparseMatrix")) # make it sparse, not necessary +dd <- length(delay) +# one line, optimized version +delay_mat <- drop0(toeplitz2(c(rep(0, n - 1), delay, rep(0, n - dd)), n, n)) + +# illustration of doing the same with a loop +# backer_delay <- c(delay, rep(0, n - dd)) +# delay_mat <- matrix(0, n, n) +# for (iter in 1:n) { +# delay_mat[iter, 1:iter] <- rev(backer_delay[1:iter]) +# } +# delay_mat <- drop0(as(delay_mat, "CsparseMatrix")) # sparseify, unnecessary + delay_mat <- delay_mat / rowSums(delay_mat) # renormalize ``` From a2b610af7687345f016504eee7394a4f17c597ed Mon Sep 17 00:00:00 2001 From: "Daniel J. McDonald" Date: Thu, 14 May 2026 11:59:04 -0700 Subject: [PATCH 5/6] fix: annoying line wrapping --- vignettes/articles/delay-distributions.Rmd | 100 ++++++++++++++++----- 1 file changed, 79 insertions(+), 21 deletions(-) diff --git a/vignettes/articles/delay-distributions.Rmd b/vignettes/articles/delay-distributions.Rmd index 1a7cc96..d3edd1c 100644 --- a/vignettes/articles/delay-distributions.Rmd +++ b/vignettes/articles/delay-distributions.Rmd @@ -33,17 +33,33 @@ knitr::opts_chunk$set( ) ``` -This package accommodates 3 different ways of specifying a delay distribution in the renewal equation, and this vignette illustrates these in turn. The renewal equation specifies that the expected incidence at time $t$ is a weighted sum of previous incidence, multiplied by $R_t$: \begin{equation} +This package accommodates 3 different ways of specifying a delay distribution +in the renewal equation, and this vignette illustrates these in turn. The +renewal equation specifies that the expected incidence at time $t$ is a +weighted sum of previous incidence, multiplied by $R_t$: +\begin{equation} \mathbb{E}[I_t \mid I_1,\ldots,I_{t-1}] = R_t \sum_{a=0}^t w_a I_{t-a}, -\end{equation} calculated by convolving the preceding $a$ days of new infections with the discretized generation (or serial) interval distribution $w$ of length $t$. It is the allowable distributions $\{w\}_{a=0}^n$ that are the focus. +\end{equation} +calculated by convolving the preceding $a$ days of new infections with the +discretized generation (or serial) interval distribution $w$ of length $t$. It +is the allowable distributions $\{w\}_{a=0}^n$ that are the focus. ## Discretization -First, it is important to recognize that the sequence $\{w\}_{a=0}^n$ is usually a discretization of a probability density function, most frequently gamma or Weibull. This is because observed incidence happens at discrete time points like days or weeks, while time is continuous. Using the default (parametric) delay distributions requires calculating a discrete approximation. +First, it is important to recognize that the sequence $\{w\}_{a=0}^n$ is +usually a discretization of a probability density function, most frequently +gamma or Weibull. This is because observed incidence happens at discrete time +points like days or weeks, while time is continuous. Using the default +(parametric) delay distributions requires calculating a discrete approximation. ## Default (parametric) delay distribution -By default, `estimate_rt()` uses a gamma distribution parameterized by the shape $k$ and scale $\theta$. This density has pdf \[ f_W(w) = \frac{1}{\Gamma(k)\theta^k} w\^{k-1} e\^{-w/\theta} I(w \> 0), \] where $I$ is the indicator function. The mean of this distribution is $k\theta$, and the variance is $k\theta^2$. Both $k$ and $\theta$ must be greater than 0. The figure below shows a few examples densities from this family. +By default, `estimate_rt()` uses a gamma distribution parameterized by the +shape $k$ and scale $\theta$. This density has pdf \[ f_W(w) = \frac{1}{\Gamma +(k)\theta^k} w\^{k-1} e\^{-w/\theta} I(w \> 0), \] where $I$ is the indicator +function. The mean of this distribution is $k\theta$, and the variance is +$k\theta^2$. Both $k$ and $\theta$ must be greater than 0. The figure below +shows a few examples densities from this family. ```{r gamma-pdf, echo = FALSE} ggplot(data.frame(x = seq(.001, 20, length.out = 100)), aes(x)) + @@ -65,11 +81,24 @@ ggplot(data.frame(x = seq(.001, 20, length.out = 100)), aes(x)) + scale_y_continuous(expand = expansion(c(0, 0.05))) ``` -The default is the orange curve shown above. Given an incidence sequence of length $n$, internally, `delay_calculator()` creates $\{w\}_{a=0}^n$ by \[ w_a = F(a + 1) - F(a - 1), \] and then renormalizing with `w / sum(w)` to ensure that it sums to 1. Here, $F$ is the cumulative distribution function for Gamma (though similar works for any continuous distribution). Note that this formula assumes that the probability of a 0 delay is *not* 0. This is related to the phenomenon of "double interval censoring": a primary infection on day $t$ may have occurred any time between 0:00 and 23:59 and the same for a secondary infection on day $t+1$. So even if the gap in reported days is only "1 day", those infections may have happened as little as a few minutes apart, or as much as 48 hours apart. - -When using this default, the delay distribution is necessarily assumed to be the same for all $t$, and we compute the convolution of `w` and `I` with the Fast Fourier Transform. - -Just to illustrate this behaviour, we show the results for the default setting on the included `cancovid` data. +The default is the orange curve shown above. Given an incidence sequence of +length $n$, internally, `delay_calculator()` creates $\{w\}_{a=0}^n$ by \[ w_a += F(a + 1) - F(a - 1), \] and then renormalizing with `w / sum(w)` to ensure +that it sums to 1. Here, $F$ is the cumulative distribution function for Gamma +(though similar works for any continuous distribution). Note that this formula +assumes that the probability of a 0 delay is *not* 0. This is related to the +phenomenon of "double interval censoring": a primary infection on day $t$ may +have occurred any time between 0:00 and 23:59 and the same for a secondary +infection on day $t+1$. So even if the gap in reported days is only "1 day", +those infections may have happened as little as a few minutes apart, or as much +as 48 hours apart. + +When using this default, the delay distribution is necessarily assumed to be +the same for all $t$, and we compute the convolution of `w` and `I` with the +Fast Fourier Transform. + +Just to illustrate this behaviour, we show the results for the default setting +on the included `cancovid` data. ```{r can-default} can_default <- estimate_rt( @@ -82,7 +111,12 @@ plot(can_default) + coord_cartesian(ylim = c(0.5, 2)) ## Constant, non-parametric delay distribution -If we don't believe that the Gamma distribution closely approximates the serial interval distribution, then we can specify our own distribution manually. For example, [Backer et al., Table S1](https://doi.org/10.2807/1560-7917.ES.2022.27.6.2200042) gives observed serial intervals for the Omicron (SGTF) and Delta (non-SGTF) COVID-19 variants during 2 weeks in 2021 in the Netherlands. For the sake of illustration, we aggregate these together and use this as our (constant) delay distribution. +If we don't believe that the Gamma distribution closely approximates the serial +interval distribution, then we can specify our own distribution manually. For +example, [Backer et al., Table S1](https://doi.org/10.2807/1560-7917.ES.2022.27.6.2200042) gives observed serial intervals for the Omicron (SGTF) and Delta +(non-SGTF) COVID-19 variants during 2 weeks in 2021 in the Netherlands. For the +sake of illustration, we aggregate these together and use this as our +(constant) delay distribution. ```{r cancov-nonpar} # Data from Backer et al. @@ -103,7 +137,8 @@ ggplot( scale_y_continuous(expand = expansion(c(0, 0.05))) ``` -This distribution looks something like a gamma, but it has finite support. We can easily use it instead. +This distribution looks something like a gamma, but it has finite support. We +can easily use it instead. ```{r can-nonpar} can_nonpar <- estimate_rt( @@ -115,7 +150,9 @@ can_nonpar <- estimate_rt( plot(can_nonpar) + coord_cartesian(ylim = c(0.5, 2)) ``` -The result is much less dramatic than the previous version. This is likely because the distribution is much more concentrated near short delays. We can see this by examining the two CDFs. +The result is much less dramatic than the previous version. This is likely +because the distribution is much more concentrated near short delays. We can +see this by examining the two CDFs. ```{r cdfs, echo = FALSE} cdfs <- data.frame( @@ -132,7 +169,10 @@ ggplot(cdfs |> pivot_longer(-x), aes(x, value, color = name)) + ## Time-varying delays -Finally, we also allow time-varying delay distributions. This is accomplished with a matrix. This requires a bit of work, but is not too challenging. For example, to create the correct matrix using the Baker et al. delays, the necessary code is the following (note that this is *not* yet time varying): +Finally, we also allow time-varying delay distributions. This is accomplished +with a matrix. This requires a bit of work, but is not too challenging. For +example, to create the correct matrix using the Baker et al. delays, the +necessary code is the following (note that this is *not* yet time varying): ```{r backer-matrix, eval=FALSE} # library(Matrix) @@ -152,13 +192,17 @@ delay_mat <- drop0(toeplitz2(c(rep(0, n - 1), delay, rep(0, n - dd)), n, n)) delay_mat <- delay_mat / rowSums(delay_mat) # renormalize ``` -We could then simply pass `delay_mat` to `estimate_rt(..., delay_distn = delay_mat)`. The result would be the same as that shown above. +We could then simply pass `delay_mat` to +`estimate_rt(..., delay_distn = delay_mat)`. The result would be the same as +that shown above. ## Variant-specific delays -For a more illustrative version of the previous example, we revisit the Canadian data, but combine it with some additional information. +For a more illustrative version of the previous example, we revisit the +Canadian data, but combine it with some additional information. -We use variant circulation information from [CoVaRR-Net's Duotang notebook](https://covarr-net.github.io/duotang/duotang.html). For simplicity, the code for this processing is hidden in the fold. +We use variant circulation information from [CoVaRR-Net's Duotang notebook](https://covarr-net.github.io/duotang/duotang.html). For simplicity, the code +for this processing is hidden in the fold.
@@ -253,7 +297,9 @@ saveRDS(counts, "duotang-counts.rds")
-We smooth the raw data using multinomial logistic regression on a third order orthogonal polynomial to produce the following estimated variant proportions in Canada. +We smooth the raw data using multinomial logistic regression on a third order +orthogonal polynomial to produce the following estimated variant proportions in +Canada. ```{r duotang-counts, echo=FALSE, message=FALSE} props <- readRDS("duotang-counts.rds") |> @@ -314,7 +360,9 @@ can_pred_class <- can_props_smoothed |> fill(variant, .direction = "downup") ``` -Using the estimated proportions, we label each date with the dominant variant at the time (and restrict ourselves to the time period of the included case data). +Using the estimated proportions, we label each date with the dominant variant +at the time (and restrict ourselves to the time period of the included case +data). ```{r, echo=FALSE, fig.height=1, dev='png'} ggplot( @@ -329,7 +377,14 @@ ggplot( ## Creating delay distributions -We use the delay distribution data from the meta analysis in [Xu, X., Wu, Y., Kummer, A.G. et al.](https://doi.org/10.1186/s12916-023-03070-8). The raw data for their analysis is available on [GitHub](https://github.com/xxyy0574/COVID-19-transmission-parameters). We use a simple version of their procedure, taking the median calculations across studies to find the mean and standard deviation of the delay for each variant separately. Then we convert these to the $k$ and $\theta$ parameters for the gamma distribution. This analysis is hidden below the fold. +We use the delay distribution data from the meta analysis in +[Xu, X., Wu, Y., Kummer, A.G. et al.](https://doi.org/10.1186/s12916-023-03070-8). +The raw data for their analysis is available on +[GitHub](https://github.com/xxyy0574/COVID-19-transmission-parameters). +We use a simple version of their procedure, taking the median calculations +across studies to find the mean and standard deviation of the delay for each +variant separately. Then we convert these to the $k$ and $\theta$ parameters +for the gamma distribution. This analysis is hidden below the fold.
@@ -381,7 +436,9 @@ saveRDS(delay_dstns_byvar, "delay-distns-byvar.rds")
-We'll use just the Serial Intervals (SI), though the incubation periods (IP), and generation times (GT) are also available. Below, we visualize the estimated delays for the 4 variants that were most prevalent in Canada. +We'll use just the Serial Intervals (SI), though the incubation periods (IP), +and generation times (GT) are also available. Below, we visualize the estimated +delays for the 4 variants that were most prevalent in Canada. ```{r variant-si-plot, echo=FALSE} delay_dstns_can <- readRDS("delay-distns-byvar.rds") |> @@ -401,7 +458,8 @@ delay_dstns_can |> scale_y_continuous(expand = expansion(c(0, 0.05))) ``` -We use these delays weighted by the probabilities of the circulating variants to estimate $R_t$. First, we build the delay matrix. +We use these delays weighted by the probabilities of the circulating variants +to estimate $R_t$. First, we build the delay matrix. ```{r tvar-matrix, message=FALSE} can_props_voc <- can_props_smoothed |> From 19b21a20a9ae008c06b15488ac13e36f24998baa Mon Sep 17 00:00:00 2001 From: "Daniel J. McDonald" Date: Thu, 14 May 2026 12:03:22 -0700 Subject: [PATCH 6/6] bump checkout action version --- .github/workflows/R-CMD-check.yaml | 2 +- .github/workflows/pkgdown.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 1dc2cbd..9156514 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -30,7 +30,7 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index fb6645f..7b1ab01 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -20,7 +20,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - uses: r-lib/actions/setup-pandoc@v2