README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,

  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```
# serosolver

[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)


`serosolver` is a modelling and inference package that uses a dynamic model to infer antibody dynamics and infection histories from cross-sectional or longitudinal serological data. The model infers individual-level infection histories, historical attack rates, and patterns of antibody dynamics by accounting for cross-reactive antibody responses and measurement error.

## Installation
1. Install [R][r-project]

2. Install the development version of serosolver from [GitHub](https://github.com/seroanalytics/serosolver):

```{r installation,eval=FALSE}
#devtools::install_github("seroanalytics/serosolver")
#install.packages("ggpubr")
#library(serosolver)
```
## Quick start and vignettes

Read the [quick start vignette][vignette-doc] to set up and run a simple implementation with a simulation model.

There are additional Rmarkdown vignettes for [Case Study 1][c1-vignette-doc] (longitudinal analysis of influenza A/H1N1p in Hong Kong) and [Case Study 2][c2-vignette-doc] (cross-sectional analysis of influenza A/H3N2 in Guangzhou Province, China), to accompany the analysis in the serosolver paper.

[r-project]: http://cran.r-project.org
[vignette-doc]: https://seroanalytics.github.io/serosolver/articles/serosolver-quick_start_guide.html
[c1-vignette-doc]: https://seroanalytics.github.io/serosolver/articles/cs1_vignette.html
[c2-vignette-doc]: https://seroanalytics.github.io/serosolver/articles/cs2_vignette.html

## Example

This is a basic example of simulating some serological data and fitting the model using the MCMC framework.

```{r example, message=FALSE,warning=FALSE,eval=FALSE}
library(serosolver)
library(plyr)
library(data.table)
library(ggplot2)
library(ggpubr)


## Load in example parameter values and antigenic map
data(example_par_tab)       # list of parameter values in model
data(example_antigenic_map) # antigenic map for each strain each year (bedford version?)

## Get all possible infection times
strain_isolation_times <- unique(example_antigenic_map$inf_times) # bunch of years

## Vector of strains that have titres (note only one representative strain per time)
sampled_viruses <- seq(min(strain_isolation_times), max(strain_isolation_times), by=2) # every two years?

## Times at which serum samples can be taken
sampling_times <- 2010:2015 # sampling times when the samples were taken 

## Number of serum samples taken
n_samps <- 2 

## Simulate some random attack rates
attack_rates <- runif(length(strain_isolation_times), 0.05, 0.15) # random attack rates

## Simulate a full serosurvey with these parameters
## That is simulate infection history, titre, and attack rates using the imported data 
all_simulated_data <- simulate_data(par_tab=example_par_tab, group=1, n_indiv=50,
                                  strain_isolation_times=strain_isolation_times,
                                  measured_strains=sampled_viruses,
                                  sampling_times=2010:2015, nsamps=n_samps,
                                  antigenic_map=example_antigenic_map,
                                  age_min=10,age_max=75,
                                  attack_rates=attack_rates, repeats=2)


## Pull out the simulated titre data and infection histories
titre_dat <- all_simulated_data$data
ages <- all_simulated_data$ages
example_inf_hist <- all_simulated_data$infection_histories
example_titre_dat <- merge(titre_dat, ages)

## Run the MCMC
# This example uses prior version 2 (i.e. beta prior on phi with parameters alpha, beta)
# We have to remove the explicit specification of phi in the parameter table
par_tab <- example_par_tab[example_par_tab$names != "phi",]
res <- run_MCMC(par_tab, example_titre_dat, example_antigenic_map, 
                filename="test", version=2,
                mcmc_pars=c(adaptive_period=20000, iterations=80000, 
                            inf_propn=0.5,hist_sample_prob=0.5,
                            save_block=10000,thin=10,thin_hist=100))

## Read in the MCMC chains and plot posteriors
chain <- read.csv(res$chain_file) # test_chain.csv in folder
inf_chain <- data.table::fread(res$history_file) # test_infection_histories.csv in folder
plot(coda::as.mcmc(chain[chain$sampno > 20000, c("mu","wane","lnlike")]))

# Plot model predicted titres for a subset of individuals
plot_infection_histories_long(chain = chain,infection_histories = inf_chain,
                         titre_dat = example_titre_dat,individuals=c(1:4),
                         antigenic_map=example_antigenic_map,par_tab=par_tab)

```