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, , 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.

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

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
}