使用blupf90运行测定日模型

使用blupf90运行测定日模型

测定日模型

重复力模型的一个基本假设是重复度量的记录值之间的遗传相关是 1 ，但是重复力模型用在奶牛的产奶性状的缺点在于，奶牛的产奶性状在整个产奶周期其均值和方差不断发生变化，也有文献指出产奶量的遗传力随着产奶日期而变化，并且重复度量的性状之间的遗传相关会因为两次度量的时间间隔增长而趋于降低，而重复力模型无法解释这种随着时间变化的协方差结构。

多性状模型假设不同性状的遗传机理不同但是彼此相关，因此当需要考虑不同记录值的遗传相关时，多性状模型相比于重复力模型更佳。但是当记录值的重复次数很多时，此时使用多性状模型需要估计的参数过多，难以得到准确的遗传参数。

分析这种随着时间或年龄变化的重复度量性状，其合适的模型应该能够解释其均值和协方差结构随着时间变化，并且可以做到有效估计需要的遗传参数。在1994年，Schaeffer 和 Dekkers 引入了随机回归(RR)模型的概念，用于分析奶牛的测定日数据。同时， Kirkpatrick et al. (1990, 1994) 引入了协方差函数 (covariance functions, CFs) 来实现纵向数据的分析。Schaeffer (2004) 概述了RR模型在动物育种中各类数据中的应用情况。在RR模型发展之前， Ptak and Schaefer (1993) 通过使用固定回归模型（fixed regression model）来分析产奶量性状。

勒让德多项式详解

看不懂。。。

勒让德多项式（Legendre polynomials）是勒让德微分方程在特定条件下的多项式解，其标准形式为：

当参数 $ n $ 为非负整数时，方程在区间 内存在有界解，这些解构成一组正交多项式序列。

其显式表达式可通过**罗德里格公式（Rodrigues’ formula）**表示：

此公式直接生成 $ n $ 阶多项式，例如 $ P_0(x) = 1 P_1(x) = x P_2(x) = \frac{1}{2}(3x^2 - 1) $ 等。

其递推公式如下

此关系用于高效计算高阶多项式，例如由 $ P_0 $ 和 $ P_1 $ 推导 $ P_2 $。
因此，我们可以通过递推公式得到下面的表达式

实际计算公式

参考《Linear models for the prediction of animal breeding values》附录G，写出R代码如下（这里代码中命名有点问题，归一化后的 DIM 不能叫做 HTD ，这是两个概念）

### Legendre covariable
Legendre <-function(x,Xmax,Xmin){
  HTD<-(2*(x-Xmin)/(Xmax-Xmin))-1
  mat<-c(0.7071,0,0,0,0)
  mat[2]<-((3/2)**0.5)*HTD
  mat[3]<-((45/8)**0.5)*HTD^2-(5/8)^0.5
  mat[4]<-((175/8)**0.5)*HTD^3-(63/8)^0.5*HTD
  mat[5]<-9.2808*HTD^4-7.9550*HTD^2+0.7955
  return(mat)
}

具体解释如下

1. 变量归一化处理

代码中的 HTD = (2*(x - Xmin)/(Xmax - Xmin)) - 1 将变量从区间 [Xmin, Xmax] 线性映射到 [-1, 1]。这是Legendre多项式的标准定义域。

---

2. 归一化Legendre多项式公式

Legendre多项式通常通过 罗德里格公式（Rodrigues’ Formula） 定义：

但在实际应用中，常对其进行归一化以满足正交性条件（大概理解了，参考勒让德多项式的正交性和归一化, 这个公式的意义就是使得勒让德多项式的模为1）：

代码中的系数与归一化后的形式一致。

---

3. 代码中各阶多项式解析

(1) 0阶项 mat[1]

• 代码值：0.7071
• 对应公式：

出处：归一化条件。

(2) 1阶项 mat[2]

• 代码值：sqrt(3/2) * HTD
• 对应公式：

出处：Legendre多项式的一次项展开。

(3) 2阶项 mat[3]

• 代码值：sqrt(45/8) * HTD^2 - sqrt(5/8)
• 对应公式：

出处：二次项展开及归一化。

(4) 3阶项 mat[4]

• 代码值：sqrt(175/8) * HTD^3 - sqrt(63/8) * HTD
• 对应公式：

出处：三次项展开及归一化。

(5) 4阶项 mat[5]

