Nothing
## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
library(sequoia)
knitr::opts_chunk$set(echo = TRUE,
fig.height=4, fig.width=6,
fig.pos="ht!")
## ----DeerAP, echo=FALSE, fig.cap="Example of 'ageprior' distribution in a red deer population", out.width="90%"----
knitr::include_graphics("Deer-AgePriorExtra.png")
## ----pipeline, echo=FALSE, fig.cap="Pipeline overview", out.width="90%"-------
knitr::include_graphics("agepriors-pipeline.png")
## ----Glossary, echo=FALSE, results="asis"-------------------------------------
Glossary <- matrix(c("$BY_i$", "Birth year", "Time unit (year, decade, month) of $i$'s birth/ hatching",
"$A_{i,j}$", "Age difference", "$BY_i - BY_j$, i.e. if $j$ is older than $i$, than $A_{i,j}$ is positive",
"$R_{i,j}$", "Relationship", "Relationship between $i$ and $j$, e.g. parent-offspring",
"$\\alpha_{A,R}$", "Ageprior", "Probability ratio of relationship $R$ given $A$, versus $R$ for a random pair"),
ncol=3, byrow=TRUE)
knitr::kable(Glossary, col.names = c("Symbol", "Term", "Definition"),
caption="Terms and abbreviations", escape=FALSE, booktabs=TRUE)
## ----AP-grif1, out.width="70%"------------------------------------------------
MakeAgePrior(Pedigree = Ped_griffin, Smooth=FALSE, Flatten=FALSE) # Note: Ped_griffin includes a birth year column
## ----RelAbbr, echo=FALSE------------------------------------------------------
as.DF <- function(V, cn=c("name", "x")) { # named vector --> data.frame
setNames(data.frame(c1 = names(V), c2=V), cn)
}
Rels <- c(M = "Mother - offspring",
P = "Father - offspring",
FS = "Full siblings", # From version XX
MS = "Maternal siblings (full + half)",
PS = "Paternal siblings (full + half)",
MGM = "Maternal grandmother",
PGF = "Paternal grandfather",
MGF = "Maternal grandfather (+ paternal grandmother)",
UA = "avuncular (niece/nephew -- aunt/uncle, full + half)",
"(M/P)(M/P/F)A" = "Avuncular; Mother's/Father's Maternal/Paternal/Full sibling")
Rels <- cbind(as.DF(Rels, c("Column (R)", "Meaning")),
"Version" = c(rep("all", 2), "from 1.0", "all", "all",
rep("up to 1.3 *", 3), "up to 1.3",
"from 2.0 *"))
knitr::kable(Rels, row.names=FALSE, booktabs=TRUE,
caption = "AgePriors column names per sequoia version. *: From sequoia v2.0 only in AgePriorExtra")
## ----AP-grif2-----------------------------------------------------------------
# sequoia() output for the griffin data is included in the package:
# see ?SeqOUT_griffin on how it was generated
round(SeqOUT_griffin$AgePriorExtra, 2)
PlotAgePrior(SeqOUT_griffin$AgePriorExtra)
## ----AP-default, fig.cap="Default minimally informative ageprior", out.width = "70%"----
MakeAgePrior(LifeHistData = LH_HSg5)
## ----AP-default-2-------------------------------------------------------------
table(LH_HSg5$BirthYear)
## ----APdiscrete, fig.cap="Discrete generations"-------------------------------
# running sequoia to also get the extra ageprior columns:
Seq_HSg5 <- sequoia(SimGeno_example, LH_HSg5, Module="par",
args.AP=list(Discrete = TRUE),
CalcLLR=FALSE, Plot=FALSE, quiet=TRUE)
PlotAgePrior(Seq_HSg5$AgePriorExtra)
## ----APDiscrete2--------------------------------------------------------------
MakeAgePrior(Discrete = TRUE, MaxAgeParent = 2, Plot=FALSE)
## ----pedHead------------------------------------------------------------------
tail(Ped_griffin) # the pedigree
## ----APall, fig.cap="Ageprior for overlapping generations (griffin example)", fig.pos="h!", out.width = "70%"----
AP.griffin <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=FALSE,
Return="all")
names(AP.griffin)
## ----BYrange------------------------------------------------------------------
AP.griffin$BirthYearRange
## ----MAP----------------------------------------------------------------------
AP.griffin$MaxAgeParent
## ----IncrMAP, out.width = "70%"-----------------------------------------------
# Specify maximum age for fathers, but estimate max mother age from pedigree
APX <- MakeAgePrior(Ped_griffin, Smooth=FALSE, MaxAgeParent = c(NA, 5))
## ----tblAR--------------------------------------------------------------------
AP.griffin[["tblA.R"]]
## ----tblAR-2------------------------------------------------------------------
AgeDifM <- outer(Ped_griffin$birthyear, Ped_griffin$birthyear, "-")
# closest relationship for each pair:
RelM <- GetRelM(Ped_griffin, patmat=TRUE, GenBack=1)
table(c(AgeDifM), c(RelM)) # vectorising first speeds up table() a lot
## ----tblAR-2a-----------------------------------------------------------------
# - 'MaxT' ensures sufficient, constant number of rows across intermediary tables
MaxAgeParent <- c(NA, NA)
MaxAgePO <- ifelse(!is.na(MaxAgeParent), MaxAgeParent, diff(AP.griffin$BirthYearRange)+1) # in actual code with thorough input checks
MaxT <- max(MaxAgePO+1, diff(AP.griffin$BirthYearRange))
# - factor() ensures levels with 0 count are included
# - drop negative age differences: all pairs are included 2x
tblA.R <- table(factor(AgeDifM, levels = 0:MaxT),
factor(RelM, levels = c("M", "P", "FS", "MHS", "PHS")))
## ----tblAR-2b-----------------------------------------------------------------
RelA <- GetRelM(Ped_griffin, patmat=TRUE, GenBack=1, Return="Array")
tblA.R <- sapply(c("M", "P", "FS", "MHS", "PHS"),
function(r) table(factor(AgeDifM[RelA[,,r]==1],
levels = 0:MaxT)))
head(tblA.R)
## ----tblAR-2c-----------------------------------------------------------------
tblA.R <- cbind(tblA.R,
"MS" = tblA.R[,"MHS"] + tblA.R[,"FS"],
"PS" = tblA.R[,"PHS"] + tblA.R[,"FS"])
tblA.R <- tblA.R[, c("M", "P", "FS", "MS", "PS")] # drop MHS & PHS columns
# add column for age diff. across all relationships (except 'Self' pairs)
tblA.R <- cbind(tblA.R,
"X" = table(factor(AgeDifM[c(RelM)!="S"], levels = 0:MaxT)))
# All pairs are included 2x in AgeDifM/RelM/RelA, to fix this:
# - drop the negative age differences (done with factor() above)
# - divide counts for row A=0 by 2 (including 'X' column)
tblA.R["0", ] <- tblA.R["0", ] / 2
# check that this is really the output from MakeAgePrior():
all(tblA.R == AP.griffin[["tblA.R"]])
## ----calcPAR------------------------------------------------------------------
NAK.R <- colSums(tblA.R, na.rm=TRUE) # total count per relationship
NAK.R
PA.R <- sweep(tblA.R, 2, STATS=NAK.R, FUN="/")
head(round(PA.R, 2))
## ----FSuseHS------------------------------------------------------------------
MinPairs.AgeKnown <- 20 # not user-setable.
# Use the product of the distributions if there are many MS & PS,
# and the average across the two distributions if there are a medium number:
if (NAK.R["FS"] / min(NAK.R[c("MS", "PS")]) < 0.5 &
all(NAK.R[c("MS", "PS")] > MinPairs.AgeKnown)) {
if (all(NAK.R[c("MS", "PS")] > 5*MinPairs.AgeKnown)) {
FS.tmp <- PA.R[,"MS"] * PA.R[,"PS"]
FS.tmp <- FS.tmp/sum(FS.tmp)
} else {
FS.tmp <- apply(PA.R[, c("MS", "PS")], 1, mean)
}
PA.R[,"FS"] <- (PA.R[,"FS"] + FS.tmp)/2
}
round(PA.R[1:5, ], 2) # note change in 'FS' column
## ----calcLRAR-----------------------------------------------------------------
LR.RU.A <- sweep(PA.R[, 1:5], 1, STATS=PA.R[,"X"], FUN="/")
head(round(LR.RU.A, 2))
# check that this is identical to output from MakeAgePrior():
LR.RU.A[!is.finite(LR.RU.A)] <- 0
all(abs(LR.RU.A - AP.griffin[["LR.RU.A.raw"]]) < 0.001) # allow for rounding errors
## ----griffinSmoothFlatten, echo=FALSE, fig.cap="Smoothed & Flattened agepriors for overlapping generations (griffin example). Note the different y-axes; non-displayed age differences have an ageprior of zero for all R.", fig.height=8, fig.width=8, message=FALSE----
AP_list <- list()
AP_list[["Original"]] <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=FALSE, Plot=FALSE)
AP_list[["Smooth"]] <- MakeAgePrior(Ped_griffin, Smooth=TRUE, Flatten=FALSE, Plot=FALSE)
AP_list[["Flatten"]] <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=TRUE, Plot=FALSE)
AP_list[["Flatten + Smooth"]] <- MakeAgePrior(Ped_griffin, Smooth=TRUE, Flatten=TRUE, Plot=FALSE)
par(mfcol=c(2,2), mai=c(.7,.8,.4,.1))
for (x in names(AP_list)) {
PlotAgePrior(AP_list[[x]], legend=FALSE)
mtext(x, side=3, line=0.5, cex=1.5)
}
## ----tailstretch--------------------------------------------------------------
SmoothAP <- function(V, MinP) {
Front <- max(1, min(which(V > MinP)), na.rm=TRUE)
End <- min(max(which(V > MinP)), length(V), na.rm=TRUE)
if (Front > 1 & V[Front] > 3*MinP) V[Front -1] <- MinP
if (End < length(V)) V[End +1] <- V[End]/2
if ((End+1) < length(V) & V[End+1] > 3*MinP) V[End+2] <- MinP
# ... dip-fixing ...
V
}
LR.RU.A <- apply(LR.RU.A, 2, SmoothAP, MinP = 0.001)
## ----PedComp-ex1, echo=TRUE, eval=TRUE----------------------------------------
AP.pedigree <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=FALSE,
quiet=TRUE, Plot=FALSE)
AP.default <- MakeAgePrior(MaxAgeParent = 3, quiet=TRUE, Plot=FALSE)
knitr::kable(list(AP.pedigree, AP.default),
caption = "Pedigree-based ageprior (left) and default, flat ageprior (right).",
booktabs=TRUE)
## ----PedComp-ex1b, echo=TRUE, eval=TRUE---------------------------------------
AP.flattened <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=TRUE, quiet=TRUE, Plot=FALSE)
## ----PedComp-ex1c, echo=FALSE, eval=TRUE--------------------------------------
if (is_html_output()) {
knitr::kable(AP.flattened, caption = "Flattened ageprior", format="html",
table.attr = "style='width:50%;'")
} else {
knitr::kable(AP.flattened, caption = "Flattened ageprior", booktabs=TRUE)
}
## ----plotFlattenNW, echo=FALSE, fig.cap="Weights versus number of pairs with relationship R, for weighed average between pedigree-derived ageprior and flat 0/1 ageprior", fig.pos="hbt", fig.height=5, , out.width="60%"----
N <- 0:310
lambdaNW <- c(1e-4, -log(0.5) / c(50,100,200), 0.2)
W <- matrix(NA, 5, length(N))
for (i in 1:5) {
W[i,] <- 1 - exp(-lambdaNW[i] * N)
}
par(mai=c(.8,.8,.1,.1))
plot(N, W[1,], type="n", ylim=c(0,1), las=1,
xlab="Number of pairs (N_R)", ylab = "Weight (W_R)")
abline(h=seq(0,1,.2), col="lightgrey")
abline(h=0.5, col="darkgrey")
segments(x0=c(100,50,200), y0=-1, y1=0.5, col="darkgrey",
lty=1:3, lwd=c(1,1,1.5))
axis(side=2, at=0.5, labels=0.5, las=1)
for (i in 1:3) {
lines(N, W[i+1,], lwd=3, col=c(4,1,4)[i], lty=c(2,1,3)[i])
}
for (i in c(1,5)) {
lines(N, W[i,], lwd=2, col="darkgreen", lty=ifelse(i==1, 3, 2))
}
legend(210,0.45, legend=c(0.2, paste("-log(0.5)/", c(50,100,200)), 0.0001),
col=c("darkgreen", 4,1,4, "darkgreen"), lty=c(2,2,1,3,3), lwd=c(2,3,3,3,2),
title="lambdaNW", bg="white", inset=.03)
## ----griffinweights-----------------------------------------------------------
AP.griffin$Weights
## ----griffinweights2----------------------------------------------------------
NAK.R # counts per relationship (column sums of tblA.R)
lambdaNW <- -log(0.5)/100 # default input value
W.R <- 1 - exp(-lambdaNW * NAK.R[1:5])
round(W.R, 4)
## ----timelineAU, echo=FALSE, fig.cap="Age difference between individual i and aunt/uncle j, and possible 'birth years' for i's parent k.", out.width="60%"----
knitr::include_graphics("agedif-AU.png")
## ----skip, fig.cap="", out.width = "70%"--------------------------------------
AP <- MakeAgePrior(MaxAgeParent = c(12, 8), Plot=FALSE)
AP[as.character(seq(1,12,by=2)), c("MS", "FS")] <- 0
PlotAgePrior(AP)
# check if valid ageprior:
chk <- sequoia:::CheckAP(AP)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.