plot.pf {numberofalleles}

Plot method for objects of class pf

Description

Plot method for objects of class pf

Usage

## S3 method for class 'pf'
plot(
  x,
  lines = TRUE,
  points = TRUE,
  line_col = "black",
  point_col = "black",
  xlab = "Number of alleles (n)",
  ylab = expression(Pr(N == n)),
  add = FALSE,
  ...
)

Arguments

x

An object of class pf

lines

if TRUE then lines will be drawn from 0 to the probability value for each x value.

points

if TRUE then points will be plotted at the probability value for each x value.

line_col

the colour of the lines drawn for each probability mass

point_col

the colour of the points plotted for each probability mass

xlab

a label for the x axis, defaults to "Number of alleles (n)"

ylab

a label for the y axis, defaults to expression(Pr(N == n))

add

If TRUE then the plotting information will be added to the existing plot.

...

Any other arguments that should be sent to plot, arrows, or points.

Value

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")