• 代码值：9.2808 * HTD^4 - 7.9550 * HTD^2 + 0.7955
• 对应公式：

出处：四次项展开及数值近似。

固定回归模型

固定回归模型的理论框架见文献 Ptak and Schaefer (1993) 。模型中直接使用了测定日的记录值，因此避免了直接将测定日的奶产量直接扩展为 305天产奶量（我的理解是是不是之前的分析都是先计算出305天产奶量，再将305天产奶量作为性状），并且在模型中添加了测定日效应（herd-test-day, HTD）用于解释所有母牛在某一天受到的影响。这个模型还能解释在相同年龄、相同季节、相同地区的母牛的一般泌乳曲线。

和重复力模型相似，固定回归模型假定在一个泌乳周期的测定日的记录值是同一个性状的重复度量值，即测定日记录值之间的遗传相关为1。通常模型中还会添加永久环境效应，模型形式如下

其中 是在测定日 ， HTD子集 ，母牛 的记录值； 是固定回归系数； 和 是母牛 的加性效应和永久环境效应， 是对于母牛 在测定日 的测定日记录的 第 个勒让德多项式（Legendre polynomials o）， 是勒让德多项式用于拟合固定泌乳曲线的阶数， 是随机残差。

用矩阵形式，模型可以写作如下

其中 是测定日产奶量的向量， 是 HTD 和 固定回归系数的解向量（看来这里勒让德多项式应该是作为协变量）， 和 是加性效应和永久环境效应的向量。假设 ， ， 。

混合模型方程组如下（这里就是很正常的公式）

举个例子，假设我们有 5头母牛的测定日的乳脂量的数据，具有HTD信息和乳脂记录值，我们的目的是估计 HTD效应的子，4阶勒让德多项式的回归系数，永久环境效应和加性效应的子。

表型数据如下，这里有2列可以作为测定时间变量的列，第一个 DIM （Days in Milk，泌乳天数）指奶牛自产犊后开始泌乳以来经过的天数，是个体层次的连续变量。例如，一头奶牛产犊后第50天，其DIM为50。第二个 HTD（Herd–Test–Day，牛群-测试-日） ，指在特定日期对某一牛群进行生产性能测试的记录，是群体层次的分类变量。例如：“牛群A在2023-10-05的测试日数据”，一般可能作为固定效应，消除不同测试日（如季节、管理差异）对结果的干扰。

但是在这个案例中，通过附录G可以看到，这里是对归一化的DIM进行计算得到了 Legendre多项式，这里其实是计算得到了一个矩阵 ，是一个 10 × 5 的矩阵，每一行就是 DIM 从低到高的10个水平，每一列就是 0阶到 4阶的 Legendre多项式 。

在 blupf90 官方文档中，运行固定效应如下，首先表型数据如下

各列内容说明如下

1. 个体号
2. 泌乳天数（DIM）
3. 牛群-测试-日效应（HTD）
4. 乳脂量
5. 0阶Legendre多项式（由DIM计算得到）
6. 1阶Legendre多项式
7. 2阶Legendre多项式
8. 3阶Legendre多项式
9. 4阶Legendre多项式

参数卡内容如下，可以看到和常规模型差不多，只是将 Legendre多项式 作为协变量（固定因子）加入到模型中。

求解得到解向量之后，我们可以通过多项式的回归系数来画出泌乳曲线。在实际计算中，固定回归的协变量会嵌套与某个环境变量中，这样不同群体的奶牛可以画出不同的泌乳曲线，因此我们可以检验不同环境因素对泌乳曲线的影响。但是在我们的例子中，我们只能画出一条适用于所有母牛的一条泌乳曲线。一个泌乳天数从4到310天的实际乳脂量的向量 计算如下

这里 是一个 10 × 5 的矩阵，每一行就是 DIM 从低到高的10个水平，每一列就是 0阶到 4阶的 Legendre多项式 。 应该就是5个估计出的协变量的值组成的向量。二者的矩阵乘积就是计算出 DIM 每一个水平的回归项之和。

最终计算出每个水平的实际乳脂量如下（标题顺序乱了，凑合看吧，就是这么个意思）

因此我们可以基于这个结果画出一条泌乳曲线。

这里我们求解得到的 EBV 是一天的乳脂量，由于我们假设测定日记录值之间的遗传相关为1，因此305天产奶量就是这里的 EBV值乘以 305 。

