ESTIMATION  AND  ANALYSIS
OF  A  COVARIANCE  FUNCTION:

SYMMETRIC  COEFFICIENTS

Written by:
      Mark Kirkpatrick
      Address:  Department of Zoology, University of Texas, Austin TX 78712  USA
      Telephone:  (512) 471-5996, -7131
      Email:  kirkp@biff.zo.utexas.edu,  kirkp@utxvms.bitnet
      Fax:  (512) 471-9651

GENERAL  INFORMATION

Introduction

"Infinite-dimensional" traits are characters like growth trajectories, reaction norms, and morphological shapes in which individuals are described by a function rather than one or a few discrete measurements.  This notebook analyzes variation in an infinite-dimensional trait.  In particular, these programs estimate and analyze the covariance function for the trait.

The infinite-dimensional approach has been developed in a series of publications that are listed below in the "References" section.  The specific methods used in this notebook are described in Kirkpatrick et al. (1990).  Applications of the methods are illustrated in Kirkpatrick and Lofsvold (1992).  This notebook assumes that the user is familiar with the approach in general and those papers in specific.

The user inputs the covariance matrix and age vector.  If a "reduced" fit (which smooths the covariance function) is being estimated then one of the options for calculating the error covariance matrix V must be chosen from the options described in the section for data input.  

For the sake of concreteness, the notation in this notebook assumes that the trait of interest is a growth trajectory.  The methods, however, are general to a variety of other types of traits (including reaction norms and shapes).  For these other applications, the variables in this notebook should be interpreted appropriately.  When analyzing a reaction norm, for example,  the user should supply the environmental states rather than ages for the variable "ageVector".

The method of symmetric coefficients used in this notebook forces there to be a continuous first derivative along the diagonal of the covariance function (no "ridge").  (An alternative method called "asymmetric coefficients" does not make that requirement.  The method of asymmetric coefficients is described in Kirkpatrick et al. (1993 ms), and is implimented in the Mathematica notebook "Asymmetric coefficients".)

The orthogonal functions used in this notebook are the Legendre polynomials.

References

Key references are denoted with (*).

     Kirkpatrick, M., and N. Heckman.  1989.  A quantitative-genetic model for growth, shape, reaction norms, and other infinite-dimensional characters.  Journal of Mathematical Biology 27: 429-450.
           
     * Kirkpatrick, M., D. Lofsvold, and M. Bulmer.  1990.  Analysis of the inheritance, selection and evolution of growth trajectories.  Genetics 124: 979-993.
           
     * Kirkpatrick, M., and D. Lofsvold.  1992.  Measuring selection and constraint in the evolution of growth.  Evolution 46: 954-971.
           
     * Gomulkiewicz, R., and M. Kirkpatrick  1992. Quantitative genetics and the evolution of reaction norms.  Evolution 46: 390-411.
     
     Kirkpatrick, M.  1993.  The evolution of size and growth in harvested natural populations.  Pages 145 - 154 in T.K. Stokes, J.M. McGlade, and R. Law (eds.), The Exploitation of Evolving Resources.  Lecture Notes in Biomathematics 99.  Springer-Verlag, Berlin.
           
     Kirkpatrick, M., W.G. Hill, and R. Thompson.  1994.  Estimating the covariance structure of traits during growth and aging, illustrated with lactation in dairy cattle.  Genetical Research (in press).

Release notes

Release of 3 March 1994.

This notebook may be freely distributed provided that no charge is made.  I would appreciate it if people who receive a copy from a third party would notify me.  (This makes it possible for me to let you know of upgrades and bugs.)  

Development of this software was supported by grants from  NSF and NIH.

This notebook is based on the notebook "new_sym_fit" written in June 1992.  Major changes are:  (1) the error covariance matrix V can now be calculated by several methods;  (2) eigenvalues and eigenfunctions are calculated;  (3) eigenvalues and eigenfunctions are now normalized over the range of the coordinate vector (rather than the range of the orthogonal functions).

INPUT  THE  DATA  AND  PARAMETERS

