Logo Search packages:      
Sourcecode: lme4 version File versions  Download package

pedigree.c

                        /* Methods for pedigrees */
#include "pedigree.h"

/**
 * Create the left Cholesky factor of the numerator relationship
 * matrix from a pedigree.
 *
 * @param x a pedigree object
 * @param ans T stored as a non-unit, lower dtCMatrix
 *
 * @return ans with elements modified to incorporate D
 */
SEXP pedigree_chol(SEXP x, SEXP ans)
{
    SEXP Sire = GET_SLOT(x, install("sire"));
    int *ai = INTEGER(GET_SLOT(ans, lme4_iSym)),
      *ap = INTEGER(GET_SLOT(ans, lme4_pSym)),
      *dam = INTEGER(GET_SLOT(x, install("dam"))),
      *sire = INTEGER(Sire), 
      i, j, n = LENGTH(Sire);
    double *ax = REAL(GET_SLOT(ans, lme4_xSym)), *F, Di, tmp;

    setAttrib(ans, install("F"), allocVector(REALSXP, n));
    F = REAL(getAttrib(ans, install("F")));
    for (i = 0; i < n; i++) {
      int k, p = sire[i] - 1, q = dam[i] - 1;
      if (sire[i] == NA_INTEGER) {
          F[i] = 0;
          Di = (dam[i] == NA_INTEGER) ? 1 : sqrt(0.75 - 0.25 * F[q]);
      } else {
          if (dam[i] == NA_INTEGER) { /* sire only */
            F[i] = 0;
            Di = sqrt(0.75 - 0.25 * F[p]);
          } else {            /* both parents in pedigree */
            Di = sqrt(0.5 - 0.25 * (F[p] + F[q]));
            F[i] = NA_REAL;
            if ((ap[i + 1] - ap[i]) > 1) {        /* skip if no progeny */
                if (p > q) {j = p; p = q; q = j;} /* ensure p <= q */
                for (j = 0, F[i] = 0; j <= p; j++) {
                  for (k = ap[j], tmp = 0;
                       k < ap[j + 1] && ai[k] <= q; k++) {
                      int ii = ai[k];
                      if (ii == p) tmp = ax[k];
                      if (ii == q) F[i] += tmp * ax[k]/2;
                  }
                }
            }
          }
      }
      for (j = ap[i]; j < ap[i + 1]; j++) ax[j] *= Di;
    }
    return ans;
}

/* NOTE: This function requires that missing parents be coded as zero */
/**
 * Create the inbreeding coefficients according to the algorithm given
 * in "Comparison of four direct algorithms for computing inbreeding
 * coefficients" by Mehdi Sargolzaei and Hiroaki Iwaisaki, Animal
 * Science Journal (2005) 76, 401--406.  This function is a modified
 * version of the code published in an appendix to that paper.
 *
 * @param x a pedigree object
 *
 * @return a list of the inbreeding coefficients
 */

SEXP pedigree_inbreeding(SEXP x)
{
    SEXP ans, sp = GET_SLOT(x, install("sire"));
    int i, j, t, n = LENGTH(sp), S, D;
    int *sire = INTEGER(sp),
      *dam = INTEGER(GET_SLOT(x, install("dam"))),
      *Anc = Calloc(n + 1, int), /* ancestor */
      *SI, *MI;            /* start and minor */
    double *F = Calloc(n + 1, double), /* inbreeding coefficients */
      *L = Calloc(n + 1, double),
      *B = Calloc(n + 1, double); 
    int *LAP = Calloc(n + 1, int);  /* longest ancestoral path */
    
    F[0] =-1; LAP[0] =-1; /* set F and lap for unknown parents */
    for(i = 1, t = -1; i <= n; i++) {     /* evaluate LAP and its maximum */
      S = sire[i]; D = dam[i]; /* parents of animal i */
      LAP[i] = ((LAP[S] < LAP[D]) ? LAP[D] : LAP[S]) + 1;
      if (LAP[i] > t) t = LAP[i];
    }
    SI = Calloc(t + 1, int);
    MI = Calloc(t + 1, int);
    for(i = 0; i <= t ; ++i) SI[i] = MI[i] = 0; /* initialize start and minor */
    for(i = 1; i <= n; i++) {       /* evaluate F */
      S = sire[i]; D = dam[i]; /* parents of animal i */
      B[i] = 0.5 - 0.25 * (F[S] + F[D]); 
                        /* adjust start and minor */
      for (j = 0; j < LAP[i]; j++) {++SI[j]; ++MI[j];} 
      if (S == 0 || D == 0) { /* both parents unknown */
          F[i] = L[i] = 0; continue;
      }
      if(S == sire[i-1] && D == dam[i-1]) { /* full-sib with last animal */
          F[i] = F[i-1]; L[i] = L[i-1]; continue;
      }
    
      F[i] = -1; L[i] = 1; 
      t = LAP[i]; /* largest lap group number in the animal's pedigree */
      Anc[MI[t]++] = i; /* initialize Anc and increment MI[t] */
      while(t > -1) { /* from the largest lap group to zero */
          j = Anc[--MI[t]]; /* next ancestor */
          S = sire[j]; D = dam[j]; /* parents of the ancestor */
          if (S) {
            if (!L[S]) Anc[MI[LAP[S]]++] = S; 
                        /* add sire in its lap group in Anc
                         * array if it is not added yet and
                         * increment the minor index for the group */ 
            L[S] += 0.5 * L[j]; /* contribution to sire */
          }
          if (D) {
            if (!L[D]) Anc[MI[LAP[D]]++] = D;
            L[D] += 0.5 * L[j]; /* contribution to dam */
          }
          F[i] += L[j] * L[j] * B[j];
          L[j] = 0; /*clear L[j] for the evaluation of the next animal */
          if (MI[t] == SI[t]) --t; /* move to the next lap group when
                              * all ancestors in group t have been
                              * evaluated */
      } 
    }
    ans = PROTECT(allocVector(REALSXP, n));
    Memcpy(REAL(ans), F + 1, n);
    Free(Anc); Free(F); Free(L); Free(B); Free(LAP); Free(SI); Free(MI);
    UNPROTECT(1);
    return ans;
}


Generated by  Doxygen 1.6.0   Back to index