%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Source Attribution: The R Package sourceR}
%devtools::build_vignettes()
\documentclass[article, nojss]{jss}

% library(knitr)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %% almost as usual
 \author{Poppy Miller\\Massey University \And 
         Jonathan Marshall\\Massey University \And
         Nigel French\\Massey University \And
         Chris Jewell\\Massey University}
 \title{Source Attribution: The R Package \pkg{sourceR}}

%for pretty printing and a nice hypersummary also set:
\Plainauthor{Poppy Miller, First Author, 
Nigel French, Second Author,
Jonathan Marshall, Third Author,
Chris Jewell, Fourth Author} %% comma-separated
\Plaintitle{Source Attribution: The R Package SourceR} %% without formatting
\Shorttitle{\pkg{sourceR}: A Source Attribution Model in R} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{
Numerous zoonotic diseases cause morbidity, mortality and productivity losses in both humans and animal populations. For many zoonotic diseases that are important to human health (such as campylobacteriosis), it is difficult to attribute human cases to sources because there is little epidemiological information on the cases. Genotyping systems allow the zoonotic pathogens to be categorised, and 
the relative distribution of the genotypes among the sources (food sources or reservoirs of bacteria) and in human cases allows inference on the likely source of each genotype. Current source 
attribution models, specifically the Island model \citep{WilGabLea08}, Hald \citep{HaldVosWed04} and modified Hald models \citep{MulJonNob09} are not fully joint and have many (often unverifiable) 
assumptions. 
  Identifiability of the parameters in this model is an issue because a large number of parameters need to be estimated, the data is imbalanced, and many of the combinations of source and type and 
have very low counts. We present techniques to overcome these issues within a Bayesian framework by developing a fully joint model which non-parametrically clusters the type effects (using a 
Dirichlet Process) allowing identification of groups of bacterial subtypes with similar pathogenicity, survival and/ or virulence mechanisms. This model is applied to \emph{Campylobacter} data 
from the Manawatu area of New Zealand (previously analysed by \cite{MulJonNob09}) using the sourceR package, and compared to current source attribution models.
}
\Keywords{bayesian, source attribution, food-bourne zoonoses, non-parametric, \proglang{R}}
\Plainkeywords{bayesian, source attribution, food-bourne zoonoses, non-parametric, R} 

%% The address of (at least) one author should be given
%% in the following format:\textmd{}
\Address{
Poppy Miller\\
  CHICAS, Faculty of Health and Medicine,\\
  Furness College,\\
  Lancaster University,\\
  Lancaster, LA1 4YG,\\
  United Kingdom\\
  E-mail: \email{p.miller@lancaster.ac.uk}\\
}

%% for those who use Sweave please include the following line (with % symbols):
%% need no \usepackage{Sweave.sty}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% My preamble
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[english]{babel}
\usepackage{tikz}
%library(tikzDevice)
% options(tikzDefaultEngine='xetex')
\usepackage{todonotes}
\usepackage{algpseudocode}
\usepackage{algorithm}
\algnewcommand\algorithmicforeach{\textbf{for each}}
\algdef{S}[FOR]{ForEach}[1]{\algorithmicforeach\ #1\ \algorithmicdo}
%\usepackage{algorithmic}

\usepackage{alltt}
\usepackage{bm}
\usepackage{listings} %% to put code chunks in
%\usepackage{natbib}

%% double space lines
%\usepackage{setspace}
%\doublespacing

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}

\begin{document}
% \SweaveOpts{concordance=TRUE}
%opts_chunk$set(concordance=TRUE)
library(tikzDevice)

% options(tikzDefaultEngine='xetex')
\input{child_docs/Introduction}
\input{child_docs/Motivating_ex}
\input{child_docs/Review}
\input{child_docs/Methods}





\section{Case Studies}

In this section we demonstrate the use of the \pkg{sourceR} package using two example data sets. The first uses simulated data to demonstrate the general use case, with both time and location extensions for the model. The second is a specific case study using the the data described in the motivation section (Section \ref{sec:motivation}). The results (proportion of cases attributed to each source) are compared to the results from the Dutch, Modified Hald and Island models.

