\documentclass[a4paper]{article}

\def\likeLTD{\textsc{LikeLTD}}
\usepackage{fullpage} 

%\VignetteIndexEntry{Using the likeLTD peak heights model}
\title{Using the \likeLTD\ peak heights model}
\author{Christopher Steele, Adrian Timpson, Mayeul d'Avezac, \\James Hetherington, David Balding}

\begin{document}
\maketitle
<<setthreads,echo=FALSE,print=FALSE,include=FALSE>>=
  if(.Call(likeLTD::.cpp.nbthreads) > 2) {
    .Call(likeLTD::.cpp.set_nbthreads, as.integer(2))
  }
@

This vignette assumes familiarity with both programming in R and forensic DNA analysis, only demonstrating how to run an analysis. For a more in depth discussion of \likeLTD\ see {\tt likeLTDguide.pdf} provided with the package.

\section{Inputs}
\likeLTD\ evaluates the weight-of-evidence (WoE) against a queried contributor (Q) in the form of a ratio of the likelihood of the evidence given the prosecution hypothesis ($H_p$) against that given the defence hypothesis ($H_d$) in which a random individual, X, replaces Q. To evaluate the WoE the user must supply a crime-stain profile (CSP), at least one reference profile (that of Q) and an allele frequency database for at least one population of interest. The reference profiles of any assumed contributors other than Q can also be provided. A number of allele frequency databases are supplied with \likeLTD, spanning the NGMSelect ({\tt DNA17}), SGM+ ({\tt SGMplus}) and Identifiler ({\tt Identifiler}) STR profiling kits, but the user may specify their own database. Similarly, the case shown here uses example CSP and reference profile files supplied with \likeLTD, but the user may specify their own.  

<<inputs, eval=FALSE, tidy=FALSE, print=FALSE>>=
  require(likeLTD)

  # Case we are going to be evaluating
  caseName = "Laboratory"
  datapath = file.path(system.file("extdata", package="likeLTD"),
      'laboratory')

  # File paths and case name for allele report
  admin = pack.admin.input.peaks(
      peaksFile = file.path(datapath, 'laboratory-CSP.csv'),
      refFile = file.path(datapath, 'laboratory-reference.csv'),
      caseName = caseName,
      detectionThresh = 20
      )
@
<<alleleReport, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Generate allele report
  allele.report.peaks(admin)
@

The important arguments to {\tt pack.admin.input.peaks} are:

\begin{description}
\item[{\tt peaksFile}:] Path to CSP file in the same format as {\tt laboratory-CSP.csv}. May contain multiple replicates.
\item[{\tt refFile}:] Path to reference file in the same format as {\tt laboratory-reference.csv}. May contain multiple reference profiles.
\item[{\tt databaseFile}:] Path to database file. Only used if {\tt kit=NULL}. The database file specified must contain a LUS column. 
\item[{\tt kit}:] One of {\tt "DNA17"}, {\tt "SGMplus"}, {\tt "Identifiler"} or {\tt NULL}. If both {\tt databaseFile} and {\tt kit} are {\tt NULL}, the DNA17 database is used as default. 
\item[{\tt linkageFile}:] File of pairwise recombination rates between loci. If {\tt NULL} the default supplied with \likeLTD\ is used. This is only used when Q and X are assumed related. 
\item[{\tt detectionThresh}:] The detection threshold used to analyse the electrophoresis results. This can either be a single value that covers all loci, or a named list with a separate detection threshold for each locus. Defaults to 20 RFU.
\end{description}

The allele report generated by \likeLTD\ gives suggestions on the the number of unprofiled contributors to include in the analysis, and on whether or not to model dropin, as well as summarising the data. Note the suggestions of the number of unprofiled contributors and dropin are based on designations of allelic, uncertain and non-allelic peaks according to the criteria given in Table \ref{tab:PHtoDiscrete}.

\begin{table}
\begin{center} 
\begin{tabular}{c|cc}
Designation & S & DS and OS \\
\hline
Non-alleleic & $x<0.05$ & $x<0.05$ \\
Uncertain & $0.05\le x<0.15$ & $0.05\le x<0.1$ \\
Allelic & $x\ge 0.15$ & $x\ge 0.1$ \\
\end{tabular}
\caption{Determination of which peaks in stutter (S), over-stutter (OS) or double-stutter (DS) positions are treated as non-allelic, uncertain or allelic when estimating {\tt nUnknowns}.  $x$ indicates the ratio of the stutter position peak height to the parent peak height.}
\label{tab:PHtoDiscrete}
\end{center}
\end{table}

If instead you wish to specify a separate detection threshold at each locus, usually for a different threshold at each epg lane, the arguments to {\tt pack.admin.input.peaks} should be:

<<detectThresh, eval=FALSE, tidy=FALSE, print=FALSE>>=
# Specifying locus specific detection thresholds
admin = pack.admin.input.peaks(
            peaksFile = file.path(datapath, 'laboratory-CSP.csv'),
            refFile = file.path(datapath, 'laboratory-reference.csv'),
            caseName = "Laboratory",
            detectionThresh = list(D10S1248=20,vWA=20,D16S539=20,D2S1338=20, # blue  
                                   D8S1179=30,D21S11=30,D18S51=30, # green
                                   D22S1045=40,D19S433=40,TH01=40,FGA=40, # black  
                                   D2S441=50,D3S1358=50,D1S1656=50,D12S391=50,SE33=50) # red
            )
@

\section{Hypothesis generation}

Based on the suggestions of the allele report, or through manual inspection, the user must decide on a set of arguments to pass to \likeLTD\ in order for it to generate the prosecution and defence hypotheses.

<<hyps, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Enter arguments
  args = list(
        nUnknowns = 1,
        doDropin = FALSE,
        ethnic = "NDU1",
        adj = 1,
        fst = 0.03,
        relationship = 0
        )

  # Create hypotheses
  hypP = do.call(prosecution.hypothesis.peaks, append(admin,args))
  hypD = do.call(defence.hypothesis.peaks, append(admin,args))
@

The arguments that may be handed to either {\tt prosecution.hypothesis.peaks} or {\tt defence.hypothesis.peaks} are:

\begin{description}
\item[{\tt nUnknowns}:] The number of unprofiled contributors to include under $H_p$. $H_d$ automatically adds an unprofiled contributor, X, that replaces Q. The allele report suggests a value which must be either 0, 1 or 2. Defaults to 0. 
\item[{\tt doDropin}:] Logical, whether or not to model dropin. Suggested by allele report. Defaults to {\tt FALSE}. 
\item[{\tt ethnic}:] Which database to use. Must match a column heading in the chosen database file. Defaults to {\tt "NDU1"} (Caucasian). Other populations included in the DNA17 database are {\tt "NDU1"} (Caucasian), {\tt "NDU2"} (African + Afro-Caribbean), {\tt "NDU3"} (South Asian), {\tt "NDU4"} (East Asian), {\tt "NDU6"} (African) and {\tt "NDU7"} (Afro-Caribbean).
\item[{\tt adj}:] Sampling adjustment. Defaults to 1. 
\item[{\tt fst}:] Fixation index, accounts for distant relatedness between Q and X. Defaults to 0.03.
\item[{\tt relationship}:] Assumed relationship between Q and X. Can take values between 0 and 7:
	\begin{enumerate}
	\item[0] Unrelated.
	\item[1] Parent/offspring.
	\item[2] Siblings.
	\item[3] Uncle (or aunt)/nephew (or niece).
	\item[4] Half-uncle (or half-aunt)/half-nephew (or half-niece).
	\item[5] Cousins.
	\item[6] Grandparent/grandchild.
	\item[7] Half-siblings.
	\end{enumerate}
The ordering of relationships does not alter the computation, so is unspecified e.g. if {\tt relationship=1}, the relationship may be Q as parent and X as offspring, or Q as offspring and X as parent.
\item[{\tt doDoubleStutter}:] Logical, whether or not to model double stutter. Defaults to {\tt TRUE}. 
\item[{\tt doOverStutter}:] Logical, whether or not to model over stutter. Defaults to {\tt TRUE}. 
\item[{\tt combineRare}:] Logical, whether or not to combine unobserved alleles into a single allele. Defaults to {\tt TRUE}. 
\item[{\tt rareThreshold}:] If {\tt combineRare=TRUE}, this parameter gives the population allele probability below which all unobserved alleles are combined. Defaults to 1 (all unobserved alleles). 
\end{description}


\section{Optimisation}

Once the hypotheses have been generated, the next step is to create the likelihood functions along with their associated parameters for optimisation, and then to perform the optimisation.

<<params, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Generate likelihood functions and optimisation parameters
  paramsP = optimisation.params.peaks(hypP,verbose=FALSE)
  paramsD = optimisation.params.peaks(hypD,verbose=FALSE)

  # reduce number of iterations for demonstration purposes
  paramsP$control$itermax=25
  paramsD$control$itermax=25
@

<<optimise, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Run optimisation
  results = evaluate.peaks(paramsP, paramsD, n.steps=1, 
      converge=FALSE)
@

Here {\tt n.steps} and {\tt converge} are specified to reduce run-time for demonstration purposes; if unspecified {\tt evaluate.peaks} determines the number of steps after the first step has completed. Other parameters the user may wish to pass to {\tt evaluate.peaks} include:

\begin{description}
\item[interim:] Logical, whether or not to generate interim reports after every step of optimisation. Defaults to TRUE. 
\item[seed.input:] Integer specification of a seed to use before optimisation. If unspecified, \likeLTD\ uses a numeric representation of the current time, date and process ID at the time of running {\tt evaluate}. 
\end{description}

The behaviour of the optimisation algorithm may be modified by specifying {\tt tolerance}, {\tt CR.start}, {\tt CR.end} and {\tt nConverged}, although we do not recommend altering these parameters unless you understand what effect each of the parameters has. 

The {\tt results} object contains the WoE at each step of optimisation ({\tt WoE}), the prosecution and defence $-\log_{10}$ likelihoods at each step of optimisation ({\tt Lp} and {\tt Ld}), the prosecution and defence objects returned by DEoptim after optimisation ({\tt Pros} and {\tt Def}) and information about the seed that was used ({\tt seed.used} and {\tt seed.input}).

\section{Generate output report}

Once optimisation has been completed, an output report summarising the results can be generated. This outputs the WoE against Q, likelihoods at each locus, the inverse match probability (IMP, the theoretical maximum WoE), DNA contribution estimates for each contributor in each replicate, degradation estimates for each contributor and dropin estimates if dropin is modelled. 

<<output, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Generate output report
  output.report.peaks(hypP,hypD,results)
@

<<optimisedResults, eval=FALSE, tidy=FALSE, print=FALSE,include=FALSE,echo=FALSE>>=
  results = list(Def=list(optim=list(bestmem = c(
    -3.459508761,     -2.599777111,     -3.385884129,    902.540872861, 
   171.908479331,    993.204056721,     49.207185386,      0.005254712, 
     1.016233455,      1.368054382,      0.984844585,      1.371468905, 
     0.535551035,      1.230762849,      1.358293511,      1.264302595, 
     0.747843015,      0.775721009,      1.000498903,      0.908954169, 
     0.822063630,      0.946761244,      0.985732220,      1.135499511, 
     0.001433327,      0.001742212))))
  names(results$Def$optim$bestmem) = c(
    "degradation1",     "degradation2",     "degradation3",         "DNAcont1", 
        "DNAcont2",         "DNAcont3",            "scale",        "gradientS", 
 "gradientAdjust1",  "gradientAdjust2",  "gradientAdjust3",  "gradientAdjust4", 
 "gradientAdjust5",  "gradientAdjust6",  "gradientAdjust7",  "gradientAdjust8", 
 "gradientAdjust9", "gradientAdjust10", "gradientAdjust11", "gradientAdjust12", 
"gradientAdjust13", "gradientAdjust14", "gradientAdjust15", "gradientAdjust16", 
           "meanD",            "meanO") 
@


\section{Likely genotypes of the unknown contributors}

Following optimisation it is possible to return the most likely marginal genotypes for each contributor or joint genotypes of all contributors under either $H_p$ or $H_d$. These can be useful for searching a national database with one of the unknown genotypes.

<<likely, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Get the most likely single-contributor genotypes
  gensMarginal = get.likely.genotypes.peaks(hypD,paramsD,
      results$Def)
  # Return joint genotypes and probabilities
  gensJoint = get.likely.genotypes.peaks(hypD,paramsD,
      results$Def,joint=TRUE)
@

It is also possible to return the probabilities for all of the genotype combinations under the specified hypothesis. This has been used to compute the WoE against a degraded incomplete reference profile, drawing incomplete genotypes from the posterior probability for which any known alleles match, but may have other uses.

<<posterior, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Get the posterior likelihoods for all genotype combinations
  gensPosterior = get.likely.genotypes.peaks(hypD,paramsD,
      results$Def,posterior=TRUE)
@

<<run,echo=FALSE,include=FALSE>>=
<<inputs>>
<<hyps>>
<<params>>
<<optimisedResults>>
@

\section{Diagnostic}

A visual inspection of the fit of the optimised parameters to the supplied CSP data can be performed using the {\tt peaks.results.plot} function, which plots the 95\% equal-tailed probability interval of the gamma distribution given the optimised parameters, assuming the joint genotype combination that is most likely (see Figure \ref{fig:diagnose}). 

%, fig.cap="95 CI of the gamma distribution given optimised parameters (black boxplots), with observed peak heights in the CSP (red bars). Y-axis displays RFU, x-axis displays allele."

<<diagnose, eval=FALSE, tidy=FALSE, print=FALSE>>=
  # Plot CSP with most likely genotypes
  peaks.results.plot(hypD,results$Def,replicate=1)
@

\begin{figure}[!h]
\begin{center}
<<plot, fig=TRUE,echo=FALSE>>=
  # Plot CSP with most likely genotypes
  peaks.results.plot(hypD,results$Def,replicate=1)
@
\caption{Boxes show the central 50\% (inter-quartile range) of the gamma distribution for each hypothesised peak, whiskers represent the 95\% equal-tailed probability interval and red bars show observed peak heights. RFU is displayed on the y-axes while allele labels corresponding to boxplots are displayed on the x-axes.}
\label{fig:diagnose}
\end{center}
\end{figure}

\end{document}