随机回归模型

1. 背景：固定回归模型的局限性

问题1：仅描述群体平均泌乳曲线

• 固定回归模型：在模型中引入泌乳天数（DIM）作为固定效应，可以拟合不同牛群的平均泌乳曲线（例如高产群 vs. 低产群的整体趋势）。
• 局限性：
  • 它假设所有奶牛的泌乳曲线形状完全由群体决定，忽略了个体间遗传差异对曲线形状（如高峰期、下降速度）的影响。
  • 因此，估计的**育种值（遗传优势）**只能反映奶牛在泌乳高峰期（曲线顶点）的遗传差异，而无法捕捉整个泌乳周期的遗传特性（如持久性）。
  • 不同泌乳天数的性状之间的遗传力假设为1，即所有测定日的性状均是同一个性状，这可能不符合实际情况。

问题2：忽略遗传协方差结构

• 固定回归模型可以调整不同泌乳阶段的残差异质性（例如泌乳初期误差大、后期误差小），但未考虑遗传效应随泌乳阶段的变化。
• 例如，一头奶牛可能在泌乳初期遗传优势显著，但另一头在后期表现更优，这种动态差异无法通过固定效应建模。

2. 解决方案：随机回归模型（RR Model）

Schaeffer and Dekkers (1994) 将固定回归模型推广位随机回归模型，将这些回归系数从固定因子改为随机因子，因此就可以看出不同动物个体的泌乳曲线的形状的差别，并且为评估泌乳持久性（即高峰期后产奶量维持能力）的遗传潜力提供了可能性。

核心思想

• 将泌乳天数（DIM）的回归系数设为随机效应，而非固定效应。
• 具体操作：
  • 固定效应部分：仍保留群体平均泌乳曲线（例如不同品种、管理组的基线）。
  • 随机效应部分：允许每头奶牛有自己的泌乳曲线（通过随机回归系数表示），这些曲线可以偏离群体平均曲线，且偏离程度受遗传控制。

数学工具实现

• Legendre多项式：最常用，因其无需预设曲线形状，且计算方便（正交性降低多重共线性）。
• 自然三次样条（Non-parametric）：更灵活，但计算复杂。
• 其他参数化方法：如随机参数曲线（例如Wood模型）。

3. 模型

经典的随机回归模型如下

其中， 是奶牛 在泌乳天数 ，HTD水平 的测定记录值， 是固定回归系数， 和 是奶牛 加性效应和永久环境效应的第 个随机回归系数。 是对于奶牛 在在泌乳天数 的 阶Legendre多项式的值。 是固定回归使用的多项式的阶数， 是加性效应和永久环境效应使用的多项式的阶数。 是残差。

写成矩阵形式如下

这里 和固定回归模型相同，但是此时 和 是加性效应和永久环境效应的随机回归的向量，而 和 是协变量矩阵，每一行的内容是测定日表型该行的DIM的Legendre多项式。如果加性效应和永久环境效应使用的阶数相同，此时 。一般来说， 和 的维度是 ，这里 就是表型行数， 等于 乘以表型中的动物个体数目（因为每一个个体都有 个多项式的随机回归系数需要估计）。

我们假设 ， ， （这里应该是按照因子排序）。这里 是亲缘关系矩阵， 和 是加性和永久环境效应拟合的多项式的方差组分矩阵。

构建 MME 如下，我感觉就是和母体效应模型差不多。

4. 举例

还是使用固定回归模型例子中的数据，此时使用随机回归模型，采用4阶多项式拟合固定泌乳曲线，采用2阶多项式拟合加性和永久环境效应（包括截距有3个回归系数）。此时每一个个体均有3个加性效应的随机系数，也有3个永久环境效应的回归系数。

加性和永久环境效应的协方差矩阵如下（类似于母体效应模型），残差设定所有泌乳阶段均为 3.710 。

这里的 和 是整个泌乳周期的协方差结构，因此 DIM 和 的遗传协方差可以通过 矩阵来计算。

举个例子，DIM 的遗传方差 （ ）可以按照下面的公式计算（因为育种值计算公式是 ， 我们有公式 , ）

其中 ，就是 的第 行（对应 DIM ）， 是阶数 。

DIM 和 DIM 的遗传协方差 （ ）的计算公式如下