\subsection{Simulated data} \label{sim_study_section}

In this section, we provide a worked example using simulated data with multiple times and locations for source attribution data generated from the model in Section \ref{sec:model} (available in the \pkg{sourceR} data sets). There are two times (1, 2) and two locations (A, B) over which the human cases vary. The data expected by \code{saBayes} is in long format, with a column for the number of human cases, the number of positive samples for each source, and columns identifying the type, time and locations. Note, the source data is the same for all locations within a time. For this data, the source prevalences ($\pi_j$) are all set to be 1. If source prevalences are not provided, \code{saBayes} will automatically set them all to 1 (with a warning). The data must be in long format, with columns giving the number of human cases for each type, a column for each of the sources giving the number of positive samples for each type, and columns giving the time, location and type id's for each observation. 

\begin{table}
\caption{Minimum parameters required for the \code{saBayes} function.}
\label{table:HaldModelMinParams}
\begin{tabular}{p{4cm}p{11cm}}
\code{formula} & A formula object of the form $y\sim x1 + x2 + ... + xJ$, where $y$ is the name of the human cases column, and $x1, ...,xJ$ are the names of the source count columns in the data.\\
\code{time} & A formula object of the form $\sim t$, where $t$ is the name of the column containing the times in the data.\\
\code{location} & A formula object of the form $\sim l$, where $l$ is the name of the column containing the locations in the data.\\
\code{type} & A formula object of the form $\sim s$, where $s$ is the name of the column containing the (sub) types in the data.\\
\code{data} & Correctly formatted data.\\
\code{priors} & List with parameters for the prior distributions for each of the model parameters.\\
\code{n\_iter} & Specifies the number of iterations to run the algorithm for.\\
\code{likelihood_dist} & Specifies the likelihood distribution to be used for the human cases from each type. Must be one of \code{nbinom} or \code{pois}.\\
\end{tabular}
\end{table}

<<run_doMCMC_sim_pois, results='hide', message=FALSE, eval = FALSE, cache = T>>=
require(sourceR)
set.seed(63164)
data(sim_SA)
data(sim_SA_true)

# Set priors
priors <- list(a = 1, r = 1, theta = c(0.01, 0.00001))

# Run model
res_sim <- saBayes(formula = Human ~ Source1 + Source2 + Source3 + Source4 + Source5, 
               time = ~Time, location = ~Location, type = ~Type,
               data = sim_SA$data, priors = priors,
               alpha_conc = 1, prev = sim_SA$prev,
               likelihood_dist = "pois", n_iter = 1010,
               mcmc_params = list(burn_in = 20, thin = 1))
@

The algorithm is run for 102,000 iterations using the \pkg{sourceR} command, with an initial burn in of 2000 iterations, followed by a further 100,000 iterations, of which every 100\textsuperscript{th} sample is saved. The acceptance rates for all parameters (except those updated using a Gibbs sampler) can be found in a list called \code{acceptance} in the output from \code{saBayes}. Trace and autocorrelation plots for the parameters (Figure~\ref{fig:trace_acf_sim_data_plots}) indicate that the Markov chain is mixing well and has converged, and that thinning by 100 is adequate. The posteriors are returned as nested lists for each parameter. The following R code demonstrates how to access the posteriors for a given source $a_{jtl}$ or type $q_i$ effects and a relative prevalence $r_{ijt}$.

<<trace_acf_sim_data_code, dev='tikz',results='hide', eval=FALSE, cache = T>>=
## Plot the marginal posterior for the source effect 2, at time 1, location A
plot(res_sim$posterior$a$time1$locationA[,"Source3"], type="l")
## Plot the marginal posterior for the type effect 21
plot(res_sim$posterior$q[,"type21"], type="l")
## Plot the marginal posterior for the relative prevalence of source effect 5, 
## type 17, at time 2
plot(res_sim$posterior$r$time2["type17","Source5",], type="l")
@


%\noindent\begin{minipage}{\linewidth}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{figure}[H]

