---
title: "nncc: nearest-neighbors matching for case-control data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{nncc: nearest-neighbors matching for case-control data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=6, 
  fig.height=4
)
```

```{r, warning=FALSE, message=FALSE}
library(nncc)
library(survival)
library(dplyr)
library(ggplot2)
data(anifood)
```

## 1. Prepare the data

Case-control studies often collect many exposures but have relatively small 
sample sizes in comparison to the number of exposures. Conventional analysis
techniques like multivariate logistic regression can only accommodate a limited 
number of exposures. The `nncc` package was designed to help address this 
limitation. `nncc` matches a case to its nearest available neighbors based 
on a user-defined measure of distance (by default Gower distance) calculated 
using all collected confounders for a given exposure. In the end, `nncc` creates 
a matched data set for each exposure of interest. 

Ideally, one prefers to follow the analysis plan pre-defined at study design.
Decisions to change your analysis plan to this approach should be considered
carefully in light of how they may affect bias and uncertainty in model estimates. 
In addition, constructing variables are a routine and sometimes necessary part of 
multifactor studies. However, constructing new outcome variables from primary 
variables will increase dependencies and the difficulty of interpretation, 
especially where there is missingness of data. 

For illustration, we created a toy example case-control data set called `anifood` which only contains 
`r table(anifood$case)["1"]` cases, `r table(anifood$case)["0"]` controls, and `r NCOL(anifood)` variables including the 
one indicating case status, `case`, which equals 0 for controls, 1 for cases.

```{r}
dim(anifood)
```

In addition to the data set, `nncc` needs to know the exposures of interest, 
the variables used for matching, and a data set that contains 
variables to be excluded from matching for each exposures (e.g., possible 
intermediates in the pathway between an exposure and case status). 

```{r}
# exposures of interest. In a real study, this list can be much longer
exp_interest <- c("exp01","exp09", "exp27")

# exposures to be controlled for any exposure of interest
exp_match <- setdiff(names(anifood), "case")

# variables to be excluded from matching for each exposure
# both exp_var and rm_vars are character variables
excl_vars %>% head
```

## 2. Establish a threshold for distance

The first step is to establish a proper distance threshold for matching. 
If too loose, there will be poor confounding control; if too tight, it will
be difficult to find enough controls for the cases due to overmatching.

The `get_threshold` function uses the maximum bipartite graph algorithm to find
each case's closest control while ensuring each control is used no more than 
once. Then, for each case, a control is selected randomly. 

A logistic regression on whether the control is the closest or the randomly
selected one based on the distance is created. By default, the threshold used
is the distance at which the probability is at least 50% that the control is
the closest vs. a randomly selected one (see the plot below generated by the 
function `threshold_model_plot`). You can choose a different probability with 
the `p_threshold` argument. 

```{r}
threshold_results <- get_threshold(anifood, exp_match, p_threshold = 0.50)
```

The `distance_density_plot` function plots the density distribution of Gower distances 
for the maximum bipartite graph algorithm matches (solid line), the density
distribution of distances for the random matches (dashed line), and the threshold (red dashed line).
As expected, a maximum bipartite graph algorithm matching results in closer 
matching than a random matching. 

```{r, warning=FALSE}
distance_density_plot(threshold_results) + ggtitle("Example of distance_density_plot")
```

The function `threshold_model_plot` plots the probability of being the maximum 
bipartite graph algorithm match vs. a random one by distance. It also plots the 
threshold for distance corresponding to a probability specified through 
the `p_threshold` argument.

```{r}
threshold_model_plot(threshold_results, p_threshold = 0.50) + ggtitle("Example of threshold_model_plot")
```

This package can also be used to analyze data from a case-control study that
was originally matched. The `original_compare_plot` function makes a density plot of 
the distance for originally matched pairs. It also gives the proportion 
of originally matched pairs that have a distance greater than the threshold.

```{r}
# create a variable (i.e., pair) indicating originally matched pairs
anifood_matched <- anifood %>% group_by(case) %>% mutate(pair = seq_along(case)) %>% ungroup
```

In this example, `r p <- original_compare_plot(anifood_matched, case, pair, threshold_results); sprintf("%.0f%%", p$prop_distance_gt_threshold[2] * 100)` of original matched case-control pairs have a distance 
greater than the threshold.

```{r}
p <- original_compare_plot(anifood_matched, case, pair, threshold_results)


