added a Gaussian noise around a sinusoidal mean

.Rbuildignore

@@ -5,3 +5,4 @@
 ^docs$
 ^pkgdown$
 ^\.github$
+^vignettes/articles$

DESCRIPTION

@@ -13,6 +13,7 @@ Biarch: true
 Depends: 
     R (>= 3.4.0)
 Imports: 
+    extraDistr,
     methods,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),

R/simulation.R

@@ -59,6 +59,7 @@ generate_vector_serial <- function(nt, mean_si, sd_si) {
 #' @param mean_si A numeric value representing the mean of the serial interval distribution.
 #' @param sd_si A numeric value representing the standard deviation of the serial interval distribution.
 #' @param i_0 An integer specifying the initial number of cases at time point 1.
+#' @param X An optional matrix of covariates with `nt` rows.
 #'
 #' @return A list containing the following components:
 #' \describe{
@@ -81,7 +82,7 @@ generate_vector_serial <- function(nt, mean_si, sd_si) {
 #' rt_fun <- function(t) { 1.5 * exp(-0.05 * t) }
 #' simulate_renewal_epidemic(rt_fun, 100, 5, 2, 10)
 #' }
-simulate_renewal_epidemic <- function(Rt_fun, nt, mean_si, sd_si, i_0){
+simulate_renewal_epidemic <- function(Rt_fun, nt, mean_si, sd_si, i_0, X=NULL){
 
   # Input validation
   if (!is.numeric(nt) || nt <= 0 || nt != as.integer(nt)) {
@@ -99,12 +100,18 @@ simulate_renewal_epidemic <- function(Rt_fun, nt, mean_si, sd_si, i_0){
   if (!is.function(Rt_fun)) {
     stop("Parameter 'Rt_fun' should be a function.")
   }
+  if (!is.null(X) && (!is.matrix(X) || nrow(X) != nt)) {
+    stop("Parameter 'X' should be a matrix with 'nt' rows.")
+  }
 
   # Time series and Rt
   t = 1:nt
   Rt <- vector(length = nt)
   for(i in seq_along(Rt)) {
-    Rt[i] = rt_fun(t[i])
+    if(is.null(X))
+      Rt[i] = rt_fun(t[i])
+    else
+      Rt[i] <- rt_fun(t[i], X[i, ])
   }
 
   # Total infectiousness and incidence with initial imports

vignettes/articles/fitting_synthetic_data_including_covariates.Rmd

@@ -0,0 +1,95 @@
+---
+title: "Fitting synthetic data including covariates"
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(epidp)
+```
+
+In this document, we explore how the incorporation of covariate information affects
+estimation of $R_t$.
+
+## Added Gaussian noise around a sinusoidal covariate-driven mean $R_t$
+```{r}
+# define sinusoidal Rt with noise
+rt_fun <- function(t, x) {
+  x[1] * exp(x[2])
+}
+nt <- 200
+t <- 1:nt
+f <- 1.3 + 1.2 * sin(4 * (pi / 180) * t)
+g <- vector(length = nt)
+g[1] <- rnorm(1, 0, 1)
+rho <- 0.8
+for(i in 2:length(g)) {
+  g[i] <- rho * g[i - 1] + rnorm(1, 0, 0.1)
+}
+X <- matrix(c(f, g), nrow = length(x), ncol = 2)
+
+# simulation parameters
+mean_si <- 6.5
+sd_si <- 4.03
+i_0 <- 10
+
+# data frame of outputs
+epidemic_df <- simulate_renewal_epidemic(rt_fun, nt, mean_si, sd_si, i_0, X)
+
+# plot
+epidemic_df %>%
+  mutate(f=f) %>% 
+  select(-c(w_dist, lambda_t)) %>%
+  pivot_longer(c(R_t, f)) %>%
+  ggplot(aes(x = t, y = value, colour = name)) +
+  geom_line() +
+  scale_y_continuous(
+    name = "Incidence (i_t)",
+    sec.axis = sec_axis(~ . / 1000, name = "Reproduction Number (R_t)")
+  ) +
+  labs(x = "Time", colour = "Variable") +
+  theme_minimal()
+```
+We first try to estimate $R_t$ without covariate information.
+We now use a Stan version of EpiFilter to estimate the $R_t$ profile.
+```{r}
+# fit model
+fit <- fit_epifilter(
+  N=length(epidemic_df$i_t),
+  C=epidemic_df$i_t,
+  w=epidemic_df$w_dist,
+  iter=200
+)
+
+# look at MCMC summaries
+print(fit, c("sigma", "R"))
+```
+
+We now overlay the estimated $R_t$ versus the actual values. The estimated $R_t$ values
+coincide reasonably with the true values.
+```{r}
+# extract posterior quantiles
+R_draws <- rstan::extract(fit, "R")[[1]]
+lower <- apply(R_draws, 2, function(x) quantile(x, 0.025))
+middle <- apply(R_draws, 2, function(x) quantile(x, 0.5))
+upper <- apply(R_draws, 2, function(x) quantile(x, 0.975))
+
+# plot
+epidemic_df %>% 
+  mutate(
+    lower=lower,
+    middle=middle,
+    upper=upper
+  ) %>% 
+  select(t, R_t, lower, middle, upper) %>% 
+  ggplot(aes(x=t, y=R_t)) +
+  geom_line(colour="red") +
+  geom_ribbon(aes(ymin=lower, ymax=upper), fill="blue", alpha=0.4) +
+  geom_line(aes(y=middle))
+```
+

vignettes/fitting_synthetic_data_using_epidp.Rmd

@@ -135,7 +135,8 @@ We now use a Stan version of EpiFilter to estimate the $R_t$ profile.
 fit <- fit_epifilter(
   N=length(epidemic_df$i_t),
   C=epidemic_df$i_t,
-  w=epidemic_df$w_dist
+  w=epidemic_df$w_dist,
+  iter=200
 )
 
 # look at MCMC summaries