一個計算我的妻子是否懷孕的貝葉斯模型
點擊上方「Datartisan數據工匠」可訂閱哦！
在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
用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_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
TAG:Datartisan數據工匠 |
※技術：一個計算我的妻子是否懷孕的貝葉斯模型
※您的孩子是否也愛看這個？
※你的童年裡是否也有這樣一隻狗？
※分辨是否懷男的那些方法，總有一種適合孕媽你！
※父親節專屬，那些年你是否成為你家娃的「模特」，老爹的溫柔和耐心更留給孩子
※孩子是否出色，跟媽媽的性格有很大關係！尤其是這8個媽媽一定記住
※雙胎媽媽，你是否和我一樣？
※孕婦可從這8個特徵中猜測自己懷的是否是男寶，其中有一條超准
※一個家庭的教育好不好，看孩子平時是否有這種「特別表情」就知
※你的老宅中是否有這樣的銅鏡？
※我們是否處在一個計算機模擬的宇宙之中？
※揭秘：薩達姆是否還有一個倖存在世的孫子
※怎麼判斷一個女孩子是否喜歡你
※那時候的你是否和我一樣？
※怎樣判斷一個男人是否真的愛你
※艾爾莎霍斯卡，不知這個超模是否符合你的審美
※你是否和我一樣 有個這樣的媽媽 任勞任怨周圍都是她的身影
※這才是判斷你是否愛一個人的唯一標準
※如果當初你牽著我的手，我是否會是另一個樣子？