一個計算我的妻子是否懷孕的貝葉斯模型

在2015年的二月21日，我的妻子已經33天沒有來月經了，她懷孕了，這真是天大的好訊息！通常月經的週期是大約一個月，如果你們夫婦打算懷孕，那麼月經沒來或許是一個好訊息。但是33天，這還無法確定這是一個消失的月經週期，或許只是來晚了，那麼它是否真的是一個好訊息？

為了能獲得結論我建立了一個簡單的貝葉斯模型，基於這個模型，可以根據你當前距離上一次經期的天數、你歷史經期的起點資料來計算在當前經期週期中你懷孕的可能性。在此篇文章中我將闡述我所使用的資料、先驗思想、模型假設以及如何使用重點抽樣法獲取資料並用R語言運算出結果。在最後，我將解釋為什麼模型的運算結果最終並不重要。另外，我將附上簡便的指令碼以供讀者自行計算.

資料

非常幸運的是，在2014年的下半年間我的妻子一直在記錄她經期起始日期，否則我只能以僅擁有小量資料而告終。總體上我們擁有8個經期的起始日期資料，但是我採用的資料不是日期而是相鄰經期起始日間相隔的天數。 已經有33天。

所以日期發生得相對規律，以28天為一個週期迴圈。最後一次月經開始日期是在1月19日，所以在2月21日，距離最後一次經期發生日。

模型的建立

我要建立一個涵蓋生理週期的模型，包括受孕期和不受孕期，這顯然需要做大量的簡化。我做了一些總體假設如下：

• 一對情侶受孕與否不受其他因素的影響。
• 女方擁有固定的經期。
• 該對想要受孕的夫妻正在積極地嘗試受孕。換言之，如Wilcox et al. (2000) 推薦的每週兩次到三次受精。
• 一旦懷孕，期間將不會有經期。

接下來是我所做的具體假設：

• 假設兩個相鄰經期間隔的天數（days_between_periods）服從的正態分佈，其中均值（mean_period）和標準差（sd_period）未知。
• 假設在一對生育能力的夫妻（is_fertile為真 ）受孕時一個週期內懷上孕的概率是0.19（更多關於選定該值的由來見參考文獻）。不幸的是，並非所有的夫妻都具備生育能力，沒有生育能力則懷孕的機率為零。如果生育率被編碼為 0-1，那麼可生育率可以被簡潔的寫為 0.19* is_fertile.
• 在某一些不能受孕的時期（n_non_pregnant_periods）的懷孕失敗率則為(1 - 0.19 * is_fertile)^n_non_pregnant_periods
• 最後，如果你在這一個週期內（從上一次生理期至這一次生理期為一個週期）將不會懷孕；那麼最新一次經期距離下一個經期的天數（next_period）將必然會大於最新一次經期距離當前日期的天數（days_since_last_perio）。即，next_period < days_since_last_period的概率為零。這麼做看上去很奇怪因為這個事件是顯然的，但是我們在模型中將會要用到它。

基本的假設就是這樣了。但是為了使其更加實際，需要考慮使用一個似然函式，一個給定了引數和一些資料、計算在給定引數下資料的概率，通常而言是一個與概率成正比例的數值——似然值。因為這個似然值可能極小所以我需要對其取對數，從而避免引起數值問題。當用R語言設計似然函式時，總體上的模式如下：

• 方程將資料和引數作為選項。
• 通過預處理，將似然值的初始值設為1.0，相應的對數為0.0。（log_like <- 0.0）
• 用R語言呼叫概率密度分佈函式（比如dnorm, dbinom and dpois），用該函式計算模型中不同部分的似然值。然後將這些似然值相乘。對應地，將取對數後的似然值log_like相加。
• 為了讓d*函式返回對數似然值，只需新增引數log=TRUE。並且注意似然值0.0對應的取對值為-inf

所以，綜上所述該模型的對數似然函式如下：

這裡資料有標量days_since_last_period以及向量days_between_periods，而其他的引數將會被被估計出來。使用這個函式，我能從任意一個資料+引數的組合中得出對數似然函式值。但是，到這裡我只完成了建模的一半工作，我還需要先驗資訊！

關於經期，受孕和生育的先驗資訊

為了完善這個模型，我需要所有引數的先驗資訊。換言之，我需要明確在獲取資料之前這個模型包含了哪些資訊。具體上，我需要實驗開始前mean_period, sd_period, is_fertile, and is_pregnant的初始值。（雖然next_period也是一個引數，我不需要給出一個它的確切初始值，因為它的分佈完全由mean_period 和sd_period確定。另外，我還需要找到在一個週期內能受孕的可能值（上文中我設定為0.19）。這裡我使用了模糊、主觀的資料嗎？不！我到生育文獻中去尋找了更加有資訊價值的依據！

對於days_between_periods的分佈，其引數為mean_period和sd_period。這裡我使用了來自文章The normal variabilities of the menstrual cycle Cole et al, 2009 中的估計值，該文測量了184個年齡來自18-36歲的女性的經期規律。相鄰經期間天數的總平均值為27.7天。每一個參與實驗者的標準差的平均值為2.4。總體樣本的間隔天數的標準差為1.6。給定了這些估計值以後我令mean_period服從（27.7,2.4）的正態分佈，令sd_period服從均值為1.6，標準差為2.05的半正態分佈。如下：