The following data are provided as an example.  They are estimates of the additive genetic components of lactation rate in British Friesian-Holstein dairy cattle during the first 10 test dates.  The data are described fully in:
    
      Pander, B.L., W.G. Hill, and R. Thompson.  1993.   Genetic parameters of test day records of British Holstein-Friesian heifers.  Anim. Prod. (in press).
      
Results from the analysis of this data set using this notebook are described in Kirkpatrick et al. (1993) (full reference given above).

Input the age vector

ageVector = {1., 2., 3., 4., 5., 6., 7., 8., 9., 10.} ;

nAges = Dimensions[ageVector][[1]];

Input the covariance matrix

dataMatrix = 
  {{188.76, 210.8, 228.64, 197.08, 196.68, 193.32, 176.24, 164.44, 150., 151.8}, 
   {210.8, 289.72, 282., 273.28, 265.64, 243.08, 245.76, 244.2, 205.4, 185.6}, 
   {228.64, 282., 325.4, 307., 281.76, 293.84, 279.16, 274.88, 231.04, 217.28}, 
   {197.08, 273.28, 307., 335.24, 317.24, 315.24, 316.2, 323.84, 276.8, 244.36}, 
   {196.68, 265.64, 281.76, 317.24, 349.8, 307.24, 315.32, 331.36, 306.04, 299.68}, 
   {193.32, 243.08, 293.84, 315.24, 307.24, 313.64, 302.92, 306.44, 274.96, 253.4}, 
   {176.24, 245.76, 279.16, 316.2, 315.32, 302.92, 307.04, 320.2, 281.48, 258.44}, 
   {164.44, 244.2, 274.88, 323.84, 331.36, 306.44, 320.2, 343.16, 298.88, 280.8}, 
   {150., 205.4, 231.04, 276.8, 306.04, 274.96, 281.48, 298.88, 279.48, 270.76}, 
   {151.8, 185.6, 217.28, 244.36, 299.68, 253.4, 258.44, 280.8, 270.76, 296.4}};

Choose the number of orthogonal functions to be fit 
(for a "full" or "reduced" fit)

The variable kPolys determines the number of orthogonal polynomials that will be fit to the data matrix.  This number cannot be larger than nAges, which is the dimensionality of the data matrix abovel.  
   A "full" fit results when kPolys = nAges.  In this case, the estimated covariance function exactly fits the original data.  A "reduced" fit occurs when kPolys < nAges.  In this case, the estimated covariance function is smoothed.

kPolys = 9;

IF YOU ARE DOING A REDUCED FIT:
Choose ONE of the following methods for estimating the error covariance matrix,  V

Methods for calculated the error covariance matrix V are described in Appendix B of Kirkpatrick et al. (1990) (a full reference for which is given at the beginning of the notebook).  The user may choose one of the following options, or supply V as calculated by his/her own method.  (Methods for estimating V when fitting an additive genetic covariance function that was estimated by parent-offspring or by full-sib breeding designs are not included here;  those are described in Kirkpatrick et al. (1990) Appendix B and could easily be added.)

Calculate the error covariance matrix V, Option (1):
    Using the approximation   
             V[ij,kl]  =  alpha (  M[ik] M[jl]  +  M[il] M[jk]  ) ,
             where alpha is an unknown constant of proportionality

This method gives an approximation to the the V matrix that is proportional to a proper estimate of V.  The unknown proportionality factor does not affect the fit.  This method can be used when there is not sufficient information about the experimental design to use another method to estimate V (such as those given below).
     NOTE:  Since we don't know the proportionality factor, the chi-squared value for the fit of the estimate to the data calculated below is incorrect.  The significance of a reduced fit cannot be tested statistically.

Clear[vMat, i, j, index3, index4];
     index3 = index3Func[nAges];
     index4 = index4Func[nAges];
     vMat = Table[0, {i, nAges(nAges+1)/2}, {j, nAges(nAges+1)/2}];
     Do[Do[
        vMat[[i,j]] = (  dataMatrix[[ index3[[i]], index3[[j]] ]]
                         * dataMatrix[[ index4[[i]], index4[[j]] ]]
                       + dataMatrix[[ index3[[i]], index4[[j]] ]]
                         * dataMatrix[[ index4[[i]], index3[[j]] ]]  );
        vMat[[j,i]] = vMat[[i,j]],
        {i, nAges(nAges+1)/2} ],
        {j, nAges(nAges+1)/2}
       ]