{\centering \includegraphics[width=\maxwidth]{trace_acf_sim_data_plots-1} 

}

\caption{Trace and acf plots for a sample of the model parameters. True values of the parameters are shown in red.}
\label{fig:trace_acf_sim_data_plots}
\end{figure}
\end{knitrout}


Medians and Chen-Shao highest posterior density credible intervals \citep{ChenShao99} can be obtained for each parameter using the \code{summary} command. 
%The credible intervals are Chen-Shao highest posterior density intervals which are the shortest possible Bayesian credible intervals, though 
%they can have high Monte Carlo error when  estimated directly from simulations via empirical shortest intervals \citep{LiuGelZhe15}. 

<<summary_sim_pois, results='hide', message=FALSE, eval = FALSE, cache = T>>=
summary(res_sim, alpha = 0.05, thin = 1, burn_in = 0)
@

The data can be subsetted using the \code{subset_posterior} command.

<<subset_sim_pois, results='hide', message=FALSE, eval = FALSE, cache = T>>=
subset_posterior(res_sim, params = c("a", "li", "q"), 
                 t = "1", l = "B", j = c("Source2", "Source1"), 
                 i = c("47", "10"), iters = c(3:10))
@

Both the full posterior and a subset of the posterior (generated using \code{subset_posterior}) can be flattened into a data frame.

<<flatten_sim_pois, results='hide', message=FALSE, eval = FALSE, cache = T>>=
flatten(res_sim)
@

The marginal density plots of the proportion of cases attributed to each source at each time and location ($\lambda_{jtl}$) show that the true values (shown by a cross on the graph) are within the 
credible intervals (Figure~\ref{fig:sim_pois_lambdaj_plots}). The residual plots for $\lambda_{i}$ (Figure~
\ref{fig:lambda_i_residuals}) show that the model is fitting well. The heatmap shows the grouping of the type effects (Figure~\ref{fig:type_effect_heatmap}) computed using a dissimilarity matrix from the clustering output of the mcmc. The coloured bar under the dendrogram gives the correct grouping from the simulated data. This shows that 
all the types have been grouped correctly if the dendrogram is cut at the true number of groups (5).

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{sim_pois_lambdaj_plots-1} \caption{Proportion of cases attributable to each source for each time (1, 2) and location (A, B) for simulated data. Error bars represent 95\% Chen-Shao credible intervals. True $\lambda_j$ values 
are shown as crosses.}\label{fig:sim_pois_lambdaj_plots}}
\end{figure}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{lambda_i_residuals-1} \caption{Residual plots (simulated data).}\label{fig:lambda_i_residuals}}
\end{figure}


\end{knitrout}



\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}


\end{kframe}

\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{type_effect_heatmap-1} \caption{Heatmap showing the grouping of the type effects (q) using simulated data (true groupings given by the 5 colours in the bar under the 
dendrogram).}\label{fig:type_effect_heatmap}}
\end{figure}


\end{knitrout}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Campylobacteriosis cases in the Manawatu (2005-2008)} \label{real_data_case_study_section}

In this section, we apply \pkg{sourceR} to the \code{campy} (campylobacteriosis) dataset from Manawatu, New Zealand, described in Section \ref{sec:motivation}.  We compare the results of our Bayesian non-parametric approach with results from the Modified Hald and Island models. Types which do not have any source cases need to be removed from the data set before running the analysis because there is no information to attribute human cases to a source if the subtype only occurs in humans.

<<run_doMCMC_real, results='hide', eval = FALSE, message=FALSE, cache = T>>=
data(campy)
set.seed(59623)
# remove rows with no source cases as there is no information for
# source attribution of human cases for these sources
zero_rows <- which(apply(campy[,c(2 : 7)], 1, sum) == 0)
campy <- campy[-zero_rows,]

# Set priors
priors <- list(a = 1, r = 1, theta = c(0.01, 0.00001))