對於引數is_fertile a以及引數is_pregnant我考慮了受精頻率作為先驗。想要確定可育的夫妻的比例幾乎是不可能的事情，因為這裡對於不育有各種不同的定義。 Van Geloven et al. (2013)做了一個小範圍的文獻回顧然後得出結論所有夫妻中有2%至5%的人被認為是不孕的。因為曾看到高達10%的情況，我決定取該範圍的上限。設定初始資料100%-5%=95%的夫妻是可孕的。

is_pregnant 是 0 1變數表示這對夫妻在最近的一輪週期中是否將要（或者說已經）受孕。在這裡我使用的先驗值是在一個週期內成功受孕的概率。當這對夫婦沒有生育能力時這個概率值顯然為0.0，但是積極地嘗試、可育的夫婦在一個週期內成功受孕的比例有多大呢？不幸的是我並沒有找到明確說明這一資料的文獻，但是我找到了比較接近的參照依據。在Increased Infertility With Age in Men and Women Dunson et al. (2004) 一書的第53頁，給出了在12個月中一直嘗試受孕但是沒有懷上的夫妻的比例，同時該資料也提供了女性不同年齡段的資料。

故上述即為對數似然函式中19%的懷孕概率值的由來，19%亦作為is_pregnant的先驗值。現在我有了所有引數的先驗，可以建立一個由先驗函式的抽樣函式了。

這裡使用了一個引數(n)，它輸出了一個n行的資料框，每一行是基於先驗數值得出的樣本資料。輸出結果如下：

使用重要性抽樣來擬合模型

現在，我已經收集了貝葉斯統計分析的三大要素：先驗資訊，似然函式以及資料。為了擬合模型我有很多方法，但是這裡有一個非常方便的方法——重要性抽樣。我之前曾寫文提及過重要性抽樣法，這裡我們來回顧一下：重要性抽樣法是一種蒙特卡洛實驗法，它建立起來非常簡單並且適用於以下情況：（1）引數空間非常小（2）先驗分佈與後驗分佈的形式區別不大。因為我的引數空間比較小，加之我使用了資訊量包含得比較豐富的先驗資料。因此，認為重點抽樣法在此例中是可用的。在重要性抽樣法中三個基本的步驟為:

1. 由先驗分佈產生大樣本（這裡可以通過sample_from_prior得到）
2. 給定了引數時，對每一個與似然值成比例的先驗資料進行賦權。（這裡可以通過 calc_log_like 得到）
3. 將權重歸一化，從而在先驗分佈的情況下形成了新的概率分佈。最終，根據此概率分佈對先驗分佈的樣本進行重新抽樣。（這裡可以用R函式抽樣）

（ 注意存在與該過程不同的多種方法，但是在用來擬合貝葉斯模型時，這是重要性抽樣法的常用版本）

因為我已經定義過 sample_from_prior 和 calc_log_like，因此需要定義一個新的方程來做後驗抽樣：

結果：懷孕的可能性

因此，在2月21日，2015，我的妻子已經沒有來月經33天了。這是一個好休息嗎？讓我們執行這個模型看看結果吧！

post這裡是一個長資料框，其中數值的表示基於這些引數得出的後驗分佈資訊。

讓我們來看看各個週期中間隔天數的均值和方差的變化吧。

像期望的那樣，後驗分佈的影像比先驗資料更狹長；並且觀察後驗資料，大致得出平均的經期週期天數在29天左右，其標準差在2-3天左右。那麼重要問題來了：我們是可育夫妻的概率為多少，以及我們在2月21日確定已經懷孕的概率為多少？為了計算這個我們取 postisfertile 與post is_pregnant，並計算眾數。當然一個捷徑是直接採用均值。

因此這是一個相當好的訊息：我們極有可能是可孕的夫妻，並且我們已經受孕的可能性高達84%！用這個模型我可以瞭解到：當經期的來臨再多延遲幾天，我們確定懷孕的概率是如何隨之而變化的。

是的，既然我們已經得到了好訊息，為何不看看在我們之前嘗試受孕的數月中我們可孕和受孕成功的可能性是如何變動的呢？

因此，這能夠解釋得通了。在距離“最後一次經期”的時間變長時，我們在當前週期懷孕成功的機率增加了，但是一旦這裡有經期發生時可能性會跌回至基線。我們看到的可孕率曲線是幾乎相同的，但在我們沒有懷孕成功的每一個週期裡，可孕率曲線稍稍有所下降。這兩個指標的影像均呈輕微的鋸齒狀，但這只是由重要性抽樣演算法的偏差導致的。另外需指出的是，雖然上述的影像非常美觀，但是檢視未懷孕期間的概率是無效的，唯一具重要性和資訊價值的是當前的概率。

一些關於這個模型的批評 但其實並不重要

• 當然，比起我相對簡單粗糙的計算，別人有可能能夠得到更優越的先驗值。還有很多可以加入考慮的預測因子，如男性的年齡，健康因子等等。
• 每個月受孕的概率本應被視作一個不確定的值而不是一個固定值，而我把它設為了固定值。但是在擁有的給定資料很少的情況下，我將其視作一個適用於多個引數的引數值。
• 沒有事物是完全符合正態分佈的，兩個經期間的天數亦然。這裡我認為假設是適用的，但是還有比我的假設遠要複雜得多經期間隔天數模型，比如 Bortot et al (2010)建立的模型。

---

---

原文網站：http://sumsar.net/blog/2015/11/a-bayesian-model-to-calculate-whether-my-wife-is-pregnant/
翻譯：lan
作者：Rasmus Bååth