健康月经周期的贝叶斯分层功能数据分析

摘要

在几乎每一个研究领域，研究随机曲线的分布可以提供所需的信息。在体外受精领域，特别感兴趣的是对整个月经周期的女性温度曲线建模。拟合曲线可以将周期划分为健康或不健康，并有助于确定一个周期中最容易受孕的日子。

关键字

层次模型，吉布斯抽样，模拟，贝叶斯规范，月经周期。

介绍

在几乎每一个研究领域，研究随机曲线的分布可以提供所需的信息。在体外受精领域，特别感兴趣的是对整个月经周期的女性温度曲线建模。拟合曲线可以将周期划分为健康或不健康，并有助于确定一个周期中最容易受孕的日子。健康的月经周期包括三个阶段，卵泡期、排卵期和黄体期。在卵泡期，基础体温(休息期间达到的最低温度，本文其余部分称为bbt)保持在一个低平台上不变。在排卵期间，bbt会增加，直到达到最高点。然后，bbt在黄体期保持恒定。在周期结束时，bbt回落到月经前的低平台，循环重复[1]。图1给出了几个健康月经周期的BBT水平的例子[2]。第2节概述了一般的层次模型，然后给出了bbt曲线的参数化建模。第2节的末尾完成了贝叶斯规范，并指定了先前分布所需的超参数。第3节简要介绍了吉布斯抽样算法如何工作，然后将其应用于模型。第4节给出了变更点问题的简短描述，并定义了识别模型中的两个变更点所需的停止规则。第5节详细介绍了数据的模拟，然后应用有意倾斜超参数的吉布斯采样算法。仿真结果将在第5节的末尾进行讨论。最后，第6节对模型和结论进行了简短的讨论。

层次模型

首先，我们将定义描述模型所必需的术语。未知数量X的先验分布是在考虑数据之前表达对X的确定性的概率分布。类似地，未知数量X的后验分布是该数量以观测数据为条件的概率分布。超参数是先验分布的参数，它不是随机分布。例如，考虑来自a的随机变量X分布，其中Σ是已知的，我们可以使用分布到模型μ。μ是一个参数，是先验分布，μ₀和Σ₀hyperparameters。

一般模型

让参考女子i,jε {1,2，.....n_我}表示女人i的周期j，并且参考妇女i的周期j中的第t天。这些指标将在本文的其余部分保留这些定义。我们可以建立以下bbt曲线模型[3.]：

在ε_ij是测量误差和η_ij平滑函数是什么来.由于预期同一女性的周期是相关的，我们可以为女性的功能分布建立先验分布叫它G_我。然而，健康的月经周期在女性中遵循类似的模式，因此我们可以为不同女性的分布集合建立一个先验分布因此，如果我们假设我们的误差是正态分布的，我们有以下模型[3.]：

基础表示

现在我们有了一般的模型，有很多不同的表示η的方法_ij。一个常见的策略是考虑一个基表示[3.]：

在哪里预先指定的基函数和为周期特定基系数。自是预先设定的，我们的主要目标是建模θ_ij我们使用的方法是给出θ_ijS为分层正态模型[3.]：

其中N_k(μ，Σ)为均值向量μ，协方差矩阵Σ的多元正态分布。女性特有的基系数由αi给出，α给出全球平均基系数，Ω和Ω₁分别为女性内部变异和女性之间变异的精度矩阵。为了完成这个模型，我们需要指定基函数建立θ的分量_ij。我们将在下一节考虑参数化层次结构时这样做。

参数化层次模型

正如我们在第1节中所讨论的，健康周期的bbt水平遵循双相模式，其中bbt水平从卵泡期的低平台开始，然后在排卵期间迅速上升到黄体期的高平台。这样我们就可以建立η模型_ij(t)如下[3.]：

在哪里在θ_ij1和θ_ij2分别为女性i的第j周期的卵泡期温度和排卵时bbt的增加。我们可以进一步定义为卵泡期和r的最后一天_ij为排卵期间BBT上升的天数(图2)。因此，如果我们考虑y_ij(t)表示Y_ij一个n的向量_ijBBT值，我们可以用矩阵形式表示我们的层次模型[3.]：