# the density plot of distance between originally matched cases and controls
p$plot_density + ggtitle("Example of original_compare_plot")

# proportion of originally matched cases and controls with a distance greater than the threshold
p$prop_distance_gt_threshold
```

## 3. Create strata: candidate matches for final analysis

The `make_knn_strata` function calculates Gower distances, by default, for pairs 
formed by every case and every control for a given exposure and creates strata
by selecting the closest 250 controls for each case (by default, can be modified 
through the `ncntls` argument). Thus, for that exposure, the number of rows in 
data frame equals the number of cases * 251. The data frames for all exposures 
are stored in a list with a length equal to `length(exp_interest)`. 

This step can be computationally intensive depending on the amount of data. 
`furrr::future_map()` can be used so that the computation can be easily adjusted
to use multiple cores on your local computer or a high-performance computing 
(HPC). More information is provided [below][9. Run the analysis on your local computer or HPC].

```{r, warning=FALSE, message=FALSE}
library(furrr)
strata1 <- future_map(exp_interest, make_knn_strata,  rmvars = excl_vars, matchvars = exp_match, df = anifood) %>%
  setNames(exp_interest)
```

```{r, warning=FALSE}
length(strata1) == length(exp_interest)

# rows in a matched data set
all.equal(anifood %>% filter(case == 1) %>% NROW %>% `*`(250 + 1), NROW(strata1[[1]]))
```

In the `make_knn_strata` function, a data frame should be supplied to the `rmvars` 
argument (see `excl_vars` above). That data frame should include two variables:
*exp_var: a character variable (only a single exposure in a row), and
*rm_vars: a character variable that contains the names of the variables that are 
not adjusted for. Names are separated by a space.


## 4. Prepare final matched sets for analysis

For a given exposure, the `make_analysis_set` function selects the closest 
20 controls for each case; subset these 20 to those that fall within the threshold; 
collapse strata that share the same controls; and remove strata without any control.

```{r, warning = FALSE}
strata2 <- future_map(exp_interest, make_analysis_set, stratified_data = strata1, data = anifood, maxdist = threshold_results$threshold) %>% setNames(exp_interest)
```

The `finalize_data` function ensures that a control retained in a data frame is 
used once; removes strata without any case or any control. In this process, 
priority is given to the smallest strata then smallest distance if a control 
is matched to multiple cases (i.e., that control exists in multiple strata). 
```{r, message=FALSE}
strata3 <- finalize_data(strata2)
```

The last step is to exclude exposures, if any, with odds ratios that are not 
estimable (e.g., none of cases and controls was exposed). 
For example, `exp09` in this example.

```{r}
# exposures to which neither cases nor controls were exposed
strata3 %>%
  lapply(function(dfm) dfm %>%
                            mutate(case = as.character(case), exp = as.character(exp)) %>%
                            filter(case != "" & exp != "") %>% 
                            with(table(case, exp))) %>%
  lapply(function(x) x==0) %>%
  lapply(function(x) {sum = sum(x); length = length(x); cbind(sum, length)}) %>%
  do.call(rbind,.) %>%
  as.data.frame %>%
  mutate(var = names(strata3)) %>%
  filter(sum >= 2 | length < 4) %>% select(var) %>%
  unclass %>%
  unlist -> expvars_invalid

expvars_invalid

# none exposed and odds ratio cannot be estimated
strata3[["exp09"]] %>% with(table(case, exp))

data_final <- strata3[setdiff(names(strata3), expvars_invalid)]
```


## 5. Calculate odds ratios

### 5.1 Mantel–Haenszel test and Fisher's exact test
The `test_mh` function performs the Mantel–Haenszel test for data frame with 
more than one stratum and the Fisher's exact test for a data frame with only 
one stratum.

```{r}
or_mh <- future_map(data_final, function(dfm) with(dfm, test_mh(case = case, exp = exp, strata = strata)))

or_mh[["exp01"]]
```

However, when the number of exposed is small, this method may yield an
infinite result.

```{r}
data_final[["exp27"]] %>% select(case, exp) %>% table

or_mh[["exp27"]]
```

### 5.2 Conditional logistic regression 
Conditional logistic regression gives results similar to those from 
Mantel–Haenszel test, but may also give less than ideal
results when few participants are exposed.
```{r, warning=FALSE}
results_clogit <- future_map(data_final, function(dfm){clogit(case ~ exp + strata(strata) , data = dfm)}) 