Calculate the error covariance matrix V, Option (2):
    Assuming we are fitting an additive genetic covariance function estimated from a balanced half-sib design

In this case, "dataMatrix" input at the top of this notebook is the additive genetic covariance matrix G.  To calculate V, we also need the phenotypic covariance matrix P.
     The data below are again for dairy cattle, from Pander et al. (1993).

Input the data and parameters

pMatrix = {{971.02, 588.38, 521.27, 458.97, 428.53, 396.37, 372.6, 353.87, 321.6, 278.82}, 
  {588.38, 934.81, 654.96, 588.46, 544.08, 500.77, 477.9, 453.41, 414.35, 360.51}, 
  {521.27, 654.96, 936.04, 655.75, 599.29, 563.51, 532.45, 499.94, 459.22, 405.76}, 
  {458.97, 588.46, 655.75, 878.96, 634.21, 586.86, 560.1, 535.35, 490.68, 429.88}, 
  {428.53, 544.08, 599.29, 634.21, 829.65, 609.58, 580.53, 551.52, 509.24, 457.48}, 
  {396.37, 500.77, 563.51, 586.86, 609.58, 794.42, 596.57, 559.73, 509.95, 456.28}, 
  {372.6, 477.9, 532.45, 560.1, 580.53, 596.57, 781.14, 591.02, 541., 483.47}, 
  {353.87, 453.41, 499.94, 535.36, 551.52, 559.73, 591.02, 785.97, 581.05, 524.59}, 
  {321.6, 414.35, 459.22, 490.68, 509.24, 509.95, 541., 581.05, 779.63, 587.07}, 
  {278.82, 360.51, 405.76, 429.88, 457.48, 456.28, 483.47, 524.59, 587.07, 890.93}};

nDamsPerSire = 10;

sSires = 700;

rResidualDF  =  16000;

Preliminary calculations:  Estimates of the mean squares

mE  =  pMatrix  -  (1/4) dataMatrix;

mA  =  (nDamsPerSire - 1) dataMatrix / 4  +  pMatrix;

Calculate V

Clear[index2Func];
index2Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexVal = 0;
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = indexVal;
              indexArg = indexArg+1,
              {i1, 1, k-i2}
             ];
           indexVal = indexVal+1,
           {i2, 0, k-1}
          ];
        vector
      )

Clear[index3Func];
index3Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = i1;
              indexArg = indexArg+1,
              {i1, i2, k}
             ],
           {i2, 1, k}
          ];
        vector )

Clear[index4Func];
index4Func[k_] := index2Func[k] + 1;

Clear[vMat, i, j, index3, index4];
     index3 = index3Func[nAges];
     index4 = index4Func[nAges];
     vMat = Table[0, {i, nAges(nAges+1)/2}, {j, nAges(nAges+1)/2}];
     Do[Do[
        covMaMa  =  ( mA[[ index3[[i]], index3[[j]] ]] * mA[[ index4[[i]], index4[[j]] ]] 
                     + mA[[ index3[[i]], index4[[j]] ]] * mA[[ index4[[i]], index3[[j]] ]] ) / (sSires + 2);
        covMeMe  =  ( mE[[ index3[[i]], index3[[j]] ]] * mE[[ index4[[i]], index4[[j]] ]] 
                     + mE[[ index3[[i]], index4[[j]] ]] * mE[[ index4[[i]], index3[[j]] ]] ) / (rResidualDF + 2);
        vMat[[i,j]] = vMat[[j,i]] = (16. / nDamsPerSire^2) (covMaMa + covMeMe),
        {i, j, nAges(nAges+1)/2} ],
        {j, 1, nAges(nAges+1)/2}
       ]

Calculate the error covariance matrix V, Option (3):
    Assuming we are fitting a phenotypic covariance function estimated from a balanced half-sib design

In this case, "dataMatrix" input at the top of this notebook is the phenotypic covariance matrix P.  To calculate V, we also need the additive genetic covariance matrix G.  The following algorithm is from the notes of MK for 1-V-1993.
     The data below are again for dairy cattle, from Pander et al. (1993).