采用例子中的 矩阵，我们可以得到 DIM 106 等于 2.6433 kg2 （ ）。DIM 106 和 DIM 140 的加性协方差为 3.0219kg 。每个 DIM 的加性和永久环境效应方差画图如下，用来指示整个泌乳周期遗传参数的变化情况。

对于所有有表型的个体， (因为 和 矩阵相同)，并且均是分块对角矩阵（因子顺序应该就是和母体效应模型一样，按照方差组分的顺序排列，以加性为例，就是个体1的0阶多项式，个体1的1阶多项式，个体1的2阶多项式，个体2的0阶多项式…）。举个例子， 的表型中的前3头母牛（4,5,6）的内容如下

使用 blupf90 的参数卡如下（这里有一点和理论对不上，加性效应是 8个水平没错，理论上永久环境效应应该是 5个水平，因为表型中只有 5 个个体，而这里永久环境效应也是 8 个水平。我知道这里设置随机效应后面多了一个数字，例如 5 8 cov 1 ，这个 1 应该是被嵌套的因子，也就是这个因子是嵌套在第一列个体号中，因为这里加性效应和永久环境效应都是用的第一列个体号，所以水平数目相同）。

这里加性效应是第 7，8，9个效应，而永久环境效应定义为第10，11，12个效应。

这里我们看第 7 个效应

5 8 cov 1

这个含义是将第 5 列作为一个回归的截距，因为我们总共有 8 个个体，所以我们需要估计 8 个截距。上面字段的内容是说，我们将第5列视为一个协变量并且嵌套于第一列个体号中（我大概知道怎么计算嵌套回归了，之后再看看DMU的嵌套回归）。以加性效应为例，由于我们这里是3个因子，因此其方差组分部分需要一个3 × 3 的矩阵。

我们得到解之后，这里固定因子部分和固定回归模型类似，我们同样可以通过固定回归系数画出一般的泌乳曲线。对于随机因子，如加性效应，每个个体都可以得到 个回归系数（这里就是3个回归系数），但是这个结果不能直接用于对动物个体进行排序，我们需要将其转化为某一个特定的泌乳天的育种值。我们通常会计算奶牛305天每一天的育种值，举个例子，个体 在泌乳天从第6天到第 天的育种值（我感觉应该是说从第6天到第m天的育种值之和）计算如下

这里 是一个维度是 的行向量，其第 个元素等于从第6天到第 天的第 个多项式的值（说白了， 就是 中从第6天到第 天的行的和）；而 是个体 的回归系数的估计值的列向量。

假设我们需要计算从第 6天到第 310天的 305天育种值（额，这不是306天吗），那么可以得到其 如下

举个例子，个体4的 305天乳脂量育种值为

同理，我们可以计算出每一个个体每一天的育种值，因此我们可以画出每一个个体的遗传泌乳曲线（横坐标是泌乳天数，纵坐标是当天的育种值）。

假设 是包含从第6天到第310天的育种值向量，那么 的计算公式如下 ( 矩阵就是 中从第6天到第 310 天的行组成的矩阵)

如果性状是产奶量，持久性育种值可以通过每天的育种值获得，举个例子，产奶量的持久性性状的预测传递力（）可以通过下面的公式计算（Schaeffer 2000）

这里 和 是某个个体在 280天和 60天的预测传递力（不清楚咋算）， 和 是280天和60天遗传水平的平均产奶量（还是不清楚咋算）。

育种值的可靠性

育种值可靠性决定于于加性方差的 PEV ，可以理解为计算 EBV 时的可用信息的丰富程度。

我们设某个个体 的育种值定义为 ，这里 。这里 可能是对于第 个年龄或胎次的权重因子（如果不同胎次的性状视为不同的性状），举个例子，如果第1胎次和第2胎次的乳脂量视为不同的性状， 可能就是 ，表示两个胎次的权重分别是0.7和0.3。 这里 向量定义为单个胎次内的育种值如何计算（公式同上）。这里加性和永久环境效应的协方差矩阵分别为 和 ，那么此时的加性方差为 ，而永久环境效应的方差为 ，而加性效应 的遗传力可以按照以下公式计算 。

假设个体 的加性效应在 MME 逆矩阵的子集为 ，那么对于个体 ，其预测误差方差 ，因此 的可靠性为 。

按照上面的例子， 是标量 1，此时

对于个体1，其矩阵 为

此时，计算 PEV 如下

因此个体1的可靠性等于 。

