使用blupf90运行阈值模型

阈值模型理论

有些性状属于有序多分类的性状，例如产犊难易程度性状，产仔数等。这些性状本身不服从正态分布，但是育种家们通常会将这类性状的表型值归因于某个潜在的无法观测到的服从正态分布的连续性状，称为liability （有的人称为易感性，我下面还是称呼为潜在变量）。因此观测到的分类表型可以解释为其潜在变量超过了某个特定阈值后的表型。因此，对于有 个分类的分类变量，这里总共有 个阈值 () 。

阈值模型需要用到正态分布的一些函数，假设我们的性状是产羔数，其有4个分类，此时潜在变量具有3个阈值的分布如下图所示，其中 表示当潜在变量超过 之后的表示产羔数为 的表型。

假设潜在变量 服从 的分布，其正态分布曲线在点 的高度为 为

这里我们用 表示正态分布的累积分布函数，因此 或 会计算出正态分布曲线左侧截止到第 个分类的面积。假设我们有 个分类，因此当 时， 。

我们定义 为表型观测值为第 个分类的概率，因此 或 可以计算为 ，其中 。或者我们用阈值模型的表示方式，每个阈值用 表示，此时 。

我们对潜在变量的分析模型如下

其中 是潜在变量的向量，其它和正常的模型一样。

因为 是没有观测到的，因此我们不可能通过构建MME来求解。

假设 ， Gianola 和 Foulley (1983) 给出了最大化后验分布 的 log 函数的估计值 的式子，其中用到了一阶导和二阶导，公式如下（看不太懂，但是这个方法只有在给定方差组分时进行遗传评估的方法，没有估计方差组分的内容）

其中 ，或 ，其中 或 是父性效应方差或加性效应方差。

很多式子不是很清楚，就先这样吧。

blupf90 实操

吉布斯抽样

在blupf90官方文档的估计方差组分章节中有吉布斯抽样的内容。吉布斯抽样（Gibbs Sampling）与REML方法完全不同。此处我们不探讨该技术的理论背景，细节可参考Sorensen和Gianola（2002）或Misztal（2008）的著作。简而言之，吉布斯抽样是一种基于当前可用信息，通过随机数生成样本，从中提取后验分布的迭代过程。在单轮抽样中，吉布斯采样器会执行以下操作：

1. 参数求解：利用当前方差组分（variance components）求解混合模型方程组；
2. 添加噪声：为每个解的估计值添加一个微小随机数（噪声）；
3. 协方差估计：基于更新后的解，通过随机数生成协方差估计值。

每一轮生成的解和方差组分样本本身不具统计意义，但通过大量重复抽样后，可通过样本直方图构建方差组分的后验分布。例如：

• 将样本均值（后验均值）作为方差组分的点估计值；
• 样本标准差（后验标准差）对应频率学派方法（如REML）中的标准误（standard error）。

基础的吉布斯抽样分析包含下面这2步：

1. 使用 THRGIBBSF90 来生成对方差组分的抽样（注意不能在同一个文件夹同时跑2个吉布斯抽样，因此会彼此覆盖结果文件 gibbs_samples）
2. 使用 POSTGIBBSF90 来整理抽样结果（确认设置的抽样次数和 burn-in 用于估计方差组分已经足够大了）

这两个模块使用的参数卡和 BLUPF90 相同，但是可能会增加一些选项来控制吉布斯抽样过程。

运行 THRGIBBSF90 时，除了参数卡名称，你还需要输入3个参数，总的迭代次数，丢弃前多少轮的结果（burn-in），每多少次迭代保存1次结果。

对方差组分的限制

当某些方差组分不存在时，这里如同 AIREML ，也支持将相应的方差组分设为 0 。但是与 AIREMLF90 不同的是，AIREMLF90 不会再估计不存在的方差组分，但是 THRGIBBSF90 仍然会估计所有的方差组分，然后将相应的不存在的方差组分设为 0 。

有一点与 AIREMLF90 相同，当估计的方差组分超过参数空间时，THRGIBBSF90 会尝试固定相应的需要估计的方差组分（我的理解，是不是说某个遗传参数的估计值超出了参数空间，就仍然使用之前的值）。在这个过程中，可能会出现某个方差组分为 0 的情况，此时这个方差组分之后的估计值会一直为 0 （因为没有信息）。因此如果出现这种某个方差组分碰巧为 0 的情况，你就要考虑可能模型欠佳，或者数据量不足以估计出这些方差组分。

THRGIBBSF90 结果文件

日志信息

打印在屏幕上的日志信息如下图，首先是每一轮抽样的方差组分结果（基于每多少次迭代保存1次结果这个参数抽取结果），这个程序本身也会计算每一次抽样的育种值结果，在抽样结束后也会打印出一些统计指标如下

其中两个值说明如下，同理另外两个就是残差。

• ave G ： 方差组分的后验均值
• SD G : 方差组分的后验标准差
• ave R ： 残差的后验均值
• SD R : 残差的后验标准差