results_clogit[["exp27"]] %>% summary() %>% `$`(conf.int)
```

When combined with multiple imputation, we recommend a Bayesian approach.

```{r, warning=FALSE, eval=FALSE}
library(rstanarm)

results_stan_clogit <- future_map(data_final, function(dfm) {
  dfm$strata <- factor(dfm$strata)
  stan_clogit(case ~ exp, 
              data = dfm, 
              strata = strata,
              # rstanarm suggests us specifying priors even if they are the 
              # defaults because they might change in the future
              prior = normal(0, 4),
              #default: prior = normal(autoscale = TRUE),
              prior_covariance = decov(),
              iter = 100,
              chains = 2,
              prior_PD = FALSE)}) # for speed only
```

### 5.3 Penalized logistic regression
When the number of participants being exposed is small, [Firth's bias-Reduced penalized-likelihood logistic regression](https://CRAN.R-project.org/package=logistf/index.html) can also be used (1). 

```{r, warning=FALSE}
library(logistf)
# regression coefficients
coef_logistf <- future_map(data_final, function(dfm){

  if(length(unique(dfm[["strata"]])) == 1) {
        x <- model.matrix(case ~ exp, data = dfm)
        plconf <- grep("exp", colnames(x))
        o <- logistf(case ~ exp, data = dfm, plconf = plconf) %>%
                 `[`(c("terms", "coefficients", "ci.lower", "ci.upper", "prob", "call", "loglik", "model"))
    } else {
        x <- model.matrix(case ~ exp + strata, data = dfm)
        plconf <- grep("exp", colnames(x))
        o <- logistf(case ~ exp + strata, data = dfm, plconf = plconf) %>%
                 `[`(c("terms", "coefficients", "ci.lower", "ci.upper", "prob", "call", "loglik", "model"))
    }
},
.progress = TRUE)


# odds ratios
or_logistf <- lapply(names(coef_logistf), function(var_name){
  coef_logistf[[var_name]] %>%
    `[`(c("terms", "coefficients", "ci.lower", "ci.upper", "prob")) %>%
    bind_cols %>%
    filter(grepl("exp", terms)) %>%
    mutate(variable = var_name, or = exp(coefficients), ci.lower = exp(ci.lower), ci.upper = exp(ci.upper)) %>%
  select(variable, terms, or, ci.lower, ci.upper, prob)}) %>%
  setNames(names(coef_logistf))

# odds ratio for exp27
or_logistf[["exp27"]] 
```


## 6. Calculate population attributable fraction (PAF)
PAF can be estimated using the method described by Bruzzi and colleagues (3). 
We implemented this method in `get_paf`.
```{r}

# prepare a data frame for calculating PAF
df_or <- bind_rows(or_logistf)

df_or
```


```{r, message=FALSE}
# point estimate of PAF
paf <- get_paf(df_or = df_or, which_or = or, exp_var = variable, exp_level = terms, df_matched = data_final)

paf

# lower confidence limit of PAF
paf_ci.lower <- get_paf(df_or = df_or, which_or = ci.lower, exp_var = variable, exp_level = terms, df_matched = data_final)
```


## 7. Multiple imputation of missingness
When the case-control data contains a large proportion of missingness, multiple 
imputation may be useful. 

* For regression coefficients generated by Mantel–Haenszel test or conditional 
logistic regression, Rubin's rule may be proper for  combining the results (4).

* Originally in our paper (2), we used Firth's bias-Reduced penalized-likelihood logistic regression, but 
the regression coefficients are not normally distributed. The results from 
imputed data sets should therefore be combined using penalized likelihood profiles (5). 
Penalized likelihood profiles are implemented in `logistf::CLIP.confint()`, but
the function was developed to combine results from imputed data sets that 
have the same structure. However, if we impute the original case-control data, 
and then apply nearest-neighbors matching to the imputed data, the final 
analytic data sets for a given exposure could have inconsistent structures 
because they could have different numbers of strata. 

Thus, we modified the `logistf::CLIP.confint()` to accommodate our method. 
The modified function is included in version 1.0.1 of the  `nncc` package along with a compatible, bundled version of the 
`logistf` package version 1.24, and the function is called 
`CLIP.confint.difflevel()`. Please cite the original paper (4) and/or the [logistf package](https://CRAN.R-project.org/package=logistf/index.html) if you use 
this modified function for publication.

* Now we recommend a Bayesian approach, based on drawing samples from the stacked posteriors of the model run on each imputed data set.


## 8. Helper function `cacheit()`
When the number of exposures of interest or the sample size of the original 
case-control study is large, it may be time-consuming for some analyses 
(e.g., strata creation). You may not want to have the same code evaluated again 
when you run your analysis later. The helper function `cacheit()` caches the results for you. Specifically, before your run `cacheit()`, create a folder named "cache" in the working directory of your project. When you run `cacheit("abc", code)` for the first time, it saves the results of your code as abc.rds in the `cache` folder; next time, 
when you run the code, it directly read the saved results in without 
evaluating your code.

```{r, eval = FALSE}
strata1 <- cacheit("abc",
                   future_map(exp_interest, make_knn_strata,  rmvars = excl_vars, matchvars = exp_match, df = anifood) %>%
  setNames(exp_interest),
  clearcache = FALSE)