下面可以证明，对于计算育种值是所有天数育种值的和，和所有天数育种值的均值 ，这两种情况下的可靠性相同。证明如下，此时所有天数育种值的和的 其实就是所有天数育种值的均值的 的 倍（ 为天数总数），则有

5. blupf90 实操

官方例子

随机回归模型-单性状模型

表型数据 datarrm 各列如下，这里最后三列应该是多项式

ID011  A  2  9.2   40.87  29  1  0.98 5 -.19245 -1.07663 .43189
ID011  A  1 10.3   39.13  30  1  0.98 6 .19245 -1.07663 -.43189
ID011  A  5 11.4   41.20  31  2  0.98 7 .57735 -.74536 -1.07790

单形状的随机回归模型的 renumf90 的参数卡如下，这里分析一下

固定因子：

• 第一个固定因子： 第2列，分类变量。
• 第一个协变量：第8列， 嵌套在第3列分类因子中的协变量。

随机因子：

• 加性效应的随机回归变量（阶数2）
• 永久环境效应的随机回归变量（阶数2）：这里使用的写法是 OPTIONAL pe ，而没有单独将永久环境效应提取出来。

# parameter file for renumf90
# animal model with 1 trait
DATAFILE
 datarrm
TRAITS
 4
FIELDS_PASSED TO OUTPUT
 
WEIGHT(S)
 
RESIDUAL_VARIANCE
  1.0
EFFECT
 2 cross alpha
#
# fixed regression:
#  covariates (position 8) nested within the class (position 3)
#
EFFECT
 8 cov
NESTED
 3 alpha
EFFECT
 1 cross alpha
RANDOM
 animal 
OPTIONAL
 pe
FILE
 ped1 
FILE_POS
 1 2 3 0 0 
PED_DEPTH
 3
RANDOM_REGRESSION
 data
RR_POSITION
 10 11 12
#
# Order:
#          a0 a1 a2
#       a0  *  *  *
#       a1  *  *  *
#       a2  *  *  *
#
(CO)VARIANCES
 0.80325  0.0702  0.0351
 0.0702   0.5751  0.1404
 0.0351   0.0702  0.6102
(CO)VARIANCES_PE
 0.3575  0.042  0.021
 0.042   0.221  0.084
 0.021   0.042  0.242

运行一下

renumf90 renum.rrm.1.txt > renumf90.log

blupf90 renf90.par > blupf90.log

输出参数卡 renf90.par 如下，可以查看 renf90.fields 可以到此时的 renf90.dat 各列内容分别位性状，L0多项式，L1多项式，L2多项式，固定因子1，协变量1，协变量1被嵌套的分类变量（原来的第三列），多项式被嵌套的分类变量（个体号）。

注意这里加性效应和永久环境效应被嵌套的列是同一列，因此水平数目相同（这和上面的官方文档一致，按照理论永久环境效应的水平数目要小于加性效应）。

运行 blupf90 报错 G Matrix not symmetric ，看来给的这个例子不能实际运行

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           8
OBSERVATION(S)
    1
WEIGHT(S)
 
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 5         2 cross 
 6         6 cov  7
 2         15 cov 8
 3         15 cov 8
 4         15 cov 8
 2         15 cov 8
 3         15 cov 8
 4         15 cov 8
RANDOM_RESIDUAL VALUES
   1.0000    
 RANDOM_GROUP
     3     4     5
 RANDOM_TYPE
 add_animal   
 FILE
renadd03.ped                                                                    
(CO)VARIANCES
  0.80325      0.70200E-01  0.35100E-01
  0.70200E-01  0.57510      0.70200E-01
  0.35100E-01  0.14040      0.61020    
 RANDOM_GROUP
     6     7     8
 RANDOM_TYPE
 diagonal     
 FILE
                                                                                
(CO)VARIANCES
  0.35750      0.42000E-01  0.21000E-01
  0.42000E-01  0.22100      0.42000E-01
  0.21000E-01  0.84000E-01  0.24200    

随机回归模型-多性状模型

多性状模型的 renumf90 的参数卡如下，这里基本和常规的单性状模型改为多性状模型一致，只需要注意一下方差组分的元素顺序（这里和母体效应模型一致，均为效应在外面，性状在里面）。

# parameter file for renumf90
# animal model with 4 traits
DATAFILE
 datarrm
TRAITS
 4 5 6 7
FIELDS_PASSED TO OUTPUT
 