接下来的内容如下，这里是一些不重要的信息，Last seeds 是随机数生成器的 Seeds 信息（不重要）， Number of samples kept 是剔除 burn-in 之后的抽样次数，solutions stored in a binary file: "last_solutions" 表示最后一次抽样结果保存为 “last_solutions” （这个文件不重要）

接下来是 DIC 和相应的信息， DIC 是 AIC 一样用于比较2个模型的优劣的指标（越小越好）（下面的公式推导没有看懂）。

其中信息如下

• detR ： 残差协方差矩阵的后验均值的行列式（）。
• # stored samples ： 有多少次抽样结果存储在文件中（等于 (总抽样次数-burn-in次数)/保存抽样结果的频数）
• D-bar :
• D(theta-bar) :
• DIC = 2*D-bar - D(theta-bar) :
• Effective number of parameters :
• solutions stored in file: "binary_final_solutions" : 抽样的后验均值结果存储在 binary_final_solutions 文件中

这里我们以 作为参数 的似然换数， 是表型数目， 表示第 次抽样（在 burn-in 之后）的残差方差。

我定义第 次的偏差 为

这里 为 ， 为某个观测值 的残差（ ）。

通过 burn-in 之后的 次抽样后， （ 的期望值的估计值）定义为

我们以同样的方式计算 后验均值的估计值。我们使用上面计算 的公式来计算后验偏差 ，只是 替换为了其后验均值

这些统计量都是使用 burn-in 之后的抽样结果，之后按照下面的公式计算统计量 DIC

THRGIBBSF90 会将每一个 值存储在文件 fort.99 中。

生成的文件

THRGIBBSF90 会生成下面的文件

• gibbs_samples ：burn-in后，特定轮次中方差组分（variance components）的抽样样本。
  （示例：用于绘制后验分布直方图或计算方差组分的后验均值）
• fort.99： burn-in后，特定轮次中模型的偏差值（Deviance）。
  （说明：常用于监控吉布斯抽样收敛性，偏差值波动趋稳表明链稳定）
• last_solutions ：最后一次抽样得到的解（如固定效应和随机效应估计值），以二进制形式保存。
  （用途：可作为后续分析的初始值或检查抽样终点状态）
• binary_final_solutions ： burn-in后 特定轮次中抽样解的后验均值，以二进制形式保存。
  （示例：直接用于基因组预测或遗传评估报告）

这里所说的特定轮次，就是基于我们所设置的每 轮抽样保存一次结果（总次数等于 (总抽样次数-burn-in次数)/n）。

这里 gibbs_samples 文件和 fort.99 文件 会被 POSTGIBBSF90 软件用于后续分析。

Post Gibbs 分析

POSTGIBBSF90 模块会执行 post-Gibbs 分析，其输入文件为参数卡名称，gibbs_samples 文件和 fort.99 文件 。

运行时第一个参数是参数卡名称，之后脚本会问你 Burn-in 的数目 （这里实际是问你需要额外丢弃的抽样数目）

这是一个需要谨慎的问题，你不能直接使用 THRGIBBSF90 使用的 burn-in 数目，因为这些抽样结果早已经被丢弃了（没有保存在文件中）。这里问的是你需要额外丢弃的抽样结果数目，因此你的回答应该是

• 0 ： 如果你不需要丢弃额外的抽样结果
• 一个整数，表示你想要丢弃额外的抽样结果

之后，脚本会问你下面的问题

这里你不能直接输入1 ，这里问的是你之前分析的 interval 参数，因此如果你要保留全部抽样结果，你就要输入 THRGIBBSF90 分析时使用的参数（或者之前interval 参数的整数倍）。

之后会输出下面的2个表格，每一个表格对应方差组分的统计值。

方差组分可以通过 Pos. eff1 eff2 trt1 trt2 这些列来确认：

• Pos. ： 一个参数的位置索引（应该就和行号一样，没用）
• eff1 和 eff2 是方差组分被定义的因子，如果没有母体加性效应或随机回归效应，这两个值相等。
• trt1 和 trt2 是相应的性状

在这个例子中，我们有2个方差组分，第一个加性方差是第4个因子，第二个就是残差。因此这里 position 1 就是加性方差，position 2 就是残差。

这里我们提取一些重要信息说明如下（我测试了一下，后验均值和标准差和日志信息一致）

表格里面有好几个方差组分的点估计值，包括mean， median ， mode （众数）。具体用哪个值更合理取决于后验分布的偏态情况。如果后验分布高度偏态，并且 mean 和 mode 结果相差很大，此时使用 mode 或 median 结果可能更好。注意这里

高后验密度区间（HPD95）是评估参数精度的常用统计量，属于区间估计范畴。其计算方式为：将排序后的样本按截断前2.5%最大值和后2.5%最小值，得到区间的上下界。

你可以通过**独立链大小（Independent chain size）和自相关系数（Autocorrelations）**可判断保存的样本是否足够：

