corrected covariate-forced model

.Rbuildignore

@@ -7,3 +7,4 @@
 ^\.github$
 ^vignettes/articles$
 ^README\.Rmd$
+^CONTRIBUTING\.md$

.github/workflows/R-CMD-check.yaml

@@ -0,0 +1,52 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+name: R-CMD-check.yaml
+
+permissions: read-all
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: macos-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
+
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2
+        with:
+          upload-snapshots: true
+          build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")'

DESCRIPTION

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

NAMESPACE

@@ -3,8 +3,9 @@
 export(fit_epifilter)
 export(fit_epifilter_covariates)
 export(simulate_renewal_epidemic)
-exportPattern("^[[:alpha:]]+")
 import(Rcpp)
 import(methods)
+importFrom(RcppParallel,RcppParallelLibs)
 importFrom(rstan,sampling)
+importFrom(rstantools,rstan_config)
 useDynLib(epidp)

R/epidp-package.R

@@ -0,0 +1,8 @@
+#' @keywords internal
+"_PACKAGE"
+
+## usethis namespace: start
+#' @importFrom rstantools rstan_config
+#' @importFrom RcppParallel  RcppParallelLibs
+## usethis namespace: end
+NULL

R/extra.R

@@ -2,5 +2,4 @@
 #' @import Rcpp
 #' @import methods
 #' @importFrom rstan sampling
-#' @exportPattern "^[[:alpha:]]+"
 NULL

R/fit_models.R

@@ -1,52 +1,55 @@
-
 #' Estimate the time-varying reproduction number from incidence data
 #'
 #' @param N number of data points
 #' @param C case counts
 #' @param w discretised serial interval
 #' @param is_sampling Boolean indicating if sampling is to be used (if TRUE) or
 #' optimisation (if FALSE). Defaults to TRUE
-#' @param ...
+#' @param ... Arguments passed to [rstan::sampling()] (if `is_sampling=TRUE`) or
+#' [rstan::optimizing()] (if `is_sampling=FALSE`)
 #'
 #' @return a stanfit object
 #' @export