WEIGHT(S)
 
RESIDUAL_VARIANCE
  1.0  0.6 -0.4 -0.1
  0.6  2.1  0.7  0.5
 -0.4  0.7  1.5  0.8
 -0.1  0.5  0.8  2.0
EFFECT
 2 2 2 2 cross alpha
#
# fixed regression:
#  covariates (position 8) nested within the class (position 3)
#
EFFECT
 8 8 8 8 cov
NESTED
 3 3 3 3 alpha
EFFECT
 1 1 1 1 cross alpha
RANDOM
 animal 
OPTIONAL
 pe
FILE
 ped1 
FILE_POS
 1 2 3 0 0 
PED_DEPTH
 3
RANDOM_REGRESSION
 data
RR_POSITION
 10 11 12
#
# Order:
#          a0          a1          a2
#          t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t3 t4
#    a0 t1  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *
#    a1 t1  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *
#    a2 t1  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *
#
(CO)VARIANCES
  0.80325   0.15795   -0.1053    -0.026325   0.0702   0.04212  -0.02808  -0.00702   0.0351   0.02106  -0.01404  -0.00351
  0.15795   1.09282    0.184275   0.131625   0.04212  0.14742   0.04914   0.0351    0.02106  0.07371   0.02457   0.01755
 -0.1053    0.184275   0.934875   0.2106    -0.02808  0.04914   0.1053    0.05616  -0.01404  0.02457   0.05265   0.02808
 -0.026325  0.131625   0.2106     1.0665    -0.00702  0.0351    0.05616   0.1404   -0.00351  0.01755   0.02808   0.0702 
  0.0702    0.04212   -0.02808   -0.00702    0.5751   0.02106  -0.01404  -0.00351   0.1404   0.08424  -0.05616  -0.01404
  0.04212   0.14742    0.04914    0.0351     0.02106  0.61371   0.02457   0.01755   0.08424  0.29484   0.09828   0.0702 
 -0.02808   0.04914    0.1053     0.05616   -0.01404  0.02457   0.59265   0.02808  -0.05616  0.09828   0.2106    0.11232
 -0.00702   0.0351     0.05616    0.1404    -0.00351  0.01755   0.02808   0.6102   -0.01404  0.0702    0.11232   0.2808 
  0.0351    0.02106   -0.01404   -0.00351    0.0702   0.04212  -0.02808  -0.00702   0.6102   0.04212  -0.02808  -0.00702
  0.02106   0.07371    0.02457    0.01755    0.04212  0.14742   0.04914   0.0351    0.04212  0.68742   0.04914   0.0351 
 -0.01404   0.02457    0.05265    0.02808   -0.02808  0.04914   0.1053    0.05616  -0.02808  0.04914   0.6453    0.05616
 -0.00351   0.01755    0.02808    0.0702    -0.00702  0.0351    0.05616   0.1404   -0.00702  0.0351    0.05616   0.6804 
(CO)VARIANCES_PE
  0.3575   0.0945   -0.063    -0.01575   0.042   0.0252  -0.0168  -0.0042   0.021   0.0126  -0.0084  -0.0021
  0.0945   0.53075   0.11025   0.07875   0.0252  0.0882   0.0294   0.021    0.0126  0.0441   0.0147   0.0105
 -0.063    0.11025   0.43625   0.126    -0.0168  0.0294   0.063    0.0336  -0.0084  0.0147   0.0315   0.0168
 -0.01575  0.07875   0.126     0.515    -0.0042  0.021    0.0336   0.084   -0.0021  0.0105   0.0168   0.042 
  0.042    0.0252   -0.0168   -0.0042    0.221   0.0126  -0.0084  -0.0021   0.084   0.0504  -0.0336  -0.0084
  0.0252   0.0882    0.0294    0.021     0.0126  0.2441   0.0147   0.0105   0.0504  0.1764   0.0588   0.042 
 -0.0168   0.0294    0.063     0.0336   -0.0084  0.0147   0.2315   0.0168  -0.0336  0.0588   0.126    0.0672
 -0.0042   0.021     0.0336    0.084    -0.0021  0.0105   0.0168   0.242   -0.0084  0.042    0.0672   0.168 
  0.021    0.0126   -0.0084   -0.0021    0.042   0.0252  -0.0168  -0.0042   0.242   0.0252  -0.0168  -0.0042
  0.0126   0.0441    0.0147    0.0105    0.0252  0.0882   0.0294   0.021    0.0252  0.2882   0.0294   0.021 
 -0.0084   0.0147    0.0315    0.0168   -0.0168  0.0294   0.063    0.0336  -0.0168  0.0294   0.263    0.0336
 -0.0021   0.0105    0.0168    0.042    -0.0042  0.021    0.0336   0.084   -0.0042  0.021    0.0336   0.284 