因此X_ij(t) t^thX的第3行，对应于妇女i的周期j的第t天。所以我们已经定义了预先指定的基函数这是特定于周期的基系数。

贝叶斯规范

现在我们有了within (Ω)和between (Ω)₁)女性变异性，滤泡期的最后一天(k_ij)表示排卵期天数(r_ij)、全球平均值(α)和测量误差方差(σ₂)留给模特。我们将使用以下方程式来模拟女性内部和女性之间的变异性:

我们用以下先验完成贝叶斯规范[3.]：

在哪里为具有形状参数a和逆尺度参数b的伽马分布，为a与b之间的离散均匀分布，mij为月经后第一天。因此，我们预先指定的超参数是α₀,Σ_αc d a_hb_h,一个_h1，和b_h2(对于h = 1,2).参数、超参数、它们的描述和分布的摘要在附录中给出。

吉布斯抽样

我们的目标是在给定一组观测数据的情况下找到模型的参数。由于找到联合分布是一个复杂的过程，我们需要使用不同的方法。一种常见的方法是使用蒙特卡洛马尔可夫链作为迭代过程。我们要用到的蒙特卡洛马尔可夫链的具体类型叫做吉布斯采样器。

定义

假设有随机变量X₁, X₂X、…_n。

1.初始化:给定一个初值。

2.迭代:对于1≤I≤K，样本x_(我)^jX的条件分布_j鉴于和对于1≤j≤n

吉布斯抽样过程是一个马尔可夫链，其中平稳分布是X的联合分布₁, X₂X、…_n.因此，对于较大的k，样本近似X的联合分布₁, X₂X、…_n。一个阻塞的吉布斯采样器使用相同的步骤，只是它将一些随机变量分组在一起，以便一个块中的变量从它们的联合分布中取样，这些联合分布以所有其他变量为条件。

吉布斯抽样算法

这是我们的模型的阻塞吉布斯采样算法，我们给了参数一些初始值[2]：

1.Cycle-specific系数:

2.专门针对女性的意思是:

3.全球的意思是:

4.Ω(h = 1,2)的分量:

5.的组件

6.误差方差:

7.滤泡期最后一天:对k使用变化点停止规则_ij(见第4节)。

8.排卵期:r使用变化点停止规则_ij(见第4节)。

因此，如果我们使用阻塞吉布斯抽样步骤来更新未知参数的条件分布，那么经过大量的迭代之后，我们将从所有未知参数的联合分布中进行抽样。

变点问题

现在我们需要考虑kij和rij，因为它们是bbt曲线中曲线形状发生变化的特定点(称为变化点)。寻找改变点P的过程是直观的。给定一组观察结果我们会得到在给定ym条件下p≤m的概率。我们希望这个概率高于某个高阈值q，因此我们有一个“停止规则”[4]：

在m处终止的当且仅当

从m = 1开始，第一次概率大于Q的是我们的变点P。

发现kij

首先我们要求k_ij假设其他参数都是已知的。让X_k是矩阵X的前m行_ij上面已经建立，但是有k_ij= k，因此。让{π_k}是k的先验分布_ij，π(.)为θ的密度函数_ij， f(.)为数据yij,m的密度函数。现在我们需要找到k的后验分布_ij为了找到对于停止规则。

让成为改变点的可能性。然后实现贝叶斯定理，我们得到[4]

因为L_k不是函数θ_ij我们可以替换θ_ij用0。也不依赖于k，那么L_k和k成正比

θ的后验分布_ij是N₂(μ_k, V_k),

注意，这与θ的后验分布相同_ij在第3节中找到，但是用了改变的矩阵X_K。利用θ的后验分布_ij，我们发现

因为我们需要找到，我们会先找到.从应用贝叶斯定理到，我们得到

因为L_Kk与W成正比_K，那么我们可以用Wk代替L_K，导致

