R/plot.pf.R

In numberofalleles: Compute the Probability Distribution of the Number of Alleles in a DNA Mixture

#' Plot method for objects of class pf
#'
#' @param x An object of class pf
#' @param lines if \code{TRUE} then lines will be drawn from 0 to the
#' probability value for each x value.
#' @param points if \code{TRUE} then points will be plotted at the
#' probability value for each x value.
#' @param line_col the colour of the lines drawn for each probability mass
#' @param point_col the colour of the points plotted for each probability mass
#' @param xlab a label for the x axis, defaults to \code{"Number of alleles (n)"}
#' @param ylab a label for the y axis, defaults to \code{expression(Pr(N == n))}
#' @param ... Any other arguments that should be sent to \code{plot},
#' \code{arrows}, or \code{points}.
#' @param add If \code{TRUE} then the plotting information will be added to
#' the existing plot.
#' @returns No return value. Has the side effect of plotting to the active graphics device.
#' @examples
#'
#' ## Load allele frequencies
#' freqs <- read_allele_freqs(system.file("extdata","FBI_extended_Cauc.csv",
#'                             package = "numberofalleles"))
#'
#' gf_loci <- c("D3S1358", "vWA", "D16S539", "CSF1PO", "TPOX", "D8S1179",
#'             "D21S11",  "D18S51", "D2S441", "D19S433", "TH01", "FGA",
#'             "D22S1045", "D5S818", "D13S317", "D7S820", "SE33",
#'             "D10S1248", "D1S1656", "D12S391", "D2S1338")
#'
#' gf_freqs <- freqs[gf_loci]
#'
#' ## Compute the pedigrees for two and three siblings respectively
#' pedSibs2 = pedtools::nuclearPed(father = "F", mother = "M",
#'                                 children = c("S1", "S2"))
#' pedSibs3 = pedtools::addChildren(father = "F", mother = "M",
#'                                  pedSibs2, ids = "S3")
#'
#' ## Compute the probability functions for each of 2 Unknowns, and 2 Sibs
#' p2U = pr_total_number_of_distinct_alleles(contributors = c("U1", "U2"),
#'                                           freqs = gf_freqs)
#' p2S =  pr_total_number_of_distinct_alleles(contributors = c("S1", "S2"),
#'                                            freqs = gf_freqs,
#'                                            pedigree = pedSibs2)
#'
#' ## And plot the two probability distribution functions on top of each other
#' plot(p2U,
#'      line_col = "red",
#'      point_col = "red",
#'      pch = 18,
#'      lwd = 2,
#'      ylim = c(0, 0.15),
#'      xlab = "Total Number of Alleles (TAC)",
#'      ylab = expression(Pr(N == n~"|"~ P)))
#'
#' plot(p2S,
#'      add = TRUE,
#'      point_col = "blue",
#'      line_col = "blue",
#'      lwd = 2)
#'
#' legend("topleft", legend = c("2 U", "2 S"), fill = c("red", "blue"), bty = "n")
#'
#' ## Compute the LR for the number of peaks given the
#' #  number of contributors and the pedigrees
#' lr = p2U$pf / p2S$pf
#' data.df = data.frame(log10lr = log10(lr), noa = p2U$noa)
#'
#' ## Plot the LR and a grid so that it's easy to see where
#' #  the LR becomes greater than 1
#' plot(log10lr~noa,
#'      data = data.df,
#'      axes = FALSE,
#'      xlab = "Total number of alleles, n",
#'      ylab = expression(log[10](LR(P[1],P[2]~"|"~N==n))),
#'      xaxs = "i", yaxs = "i",
#'      xlim = c(22,93),
#'      pch = 18)
#' axis(1)
#' y.exponents = seq(-20, 15, by = 5)
#' for(i in 1:length(y.exponents)){
#'   if(y.exponents[i] == 0)
#'     axis(2, at=y.exponents[i], labels="1", las = 1)
#'   else if(y.exponents[i] == 1)
#'     axis(2, at=y.exponents[i], labels="10", las = 1)
#'   else
#'     axis(2, at = y.exponents[i],
#'          labels = eval(substitute(
#'            expression(10^y),
#'            list(y = y.exponents[i]))),
#'          las = 1)
#' }
#' grid()
#' box()
#'
#' ## Let's look at 2 sibs versus 3 sibs
#'
#' p3S = pr_total_number_of_distinct_alleles(contributors = c("S1", "S2", "S3"),
#'                                           freqs = gf_freqs,
#'                                           pedigree = pedSibs3)
#' plot(p3S,
#'      line_col = "green",
#'      point_col = "green",
#'      pch = 18,
#'      lwd = 2,
#'      ylim = c(0, 0.15),
#'      xlab = "Total Number of Alleles (TAC)",
#'      ylab = expression(Pr(N == n~"|"~ P)))
#'
#' plot(p2S,
#'      add = TRUE,
#'      pch = 18,
#'      point_col = "blue",
#'      line_col = "blue",
#'      lwd = 2)
#' legend("topleft", legend = c("3 S", "2 S"), fill = c("green", "blue"), bty = "n")
#'
#' ## And finally two sibs and one unknown versus three sibs
#' p2SU = pr_total_number_of_distinct_alleles(contributors = c("S1", "S2", "U1"),
#'                                            freqs = gf_freqs,
#'                                            pedigree = pedSibs2)
#' plot(p3S,
#'      line_col = "green",
#'      point_col = "green",
#'      pch = 18,
#'      lwd = 2,
#'      ylim = c(0, 0.15),
#'      xlab = "Total Number of Alleles (TAC)",
#'      ylab = expression(Pr(N == n~"|"~ P)))
#' plot(p2SU,
#'      add = TRUE,
#'      pch = 18,
#'      point_col = "orange",
#'      line_col = "orange",
#'      lwd = 2)
#' legend("topleft", legend = c("3 S", "2 S + U"),
#'        fill = c("green", "orange"), bty = "n")
#' @export
plot.pf = function(x,
                   lines = TRUE,
                   points = TRUE,
                   line_col = "black",
                   point_col = "black",
                   xlab = "Number of alleles (n)",
                   ylab = expression(Pr(N == n)),
                   add = FALSE,
                   ...){
  xVals = x$noa
  yVals = x$pf

  if(!add)
    plot(xVals, yVals, type = "n", xlab = xlab, ylab = ylab, ...)

  if(points)
    points(xVals, yVals,
           col = point_col,
           ...)

  if(lines)
    graphics::arrows(x0 = xVals,
                     x1 = xVals,
                     y0 = rep(0, length(yVals)),
                     y1 = yVals,
                     length = 0, ## No head
                     col = line_col,
                     ...
    )
}