```


## 9. Run the analysis on your local computer or HPC

Considering the intensive computation for this approach, 
especially when multiple imputation is conducted, a [future](https://CRAN.R-project.org/package=future/vignettes/future-1-overview.html) 
can be used to accelerate the analysis. Below we briefly introduce the use of futures on your local multi-core computer and HPC.

### 9.1 local computer
If you have a multi-core computer, the following code can assign a specified number of cores for the analysis.

```{r, eval = FALSE}
library(nncc)
library(dplyr)
library(furrr)
library(future.batchtools)
# the workers argument is used to define the number of cores for the analysis. By default, all cores will be used.
plan(multisession, workers = 3)

strata1 <- future_map(exp_interest, make_knn_strata,  rmvars = excl_vars, matchvars = exp_match, df = anifood) %>%
  setNames(exp_interest)

```


### 9.2 HPC

To run futures on HPC, you may need a template file, a job script, and a R script. The template file and 
the job script are created using a Linux application such as nano and saved as .sh files.

There are some exemplary template files [here](https://github.com/mllg/batchtools/tree/master/inst/templates). 
The template file should be saved as `batchtools.sge.tmpl` in your current 
working directory or as `.batchtools.sge.tmpl` in your home directory on the HPC. 
Below is an example of a template file for Sun Grid Engine.

```
#!/bin/bash

## name of the job
#$ -N <%= job.name %>

## tell the queue system to use the current directory as the working directory
## Or else the script may fail as it will execute in your top level home directory /home/username
#$ -cwd

## Use environment variables
#$ -V

## specify a queue
#$ -q <%= resources$queue %>

## free to add other modules if necessary
module load R/3.6.2

Rscript -e 'batchtools::doJobCollection("<%= uri %>")'
exit 0
```

Below is an example of a job script to be submitted to the scheduler through a shell:
```{r, eval=FALSE, comment = ""}
#!/bin/bash -l

# name of the job is helloR
#$ -N helloR

#$ -cwd

#$ -V

#$ -pe smp 2-12

#$ -q all.q

module load R/3.6.2

# to execute fugure-map.R
Rscript future-map.R

exit 0
```

The R code looks like:

```{r, eval=FALSE}
library(furrr)
library(future.batchtools)

plan(batchtools_sge)

future_map(<a_list_or_vector>, function(x){...}, .progress = TRUE)
```


### References:
1. Firth D. Bias reduction of maximum likelihood estimates. Biometrika 1993;80:27-38.
2. Bruzzi P, Green SB, Byar DP, et al. Estimating the population attributable risk for multiple risk factors using case-control data. Am J Epidemiol 1985;122(5):904-14.
3. Cui Z, Marder EP, Click ES, Hoekstra RM, Bruce BB. Nearest-neighbors matching for case-control study analyses: better risk factor identification from a study of sporadic campylobacteriosis in the United States. Epidemiology 2022;33(5):633-41.  
4. Rubin DB. Multiple imputation for nonresponse in surveys. New York: John Wiley and Sons; 1987.
5. Heinze G, Ploner M, Beyea J. Confidence intervals after multiple imputation: combining profile likelihood information from logistic regressions. Stat Med 2013;32(29):5062-76.