Input the data and parameters

gMatrix = 
  {{188.76, 210.8, 228.64, 197.08, 196.68, 193.32, 176.24, 164.44, 150., 151.8}, 
   {210.8, 289.72, 282., 273.28, 265.64, 243.08, 245.76, 244.2, 205.4, 185.6}, 
   {228.64, 282., 325.4, 307., 281.76, 293.84, 279.16, 274.88, 231.04, 217.28}, 
   {197.08, 273.28, 307., 335.24, 317.24, 315.24, 316.2, 323.84, 276.8, 244.36}, 
   {196.68, 265.64, 281.76, 317.24, 349.8, 307.24, 315.32, 331.36, 306.04, 299.68}, 
   {193.32, 243.08, 293.84, 315.24, 307.24, 313.64, 302.92, 306.44, 274.96, 253.4}, 
   {176.24, 245.76, 279.16, 316.2, 315.32, 302.92, 307.04, 320.2, 281.48, 258.44}, 
   {164.44, 244.2, 274.88, 323.84, 331.36, 306.44, 320.2, 343.16, 298.88, 280.8}, 
   {150., 205.4, 231.04, 276.8, 306.04, 274.96, 281.48, 298.88, 279.48, 270.76}, 
   {151.8, 185.6, 217.28, 244.36, 299.68, 253.4, 258.44, 280.8, 270.76, 296.4}};

nDamsPerSire = 10;

sSires = 700;

rResidualDF  =  16000;

Preliminary calculations:  estimates of the mean squares

mE  =  dataMatrix  -  (1/4) gMatrix;

mA  =  (nDamsPerSire - 1) gMatrix / 4  +  dataMatrix;

Calculate V

Clear[index2Func];
index2Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexVal = 0;
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = indexVal;
              indexArg = indexArg+1,
              {i1, 1, k-i2}
             ];
           indexVal = indexVal+1,
           {i2, 0, k-1}
          ];
        vector
      )

Clear[index3Func];
index3Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = i1;
              indexArg = indexArg+1,
              {i1, i2, k}
             ],
           {i2, 1, k}
          ];
        vector )

Clear[index4Func];
index4Func[k_] := index2Func[k] + 1;

Clear[vMat, i, j, index3, index4];
     index3 = index3Func[nAges];
     index4 = index4Func[nAges];
     vMat = Table[0, {i, nAges(nAges+1)/2}, {j, nAges(nAges+1)/2}];
     Do[Do[
        covMaMa  =  ( mA[[ index3[[i]], index3[[j]] ]] * mA[[ index4[[i]], index4[[j]] ]] 
                     + mA[[ index3[[i]], index4[[j]] ]] * mA[[ index4[[i]], index3[[j]] ]] ) / (sSires + 2);
        covMeMe  =  ( mE[[ index3[[i]], index3[[j]] ]] * mE[[ index4[[i]], index4[[j]] ]] 
                     + mE[[ index3[[i]], index4[[j]] ]] * mE[[ index4[[i]], index3[[j]] ]] ) / (rResidualDF + 2);
        vMat[[i,j]] = vMat[[j,i]] = ( (nDamsPerSire - 1) / nDamsPerSire )^2 covMeMe  +  (1/nDamsPerSire)^2 covMaMa,
        
        {i, j, nAges(nAges+1)/2} ],
        {j, 1, nAges(nAges+1)/2}
       ]

INTERMEDIATE  CALCULATIONS

Function definitions  (based on Legendre polynomials)

Define the normalized Legendre polynomials, legPoly[j, x]

(Under Version 2 of Mathematica, the following line can be replaced with:

     legPoly = Compile[ {j, x}, Sqrt[(2j + 1)/2] LegendreP[j, x] ;
     
This will speed computation of the X matrix.)

legPoly[j_, x_] := Sqrt[(2j + 1)/2] LegendreP[j, x]//N ;

Define the index functions

Clear[index1Func];
index1Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexVal = 0;
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = indexVal;
              indexVal = indexVal+1;
              indexArg = indexArg+1,
              {i1, 1, k-i2}
             ];
           indexVal = 0,
           {i2, 0, k-1}
          ];
        vector )

