---
title: "Sequoia -- age effects"
author: "Jisca Huisman ( jisca.huisman @ gmail.com )"
date: "22 April 2020"
output:
  bookdown::html_document2:
  bookdown::pdf_book:
    extra_dependencies:
    - amssymb
    - amsmath
    - graphicx
    - Xcolor
linkcolor: blue
config:
  toc: yes
  toc_depth: 3
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Sequoia -- age effects} 
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
# output: 
#   bookdown::pdf_book: 
#   pdf_document:
library(knitr)
library(sequoia)
knitr::opts_chunk$set(echo = TRUE, 
                      fig.path='figs/', fig.height=4, fig.width=6,
                      fig.pos="h!")

# superseded: \newcommand{\grey}[1]{\textcolor{gray}{#1}}  (in main text, pdf only)

# adapted from: https://bookdown.org/yihui/rmarkdown-cookbook/font-color.html
MkGrey <- function(x) {    
  if (knitr::is_latex_output()) {
    sprintf("\\textcolor{gray}{%s}", x)
  } else if (knitr::is_html_output()) {
    sprintf("<span style='color: gray;'>%s</span>", x)
  } else x
}
```

# In short
The relative age of individuals can be very useful in pedigree reconstruction, both to filter potential relatives and during assignment. Age difference based probability ratios for different relationships are calculated based on a skeleton parentage-only pedigree and the birth years provided in `LifeHistData` (Figure \@ref(fig:pipeline)). 

This 'ageprior' (Glossary in Table \@ref(tab:Glossary); please note that this is not an official term, nor a proper prior) was calculated quietly in the background up until \verb+sequoia+ version 1.3.3, but from version 2.0 you will see heatmaps and messages such as

\begin{verbatim}
	Ageprior: Default 0/1, overlapping generations, MaxAgeParent=5,5
\end{verbatim}
\begin{verbatim}
	Ageprior: Pedigree-based, discrete generations, MaxAgeParent = 1,1
\end{verbatim}
\begin{verbatim}
	Ageprior: Pedigree-based, overlapping generations, smoothed, MaxAgeParent = 20,16
\end{verbatim} 

These are *not* warnings, but inform you what is going on what is going on behind the scenes. This will hopefully minimise unintended age priors as a source of error. There are four parts to the message:
\begin{itemize} \itemsep0em
  \item `Default 0/1` versus `Pedigree-based` indicates whether the ageprior at that point is simply the default (see further) or based on the skeleton parentage-only pedigree;
	\item `discrete` or `overlapping` generations
	\item `smoothed` and/or `flattened` indicate small-sample corrections for the pedigree-based age prior, see further (if it isn't mentioned, it wasn't applied)
	\item Maximum age of mothers, fathers
\end{itemize}

In addition `MakeAgePrior` allows more explicit control from version 2.0 onwards, so that 'hand-tailoring' will only be required in rare cases. Arguments can be passed to `MakeAgePrior` when running `sequoia` via `args.AP`. 


```{r pipeline, echo=FALSE, fig.cap="Pipeline overview", out.width="90%"} 
if (knitr::is_latex_output()) {
  knitr::include_graphics("figs/agepriors_pipeline.pdf")
} else if (knitr::is_html_output()) {
  knitr::include_graphics("figs/agepriors_pipeline.png")
}
```


The ageprior matrix now has 5 columns (more in sequoia versions before 2.0), corresponding to parent-offspring and siblings (Table \@ref(tab:RelAbbr)). The number of rows equals the maximum age of parents plus one; the first row is for individuals born in the same year (age difference $A=0$), the second row for individuals born one year apart, etc. Some values are fixed to zero: parents and their offspring can never be exactly the same age.


## Why \& how: use of age in pedigree reconstruction
Based on the species' age at first and last reproduction, some age differences between parent and offspring or between siblings are more likely than others, and some are downright impossible. Age differences are calculated from the birth years (or other time unit) in `LifeHistData`, and are used both to filter potential relatives (and reduce the number of time-consuming computations) and during assignment, for example to distinguish between half-siblings, grandparent--grand-offspring and full avuncular pairs. 

By default, for the parentage assignment only a 0/1 coded 'flat' ageprior is used to indicate which age-differences are possible. If any parents with known age are assigned, the age prior is updated before use during sibship clustering and grandparent assignment (Figure \@ref(fig:pipeline)). The ageprior is not updated during or after the sibship clustering. 

If an ageprior is provided as argument to `sequoia` (as an element of `SeqList`), this ageprior will be used instead of calling `MakeAgePrior` (details further down).


```{r 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$ across all age differences"),
       ncol=3, byrow=TRUE)