# set prevalences
# Number of samples  tested  for c. jejuni, for each source.
tot_samples<-c(239, 196, 127, 595, 552, 192 + 332)
# Number of samples positive for c. jejuni, for each source.
pos_samples<-c(181, 113, 109, 97, 165, 24 + 62)   
prevs <- data.frame(value = pos_samples / tot_samples, 
                    source_id = colnames(campy[, 2:7]))

# Run model
# the model assumes one time and location if none are specified in saBayes
res_real <- saBayes(formula = Human ~ ChickenA + ChickenB + ChickenC + 
                 Bovine + Ovine + Environment, 
               type = ~Type, data = campy, priors = priors, alpha_conc = 1, 
               prev = prevs, likelihood_dist = "pois", n_iter = 1020,
               mcmc_params = list(burn_in = 20, thin = 1))
@

Trace and autocorrelation plots for the parameters indicate that the Markov chain is mixing well and has converged, and that thinning by 500 is adequate for most of the parameters (Figure~
\ref{fig:trace_params_real}). The residual plots for the $\lambda_{i}$s (Figure~\ref{fig:lambda_i_residuals_real}) show that the model is fitting the data well. The proportion of cases attributed to each source ($\lambda_j$) using the new model can be compared to the previous models (Figure~\ref{fig:lambda_j_real}). The new 
model has very similar medians to the modified Hald and Island models. The credible intervals are much narrower than the modified Hald model, but still relatively wide compared to the Island model.
The heatmap and dendrogram of the type effects (Figure~
\ref{fig:type_effect_heatmap_real}) shows that there are 3 main groups of type effects. The violin plots (Figure~\ref{fig:type_effect_violinplots_real}) show that the largest group of types have very small type effects. These correspond to types that are observed in source samples, but no human cases. There is a group of 5 types which have very large type effects (including type 474 which is endemic to NZ and largely associated 
with poultry). Although the clustering was determined without reference to genetic relatedness of the types, three members of this group (subtypes 38, 48 and 474) are members of the same clonal complex (CC48) and therefore genetically closely related \citep{pubmlstCampy}. Subtype 52 was frequently placed in both the groups with the largest and middling type effects (as can be seen in Figure\ref{fig:type_effect_heatmap_real}, although overall it was attributed to the group with the largest type effects.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{trace_params_real-1} \caption{Trace and acf plots for a sample of the model 
parameters (Manawatu \emph{Campylobacter} data).}\label{fig:trace_params_real}}
\end{figure}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{lambda_i_residuals_real-1} \caption{Residual plots (Manawatu \emph{Campylobacter} data).}\label{fig:lambda_i_residuals_real}}
\end{figure}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{lambda_j_real-1} \caption{Proportion of human campylobacteriosis cases 
attributable to each source (Manawatu \emph{Campylobacter} data). Error bars represent 95\% confidence or credible intervals.}\label{fig:lambda_j_real}}
\end{figure}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\end{kframe}\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{type_effect_heatmap_real-1} \caption{Heatmap showing the grouping of the type effects (q) using the Manawatu \emph{Campylobacter} data.}
\label{fig:type_effect_heatmap_real}}
\end{figure}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{figure}
\centering{
\includegraphics[width=\maxwidth]{type_effect_boxplots_real1-1} \caption{Violin plots of the type effects (q) using the Manawatu \emph{Campylobacter} data.}\label{fig:type_effect_violinplots_real}}
\end{figure}
\end{knitrout}






\input{child_docs/Discussion}
\input{child_docs/Conclusions}

\section*{Acknowledgements}
The research for this paper was financially supported by the Ministry for Primary Industries, the Institute of Fundamental Sciences (Massey University), the mEpiLab (Massey University), and CHICAS (Lancaster University). This project was also supported by the Livestock Improvement Corporation (LIC) and Soroptimist International Palmerston North. We acknowledge the following individuals and groups: mEpiLab
(Massey University), MidCentral Public Health Services and Petra Mullner (for the Manawatu data set) and Geoff Jones (for his helpful input on automatic clustering methods).

\bibliography{source_attribution_jss}
\appendix
\input{child_docs/McMC_alg}

\end{document}