Kinship_Matrix_Code_Details
In kinship2: Pedigree Functions

Introduction

The kinship matrix is foundational for random effects models with family data.
For $n$ subjects it is an $n \times n$ matrix whose $ij$ element contains the expected fraction of alleles that would be identical by descent if we sampled one from subject $i$ and another from subject $j$. Note that the diagonal elements of the matrix will be 0.5 not 1: when we randomly sample twice from the same subject (with replacement) we will get two copies of the gene inherited from the father 1/4 of the time, the maternal copy twice (1/4) or one of each 1/2 the time. The formal definition is $K(i,i) = 1/4 + 1/4 + 1/2 K(m,f)$ where $m$ and $f$ are the father and mother of subject $i$.

The algorithm used is found in K Lange, \emph{Mathematical and Statistical Methods for Genetic Analysis}, Springer 1997, page 71--72.

The key idea of the recursive algorithm for $K(i,j)$ is to condition on the gene selection for the first index $i$. Let $m(i)$ and $f(i)$ be the indices of the mother and father of subject $i$ and $g$ be the allele randomly sampled from subject $i$, which may of either maternal or paternal origin.

\begin{align} K(i,j) &= P(\mbox{$g$ maternal}) * K(m(i), j) + P(\mbox{$g$ paternal}) * K(f(i), j) \label{recur0} \ &= 1/2 K(m(i), j) + 1/2 K(f(i), j) \label{recur1} \ K(i,i) &= 1/2(1 + K(m(i), f(i))) \label{self} \end{align}

The key step in equation \eqref{recur0} is if $g$ has a maternal origin, then it is a random selection from the two maternal genes, and its IBD state with respect to subject $j$ is that of a random selection from m(i) to a random selection from $j$. This is precisely the definition of $K(m(i), j)$. The recursion does not work for $K(i,i)$ in equation \eqref{self} since once we select a maternal gene the second choice from $j$ cannot use a different maternal gene.

For the recurrence algorithm to work properly we need to compute the values of $K$ for any parent before the calculations for their children. Pedigree founders (those with no parents) are assumed to be unassociated, so for these subjects we have \begin{align} K(i,i) &= 1/2 K(i,j) &=0 \; i\ne j \end{align}

The final formula slightly different for the $X$ chromosome. Equation \ref{recur0} still holds, but for males the probability that a selected $X$ chromosome is maternal is 1, so when $i$ a male the recurrence formula becomes $K(i,j) = K(m(i),j)$.
For females it is unchanged. All males will have $K(i,i) = 1$ for the $X$ chromosome.

In order to have already-defined terms on the right hand side of the recurrence formula for each element, subjects need to be processed in the following order \begin{enumerate} \item Generation 0 (founders) \item $K(i,j)$ where $i$ is from generation 1 and $j$ from generation 0. \item $K(i,j)$ with $i$ and $j$ from generation 1 \item $K(i,j)$ with $i$ from generation 2 and $j$ from generation 0 or 1 \item $K(i,j)$ with $i$ and $j$ from generation 2. \item \ldots \end{enumerate}

The kindepth routine assigns a plotting depth to each subject in such a way that parents are always above children.
For each depth we need to do the compuations of formula \eqref{recur} twice. The first time it will get the relationship between each subject and prior generations correct, the second will correctly compute the values between subjects on the same level. The computations within any stage of the above list can be vectorized, but not those between stages.

Let [[indx]] be the index of the rows for the generation currently being processed, say generation $g$.
We add correct computations to the matrix one row at a time; all of the calculations depend only on the prior rows with the exception of the [i,i] element. This approach leads to a for loop containing operations on single rows/columns.

At one point below we use a vectorized version. It looks like the snippet below

for (g in 1:max(depth)) {
    indx <- which(depth==g)
    kmat[indx,] <- (kmat[mother[indx],] + kmat[father[indx], ])/2
    kmat[,indx] <- (kmat[,mother[indx]] + kmat[,father[indx],])/2
    for (j in indx) kmat[j,j] <- (1 + kmat[mother[j], father[j]])/2
}

The first line computes all the values for a horizontal stripe of the matrix. It will be correct for columns in generations $<g$, unreliable for generation $g$ with itself because of incomplete parental relationships, and zero for higher generations. The second line does the vertical stripe, and because of the line before it does have the data it needs and so gets all the stripe correct. Except of course for the diagonal elements, for which formula \eqref{recur1} does not hold. We fill those in last. We know that vectorized calculations are always faster in R and I was excited to figure this out. The unfortunate truth is that for this code it hardly makes a difference, and for the X chromosome calculation leads to impenetrable if-then-else logic.