-fit_epifilter <- function(N, C, w, is_sampling=TRUE, ...) {
+fit_epifilter <- function(N, C, w, is_sampling = TRUE, ...) {
   wmax <- length(w)
-  standata <- list(N=N,
-                   C=C,
-                   wmax=wmax,
-                   w=w)
-  if(is_sampling)
+  standata <- list(
+    N = N,
+    C = C,
+    wmax = wmax,
+    w = w
+  )
+  if (is_sampling) {
     out <- rstan::sampling(stanmodels$epifilter, data = standata, ...)
-  else
+  } else {
     out <- rstan::optimizing(stanmodels$epifilter, data = standata, ...)
+  }
   return(out)
 }
 
 #' Estimate the time-varying reproduction number from incidence data using covariate
 #' data
 #'
-#' @param N number of data points
-#' @param C case counts
-#' @param w discretised serial interval
+#' @inheritParams fit_epifilter
 #' @param X matrix of covariates (which likely should include a first column of 1s)
 #' of dimensions N x (N_covariates + 1)
-#' @param ...
 #'
 #' @return a stanfit object
 #' @export
-fit_epifilter_covariates <- function(N, C, w, X, is_sampling=TRUE, ...) {
+fit_epifilter_covariates <- function(N, C, w, X, is_sampling = TRUE, ...) {
   wmax <- length(w)
   N_covariates <- dim(X)[2]
-  standata <- list(N=N,
-                   C=C,
-                   w=w,
-                   wmax=wmax,
-                   N_covariates=N_covariates,
-                   X=X)
-  if(is_sampling)
+  standata <- list(
+    N = N,
+    C = C,
+    w = w,
+    wmax = wmax,
+    N_covariates = N_covariates,
+    X = X
+  )
+  if (is_sampling) {
     out <- rstan::sampling(stanmodels$epifilter_covariates, data = standata, ...)
-  else
+  } else {
     out <- rstan::optimizing(stanmodels$epifilter_covariates, data = standata, ...)
+  }
   return(out)
 }

R/simulation.R

@@ -1,5 +1,3 @@
-
-
 #' Generate a Discretized Serial Interval Distribution
 #'
 #' This function generates a discretized serial interval distribution based on
@@ -19,9 +17,8 @@
 #'
 #' @examples
 #' generate_vector_serial(10, 5, 2)
-#' }
+#' @export
 generate_vector_serial <- function(nt, mean_si, sd_si) {
-
   if (!is.numeric(nt) || nt <= 0 || nt != as.integer(nt)) {
     stop("Parameter 'nt' should be a positive integer.")
   }
@@ -77,11 +74,12 @@ generate_vector_serial <- function(nt, mean_si, sd_si) {
 #' the time-varying reproduction number.
 #'
 #' @examples
-#' rt_fun <- function(t) { 1.5 * exp(-0.05 * t) }
+#' rt_fun <- function(t) {
+#'   1.5 * exp(-0.05 * t)
+#' }
 #' simulate_renewal_epidemic(rt_fun, 100, 5, 2, 10)
 #' @export
-simulate_renewal_epidemic <- function(Rt_fun, nt, mean_si, sd_si, i_0, X=NULL){
-
+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)) {
     stop("Parameter 'nt' should be a positive integer.")
@@ -103,27 +101,31 @@ simulate_renewal_epidemic <- function(Rt_fun, nt, mean_si, sd_si, i_0, X=NULL){
   }
 
   # Time series and Rt
-  t = 1:nt
+  t <- 1:nt
   Rt <- vector(length = nt)
-  for(i in seq_along(Rt)) {
-    if(is.null(X))
-      Rt[i] = Rt_fun(t[i])
-    else
+  for (i in seq_along(Rt)) {
+    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
-  Lt = rep(0, nt); It = Lt; It[1] = i_0; Lt[1] = It[1]
+  Lt <- rep(0, nt)
+  It <- Lt
+  It[1] <- i_0
+  Lt[1] <- It[1]
 
   w_dist <- generate_vector_serial(nt, mean_si, sd_si)
 
   # Simulate from standard renewal model
-  for(i in 2:nt){
+  for (i in 2:nt) {
     # Total infectiousness is a convolution
-    Lt[i] = sum(It[seq(i-1, 1, -1)] * w_dist[1:(i-1)])
+    Lt[i] <- sum(It[seq(i - 1, 1, -1)] * w_dist[1:(i - 1)])
     # Poisson renewal model
-    It[i] = stats::rpois(1, Lt[i] * Rt[i])
+    It[i] <- stats::rpois(1, Lt[i] * Rt[i])
   }
 
-  data.frame(t=t, i_t=It, R_t=Rt, lambda_t=Lt, w_dist=w_dist)
+  data.frame(t = t, i_t = It, R_t = Rt, lambda_t = Lt, w_dist = w_dist)
 }

README.Rmd

@@ -17,6 +17,7 @@ knitr::opts_chunk$set(
 
 <!-- badges: start -->
 [![codecov](https://codecov.io/github/ben18785/epidp/graph/badge.svg?token=STG0INT235)](https://codecov.io/github/ben18785/epidp)
+[![R-CMD-check](https://github.com/ben18785/epidp/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/ben18785/epidp/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
 The goal of `epidp` is to allow covariate information to inform estimates of the time-varying reproduction number, $R_t$.
@@ -41,13 +42,14 @@ library(magrittr)
 library(purrr)
 library(tidyr)
 
-rt_fun = function(t){
-  if(t <= 60)
-    R = 2
-  else if (t <= 90)
-    R = 0.5
-  else
-    R = 1
+rt_fun <- function(t) {
+  if (t <= 60) {
+    R <- 2
+  } else if (t <= 90) {
+    R <- 0.5
+  } else {
+    R <- 1
+  }
   R
 }
 
@@ -64,7 +66,7 @@ epidemic_df <- simulate_renewal_epidemic(rt_fun, nt, mean_si, sd_si, i_0)
 transform_factor <- 300
 epidemic_df %>%
   select(-c(w_dist, lambda_t)) %>%
-  mutate(R_t = R_t * transform_factor) %>% 
+  mutate(R_t = R_t * transform_factor) %>%
   pivot_longer(c(i_t, R_t)) %>%
   ggplot(aes(x = t, y = value, colour = name)) +
   geom_line() +
@@ -85,22 +87,22 @@ period when case counts are low.
 ```{r}
 # fit model
 fit <- fit_epifilter(
-  N=length(epidemic_df$i_t),
-  C=epidemic_df$i_t,
-  w=epidemic_df$w_dist,
-  is_sampling=FALSE,
-  as_vector=FALSE
+  N = length(epidemic_df$i_t),
+  C = epidemic_df$i_t,
+  w = epidemic_df$w_dist,
+  is_sampling = FALSE,
+  as_vector = FALSE
 )
 
 # plot
 R <- fit$par$R
-epidemic_df %>% 
-  mutate(estimated=R) %>% 
-  rename(true=R_t) %>% 
-  select(t, estimated, true) %>% 
-  pivot_longer(c(estimated, true)) %>% 
-  ggplot(aes(x=t, y=value)) +
-  geom_line(aes(colour=name)) +
+epidemic_df %>%
+  mutate(estimated = R) %>%
+  rename(true = R_t) %>%
+  select(t, estimated, true) %>%
+  pivot_longer(c(estimated, true)) %>%
+  ggplot(aes(x = t, y = value)) +
+  geom_line(aes(colour = name)) +
   scale_color_brewer("R_t", palette = "Dark2") +
   ylab("R_t") +
   xlab("Time") +

README.md

@@ -6,6 +6,7 @@
 <!-- badges: start -->
 
 [![codecov](https://codecov.io/github/ben18785/epidp/graph/badge.svg?token=STG0INT235)](https://codecov.io/github/ben18785/epidp)
+[![R-CMD-check](https://github.com/ben18785/epidp/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/ben18785/epidp/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
 The goal of `epidp` is to allow covariate information to inform
@@ -53,13 +54,14 @@ library(tidyr)
 #> 
 #>     extract
 
-rt_fun = function(t){
-  if(t <= 60)
-    R = 2
-  else if (t <= 90)
-    R = 0.5
-  else
-    R = 1
+rt_fun <- function(t) {
+  if (t <= 60) {
+    R <- 2
+  } else if (t <= 90) {
+    R <- 0.5
+  } else {
+    R <- 1
+  }
   R
 }
 
@@ -76,7 +78,7 @@ epidemic_df <- simulate_renewal_epidemic(rt_fun, nt, mean_si, sd_si, i_0)
 transform_factor <- 300
 epidemic_df %>%
   select(-c(w_dist, lambda_t)) %>%
-  mutate(R_t = R_t * transform_factor) %>% 
+  mutate(R_t = R_t * transform_factor) %>%
   pivot_longer(c(i_t, R_t)) %>%
   ggplot(aes(x = t, y = value, colour = name)) +
   geom_line() +
@@ -101,22 +103,22 @@ when case counts are low.
 ``` r
 # fit model
 fit <- fit_epifilter(
-  N=length(epidemic_df$i_t),
-  C=epidemic_df$i_t,
-  w=epidemic_df$w_dist,
-  is_sampling=FALSE,
-  as_vector=FALSE
+  N = length(epidemic_df$i_t),
+  C = epidemic_df$i_t,
+  w = epidemic_df$w_dist,
+  is_sampling = FALSE,
+  as_vector = FALSE
 )
 
 # plot
 R <- fit$par$R
-epidemic_df %>% 
-  mutate(estimated=R) %>% 
-  rename(true=R_t) %>% 
-  select(t, estimated, true) %>% 
-  pivot_longer(c(estimated, true)) %>% 
-  ggplot(aes(x=t, y=value)) +
-  geom_line(aes(colour=name)) +
+epidemic_df %>%
+  mutate(estimated = R) %>%
+  rename(true = R_t) %>%
+  select(t, estimated, true) %>%
+  pivot_longer(c(estimated, true)) %>%
+  ggplot(aes(x = t, y = value)) +
+  geom_line(aes(colour = name)) +
   scale_color_brewer("R_t", palette = "Dark2") +
   ylab("R_t") +
   xlab("Time") +