Clear[index2Func];
index2Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexVal = 0;
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = indexVal;
              indexArg = indexArg+1,
              {i1, 1, k-i2}
             ];
           indexVal = indexVal+1,
           {i2, 0, k-1}
          ];
        vector
      )

Clear[index3Func];
index3Func[k_] := 
      ( Clear[junk, vector, i1, i2];
        vector = Array[junk, k(k+1)/2];
        indexArg = 1;
        Do[
           Do[
              vector[[indexArg]] = i1;
              indexArg = indexArg+1,
              {i1, i2, k}
             ],
           {i2, 1, k}
          ];
        vector )

Clear[index4Func];
index4Func[k_] := index2Func[k] + 1;

Clear[index5Func];
index5Func[k_] := index3Func[k] - 1;

Preliminary manipulations of the data

Calculate the vector a* of normalized ages

smallestAge = ageVector[[1]];
largestAge = ageVector[[nAges]];
ageRange = largestAge - smallestAge;

aStarVec = Array[0, nAges];

Do[
    aStarVec[[i]] = -1 + 2(ageVector[[i]] - smallestAge)/ageRange,
    {i, nAges}
  ]

Form the data vector from the data matrix:

Clear[i1, i2, i3];
mVec = Array[0, nAges(nAges + 1)/2];

i1 = 1;
Do[
    Do[
        mVec[[i1]] = dataMatrix[[i2, i3]];
        i1 = i1+1,
        {i2, i3, nAges}
      ],
    {i3, 1, nAges}
  ]

Calculate the X matrix from the orthogonal functions
(Legendre polynomials)

Clear[phiMat, xMat, i, j];
phiMat = Table[ legPoly[ j - 1, aStarVec[[i]] ] , {i, nAges}, {j, kPolys} ];
xMat = Table[ 0, {i, nAges(nAges+1)/2}, {j, kPolys(kPolys+1)/2}];

Clear[index2, index3, index4, index5, kronDelta];
index2 = index2Func[kPolys];
index3 = index3Func[nAges];
index4 = index4Func[nAges];
index5 = index5Func[kPolys];
kronDelta[i_, j_] := If[(i == j),
                         1,
                         0]

xMat = Table[ ( ( phiMat[[ index3[[i]], (index5[[j]] + 1) ]] *
                  phiMat[[ index4[[i]], (index2[[j]] + 1) ]] +
                  phiMat[[ index3[[i]], (index2[[j]] + 1) ]] *
                  phiMat[[ index4[[i]], (index5[[j]] + 1) ]] ) /
                ( 1.0 + kronDelta[ index2[[j]], index5[[j]] ] ) ),
              {i, nAges(nAges+1)/2},
              {j, kPolys(kPolys+1)/2} ];

RESULTS

The coefficients:  the c vector and C matrix

Calculations

Clear[cVec];

If[ (kPolys == nAges),

    (cVec = Inverse[xMat] . mVec;), (* Calculation for full fit *)
     
    (vInv = Inverse[vMat];        (* Calculations for reduced fit *)
     xMatTrans = Transpose[xMat];  
     cVec = Inverse[xMatTrans . vInv . xMat] . xMatTrans . vInv . mVec) ]

Clear[cMat, i, j, k];
cMat = Table[0.0, {i, kPolys}, {j, kPolys}];
k = 1;
Do[
   Do[
      cMat[[i,j]] = cVec[[k]];
      cMat[[j,i]] = cVec[[k]];
      k = k+1,
      {i, j, kPolys}
   ],
  {j, 1, kPolys} 
];

The result

cMat

The estimated covariance function

The ages t1 and t2 are in the original units in which age was measured.

Calculations

Clear[aStar, sum, i, j, covTilda, t1, t2];
aStar[a_] := -1  +  2 (a - ageVector[[1]]) /
                      (ageVector[[nAges]] - ageVector[[1]]);

sum = 0.0;
Do[
   Do[(sum = sum + cMat[[i,j]] * 
                   legPoly[ i-1, aStar[t1] ] *
                   legPoly [ j-1, aStar[t2] ]),
      {i, 1, kPolys}],
   {j, 1, kPolys}];
covTilda[t1_, t2_] := Release[Simplify[sum]];