The program can be called with a pedigree, a pedigree list, or raw data. The first argument is [[id]] instead of the more generic [[x]] for backwards compatability with an older version of the routine. We give founders a fake parent of subject $n+1$ who is not related to anybody (even themself); it avoids some if-then-else constructions.

For S3 method dispatch to work on pedigree and pedigreeList objects, we need to define a kinship(), kinship.default(), kinship.pedigree(), and kinship.pedigreeList(). The default version takes the vectors

Pedigree Object

The method for a pedigree object is an almost trivial modification. Since the mother and father are already indexed into the id list it has two lines that are different, those that create mrow and drow. The other change is that now we potentially have information available on monozygotic twins. If there are any such, then when the second twin of a pair is added to the matrix, we need to ensure that the pairs kinship coefficient is set to the self-self value. This can be done after each level is complete, but before children for that level are computed. If there are monozygotic triples, quadruplets, etc. this computation gets more involved.

The total number of monozygotic twins is always small, so it is efficient to fix up all the monzygotic twins at each generation. A variable [[havemz]] is set to TRUE if there are any, and an index array [[mzindex]] is created for matrix subscripting.

PedigreeList Object

For the Minnesota Family Cancer Study there are 461 families and 29114 subjects. The raw kinship matrix would be 29114 by 29114 which is over 5 terabytes of memory, something that clearly will not work within S.
The solution is to store the overall matrix as a sparse Matrix object. Each family forms a single block. For this study we have [[n <- table(minnbreast$famid); sum(n*(n+1)/2)]] or 1.07 million entries; assuming that only the lower half of each matrix is stored. The actual size is actually smaller than this, since each family matrix will have zeros in it --- founders for instance are not related --- and those zeros are also not stored.

The result of each per-family call to kinship will be a symmetric matrix. We first turn each of these into a dsCMatrix object, a sparse symmetric form. The [[bdiag]] function is then used to paste all of these individual sparse matrices into a single large matrix.

Why do we note use [[(i in famlist)]] below? A numeric subscript of [[[9]]] %' selects the ninth family, not the family labeled as 9, so a numeric family id would not act as we wished. If all of the subject ids are unique, across all families, the final matrix is labeled with the subject id, otherwise it is labeled with family/subject.

MakeKinship

The older [[makekinship]] function, from before the creation of pedigreeList objects, accepts the raw identifier data, along with a special family code for unrelated subjects, as produced by the [[makefamid]] function. All the unrelated subjects are put at the front of the kinship matrix in this case rather than within the family. Because unrelateds get put into a fake family, we cannot create a rational family/subject identifier; the id must be unique across families. We include a copy of the routine for backwards compatability, but do not anticipate any new usage of it. Like most routines, this starts out with a collection of error checks.

Monozygotic Twins, for Pedigree object

Return now to the question of monzygotic sets, used specifically in the kinship for pedigree objects. Consider the following rather difficult example: \begin{verbatim} 1 2 1 3 5 6 3 7 10 9 \end{verbatim} Subjects 1, 2, 3, and 7 form a monozygotic quadruple, 5/6 and 9/10 are monzygotic pairs.
First create a vector \code{mzgrp} which contains for each subject the lowest index of a monozygotic twin for that subject.
For non-twins it can have any value.
For this example that vector is set to 1 for subjects 1, 2, 3, and 7, to 5 for 5 and 6, and to 9 for 9 and 10. Creating this requires a short while loop. Once this is in hand we can identify the sets.

Next make a matrix that has a row for every possible pair. Finally, remove the rows that are identical. The result is a set of all pairs of observations in the matrix that correspond to monozygotic pairs.

Older routines

For testing purposes we have a version of the kinship function prior to addition of the chrtype argument.

oldkinship <- function(id, ...) {
    UseMethod('oldkinship')
    }