knitr::kable(Glossary, col.names = c("Symbol", "Term", "Definition"),
            caption="Terms and abbreviations", escape=FALSE, booktabs=TRUE)
```


## Interpretation
The 'AgePrior' matrix is a set of age-difference based probability ratios which can be interpreted as

\begin{quotation}\slshape
If I were to pick two individuals with an age difference $A$, and two individuals at random, how much more likely are the first pair to have relationship $R$, compared to the second pair?
\end{quotation}

For example, in a species where females lay large clutches but only breed for a few years, two individuals $i$ and $j$ born in the same year ($BY_i=BY_j$, $A_{i,j}=0$) may be more likely to be maternal siblings ($R_{i,j} =$MS, relationships in Table \@ref(tab:RelAbbr)) than two randomly picked individuals  ($\alpha_{0,MS} > 1$). In contrast, in a species where singletons are the norm and twins are rare, individuals born in the same year are unlikely to be maternal siblings ($\alpha_{0,MS} < 1$). 

# Maths

The ageprior is defined as 
\begin{equation}
	\alpha_{A,R} = \frac{P(R|A)}{P(R)} \text{ .}
	\label{eq:AP}
\end{equation}

Using Bayes' rule that  
\begin{equation}
	P(R|A) = \frac{P(A|R) P(R)}{P(A)}
\end{equation}

it follows that
\begin{equation}
	\alpha_{A,R} = \frac{P(R|A)}{P(R)} = \frac{P(A|R)}{P(A)} \text{ ,}
	\label{eq:AP2}
\end{equation}
When a pedigree and birth years are provided, $P(A)$ and $P(A|R)$ are simply calculated as the observed proportions of individual pairs with age-difference $A$ (and relationship $R$) (for an example, see further). 

When using the ageprior for relationships with and between dummy individuals ('stand-ins' for unsampled parents of sibship clusters), the assumption is made that the effect of age on relationship probability doesn't differ (much) between sampled and unsampled individuals:

\begin{equation}
	\frac{P(R|A, \text{sampled})}{P(R| \text{sampled})} \approx \frac{P(R|A | \text{unsampled})}{P(R | \text{unsampled})} 
	\label{eq:APsampled}
\end{equation}

This does \emph{not} mean that sampling probability is assumed to be independent of $A$: any sampling period is finite, thus there are always fewer pairs with large age differences than with small age differences. It does assume that the probability to be sampled is independent of relationship $R$, so that 

\begin{equation}
	\frac{P(A|R, \text{sampled})}{P(A| \text{sampled})} \approx
	\frac{P(A|R, \text{unsampled})}{P(A| \text{unsampled})}  \text{.}
\end{equation}


During pedigree reconstruction, this ageprior can simply be multiplied by the genotypes-dependend probability,

\begin{equation}
	P(R_{i,j}| G_i, G_j, A_{i,j}) = P(R_{i,j} | G_i, G_j) \times P(R_{i,j} | A_{i,j})  \text{ ,}
\end{equation}
as it can be reasonably assumed that conditional on the relationship $R$, the genotypes $G$ and age difference $A$ are independent (at least on the time scales relevant for pedigree reconstruction). The calculation of $P(R_{i,j} | G_i, G_j)$ is described in detail in (Huisman, 2017). 

 

## When \& why small sample correction is needed
The main challenge is when there are very few or no sampled pairs with age difference $A$, while that age difference is biologically plausible for one or more relationships in Table \@ref(tab:RelAbbr). In other words, if $P(A|\text{sampled}) = 0$ in equation \@ref(eq:AP2), this leads to division by zero and undefined $\alpha_{A,R}$, while we know that $P(R|A) >0$. 

Or, as another extreme example, when no candidate fathers were sampled ($P(A | R=P, \text{sampled}) = 0$) then $\alpha_{A=0,R=P} =0$. The dummy fathers that will stand in for these unsampled fathers will only be assigned when any $\alpha_{A,R=P} >0$. Consequently, when simply using the 'raw' ratio from the observed data, no dummy fathers will be assigned. 

To mitigate these issues, the \verb+MakeAgePrior+ function uses two types of 
small sample correction, adjustable with the arguments 'Flatten' and 'Smooth'. See section [Small sample correction] for details. 
 
```{r 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" = "Mother's/Father's Maternal/Paternal/Full sibling")

Rels <- cbind(as.DF(Rels, c("Column", "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. *: From sequoia v2.0 only in AgePriorExtra")
``` 



# Implementation
The implementation of `MakeAgePrior` differs considerably depending on whether or not a pedigree is provided as input, and whether discrete or overlapping generations are assumed. The details are described for each of these situations, and a worked example is given for a pedigree with overlapping generations. 


## No pedigree, discrete generations
When generations do not overlap, it is often most straightforward to use generation number instead of birth year, such that parents and offspring are born 1 time unit apart. We can initialise this ageprior as follows:
```{r APDiscrete, fig.cap="Ageprior for discrete generations", out.width="70%"}
MakeAgePrior(Discrete = TRUE)
```

and '1' (pale green in heatmap) indicates that a age-relationship combination is possible, and 0 (black) the impossible combinations.

It is possible to use discrete generations with a larger time interval. For example, pink salmon have a strict two-year life-cycle with odd-year and even-year populations:
```{r APDiscrete2}
MakeAgePrior(Discrete = TRUE, MaxAgeParent = 2, Plot=FALSE)
```


## No pedigree, (potentially) overlapping generations
When generations (may) overlap, and there is no pedigree yet, you must either provide an estimate of the maximum age of parents, or `LifeHistData` with a representative set of birth years. 

For example, if females breed to maximum age 3, and males to maximum age 2,
```{r APMaxPO, fig.cap="Ageprior for overlapping generations, MaxAgeParent is 3 for dams and 2 for sires.", out.width="70%"}
MakeAgePrior(MaxAgeParent = c(3, 2))
```

When `MaxAgeParent` is not provided, it is set to the maximum age difference observed in `LifeHistData`. 

Note that when assuming overlapping generations, aunts and uncles may be born in the same year as their nieces/nephews ($\alpha_{A=0,R=UA} >0$). Considering these relationship alteratives (maternal/paternal full/half avuncular) may lower the assignment rate of siblings compared to a situation where discrete generations are assumed -- both correct and incorrect ones. 


## With pedigree, discrete generations
When generations do not overlap the ageprior as estimated from a pedigree will look indistinguishable from the ageprior defined without any pedigree. 

When a pedigree is provided and \verb+Discrete=TRUE+, a check is done to confirm that all parent-offspring pairs do indeed have an age difference of 1 (typically), and that all sibling pairs have an age difference of 0 (strictly); any violations result in an error. When \verb+Discrete=NULL+ (the default), and it is found that these conditions are met, \verb+Discrete+ is set to TRUE internally, and thereby \verb+Smooth+ and \verb+Flatten+ are set to FALSE, but only if there are at least 20 pairs each of mother-offspring, father-offspring, maternal siblings, and paternal siblings (full siblings count here as both maternal and as paternal siblings). 


## With pedigree, overlapping generations
When a large number of parents with known age is assigned, the age difference distribution per relationship can simply be calculated from a table of counts. When the number is small, a weighed average is taken between these observed proportions and the default flat 0/1 prior (see section [Flatten]). Even with moderate sample sizes bumps and gaps can exist in the empirical distribution purely by chance, which are 'Smoothened' out by default (see section [Smooth]). 

The process is illustrated using an imaginary population  of griffins were each year exactly 20 baby griffins are born. From 2001 to 2010, all juveniles are sampled, and parents assigned.

```{r pedHead}
data(Ped_griffin)
tail(Ped_griffin)
```

We can calculate the ageprior matrix, and return all the intermediate steps as well:
```{r APall, fig.cap="Ageprior for overlapping generations (griffin example)", fig.pos="h!"}
AP.griffin <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=FALSE, Return="all")
```

Were we see that both dams and sires have a maximum age of 3 years, and that for dams (M: Maternal) age 1 is most common (darkest green), while for sires (P:Paternal) age 2 is most common. 
```{r APnames}
names(AP.griffin)
```


### Counting the number of pairs per age \& relationship
The list element `tblA.R` contains the raw counts of each age $A$ per relationship $R$:
```{r tblAR}
AP.griffin[["tblA.R"]]
```

Maternal ages in this table are calculated as follows:
```{r matAge}
for (x in c("id", "dam", "sire")) {
  Ped_griffin[,x] <- as.character(Ped_griffin[,x])
}
Ped.g2 <- merge(Ped_griffin, 
                setNames(Ped_griffin[, c("id", "birthyear")], c("dam", "BY.dam")),
                all.x=TRUE)
Ped.g2$Age.dam <- with(Ped.g2, birthyear - BY.dam)
table(Ped.g2$Age.dam)
```
and analogous for paternal age. 

Siblings are identified using function `GetRelCat`, and all pairwise age differences are calculated using `outer`:
```{r sibAge}
# this is part of the code in 'MakeAgePrior':
Ped.R <- setNames(Ped_griffin, c("id", "dam", "sire", "BirthYear"))  
RCM <- sapply(seq_along(Ped.R$id), GetRelCat, Ped.R, GenBack=1)  
# NOTE: GetRelCat is slow for very large pedigrees
AgeDifM <- outer(Ped.R$BirthYear, Ped.R$BirthYear, "-")
diag(AgeDifM) <- NA

MaxT <- 9   # determined internally based on LifeHistData
RR <- c("M", "P", "FS", "MS", "PS")
tblA.R <- matrix(NA, MaxT+1, length(RR)+1, dimnames=list(0:MaxT, c(RR, "X")))

tblA.R[, "FS"] <- table(factor(abs(AgeDifM[RCM == "FS"]), levels=0:MaxT)) /2
tblA.R[, "MS"] <- table(factor(abs(AgeDifM[RCM %in% c("FS", "MHS")]), levels=0:MaxT)) /2
tblA.R[, "PS"] <- table(factor(abs(AgeDifM[RCM %in% c("FS", "PHS")]), levels=0:MaxT)) /2

# Reference: age difference distribution across all pairs of individuals
tblA.R[, "X"] <- table(factor(abs(AgeDifM), levels=0:MaxT))/2
  
tblA.R
```

Taking the absolute value and dividing by 2 ensures that each pair $i$, $j$ is counted exactly once, both when they are born in different years and when they are born in the same year. Using `factor` with pre-set levels ensures that the sub-tables are lined up properly, and speeds up `table` enormously. 


### Calculating ageprior P(A|R)/P(A)
From the count data $P(A|R)$ is calculated by simply dividing by the column sums: 

```{r calcPAR}
tblA.R <- AP.griffin[["tblA.R"]]    # with dam & sire columns
PA.R <- sweep(tblA.R, 2, STATS=colSums(tblA.R), FUN="/")
round(PA.R, 2)
```

and the age-difference based probability ratio is calculated by dividing each relationship column by the `X` column
```{r calcLRAR}
LR.A.R <- sweep(PA.R[, 1:5], 1, STATS=PA.R[,"X"], FUN="/")
round(LR.A.R, 2)
```

This is the non-flattened, non-smoothed ageprior, stored in list element `[["LR.RU.A.unweighed"]]`. Since in this example `Smooth=FALSE` and `Flatten=FALSE`, it is identical to list element `[["LR.RU.A"]]`, which is the default output. 

### visualisation
The ageprior (or any similar matrix) can be visualised using `PlotAgePrior`, or using for example `image` in base R or `ggplot2` 's `geom_tile()`. 



# Grandparents and avuncular relationships

## change from v1.3.3 to v2.0
The ageprior distribution for avuncular pairs (aunts/uncles) has changed considerably from sequoia v1.3.3 to v2.0. Earlier, it was assumed that the ageprior for maternal and paternal full- and half- aunts and uncles was approximately similar, and approximately symmetrical around zero. From version 2.0 both these assumptions have been dropped, and there are now 6 avuncular classes (FS, MS, or PS of dam or sire) with agepriors that are not symmetrical around zero. The latter allows for the fact that aunts/uncles are often older than their nieces/nephews, but not necessarily so. 

The agepriors for grandparents and avuncular relationships have always been estimated from parent and sibling age distributions, but this estimation has now been moved from `MakeAgePrior` in R (up to v.1.3.3) to the Fortran part. The extra columns are calculated just *before* parentage assignment or sibship clustering. This is in contrast to the 'regular' ageprior, which is updated just after parentage assignment (see also pipeline in Figure \@ref(fig:pipeline)). 


## Mathematical details
The probability that individual $i$ and say its maternal grandfather $j$ have an age difference $A_{i,j}$ is the sum over the possible ages of $i$'s 'in-between' mother $k$, weighed by the agepriors for $k$ being $i$'s mother, and $j$ being $k$'s father, given their age differences.

Or, in mathematical notation: 

\begin{equation}
P(A_{i,j} = y | j \text{ MGF of }i) \propto \sum_{x=0}^{Tmax} \sum_{z=0}^{Tmax} I(x + z = y) P(A_{i,k} = x | k \text{ dam of } i) P(A_{k,j} = z | j \text{ sire of } k) 
\end{equation}

where
\begin{align}
P(A_{i,k} = x | k \text{ dam of } i) = \frac{\alpha_{x, dam}}{\sum \alpha_{., dam}}  \\
P(A_{k,j} = z | j \text{ sire of } k)  = \frac{\alpha_{z, sire}}{\sum \alpha_{., sire}} 
\end{align}

and $I(x + z = y)$ is an indicator that equals $1$ if $x+z=y$, and $0$ otherwise, and $Tmax$ is the maximum age of parents. 

For avuncular pairs, we consider that $j$ might be either an older sibling or a younger sibling of $i$'s parent $k$:

\begin{equation}
P(A_{i,j} = y | j \text{ MFA of }i) \propto \sum_{x=0}^{Tmax} \sum_{z= -Tmax}^{Tmax} I(x + z = y) P(A_{i,k} = x | k \text{ dam of } i) P(A_{k,j} = z | j \text{ full sib of } k) 
\end{equation}

The grandparent and avuncular agepriors are then scaled to allow direct comparison with the other agepriors. The scaling factor used is the average of the sums of the mother and father agepriors; these may differ somewhat from each other if their range of possible ages differs. 




## Griffin example 
The calculation of the extra columns can (currently) only be done by running `sequoia`.
```{r griffin_sim, eval=FALSE}
GenoS <- SimGeno(Ped_griffin, nSnp=400, ParMis=0.4)

Ped_griffin$id <- as.character(Ped_griffin$id)
griffin.sex <- sapply(Ped_griffin$ID, function(x) substr(x, start=nchar(x), stop=nchar(x)))
LH.griffin <- data.frame(ID = Ped_griffin$ID,
                         Sex = ifelse(griffin.sex=="F", 1, 2),
                         BirthYear = Ped_griffin$birthyear)

SeqOUT.griffin <- sequoia(GenoS, LH.griffin,
                  MaxSibIter = 10,
                  args.AP = list(Smooth = FALSE)) 
```

```{r griffinAPxtra, fig.cap="Ageprior with extra columns for Griffin example", fig.width=8}
data(SeqOUT_griffin)   # example output included with package
PlotAgePrior(SeqOUT_griffin$AgePriorExtra)
```

This `AgePriorExtra` (Figure \@ref(fig:griffinAPxtra)) illustrates the following universal properties:

* the ageprior for siblings is symmetrical around 0; 
* grandparents are at least 2 time units older than their grand-offspring;
* aunts/uncles tend to be older than their nieces/nephews, but may be younger;

as well as subtle differences between the different avuncular classes, which may be 
more pronounced for some life histories (abbreviations: e.g. MPA = maternal-paternal avuncular, or paternal half-sibling of mother). 


# Small sample correction
With small to moderate sample sizes some age/relationship combinations may be over- or under- represented by chance, or even completely absent. In the latter case, the uncorrected ageprior  $\alpha_{A,R} = P(A|R)/P(A) = 0$ specifies for subsequent pedigree reconstruction specifies that this age/relationship combination is always impossible, including for unsampled pairs and for individuals of unknown (estimated) age. This can result of all kinds of false negatives, as well as false positives: e.g. when half-siblings $i$ and $j$ can (wrongly-presumed) not be maternal siblings, they must be paternal siblings, and subsequently more individuals may be added to this erroneous paternal sibship. Alternatively, certain age-differences may be biologically plausible but do not occur at all in the sample, causing the denominator in Equation \@ref(eq:AP2) $P(A)=0$, and $\alpha_{A,R}$ to be undefined. 

To minimise these problems, `MakeAgePrior` applies two generic, one-size-fits-most corrections: a fairly subtle `Smooth`ing of dips and padding of the tails, and a more rigorous `Flatten`ing approach. 



## Smooth
By default, `MakeAgePrior` smooths the distribution, by 'stretching' the tails and smoothing over any dips. It is set to `TRUE` by default, as in many cases only a subset of the possible parental ages is sampled, and there is only rarely a detrimental effect of smoothing when none is required. Exception is when generations are discrete (non-overlapping); when this is detected or declared by the user, `Smooth` is set to `FALSE` automatically. 

There is no way to change the behaviour of `Smooth`, except turning it on or off. If for example a longer tail is required at either end, or more thorough smoothing is needed, either consider using `Flatten=TRUE`, or edit the ageprior matrix by hand and supply it as an argument to `sequoia`:

```{r OwnAP, eval=FALSE, error=FALSE}
MyAgePrior <- MakeAgePrior(LifeHistData)   # with or without scaffold/old pedigree 
# (...)  # edit by hand
SeqOUT <- sequoia(Genotypes, LifeHistData, SeqList(AgePriors = MyAgePrior))
```


## Flatten
The ageprior is 'flattened' by taking a weighed average between the flat, default 0/1 ageprior without any pedigree information, and the ageprior based on the pedigree and birth years:

\begin{equation}
	\alpha* = W_R \times \frac{P(A|R, \text{sampled})}{P(A, \text{sampled})} +
	 (1-W_R) \times I(A,R)
\end{equation}

where $W_R$ is a weight based on $N_{R}$, the number of pairs with relationship $R$ and known age difference, and $I(A,R)$ is a 0/1 indicator whether the age-relationship combination is biologically possible (0 for $A=0$ and parent-offspring, 1 otherwise).

The weighing factor $W_R$ follows a sigmoid curve, calculated as $W_R = 1 - \exp(-\lambda_{N,W} * N_R)$. By default $\lambda_{N,W}=-\log(0.5)/100$, which corresponds to $W_R< 0.5$ if $N_R < 100$, but this can be altered with argument \verb+lambdaNW+ (see Figure \@ref(fig:plotFlattenNW)). When $\lambda_{N,W}$ is large (say $>0.2$), `Flatten` is effectively only applied for relationships with very small sample size and the 'raw' pedigree-based estimate used otherwise, while when $\lambda_{N,W}$ is small (say $<0.0001$) there is effectively no contribution of the pedigree.

```{r 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")
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)
```


### Why not use the sampling probability
In theory, the sampling probability per relationship could be estimated, and used to calculate some correction factor for equation \@ref(eq:AP2). In practice, this requires some strong assumptions about the number of assignable parents (i.e., distinction between pedigree founders and non-founders) and the number of sibling pairs (i.e. distribution of sibship sizes). This correction factor would depend on the total number of sampled individuals in a non-nonsensical way: removing individuals without assigned parents from the pedigree would 'magically' increase our faith in the ageprior estimated from the remainder of the pedigree. Therefore, the degree of correction by `Flatten` depends only on the number of individuals per relationship with known age difference, $N_R$. 



## Example: griffins
In the `smooth`ed to the original ageprior, future parents may be age 4 (with an ageprior half that of age 3), and very rarely even age 5, 
(Figure \@ref(fig:griffinSmoothFlatten)), while in the scaffold pedigree they are never older than age 3 (see also section [With pedigree, overlapping generations]). Resultingly, the age difference between siblings can be larger too. If the observed parental ages had started at say age 7, the ageprior for age 6 would be set to a small value (0.001) too, but not age 5: it is assumed that larger age differences are easier to have been missed than smaller age differences. 

```{r 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[["Flat"]] <- MakeAgePrior(Ped_griffin, Smooth=FALSE, Flatten=TRUE, Plot=FALSE)
AP_list[["Flat + Smooth"]] <- MakeAgePrior(Ped_griffin, Smooth=TRUE, Flatten=TRUE, Plot=FALSE)

par(mfrow=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)
}
```

With `Flatten=TRUE` the ageprior looks rather different, with (as is typical) negligible additional effect of `Smooth` or not (bottom row in Figure \@ref(fig:griffinSmoothFlatten)). The weight for maternal siblings (MS) is larger than for paternal siblings (PS), 
```{r griffinweights}
AP.griffin$Weights
```

as by chance there is a large maternal sibship (see `SummarySeq(Ped_griffin)`), resulting in more maternal than paternal sibling pairs:
```{r griffinNR}
N.R <- colSums(AP.griffin$tblA.R)
N.R
```

and we can use these counts of pairs to double check the weights:
```{r griffinweights2}
lambdaNW <- -log(0.5)/100   # the default
W.R <-  1 - exp(-lambdaNW * N.R)
round(W.R, 4)
```



## Example 2: No candidate fathers genotyped
In some study systems, it is very difficult or too expensive to genotype the (candidate/known) parents of one sex. In this example we assume that no candidate fathers are genotyped in a population with discrete generations, but the same principle holds when mothers are not genotyped and/or generations overlap. 

```{r SimNoSire, message=FALSE, results="hide"}
data(Ped_HSg5, LH_HSg5)
# Simulate data: 20% of mothers non-genotyped & 100% of fathers
GM <- SimGeno(Ped_HSg5, nSnp=400, ParMis=c(0.2, 1.0))  

ParOUT <- sequoia(GM, LH_HSg5, MaxSibIter = 0, Plot=FALSE)
```

```{r simNoSirePlot, fig.cap="Example: no fathers genotyped", out.width="70%"}
PlotAgePrior(ParOUT$AgePriors)
```

In the output the maximum age of parents is stated as 1 for mothers (from the assigned mother-offspring pairs), and 5 for fathers (using as best guess the observed age range in the lifehistory data). Since there are relationships with fewer than 20 pairs with known age difference (P, PS and FS), `Flatten` is automatically set to `TRUE`. 

The ageprior for fathers $\alpha_{A,R=P}$ equals $1$ for age differences $A$ between 1 and `MaxAgeParent`, inclusive. Since there is no convincing evidence that the generations are non-overlapping, the default `Smooth=TRUE` is used, which adds a 'tail' of 2 extra possible years of age differences for all relationship categories. 

If we know for sure that fathers cannot be older than 3 years, but want to estimate the maximum age of mothers from the data, we can use
```{r NoSire2, eval=FALSE}
ParOUT <- sequoia(GM, LH_HSg5, MaxSibIter = 0,
                  args.AP = list(MaxAgeParent=c(NA, 3),
                                 Smooth = FALSE))
```

Here `Flatten=FALSE` would simply be ignored, as the ageprior would be undefined. Using `Flatten=FALSE` with 5--20 pairs for some relationships is discouraged and results in a warning.