Results

covTilda[t1, t2]

A smoothed estimate of the data matrix 
(with 3-D plots of the smoothed matrix and of the discrepancy between the smoothed matrix and the original data matrix)

This step is only useful with a reduced fit;  otherwise it simply re-plots the original data.

Calculations

Clear[smoothVec];
smoothVec = xMat . cVec;

Clear[smoothMat, i, j, k];
smoothMat = Table[0.0, {i, nAges}, {j, nAges}];
k = 1;
Do[
   Do[
      smoothMat[[i,j]] = smoothVec[[k]];
      smoothMat[[j,i]] = smoothVec[[k]];
      k = k+1,
      {i, j, nAges}
   ],
  {j, 1, nAges} 
];

The smoothed matrix

MatrixForm[smoothMat]

A 3-D plot of the smoothed matrix

ListPlot3D[smoothMat, Shading->False, PlotRange->All]

A 3-D plot of the discrepancy between the smoothed matrix and the  data matrix

ListPlot3D[dataMatrix - smoothMat, 
           Shading->False, PlotRange->All]

3-D plot of the continuous covariance and correlation functions

nPlotAges = 30;

Calculations

plotAgeVec = Table[ageVector[[1]] + 
                  (i - 1) (ageVector[[nAges]] - ageVector[[1]]) / 
                  (nPlotAges - 1),
                  {i, 1, nPlotAges}];
plotMat = Table[0, {i, nPlotAges}, {j, nPlotAges}];
Do[
   plotMat[[i,j]] = plotMat[[j,i]] =
      covTilda[plotAgeVec[[i]], plotAgeVec[[j]]],
   {i, nPlotAges}, {j, nPlotAges}];

rhoMat = Table[0, {i, nPlotAges}, {j, nPlotAges}];
Do[
   rhoMat[[i,j]] = rhoMat[[j,i]] =
      plotMat[[i,j]] / Sqrt[ plotMat[[i,i]] plotMat[[j,j]] ],
   {i, nPlotAges}, {j, nPlotAges}];

Plot of the covariance function

ListPlot3D[plotMat, PlotRange->{0, Max[plotMat]}, 
           Shading->False];

Plot of the corresponding correlation function

ListPlot3D[rhoMat, PlotRange->{0.5, 1}, 
           Shading->False];

Goodness-of-fit test of the smoothed estimate of the covariance function

This test is based on the method described in Appendix B of Kirkpatrick et al. (1990) (the full reference for which is given at the beginning of the notebook).  The resulting chi squared value is a test statistic for the goodness-of-fit of the smoothed covariance function resulting (from a reduced fit) to the original data matrix.  A significant value indicates that the smoothed fit is NOT consistent with the original data.

Calculations

If[ (kPolys < nAges),

    (i1 = 1;   (* We have a reduced fit *)
     Do[
         Do[
            smoothVec[[i1]] = smoothMat[[i2, i3]];
            i1 = i1+1,
            {i2, i3, nAges}
           ],
           {i3, 1, nAges}
       ];
     chiSq = Chop[(mVec - smoothVec) . vInv . (mVec - smoothVec)] ),
    
    chiSq = 0  (* We have a full fit *)
  ];

Result

The following chi squared value is to be compared against a standard table of critical values:

chiSq

The relevant degrees of freedom for the test are:

df = nAges(nAges + 1)/2  -  kPolys(kPolys + 1)/2

Eigenfunctions and eigenvalues of the covariance function

These are calculated analytically from the coefficient matrix.  See Kirkpatrick et al. 1990.

Calculations

eigens = Eigensystem[cMat];
eVal = eigens[[1]]//N;
eVecsTemp = eigens[[2]]//N;

Sort the eigenvalues by size, and re-order the eigenvectors:

lambda = Sort[eVal, Greater];
eVec = Table[0, {i, kPolys}];
Do[ Do[ If[ eVal[[j]] == lambda[[i]],
            eVec[[i]] = eVecsTemp[[j]] ],
        {j, 1, kPolys} 
      ],
    {i, 1, kPolys} ]