当数据可用时，这些概率很容易计算。当k≥m时，

所以X_K不依赖于k(因此W_K= W_米)。因此，我们的停止规则是当且仅当在m处终止

一旦我们达到了停止规则[4，让k_ij= m。

发现r_ij

我们可以用同样的方法求r_ij通过使用转换来“向后”检查数据，即首先考虑最后一天，并以相反的顺序进行循环。所以我们新的反向数据是

让是我们的变换矩阵。然后

因此，如果我们让Xr*对应于上述系统的前m行，则

这个变换可以让我们用和求k时一样的方法_ij。让{π_r}是r(第二个变更点)的先验分布。θ*的新后验分布_ij是N₂(μ_r, V_r),

请注意，由于从而给出了在第3节中发现的θ的类似后验分布_ij。使用*的后验分布_ij我们发现

因此，我们的停止规则是当且仅当在m处终止

一旦我们达到停止规则，r =n_ij-m(因为我们转换了数据)。自r_ijk和k之间的天数是多少_ij然后是r_ij= r -k_ij

模拟

数据模拟

现在已经描述了该方法，我们将介绍一个模拟来演示该方法。模拟数据由20名女性组成，每名女性有10个月经周期，所以我们总共有200个月经周期。每个周期为30天[5]。我们任意地将全局均值和协方差矩阵赋值为

每个女性的特定平均值是根据第2节中讨论的模型生成的，因此然后根据，生成每个循环的kij和rij.最后，根据。生成循环特定系数而测量误差项根据在σ²= 0.01。现在我们生成了模型所需的所有参数，我们根据生成了bbt值

初始化

为了运行算法，我们需要指定超参数(我们最初的猜测)。为了考虑误差，我们故意扰动超参数，使其远离数据模拟部分中建立的实际参数。因此我们将超参数设置为

这些超参数在变化点问题一节中概述的吉布斯采样步骤中使用。但在我们实现吉布斯采样算法之前，我们需要使用超参数初始化未知数的先验分布，如下所示:

我们还需要从它们的先验分布初始化α，因为吉布斯采样的第一步是寻找θ的后验分布_ij条件是αi。我们会得到α_我来自先验分布的值

初始化的最后一步是设置k_ij和r_ij将所有值设置为固定值7。

运行算法

现在初始化已经完成，我们可以运行算法了。我们总共运行了6000次迭代，其中1000次迭代是旧的。为了确保参数是收敛的，我们执行了一些变量的跟踪图。图3和4显示α的轨迹图₁₁和α₁₂分别。

表1和2总结估计模型参数的结果。周期特定参数(k_ij, r_ij,θ_{j - 1}我和θ_j2)是所有周期的总结，而女性特定参数(α_i1和α_i2)是对所有女性的总结。从表中可以看出，这些参数的分布得到了很好的估计。图5和6显示两个示例周期的真实分布和估计分布的图形。

结论

由此可见，阻塞吉布斯抽样算法是一种准确估计分层模型各参数的方法。最重要的是，变化点分析准确估计了卵泡期的最后一天和排卵期间温度上升的天数(附录)。这两个数字对于确定月经周期中最容易受孕的时期至关重要。本主题的自然延伸将是使用非参数方法对不健康循环建模。将参数模型和非参数模型结合起来可以进行更广泛的分析和更多不寻常的数据集。

参考文献

1. 麦克卢尔·布朗JC。研究生妇产科第四版355。伦敦:巴特沃思出版社，1973年。
2. 卵形生育图。2015年www.ovagraph.com/chart-tags/typical-cycles.html。
3. 斯卡帕B和唐森DB。基于污染信息先验的贝叶斯分层函数数据分析。生物识别技术。2009;65: 772 - 780。
4. 卡特RL和布莱特BJN。贝叶斯变点问题及其在女性排卵预测和检测中的应用。生物识别技术。1981;371: 743 - 751。
5. Royston JP和Abrams RM。一种检测女性基础体温变化的客观方法。Biometrics.1980;36: 217 - 224