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
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
diff --git a/vignettes/articles/delay-distributions.Rmd b/vignettes/articles/delay-distributions.Rmd
index 351a04d..d3edd1c 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,
@@ -35,78 +34,89 @@ 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
+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$:
+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.
+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
+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 +129,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 +137,22 @@ 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 +162,78 @@ 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
+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)
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
```
-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 +249,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 +294,12 @@ 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 +313,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 +331,44 @@ 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
+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()) +
@@ -320,18 +379,20 @@ ggplot(can_pred_class, aes(x = Date, fill = fct_relevel(var, "Ancestral lineage"
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.
+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 +411,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 +433,16 @@ 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
+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,25 +453,43 @@ 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 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}
-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_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)
+ )
+
+# 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()
+
+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
```
@@ -416,13 +497,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))
```
-
-
-
-