prophylactic.Rmd

---
title: "Does prophylactic AMU increase the susceptibility to infections?"
output:
  html_document:
    theme: cerulean
    toc: yes
    toc_float: yes
  pdf_document:
    toc: yes
editor_options: 
  chunk_output_type: console
---

```{css horizontal scrolling, echo = FALSE}
pre, code {white-space:pre !important; overflow-x:scroll !important}

html { overflow-x: scroll; }
```

```{r general options, include = FALSE}
knitr::knit_hooks$set(
  margin = function(before, options, envir) {
    if (before) par(mgp = c(1.5, .5, 0), bty = "n", plt = c(.105, .97, .13, .97))
    else NULL
  },
  prompt = function(before, options, envir) {
    options(prompt = if (options$engine %in% c("sh", "bash")) "$ " else "> ")
  })

knitr::opts_chunk$set(cache = FALSE, autodep = TRUE, message = FALSE, warning = FALSE,
                      dev.args = list(pointsize = 11), fig.height = 3.5,
                      fig.width = 4.24725, fig.retina = 2, fig.align = "center")

options(width = 137)
```

> If you want to run this code locally, look at the last two sections

## Introduction

Here we aim at testing whether the prophylactic use of antibiotics (exposure)
has an effect on the occurence of clinical signs (outcome). The main challenge
in this analysis is the preparation of the data in such a way that AMU exposure
can be considered prophylactic AMU exposure. The overall strategy for data
preparation is described in **Figure 1**.

<div style="text-align:center"> <img src="figure1.png" width="600"/> </div>

**Figure 1: Preparation of data before analysis.** The horizontal arrow
represents the time line of a flock, divided into weeks that are represented by
rectangles, starting on week 1 (on the left). For any given week *w*
(represented here by the hashed rectangle), we compute (i) an exposure variable
based on the use or not of antibiotics on week *w* (step 1, in green) and (ii)
an outcome variable based on the occurence or not of clinical signs over an
*observation period* of *x* weeks after week *w* (step 2, in red). Statistical
analyses will consist in investigating whether AMU on week *w* (exposure)
affects the occurence of clinical signs over the observation period (outcome).
In order to make sure that AMU exposure on week *w* does correspond to AMU that
we can consider as *prophylactic*, we need to filter out from the analysis all
the weeks *w* preceded by (i) the presence of clinical signs over the *y* weeks
before week *w* (including week *w*) and (ii) AMU over the *z* weeks before week
*w* (naturally excluding week *w* since this information is used to compute the
exposure variable). This step 0 is shown in blue on the figure. Finally, the
analysis may include potential confounding factors (shown in orange) such as the
age of the chicken (i.e. week *w*) as well as AMU during the first *a* weeks of
the flock's life (in grey) and during the *x* weeks of the observation period.

Note that the duration *y* over which the filtering based on clinical signs is
done and the duration *z* over which the filtering based on AMU is done do not
have to be identical, although it probably makes sense if they are. On
**Figure 1**, CS corresponds to a given set of clinical signs and AMU
corresponds to a given set of antimicrobials. Note that the sets of clinical
signs used for the filtering process and the computation of the outcome variable
do not have to be identical, although, again, it probably makes sense if they
are. Same remark for the sets of antimicrobials used for the filtering process,
and the computations of the exposure variable and the confounding variables over
the first weeks of the flock and the observation period. In the code below, we
create the following 4 functions allowing us to prepare the data according to
the strategy outlined in **Figure 1**:

* `make_exposure(df, antimicrobials)` computes a boolean exposure variable using
  the set of `antimicrobials` antimicrobials of the data frame `df`.

* `make_outcome(df, outcome, loc)` computes a boolean outcome variable for the
  `loc` weeks that follow week *w*, using the set `outcome` of clinical signs of
  the data frame `df`. Note that this function will also be used to compute the
  AMU covariable during the observation period, in which case `outcome` will
  then contain the set of antimicrobials to consider.

* `prophylactic(df, clinical_signs, antimicrobials, lcs, lam)` computes a
  boolean variable expressing whether any given week of the data `df` can be
  included in the analysis or not. This is based on the filtering process
  explained in **figure 1** and considering the `clinical_signs` set of clinical
  signs over the `lcs` weeks before week *w* (including week *w*) and the
  `antimicrobials` set of antimicrobials over `lam` weeks before week *w*
  (excluding week *w*).

* `first_weeks(df, antimicrobials, wk)`: computes the boolean AMU covariable
  during the first `wk` weeks of the flock' life, using the `antimicrobials` set
  of antimicrobials.