oldkinship.default <- function(id, dadid, momid, ...) {
    n <- length(id)
    if (n==1) 
        return(matrix(.5,1,1, dimnames=list(id, id)))
    if (any(duplicated(id))) stop("All id values must be unique")
    kmat <- diag(n+1) /2
    kmat[n+1,n+1]    <- 0 

    pdepth <- kindepth(id, dadid, momid)
    mrow <- match(momid, id, nomatch=n+1) #row number of the mother
    drow <- match(dadid, id, nomatch=n+1) #row number of the dad 

    for (depth in 1:max(pdepth)) {
        indx <- (1:n)[pdepth==depth]
        for (i in indx) {
            mom <- mrow[i]
            dad <- drow[i]
            kmat[i,]  <- kmat[,i] <- (kmat[mom,] + kmat[dad,])/2
            kmat[i,i] <- (1+ kmat[mom,dad])/2
            }
        }

    kmat <- kmat[1:n,1:n]
    dimnames(kmat) <- list(id, id)
    kmat
    }

oldkinship.pedigree <- function(id, ...) {
    n <- length(id$id)
    if (n==1) 
        return(matrix(.5,1,1, dimnames=list(id$id, id$id)))
    if (any(duplicated(id$id))) stop("All id values must be unique")
    kmat <- diag(n+1) /2
    kmat[n+1,n+1]    <- 0 

    pdepth <- kindepth(id)
    mrow <- ifelse(id$mindex ==0, n+1, id$mindex)
    drow <- ifelse(id$findex ==0, n+1, id$findex)

    # Are there any MZ twins to worry about?
    if (!is.null(id$relation) && any(id$relation$code=="MZ twin")) {
        havemz <- TRUE
        temp <- which(id$relation$code=="MZ twin")
        ## drop=FALSE added in case only one MZ twin set
        mzmat <- as.matrix(id$relation[,c("indx1", "indx2")])[temp,,drop=FALSE]

        # any triples, quads, etc?
        if (any(table(mzmat) > 1)) { #yes there are
            # each group id will be min(member id)
            mzgrp <- 1:max(mzmat)  #each person a group
            indx <- sort(unique(as.vector(mzmat)))
            # The loop below will take k-1 iterations for a set labeled as
            #   1:2, 2:3, ...(k-1):k;  this is the worst case.
            while(1) {
                z1 <- mzgrp[mzmat[,1]]
                z2 <- mzgrp[mzmat[,2]]
                if (all(z1 == z2)) break
                mzgrp[indx] <- tapply(c(z1, z1, z2, z2), c(mzmat,mzmat), min)
            }
            # Now mzgrp = min person id for each person in a set
            matlist <- tapply(mzmat, mzgrp[mzmat], function(x) {
                x <- sort(unique(x))
                temp <- cbind(rep(x, each=length(x)), rep(x, length(x)))
                temp[temp[,1] != temp[,2],]
                })
            }
        else {  #no triples, easier case
            matlist <- tapply(mzmat, row(mzmat), function(x) 
                            matrix(x[c(1,2,2,1)],2), simplify=FALSE)
            }
        }
    else havemz <- FALSE

    for (depth in 1:max(pdepth)) {
        indx <- (1:n)[pdepth==depth]
        for (i in indx) {
            mom <- mrow[i]
            dad <- drow[i]
            kmat[i,]  <- kmat[,i] <- (kmat[mom,] + kmat[dad,])/2
            kmat[i,i] <- (1+ kmat[mom,dad])/2
            }
        if (havemz) {
            for (i in 1:length(matlist)) {
                temp <- matlist[[i]]
                kmat[temp] <- kmat[temp[1], temp[1]]
            }
        }
    }

    kmat <- kmat[1:n,1:n]
    dimnames(kmat) <- list(id$id, id$id)
    kmat
}    

oldkinship.pedigreeList <- function(id, ...) {
    famlist <- unique(id$famid)
    nfam <- length(famlist)
    matlist <- vector("list", nfam)
    idlist  <- vector("list", nfam) #the possibly reorderd list of id values

    for (i in 1:length(famlist)) {
        tped <- id[i]  #pedigree for this family
        temp <- try(oldkinship(tped, ...), silent=TRUE)
        if (class(temp)=="try-error") 
            stop(paste("In family", famlist[i], ":", temp))
        else matlist[[i]] <- as(forceSymmetric(temp), "dsCMatrix")
        idlist[[i]] <- tped$id
    }

    result <- bdiag(matlist)
    if (any(duplicated(id$id))) 
        temp <-paste(rep(famlist, sapply(idlist, length)),
                     unlist(idlist), sep='/') 
    else temp <- unlist(idlist)

    dimnames(result) <- list(temp, temp)
    result
}