多泌乳期随机回归测定日模型-单性状模型

这里指分析奶牛（或其他哺乳动物）的多个泌乳周期（如第1次、第2次泌乳等），研究不同泌乳阶段的生产性能。这里应该是将不同的胎次的表型视为不同的性状来分析。

表型数据如下，这里的案例应该是 2 个胎次的数据。

ID011  A  2  9.2   40.87  29  1  0.98 1  5 -.19245 -1.07663 .43189   0 0 0
ID011  A  1 10.3   39.13  30  1  0.98 2  6 0 0 0 .19245 -1.07663 -.43189
ID011  A  5 11.4   41.20  31  2  0.98 2  7 0 0 0 .57735 -.74536 -1.07790

参数卡如下，这里分析如下

固定因子：

• 固定因子1 ： 第2列
• 协变量1：第8列，嵌套在第3列中

随机因子：

• 加性效应的随机回归变量（阶数2）

这里是用了 6 个多项式，前3个多项式是第1胎次的3个多项式，后3个多项式是第2胎次的3个多项式（从表型数据来看，应该是如果是第1胎次，那么第2胎次的3个多项式就会设为0，）。

• 永久环境效应的随机回归变量（阶数2）

# parameter file for renumf90
# animal model with 1 traits
DATAFILE
 datamltdm
TRAITS
 4
FIELDS_PASSED TO OUTPUT
 
WEIGHT(S)
 
RESIDUAL_VARIANCE
  1.0
EFFECT
 2 cross alpha
#
# fixed regression:
#  covariates (position 8) nested within the class (position 3)
#
EFFECT
 8 cov
NESTED
 3 alpha
EFFECT
 1 cross alpha
RANDOM
 animal 
OPTIONAL
 pe
FILE
 ped1 
FILE_POS
 1 2 3 0 0 
PED_DEPTH
 3
RANDOM_REGRESSION
 data
# First lactation  = 11 12 13
# Second lactation = 14 15 16
RR_POSITION
 11 12 13 14 15 16
#
# Order:
#          First    Second
#          11 12 13 14 15 16
#          a0 a1 a2 a1 a2 a3
#       a0  *  *  *  *  *  *
#       a1  *  *  *  *  *  *
#       a2  *  *  *  *  *  *
#       a0  *  *  *  *  *  *
#       a1  *  *  *  *  *  *
#       a2  *  *  *  *  *  *
#
(CO)VARIANCES
 0.80325   0.0702    0.0351    0.200813  0.01755   0.008775
 0.0702    0.5751    0.1404    0.01755   0.143775  0.0351
 0.0351    0.0702    0.6102    0.008775  0.01755   0.15255
 0.200813  0.01755   0.008775  0.80325   0.0702    0.0351
 0.01755   0.143775  0.0351    0.0702    0.5751    0.1404
 0.008775  0.01755   0.15255   0.0351    0.0702    0.6102
(CO)VARIANCES_PE
 0.3575    0.042    0.021    0.089375  0.0105   0.00525
 0.042     0.221    0.084    0.0105    0.05525  0.021
 0.021     0.042    0.242    0.00525   0.0105   0.0605
 0.089375  0.0105   0.00525  0.3575    0.042    0.021
 0.0105    0.05525  0.021    0.042     0.221    0.084
 0.00525   0.0105   0.0605   0.021     0.042    0.242

多泌乳期随机回归测定日模型-多性状模型

同样地，这里也只需要注意提供的方差组分的元素的位置，参数卡中这部分内容的说明如下

#
# Order:
#          First lactation                     Second lactation
#          11          12          13          14          15          16
#          a0          a1          a2          a0          a1          a2
#          t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t3 t4
#    a0 t1  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#    a1 t1  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#    a2 t1  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#    a0 t1  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#    a1 t1  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#    a2 t1  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t2  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t3  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#       t4  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
#

参考文献

1. Linear models for the prediction of animal breeding values[M]. Cabi, 2014.