In these 4 functions, `df` is a data frame of the structure of viparc data, as
described [here](https://viparc.github.io/clires_data). These 4 functions are
designed so that they assemble perfectly well into pipelines. Once the data are
prepared, the analysis can be performed by different statistical tests / models:

* the simplest one is certainely a [**Fisher exact test**](#FET). However, it
appears that both prophylactic AMU and the occurence of clinical signs strongly
depends on the age of the flock (as well as, to a lesser extend, on the AMU
during the first weeks of the flock and the observation period). For that reason,
we considered more sophisticatted analyses allowing to correct for confounding
effects. These are

* [**logistic regression**](#GLM). Here, we compared correction by partial
likelihood ratio tests and **stratification**. It produces similar results.
Because the effect of age on the occurence of clinical signs may be non-linear
we also considered

* [**general additive models**](#GAM) in order to model the age - outcome
relationship in an as flexible way as possible. After the preliminary analyses
shown below, it is this option that we selected for full analysis.

We then developped a number of [pipelines](#pipelines) combining both data
preparation and data analysis. `pipeline1()` is based on the Fisher exact test,
`pipeline3()` is based on the GLM and `pipeline5()` is based on the GAM.
`pipeline2()`, `pipeline4()` and `pipeline6()` are the same as `pipeline1()`,
`pipeline3()` and `pipeline5()` respectively, with simplified interfaces. At the
[end](#end), we created 3 functions that allows to explore the exposure -
outcome relationships using various choices for the definitions of the duration
of the first week, the filtering periods and the observation period, as well as
the definition of the sets of antimicrobials and clinical signs to consider:

* `run_plan(pipeline, df, ...)` runs the pipeline `pipeline` on the data frame
`df` using the sets of parameters as defined by `...`.

* `develop()` and `simplify()` allow to manipulate the output of `run_plan()`.
  
## Packages

Installing the required packages:

```{r starts here}
installed_packages <- row.names(installed.packages())
required <- c("dplyr", "magrittr", "mgcv", "purrr", "readr")
to_install <- setdiff(required, installed_packages)
if (length(to_install)) install.packages(to_install)
if (! require("mgcv.helper")) {
  if (! require("remotes")) install.packages("remotes")
  remotes::install_github("samclifford/mgcv.helper")
}
```

Loading `magrittr` for interactive use:

```{r}
library(magrittr)
```

## Loading the data

```{r eval = FALSE}
viparc <- readr::read_csv("https://raw.githubusercontent.com/viparc/prophylactic/master/data/viparc_qualitative.csv",
                          col_types = paste(c("ciiiddd", rep("l", 51)), collapse = "")) %>% 
  dplyr::select(-nb_chicken, -chicken_disease_death, -chicken_sudden_death, -nb_chicken_sold)
```

```{r eval = TRUE, include = FALSE}
viparc <- readr::read_csv("data/viparc_qualitative.csv",
                          col_types = paste(c("ciiiddd", rep("l", 51)), collapse = "")) %>% 
  dplyr::select(-nb_chicken, -chicken_disease_death, -chicken_sudden_death, -nb_chicken_sold)
```

We have 3 types of variables in this data frame:

* week ID with the variables `farm`, `flock` and `week`;
* presence of clinical signs with the variables `respiratory`, `diarrhoea`,
  `cns`, `malaise`, `leg_lesions` and `sudden_death`;
* AMU (presence / absence too) with the `r ncol(viparc) - 9` variables
  `amoxicillin_use` to `unknown_use`.

For the purpose of the examples used below, we will used the following sets of
antimicrobials and clinical signs:

```{r}
am_set <- c("colistin_use", "oxytetracycline_use")
cs_set <- c("malaise", "diarrhoea")
```

## Filtering out weeks

The week we are interested in are the weeks for which

* there is no clinical signs during the weeks $w, \dots, w - x$;
* there is no AMU during the weeks $w - 1, \dots, w - y$.

Here, we build the function `prophylactic()` that helps performing this
filtering. This function makes use of the function `lagging()` defined below and
that lags the variables `var` of the data frame `df`. Note that this later must
have at least 3 columns named respectively `farm`, `flock` and `week`.

```{r}
lagging <- function(df, var, lag) {
  require(magrittr)                
  stopifnot(c("farm", "flock", "week", var) %in% names(df))
  df %>%
    dplyr::select(!! c("farm", "flock", "week", var)) %>% 
    dplyr::mutate(week = week + lag) %>% 
    dplyr::rename_at(var, paste0, lag)
}
```

Let's try it:

```{r}
lapply(0:2, lagging, df = viparc, var = am_set) %>% 
  purrr::reduce(dplyr::full_join, c("farm", "flock", "week")) %>% 
  dplyr::left_join(viparc, ., c("farm", "flock", "week")) %>% 
  dplyr::select(farm, flock, week, dplyr::matches("colistin|oxytetracycline"))
```

It works as expected (note that we could have used the `dplyr::lag()` function
too). Now we can define the `prophylactic()` function that uses the
above-defined `lagging()` function:

```{r}
prophylactic <- function(df, clinical_signs, antimicrobials, lcs, lam) {
  require(magrittr)
  week_id <- c("farm", "flock", "week")
  purrr::map2(list(c(0, seq_len(lcs)), seq_len(lam)), # doing the lagging on clinical signs and AMU
              list(clinical_signs, antimicrobials),
              ~ purrr::map(.x, lagging, df = df, var = .y)) %>%
    unlist(FALSE) %>% 
    purrr::reduce(dplyr::left_join, week_id) %>% # putting together into a single data frames
# note: we don't need to join with the original data because the original data
# is already here thanks to the lag of 0 in the above lagging() call to the clinical signs.    
    dplyr::select(-week_id) %>% 
    dplyr::mutate_all(dplyr::coalesce, FALSE) %>% # we want to get rid off FALSEs and also NAs
    rowSums() %>%
    not()
}
```

Let's try it:

```{r}
viparc %>% 
  dplyr::mutate(prphlctc = prophylactic(., cs_set, am_set, 2, 2)) %>% 
  dplyr::select(farm, flock, week, am_set, cs_set, prphlctc) %>%
  head(20)
```

The filtering can be performed in one step:

```{r}
viparc %>% 
  dplyr::filter(prophylactic(., cs_set, am_set, 2, 2))
```

Seems to work as expected.

## Computing clinical signs outcomes

The `make_outcome()` function is another use of the `lagging()` function:

```{r}
make_outcome <- function(df, outcome, loc) {
  require(magrittr)
  week_id <- c("farm", "flock", "week")
  lapply(-seq_len(loc), lagging, df = df, var = outcome) %>% # this is a list of loc slots
    purrr::reduce(dplyr::left_join, week_id) %>% # puts together into a single data frame
    dplyr::left_join(df, ., week_id) %>% # merges with original data frame
    dplyr::select(dplyr::matches("-\\d+$")) %>% 
    rowSums() %>% 
    as.logical()
}
```

where `outcome` is a vector of clinical signs names used to compute the outcome,
`loc` is the duration, in weeks, of the observation period. Let's try it:

```{r}
viparc[356:386, ] %>% 
  dplyr::mutate(disease = make_outcome(., cs_set, 2)) %>% 
  dplyr::select(farm, flock, week, cs_set, disease)
```

Seems to work as expected. Note the missing values at the end of the flocks.

## Computing AMU exposure

The `make_exposure()` function is very simple:

```{r}
make_exposure <- function(df, antimicrobials) {
  require(magrittr)
  df %>%
    dplyr::select(antimicrobials) %>%
    rowSums() %>%
    as.logical()
}
```

Where `antimicrobials` is a vector of antimicrobials names that should be
considered to compute the exposure. Let's try it:

```{r}
viparc %>% 
  dplyr::mutate(amu = make_exposure(., am_set)) %>% 
  dplyr::select(farm, flock, week, am_set, amu)
```

Seems to work as expected.

## Optional: covariables

The final statistical analysis can be performed correcting for potential
confounders. Potential confounders that we may want to consider are the AMU
during the first few weeks of the cycle, AMU during the observation period, and
the age of the flock. For the latter one, the variable is simply the week number.
The two other variables are respectively computed by the `first_weeks()` and
`make_outcome()` functions as detailed below.

### AMU during the first week(s)

The following function computes a boolean variable that specifies whether any of
the antimicrobials defined in the `antimicrobials` vector has been used during
the first `wk` weeks of the flock's life:

```{r}
first_weeks <- function(df, antimicrobials, wk = 1) {
  df %>%
    dplyr::filter(week < wk + 1) %>%
    dplyr::group_by(farm, flock) %>% 
    dplyr::group_modify(~ tibble::tibble(out = any(unlist(.x[, antimicrobials])))) %>% 
    dplyr::ungroup() %>% 
    dplyr::left_join(df, ., c("farm", "flock")) %>% 
    dplyr::select(out) %>% 
    unlist()
}
```

Let's try it:

```{r}
viparc %>% 
  dplyr::mutate(amu_first_weeks = first_weeks(., am_set)) %>%
  dplyr::select(farm, flock, week, am_set, amu_first_weeks)
```

And with the first 2 weeks:

```{r}
viparc %>% 
  dplyr::mutate(amu_first_weeks = first_weeks(., am_set, 2)) %>% 
  dplyr::select(farm, flock, week, am_set, amu_first_weeks)
```

Looking at another cycle from another farm:

```{r}
viparc %>% 
  dplyr::mutate(amu_first_weeks = first_weeks(., am_set)) %>% 
  dplyr::select(farm, flock, week, am_set, amu_first_weeks) %>% 
  dplyr::filter(farm == "75-013")
```

Seems to work as expected.

### AMU over the watching period

If one wants to compute a covariable accounting for AMU during the watching
period, this can simply be done by reusing the `make_outcome()` function defined
above:

```{r}
viparc %>%
  dplyr::mutate(amu_observation = make_outcome(., am_set, 2)) %>% 
  dplyr::select(farm, flock, week, am_set, amu_observation) %>%
  head(20)
```

## Examples of uses

Here we show how to assemble the functions `prophylactic()`, `make_exposure()`,
`make_outcome()` and `first_weeks()` into a pipeline to test a hypothesis.

**WARNING 1:** in these pipelines, the filtering of the weeks should always be
the last step of the data preparation, just before the statistical analysis.

**WARNING 2:** `na.exclude()` calls should always be done after the selection of
variables.

For the sake of the examples, we will consider the following data:

```{r}
(data1 <- viparc %>% 
  dplyr::mutate(amu = make_exposure(., am_set),
                disease = make_outcome(., cs_set, 2)) %>% 
  dplyr::filter(prophylactic(., cs_set, am_set, 2, 2)) %>%
  dplyr::select(amu, disease) %>% 
  na.exclude())
```

### Fisher exact test {#FET}

Simple use of a Fisher exact test, without any co-variable:

```{r}
data1 %$%
  table(amu, disease) %T>% print() %>% 
  fisher.test()
```

### Logistic regression {#GLM}

A logistic regression can equivalently be used:

```{r}
data1 %>% 
  glm(disease ~ amu, binomial, .) %>% 
  summary()
```

where the estimate of the coefficient (`amuTRUE`) is the log odds ratio.

### Correcting for potential confounders {#GAM}

In this section we address the problem of confounders and explore on an
example different options to control for them.

#### The presence of potential counfounders

Before we start this section, let's define another utilitary function that adds
a vector of suffixes `suff` to a character vector `x` where we want all possible
combinations of `x` and `suff`:

```{r}
add_suffix <- function(x, suff) {
  require(magrittr)
  expand.grid(x, suff) %>%
    purrr::reduce(purrr::map2, paste0) %>%
    unlist()
}
```

Let's try it:

```{r}
add_suffix(letters[1:5], 1:3)
```

Let's now consider the following example:

```{r}
(data2 <- viparc %>% 
  dplyr::mutate(amu = make_exposure(., am_set),
                disease = make_outcome(., cs_set, 2),
                amu_first_weeks = first_weeks(., am_set),
                amu_observation = make_outcome(., am_set, 2)) %>% 
  dplyr::filter(prophylactic(., cs_set, am_set, 2, 2)) %>% 
  dplyr::select(amu, disease, week, amu_first_weeks, amu_observation))
```

As we can see below, both exposure (`amu`) and outcome (`disease`) depend on
`week`, `amu_first_weeks` and `amu_observation`:

```{r}
purrr::map(add_suffix(c("amu", "disease"), paste(" ~", c("week", "amu_first_weeks", "amu_observation"))),
           ~ .x %>% 
             formula() %>% 
             glm(binomial, na.exclude(data2)) %>% 
             summary() %>% 
             coef() %>% 
             `[`(2, 4) %>% 
             data.frame(model = .x, p = ., stringsAsFactors = FALSE)) %>% 
  dplyr::bind_rows()
```

These three variables `week`, `amu_first_weeks` and `amu_observation` are thus
potential confounders that need to be corrected for.

#### Correcting for potential confounders

If one wants to consider (and correct for) covariables -- such as age, AMU on
the first week or AMU during the observation period --, the pipeline would look
as below:

```{r}
data2 %>% 
  na.exclude() %T>%
  {dplyr::group_by(., amu, disease, amu_first_weeks, amu_observation) %>%
     dplyr::tally() %>%
     dplyr::ungroup() %>%
     print()} %>%  
  glm(disease ~ amu_first_weeks + amu_observation + week + amu, binomial, .) %T>%
  {summary(.) %>% print()} %>% 
  anova(test = "LRT")
```

**NOTE:** If one wanted to correct for potential farm and / or flock effect,
he'd have to consider mixed effect models (not done here). (Not even sure we
need to do this here.)

Note that accounting for potential confounders makes the effect of the AMU
exposure not significant anymore. Furthermore the estimated odds ratio is now
below 1. Let's try to drop the `week` variable:

```{r}
data2 %>% 
  dplyr::select(-week) %>% 
  na.exclude() %>%  
  glm(disease ~ amu_first_weeks + amu_observation + amu, binomial, .) %T>% {summary(.) %>% print()} %>% 
  anova(test = "LRT")
```

which suggests that the `week` variable (i.e. the age of the flock) is really
the variable that causes the spurious relationship between prophylactic use
(`amu`) and occurence of infection (`disease`). Let's try again with `week` but
this time also including interaction with `amu`:

```{r}
data2 %>% 
  na.exclude() %>%  
  glm(disease ~ amu_first_weeks + amu_observation + week * amu, binomial, .) %>% 
  anova(test = "LRT")
```

It doesn't really help. Let's try modeling the age effect in a more flexible way
than a linear relationship, using GAM:

```{r}
summary(mgcv::gam(disease ~ amu_first_weeks + amu_observation + s(week) + amu,
                  binomial, na.exclude(data2)))
```

Not better. Let's try to stratify. For example, under the age of 6 weeks there
seems to have not much effect of `week`, `amu_first_weeks` and
`amu_observation`. Look at the significativity of the confounders for the
1-5-week age strata:

```{r}
data2 %>%
  dplyr::filter(week < 6) %>% 
  na.exclude() %>%  
  glm(disease ~ amu_first_weeks + amu_observation + week + amu, binomial, .) %>% 
  anova(test = "LRT")
```

and for the 1-6-week age strata:

```{r}
data2 %>% 
  dplyr::filter(week < 7) %>% 
  na.exclude() %>%  
  glm(disease ~ amu_first_weeks + amu_observation + week + amu, binomial, .) %>% 
  anova(test = "LRT")
```

Let's thus consider the 6-weeks threshold and investigate the relationship
between `disease` and `amu`:

```{r}
data2 %>% 
  dplyr::filter(week < 6) %>% 
  dplyr::select(amu, disease) %>% 
  na.exclude() %>%  
  glm(disease ~ amu, binomial, .) %>% 
  anova(test = "LRT")
```

No significant relationship between prophylaxy and disease occurence in the
1-to-6-week age band.

## Pipelines {#pipelines}

In this section we show how to build pipelines that combine data preparation and
data analysis.

### Using the Fisher exact test

Let's redefine a Fisher exact test function that returns the results of the test
together with the input contingency table:

```{r}
fisher_test <- function(x, ...) {
  list(x, fisher.test(x, ...))
}
```

Let's define a first pipeline using this new function:

```{r}
pipeline1 <- function(df, clinical_signs_filtering = NULL, antimicrobials_filtering = NULL, lcs,
                      lam, antimicrobials_input = NULL, clinical_signs_output = NULL, watching, ...) {
  require(magrittr)
  var_names <- setdiff(names(df), c("farm", "flock", "week"))
  cs_names <- grep("_use$", var_names, value = TRUE, invert = TRUE)
  am_names <- grep("_use$", var_names, value = TRUE)
  if (is.null(clinical_signs_filtering)) clinical_signs_filtering <- cs_names
  if (is.null(clinical_signs_output)) clinical_signs_output <- cs_names
  if (is.null(antimicrobials_filtering)) antimicrobials_filtering <- am_names
  if (is.null(antimicrobials_input)) antimicrobials_input <- am_names
  df %>% 
    dplyr::mutate(amu = make_exposure(., antimicrobials_input), # creates amu variable
                  disease = make_outcome(., clinical_signs_output, watching)) %>% # creates disease variable
    dplyr::filter(prophylactic(., clinical_signs_filtering, antimicrobials_filtering, lcs, lam)) %>% 
    dplyr::select(amu, disease) %>% 
    na.exclude() %$%
    table(amu, disease) %>% 
    fisher_test(., ...)
}
```

Let's try it:

```{r}
pipeline1(viparc,
# 1. FILTERING -----------------------------------------------------------------
          clinical_signs_filtering = cs_set,
          antimicrobials_filtering = am_set,
          lcs = 2, # the duration over which we perform the filtering on clinical signs
          lam = 2, # the duration over which we perform the filtering on antimicrobials use
# 2. INPUT ---------------------------------------------------------------------
          antimicrobials_input = am_set,
# 3. OUTPUT --------------------------------------------------------------------
          clinical_signs_output = cs_set,
          watching = 2) # the duration of the observation period
```

A simplified version of `pipeline1()` could be:

```{r}
pipeline2 <- function(df, clinical_signs = NULL, antimicrobials = NULL, filtering, watching, ...)
  pipeline1(df, clinical_signs, antimicrobials, filtering, filtering, antimicrobials, clinical_signs, watching, ...)
```

And a call of this function would look like:

```{r}
a <- pipeline2(viparc,
# 1. DEFINING THE CLINICAL SIGNS WE CONSIDER -----------------------------------
          clinical_signs = cs_set,
# 2. DEFINING THE ANTIMICROBIALS USAGES WE CONSIDER ----------------------------
          antimicrobials = am_set,
# 3. THE DURATION OF THE PERIOD BEFORE THE FOCAL WEEK OVER WHICH WE DO THE FILTERING
          filtering = 2,
# 4. THE DURATION OF THE PERIOD AFTER THE FOCAL WEEK OVER WHICH WE OBSERVE -----
          watching = 2)
```

Let's explore the structure of the produced object:

```{r}
str(a)
```

Thanks to `pipeline2()`, it's very easy to explore different options for the
analysis. For example, let's compare this

```{r}
pipeline2(viparc,
# 1. DEFINING THE CLINICAL SIGNS WE CONSIDER -----------------------------------
          clinical_signs = cs_set,
# 2. DEFINING THE ANTIMICROBIALS USAGES WE CONSIDER ----------------------------
          antimicrobials = am_set,
# 3. THE DURATION OF THE PERIOD BEFORE THE FOCAL WEEK OVER WHICH WE DO THE FILTERING
          filtering = 2,
# 4. THE DURATION OF THE PERIOD AFTER THE FOCAL WEEK OVER WHICH WE OBSERVE -----
          watching = 2)
```

with this:

```{r}
pipeline2(viparc,
# 1. DEFINING THE CLINICAL SIGNS WE CONSIDER -----------------------------------
          clinical_signs = c(cs_set, "respiratory"),
# 2. DEFINING THE ANTIMICROBIALS USAGES WE CONSIDER ----------------------------
          antimicrobials = am_set,
# 3. THE DURATION OF THE PERIOD BEFORE THE FOCAL WEEK OVER WHICH WE DO THE FILTERING
          filtering = 2,
# 4. THE DURATION OF THE PERIOD AFTER THE FOCAL WEEK OVER WHICH WE OBSERVE -----
          watching = 2)
```

If you want all the clinical signs and all the antimicrobials:

```{r}
pipeline2(viparc, filtering = 2, watching = 2)
```

### Logistic regression with covariables

We can use similar rationale in order to define pipeline functions that perform
logistic regressions with covariables:

```{r}
pipeline3 <- function(df, clinical_signs_filtering = NULL, antimicrobials_filtering = NULL, lcs,
                      lam, antimicrobials_input = NULL, clinical_signs_output = NULL, watching,
                      antimicrobials_week1 = NULL, wk, antimicrobials_watching = NULL) {
  require(magrittr)
  var_names <- setdiff(names(df), c("farm", "flock", "week"))
  cs_names <- grep("_use$", var_names, value = TRUE, invert = TRUE)
  am_names <- grep("_use$", var_names, value = TRUE)
  if (is.null(clinical_signs_filtering)) clinical_signs_filtering <- cs_names
  if (is.null(clinical_signs_output)) clinical_signs_output <- cs_names
  if (is.null(antimicrobials_filtering)) antimicrobials_filtering <- am_names
  if (is.null(antimicrobials_input)) antimicrobials_input <- am_names
  if (is.null(antimicrobials_week1)) antimicrobials_week1 <- am_names
  if (is.null(antimicrobials_watching)) antimicrobials_watching <- am_names
  df %>% 
    dplyr::mutate(amu = make_exposure(., antimicrobials_input),
                  disease = make_outcome(., clinical_signs_output, watching),
                  amu_first_weeks = first_weeks(., antimicrobials_week1, wk),
                  amu_observation = make_outcome(., antimicrobials_watching, watching)) %>% 
    dplyr::filter(prophylactic(., clinical_signs_filtering, antimicrobials_filtering, lcs, lam)) %>% 
    dplyr::select(farm, flock, week, amu_first_weeks, amu_observation, amu, disease) %>% 
    na.exclude() %>% 
    glm(disease ~ amu_first_weeks + amu_observation + week + amu, binomial, .)
}
```

A wrapper with reduced flexibility in the interface:

```{r}
pipeline4 <- function(df, clinical_signs = NULL, antimicrobials = NULL, filtering, watching, first_wks)
  pipeline3(df, clinical_signs, antimicrobials, filtering, filtering, antimicrobials, clinical_signs, watching, antimicrobials, first_wks, antimicrobials)
```

Let's try it:

```{r}
pipeline4(viparc, filtering = 2, watching = 2, first_wks = 1) %T>% 
  {summary(.) %>% print()} %>% 
  anova(test = "LRT")
```

### Using GAM with covariables

Same thing as above but using a GAM instead of a GLM.:

```{r}
pipeline5 <- function(df, clinical_signs_filtering = NULL, antimicrobials_filtering = NULL, lcs,
                      lam, antimicrobials_input = NULL, clinical_signs_output = NULL, watching,
                      antimicrobials_week1 = NULL, wk, antimicrobials_watching = NULL) {
  require(magrittr)
  var_names <- setdiff(names(df), c("farm", "flock", "week"))
  cs_names <- grep("_use$", var_names, value = TRUE, invert = TRUE)
  am_names <- grep("_use$", var_names, value = TRUE)
  if (is.null(clinical_signs_filtering)) clinical_signs_filtering <- cs_names
  if (is.null(clinical_signs_output)) clinical_signs_output <- cs_names
  if (is.null(antimicrobials_filtering)) antimicrobials_filtering <- am_names
  if (is.null(antimicrobials_input)) antimicrobials_input <- am_names
  if (is.null(antimicrobials_week1)) antimicrobials_week1 <- am_names
  if (is.null(antimicrobials_watching)) antimicrobials_watching <- am_names
  df %>% 
    dplyr::mutate(amu = make_exposure(., antimicrobials_input),
                  disease = make_outcome(., clinical_signs_output, watching),
                  amu_first_weeks = first_weeks(., antimicrobials_week1, wk),
                  amu_observation = make_outcome(., antimicrobials_watching, watching)) %>% 
    dplyr::filter(prophylactic(., clinical_signs_filtering, antimicrobials_filtering, lcs, lam)) %>% 
    dplyr::select(farm, flock, week, amu_first_weeks, amu_observation, amu, disease) %>% 
    na.exclude() %>% 
    mgcv::gam(disease ~ amu_first_weeks + amu_observation + s(week) + amu, binomial, ., method = "REML")
}
```

And the wrapper with the simple interface:

```{r}
pipeline6 <- function(df, clinical_signs = NULL, antimicrobials = NULL, filtering, watching, first_wks)
  pipeline5(df, clinical_signs, antimicrobials, filtering, filtering, antimicrobials, clinical_signs, watching, antimicrobials, first_wks, antimicrobials)
```

Let's try it:

```{r}
pipeline6(viparc, filtering = 2, watching = 2, first_wks = 1) %>% 
  summary()
```

## The effects of durations and sets {#end}

### With GLM

Let's define an S3 generic function that extracts all the information that we
are interested in from the output of an object of class `glm` or `mgcv`:

```{r}
get_outputs <- function (x, ...) {
   UseMethod("get_outputs", x)
 }
```

Let's now define a `glm` method for that generic function:

```{r}
get_outputs.glm <- function(x) {
  data <- x$model
  m <- with(data, table(amu, disease))
  conf_int <- confint(x)
  list(n    = nrow(data),
       ct   = addmargins(m),
       prob = round(sweep(m, 1, rowSums(m), "/"), 2),
       or   = exp(tail(coef(x), 1)),
       p    = tail(anova(x, test = "LRT")$`Pr(>Chi)`, 1),
       ci   = conf_int[nrow(conf_int), ])
}
```

Let's try it:

```{r}
get_outputs(pipeline4(viparc, filtering = 2, watching = 2, first_wks = 1))
```

Let's now define a batch of durations values:

```{r}
(batch <- expand.grid(filtering = 1:2, watching = 1:2, first_wks = 1:2))
```

And let's run `pipeline4()` on this batch:

```{r}
batch %>% 
  purrr::pmap(pipeline4, df = viparc) %>% 
  lapply(get_outputs) %>% 
  head(2)
```

### With GAM

Let's now define a `gam` method for the `get_outputs()` generic:

```{r}
get_outputs.gam <- function(x, level = .95) {
  data <- x$model
  m <- with(data, table(amu, disease))
  summ <- summary(x)
  conf_int <- mgcv.helper::confint.gam(x, level = level)
  list(n    = nrow(data),
       ct   = addmargins(m),
       prob = round(sweep(m, 1, rowSums(m), "/"), 2),
       or   = unname(tail(summ$p.coef, 1)),
       p    = unname(tail(summ$p.pv, 1)),
       ci   = unlist(conf_int[nrow(conf_int), 5:6]))
}
```

Let's try it:

```{r}
get_outputs(pipeline6(viparc, filtering = 2, watching = 2, first_wks = 1))
```

Let's now run it on a batch of durations values, and put the output in shape:

```{r}
batch %>% 
  tibble::as_tibble() %>%
  assign("plan", ., 1) %>% # saving intermediate value
  purrr::pmap(pipeline6, df = viparc) %>% 
  lapply(get_outputs) %>% 
  purrr::transpose() %>%
  tibble::as_tibble() %>% 
  dplyr::mutate_at(c("n", "or", "p"), unlist) %>% 
  dplyr::bind_cols(plan, .) # intermediate value re-used here
```

### `run_plan()` and Co

Let's make a function out of this above code:

```{r}
run_plan <- function(f, df, clinical_signs = NULL, antimicrobials = NULL, filtering = 1, watching = 1, first_wks = 1, nb_cores = 1) {
  require(magrittr)

  var_names <- setdiff(names(df), c("farm", "flock", "week"))  
  all_antimicrobials <- list(grep("_use$", var_names, value = TRUE))
  all_clinical_signs <- list(grep("_use$", var_names, value = TRUE, invert = TRUE))

  if (is.null(antimicrobials)) antimicrobials <- all_antimicrobials
  else if (! is.list(antimicrobials)) antimicrobials <- list(antimicrobials)
  if (is.null(clinical_signs)) clinical_signs <- all_clinical_signs
  else if (! is.list(clinical_signs)) clinical_signs <- list(clinical_signs)
  
  is_na <- function(x) {
    if (length(x) > 1) return(FALSE)
    is.na(x)
  }
  sel <- which(sapply(antimicrobials, is_na))
  if (length(sel)) antimicrobials[sel] <- all_antimicrobials
  sel <- which(sapply(clinical_signs, is_na))
  if (length(sel)) clinical_signs[sel] <- all_clinical_signs
  rm(is_na, sel, all_antimicrobials, all_clinical_signs, var_names)

  run <- function(x) {
    if (nb_cores < 2) return(purrr::pmap(x, f, df = df))
    c(FUN = f, x, MoreArgs = list(list(df = df)), SIMPLIFY = FALSE, mc.cores = nb_cores) %>% 
        do.call(parallel::mcmapply, .)
  }
  
  environment() %>% 
    as.list() %>% 
    `[<-`(c("f", "df", "run", "nb_cores"), NULL) %>% 
    expand.grid() %>% 
    assign("plan", ., 1) %>% # saving intermediate value
    run() %>% 
    lapply(get_outputs) %>%
    purrr::transpose() %>%
    tibble::as_tibble() %>%
    dplyr::mutate_at(c("n", "or", "p"), unlist) %>%
    dplyr::bind_cols(tibble::as_tibble(plan), .) %>% # reusing the intermediate value
    dplyr::mutate(lower = unlist(lapply(ci, head, 1)),
                  upper = unlist(lapply(ci, tail, 1))) %>% 
    dplyr::select(-ci)
}
```

Let's also write two functions that will help the exploration of the outputs of
the `run_plan()` function. The first one develops the list of clinical signs and
/ or antimicrobials so that they can easily be read in the data frame:

```{r}
develop <- function(df, var = NULL) {
  names_df <- names(df)
  var_ref <- c("clinical_signs", "antimicrobials")
  if (is.null(var)) var <- intersect(names_df, var_ref)
  stopifnot(var %in% names_df)
  stopifnot(var %in% var_ref)
  dplyr::mutate_at(df, var, ~ sapply(., paste, collapse = ", "))
}
```

The second one removes the columns that contain no variability:

```{r}
simplify <- function(df) {
  df[, sapply(df, function(x) length(unique(x))) > 1]
}
```

Let's now try these three functions. With default values:

```{r}
run_plan(pipeline6, viparc)
```

By default, the set of all clinical signs and the set of all antimicrobials are
considered, and the durations of the begining of the flock, the filtering period,
and the observation period are fixed to 1 week. This can be changed by assigning
values to the arguments of the `run_plan()` function. For example, considering 2
different durations of the filtering period:

```{r}
(tmp <- run_plan(pipeline6, viparc, filtering = 1:2))
```

Note that you can access the contingency tables this way:

```{r}
tmp$ct
```

Access to the probabilities would be done the same way:

```{r}
tmp$prob
```

and same thing for the list of clinical signs and antimicrobials. If you want to
see only the columns for which values change:

```{r}
simplify(tmp)
```

Let's see a more complex example:

```{r}
(tmp <- run_plan(pipeline6, viparc, filtering = 1:2, watching = 1:2,
                 antimicrobials = list("colistin_use", c("colistin_use", "cefotaxime_use")),
                 clinical_signs = "respiratory", first_wks = 1))
```

This output can be visualized with the `simplify()` function:

```{r}
simplify(tmp)
```

the `develop()` function:

```{r}
develop(tmp)
```

or the two in combination:

```{r}
simplify(develop(tmp))
```

which is equivalent to:

```{r}
develop(simplify(tmp))
```

### Generating tables

Let's download the data on the antimicrobials classes:

```{r}
ab_classes <- readr::read_csv("https://raw.githubusercontent.com/viparc/clires_data/master/data/antimicrobial_classes.csv", col_types = "cc")
ab_classes$antimicrobial <- paste0(ab_classes$antimicrobial, "_use")
```

It looks like this:

```{r}
ab_classes
```

We can use the `run_plan()` function to compute odds ratios across
antimicrobials classes, clinical signs and values 1 to 3 for the durations in
weeks of the initial, filtering and watching periods (it takes about 15' to run
on 3 cores):

```{r eval = FALSE}
all_cs <- function(x) {
  sel <- sapply(x, length)
  x[which(sel > 1)] <- "all"
  x
}

replace_empty <- function(x) {
  x[which(x == "")] <- "all"
  x
}

tmp <- run_plan(pipeline6, viparc, filtering = 1:3, watching = 1:3,
                antimicrobials = c(split(ab_classes$antimicrobial, ab_classes$class), NA),
                clinical_signs = c(as.list(setdiff(grep("_use$", names(viparc), value = TRUE, invert = TRUE), c("farm", "flock", "week"))), NA),
                first_wks = 1:3, nb_cores = 3) %>% 
  dplyr::mutate(clinical_signs = unlist(all_cs(clinical_signs)),
                antimicrobials = replace_empty(names(antimicrobials))) %>% 
  dplyr::rename(class = antimicrobials)
```

We can the split this output by the durations of the initial, filtering and
watching periods, and spread the data to show antimicrobials classes in rows and
clinical signs in columns:

```{r eval = FALSE}
tmp %>% 
  dplyr::mutate(or    = exp(or),
                lower = exp(lower),
                upper = exp(upper)) %>% 
  dplyr::group_by(filtering, watching, first_wks) %>% 
  dplyr::group_split() %>% 
  lapply(. %>%
           dplyr::select(clinical_signs, class, or) %>% 
           tidyr::pivot_wider(names_from = clinical_signs, values_from = or))
```

## Session info

The analyses were performed on the system and versions as defined below:

```{r session_info}
sessionInfo(required)
```

If you want to run the code locally, you can copy and paste code chunks to your
commands line. Make sure though that your system and packages versions are the
same or more recent than the ones listed above.

## Running the script locally

You can download the R markdown file used to generate this page as so:

```{r eval = FALSE}
download.file("https://raw.githubusercontent.com/viparc/prophylactic/master/prophylactic.Rmd", "prophylactic.Rmd")
```

You can transforme this `prophylactic.Rmd` R markdown document in an R script
named `prophylactic.R` as so:

```{r eval = FALSE}
if (!require("kintr")) install.packages("knitr")
knitr::purl("prophylactic.Rmd")
script <- readLines("prophylactic.R")
writeLines(script[-c(seq_len(which(grepl("starts here", script))[1]),
                     which(grepl("session_info", script))[1]:length(script))],
           "prophylactic.R")
```

To execute the code:

```{r eval = FALSE}
source("prophylactic.R")
```