Normalize the eigenvectors (which some versions of Mathematica  don't do automatically):

Do[ eVec[[i]] = eVec[[i]] / 
                Sqrt[ Sum[eVec[[i,j]]^2, {j, kPolys}] ],
    {i, kPolys} ]

Chop[eVec]

Compute the first three eigenfunctions:
   NOTE that the following is based on the range of the Legendre polynomials, [-1, 1].  If orthogonal functions defined over a different range are used, then the value "2" in the square roots must be replaced by (Maximum argument  -  Minimum argument) for those functions.

Clear[eigenFunc1, eigenFunc2, eigenFunc3, i, t];
eigenFunc1[t_] := Release[ Simplify[Sqrt[2.0/ageRange] *
                   Sum[ eVec[[1,i]] legPoly[i-1, aStar[t]],
                        {i, 1, kPolys} ] ] ]
eigenFunc2[t_] := Release[ Simplify[Sqrt[2.0/ageRange] *
                   Sum[ eVec[[2,i]] legPoly[i-1, aStar[t]],
                        {i, 1, kPolys} ] ] ]
eigenFunc3[t_] := Release[ Simplify[Sqrt[2.0/ageRange] *
                   Sum[ eVec[[3,i]] legPoly[i-1, aStar[t]],
                        {i, 1, kPolys} ] ] ]

Renormalize the eigenvalues to accomodate the range of the coordinate vector:
   NOTE that the following is based on the range of the Legendre polynomials, [-1, 1].  If orthogonal functions defined over a different range are used, then the value "2" must be replaced by (Maximum argument  -  Minimum argument) for those functions.

lambda = lambda ageRange / 2.0;

eigenSum = Sum[lambda[[i]], {i, kPolys}];

Calculate the "squeezed" eigenvalues:

lambdaMean = Sum[ lambda[[i]], {i, kPolys} ] / kPolys

squeezeFactor = lambdaMean / (lambdaMean - lambda[[kPolys]])

squeezedLambda = Chop[ Table[lambdaMean + (lambda[[i]]- lambdaMean) squeezeFactor,
                       {i, kPolys}] ];

squeezedEigenSum = Sum[squeezedLambda[[i]], {i, kPolys}];

Results:  The eigenvalues

NOTE that the eigenvalues and eigenvectors are normalized here using a criterion different than that in Kirkpatrick et al. (1990) and Kirkpatrick and Lofsvold (1992).  (The discrepancy involves a multiplicative constant.  In this package the eigenfunctions have norm unity when integrated over the range of the coordinate vector.  In the earlier publications, the integration was over the range of the orthogonal functions.)

The eigenvalues are:

lambda

The fraction of all variation that these eigenvalues account for:

varFrac = Chop[ Table[lambda[[i]] / eigenSum, {i,kPolys}] ]

A plot of the eigenvalues:

ListPlot[lambda, 
         PlotJoined->True, AxesOrigin->{0,0}]

Results:  The "squeezed" eigenvalues

To correct for biases that inflate the range of estimate eigenvalues, the eigenvalue range can be "squeezed".  See Kirkpatrick and Lofsvold (1992) Evolution.

The "squeezed" eigenvalues:

squeezedLambda

The fraction of all variation that these eigenvalues account for:

squeezedVarFrac = Chop[ Table[squeezedLambda[[i]] / squeezedEigenSum,
                        {i,kPolys}] ]

eigenSum

Results:  The first eigenfunction

The first eigenvalue is:

lambda[[1]]

Its associated eigenfunction is:

eigenFunc1[t]

Plot of the first eigenfunction:

Plot[ eigenFunc1[t], 
      {t, ageVector[[1]], ageVector[[nAges]]} ]

Results:  The second eigenfunction

The first eigenvalue is:

lambda[[2]]

Its associated eigenfunctions is:

eigenFunc2[t]

Plot of the first eigenfunction:

Plot[ eigenFunc2[t], 
      {t, ageVector[[1]], ageVector[[nAges]]} ]

Results:  The third eigenfunction

The first eigenvalue is:

lambda[[3]]

Its associated eigenfunctions is:

eigenFunc3[t]

Plot of the first eigenfunction:

Plot[ eigenFunc3[t], 
      {t, ageVector[[1]], ageVector[[nAges]]} ]