1. 自相关性（Autocorrelations）
   • 相邻样本通常高度相关（因下一个样本基于当前样本生成）；
   • 若远隔样本间仍存在高自相关性，表明样本间依赖性强，绝对值相似度高，此时样本无法充分反映后验分布的信息；
   • 这个参数告诉我们自相关系数指示达到0相关所需的样本间隔数。
2. 独立链大小（Independent chain size）
   • 对应可视为独立样本的数量（如独立链大小为3，则等效于仅用3个独立样本计算后验均值和标准差，显然不足）；
   • 该指标直接提示是否需增加样本量。

展示完上面的表格后，软件会提示下面的信息

POSTGIBBSF90 支持通过外部软件 Gnuplot 绘制样本的时间序列图或直方图。操作说明如下：

1. 输入1：生成时间序列图（需系统已安装 Gnuplot；安装命令为 sudo yum install gnuplot ）；
2. 输入2：生成抽样结果的直方图；
3. 输入0：退出程序。

这里我们输入1 表示画时间序列图，这里会提示 position ，再次输入参数的 position 值（就是上面表格中的位置索引）。

下面对时间序列图（Time-series plot）做个介绍

• 用途：直观判断预烧期长度或总抽样量是否充足。
• 诊断依据：
  • 初始抽样结果几乎保持为初始值不变 ，则需延长预烧期；
  • 样本长期滞留于特定值域 ，这说明后验均值无法代表整体分布，需增加抽样总量。

如果我们输入2表示画抽样直方图，之后会提示 Type position and # bins ，让你输入参数的位置索引和 画图柱子数目

直方图（Histogram） 介绍如下

• 用途：展示后验分布形态。
• 典型分布：单峰分布（可能伴有右偏长尾）；
• 异常情况：
  • 双峰分布或分布过于平坦 → 需检查参数卡文件是否存在人为错误，或增加总抽样次数。

生成的文件

POSTGIBBSF90 会生成下面7个文件

• postout
  与屏幕显示内容完全一致的输出文件（2个统计信息的表格内容）。
• postind
  方差组分（variance components）与其在模型中的位置索引的对应关系文件。
• postmean
  以用户友好格式（如表格）展示的后验均值结果文件。
• postmeanCorr
  以用户友好格式展示的相关性后验均值文件。
  （示例：遗传相关性的点估计值）
• postsd
  后验标准差（Posterior Standard Deviation）的数值记录文件。
  （用途：评估参数估计的精确性，值越小表示置信度越高）
• fort.998
  基于 fort.99 文件计算的**对数似然值（logL）**记录文件。
• postgibbs_samples
  经过预烧期（burn-in）剔除和**间隔抽样（thinning）**处理后的最终抽样结果集合。

这里 postgibbs_samples 文件是文本文件，可以用于计算一些你自己感兴趣的别的统计指标。其内容如下表，说明如下

• 存储的抽样结果的次序数字（行号）
• 抽样结果的实际的抽样次数
• 参数数目
• 第4列及其它列就是抽样的参数结果（按照 position 的顺序）

1
2
3
4
5

1 1010 2 35.60 64.26
2 1020 2 36.98 63.53
3 1030 2 30.92 66.49
4 1040 2 34.93 65.60
5 1050 2 37.38 63.26

计算方差组分任意函数的后验分布

当我们需要计算由方差组分的函数的后验分布，如遗传力或遗传相关的后验分布。我们这里有2种方式，第一种方式，通过对 postgibbs_samples 文件内容计算每一次抽样的遗传力或遗传相关，得到其后验分布（计算标准误或画个直方图），示例代码如下

1
2
3
4
5

a=read.table("postgibbs_samples",header=FALSE)
# tell what each thing is
colnames(a)=c("i","iter","varu","vare")
a$h2=a$varu/(a$varu+a$vare)
hist(a$h2)

第二种方式就是通过 option se_covar_function 来指定需要计算的方差组分的函数。

阈值模型

不估计方差组分

表型就是按照1，2，3这种方式进行编码。

参数卡和正常模型的参数卡相同，只是增加了 2 个option。这里 OPTION cat 用于定义有序多分类性状的数目，数值是该性状的分类的水平数目（这里是3个水平，所以是3；如果是0表示是连续性状）。

OPTION fixed_var mean 表示采用固定的方差组分（不估计方差组分）计算 BLUP 结果的后验均值和后验标准差。

1
2

OPTION cat 3
OPTION fixed_var mean

此时之后这里没有估计的方差组分的参数，因此不用跑 postgibbsf90

估计方差组分

如果是下面的写法，OPTION solution mean 这个选项是需要估计方差组分的。

1
2

OPTION cat 3
OPTION solution mean

后续需要跑一下 postgibbsf90 ，确认结果没有问题

参考文献

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

2. http://nce.ads.uga.edu/wiki/doku.php?id=readme.gibbsf90plus