原文鏈接：http://tecdat.cn/?p=24334

最近我們被客戶要求撰寫關(guān)于貝葉斯線性回歸的研究報告，包括一些圖形和統(tǒng)計輸出。

像任何統(tǒng)計建模一樣，貝葉斯建?？赡苄枰獮槟愕难芯繂栴}設(shè)計合適的模型，然后開發(fā)該模型，使其符合你的數(shù)據(jù)假設(shè)并運行

1. 了解?Stan

統(tǒng)計模型可以在R或其他統(tǒng)計語言的各種包中進行擬合。但有時你在概念上可以設(shè)計的完美模型，在限制了你可以使用的分布和復(fù)雜性的軟件包或程序中很難或不可能實現(xiàn)。這時你可能想轉(zhuǎn)而使用統(tǒng)計編程語言，如Stan。

Stan是一種新式的語言，它提供了一種更全面的學習和實現(xiàn)貝葉斯模型的方法，可以適應(yīng)復(fù)雜的數(shù)據(jù)結(jié)構(gòu)。Stan開發(fā)團隊的一個目標是通過清晰的語法、更好的采樣器（這里的采樣是指從貝葉斯后驗分布中抽取樣本）以及與許多平臺（包括R、RStudio、ggplot2和Shiny）的集成，使貝葉斯建模更易于使用。

在這個入門教程中，我們將從一個線性模型開始，經(jīng)歷模型建立的迭代過程。在我們的高級stan教程中，我們將探索更復(fù)雜的模型結(jié)構(gòu)。

首先，在建立模型之前，你需要定義你的問題并了解你的數(shù)據(jù)。探索它們，繪制它們，計算一些匯總統(tǒng)計。

一旦你對你的數(shù)據(jù)和你想用統(tǒng)計模型回答的問題有了了解，你就可以開始建立貝葉斯模型的迭代過程。

• 設(shè)計你的模型。

• 選擇先驗

• 對后驗分布進行采樣。

• 檢查模型收斂（traceplots、rhats ）

• 使用后驗預(yù)測批判性地評估模型并檢查它們與您的數(shù)據(jù)的比較情況

• 重復(fù)…

模擬數(shù)據(jù)也是很好的做法，以確保你的模型正確，作為測試你的模型的另一種方式。

2. 數(shù)據(jù)

首先，讓我們找到一個可以擬合簡單線性模型的數(shù)據(jù)集。?氣候變化對地球最顯著的影響之一是北半球每年海冰范圍的減少。讓我們使用 Stan 的線性模型探索海冰范圍如何隨時間變化。

通過運行setwd("your-file-path")?包含您自己的文件路徑的代碼?，將您的工作目錄設(shè)置為您保存數(shù)據(jù)的文件夾??，F(xiàn)在，讓我們加載數(shù)據(jù)：

#?添加stringsAsFactors?=?F意味著數(shù)字變量將不會被#?作為因子/分類變量讀入ece?<-?red.cv("sv",?stinsAsFators?=?F)

我們來看一下數(shù)據(jù)：

我們可以用這些數(shù)據(jù)提出什么研究問題？以下情況如何：

研究問題：?北半球的海冰范圍是否會隨著時間的推移而減少？

為了探索這個問題的答案，首先我們可以做一個數(shù)字。

plot(?th?~?yr,?data)

圖 1. 北半球海冰范圍隨時間的變化。

現(xiàn)在，讓我們使用?lm().

l1?<-?lm(exnoh?~?yer,?data?=?sie)summary(l1)

我們可以將該模型添加到我們的繪圖中：

ablne(m1,?l?=?2,?ty?=?2,?w?=?3)

圖 2. 北半球海冰范圍隨時間的變化（加上線性模型擬合）。

記住線性模型的方程：

y = α + β ? x + 誤差

在?Stan?你需要指定你想模型。? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

也許我們已經(jīng)找到了問題的答案，但本教程的重點是探索使用編程語言?Stan，所以現(xiàn)在讓我們嘗試在 Stan 中編寫相同的模型。

準備數(shù)據(jù)

讓我們重命名變量并將年份從 1 索引到 39。關(guān)于貝葉斯模型的一個關(guān)鍵是您必須使用信息分布來描述數(shù)據(jù)中的變化。因此，您希望確保您的數(shù)據(jù)符合這些分布，并且它們將適用于您的模型。在這種情況下，我們真的想知道從數(shù)據(jù)集的開始到數(shù)據(jù)集結(jié)束的海冰是否發(fā)生了變化，而不是 1979 年到 2017 年。我們不需要我們的模型估計 500 年或 600 年的海冰是什么樣的，就在我們的數(shù)據(jù)集的持續(xù)時間內(nèi)。因此，我們將年份數(shù)據(jù)設(shè)置為索引 1 到 30 年。

x?<-?I(year?-?1978)

我們可以使用新數(shù)據(jù)重新運行該線性模型。

summary(lm1)

我們還可以從我們的簡單模型中提取一些關(guān)鍵的匯總統(tǒng)計數(shù)據(jù)，以便我們Stan?稍后可以將它們與模型的輸出進行比較?。

coeff[1]?#?截距值coeff[2]?#?斜率sigma(lm1)?#?殘差

現(xiàn)在讓我們將其轉(zhuǎn)換為用于輸入Stan?模型的數(shù)據(jù)框?。傳遞給 Stan 的數(shù)據(jù)需要是命名對象列表。此處給出的名稱需要與模型中使用的變量名稱相匹配。

庫

請確保安裝了以下庫（這些是本Stan?教程和下一個教程的庫?）。?rstan?是最重要的，如果您沒有 C++ 編譯器，則需要一些額外的東西。

3.?Stan?程序

我們將首先用語言編寫一個線性模型?Stan。這可以寫在你的 R 腳本中，或者單獨保存為一個?.stan?文件并調(diào)用到?R.

一個?Stan?程序具有三個必需的“塊”：

• “數(shù)據(jù)”?塊：您可以在其中聲明數(shù)據(jù)類型、它們的維度、任何限制（即 upper = 或 lower = ，用作檢查?Stan）及其名稱。

• “參數(shù)”?塊：您可以在此處指明要建模的參數(shù)和名稱。對于線性回歸，我們希望對回歸線周圍的誤差的截距、任何斜率和標準偏差進行建模。

• “模型”?塊：這是包含任何抽樣語句的地方，包括正在使用的模型。模型塊是指明要為參數(shù)包含的任何先驗分布的地方。如果未定義Stan?先驗，則?使用默認先驗?uniform(-infinity, +infinity)。您可以在聲明參數(shù)時使用上限或下限來限制先驗（即?lower = 0> 確保參數(shù)為正）。

采樣由?~?符號表示，并且?Stan?已經(jīng)包含許多常見的分布作為矢量化函數(shù)。

還有四個可選塊：

• “功能”

• "轉(zhuǎn)化的數(shù)據(jù)"

• "轉(zhuǎn)換后的參數(shù)

• "生成的數(shù)量"

注釋//?在 Stan中用 表示?。該write("model code", "file_name")?允許我們在 R 腳本中編寫 Stan 模型并將文件輸出到工作目錄（或者您可以設(shè)置不同的文件路徑）。

write("//簡單線性回歸的模型數(shù)據(jù)?int?<?lower?=?1?>?N;?//?樣本大小vector[N]?x;//?預(yù)測vecor[N]?y;//?結(jié)果參數(shù)??real?alha;?//?截距?real?beta;?//?斜率(回歸系數(shù))?real?<?loer?=?0?>?sima;?//?誤差SD模型??y?~?nrmal(alpha?+?x?*?beta?,?siga);產(chǎn)生的數(shù)量?//?后驗預(yù)測分布"?。"md1.stan"

首先，我們應(yīng)該檢查我們的?Stan?模型以確保我們編寫了一個文件。

現(xiàn)在讓我們保存該文件路徑。

?"model1.stan"

在這里，我們隱式地將?uniform(-infinity, +infinity)?先驗用于我們的參數(shù)。

4. 運行Stan?模型

Stan?程序C++?在被使用之前被遵守?。這意味著在 R 可以使用模型之前需要運行 C++ 代碼。為此，您必須?C++?安裝編譯器。編譯后，您可以在每個會話中多次使用模型，但在開始新R?會話時必須重新編譯?。有許多?C++?編譯器，而且它們在不同系統(tǒng)中通常是不同的。如果您的模型一堆錯誤，請不要擔心。只要模型可以與該stan()?函數(shù)一起使用?，它就可以正確編譯。如果我們想使用以前編寫的?.stan?文件，我們在file?函數(shù)中使用?參數(shù)?stan_model()?。

我們通過使用stan()?函數(shù)擬合我們的模型?，并為它提供模型、數(shù)據(jù)，并指示預(yù)熱的迭代次數(shù)（這些迭代稍后不會用于后驗分布，因為它們只是模型“預(yù)熱” ”），總迭代次數(shù)，我們要運行的鏈數(shù)，我們要使用的內(nèi)核數(shù)（Stan?為并行化而設(shè)置），它表示同時運行的鏈數(shù)（即，如果您的計算機有四個內(nèi)核） ，您可以在每個鏈上運行一個鏈，同時創(chuàng)建四個鏈）和細化，這是我們想要存儲我們的預(yù)熱后迭代的頻率。“thin = 1”將保留每次迭代，“thin = 2”將保留每一秒，依此類推……

Stan?如果warmup =?未指定參數(shù)，則自動使用一半的迭代作為預(yù)熱?。

stan(fie?=moel1,?data?=?data,?wrmup?=?500,?ier?=?1000,?chais?=?4,?cres?=?2,?thn?=?1)

訪問fit?對象的內(nèi)容?

結(jié)果?stan()?保存為?stanfit?對象（S4 類）。

我們可以通過執(zhí)行對象的名稱來獲取參數(shù)估計和采樣器診斷的匯總統(tǒng)計信息：

fit

模型輸出展示了什么？你怎么知道你的模型已經(jīng)收斂了？您能看到指示您的 C++ 編譯器已運行的文本嗎？

從這個輸出中，我們可以通過查看Rhat?每個參數(shù)的值來快速評估模型收斂性?。當這些值等于或接近 1 時，鏈已經(jīng)收斂。還有許多其他診斷方法，但這對 Stan 來說很重要。

我們還可以通過從模型對象中提取參數(shù)來查看參數(shù)的完整后驗。有很多方法可以查看后驗。

poteir?<-?exrat(fit)

extract()?將每個參數(shù)的后驗估計放入一個列表中。

讓我們與我們之前使用“l(fā)m”的估計進行比較：

plot(y?~?x)

圖 3. 北半球海冰范圍隨時間的變化（比較?Stan?線性模型擬合和一般?lm?擬合）。

結(jié)果與lm?輸出相同?。這是因為我們使用了一個簡單的模型，并且在我們的參數(shù)上放置了非信息先驗。

將回歸線估計中的可變性可視化的一種方法是繪制來自后驗的多個估計。

plot(y?~?x,?pch?=?20)

圖 4. 北半球海冰范圍隨時間的變化（Stan?線性模型擬合）。

5. 改變先驗

我們再試一次，但現(xiàn)在對海冰和時間之間的關(guān)系采用更有信息量的預(yù)設(shè)。我們將使用小標準差的正態(tài)先驗。如果我們使用標準差非常大的正態(tài)先驗（比如1000，或者10000），它們的作用與均勻先驗非常相似。

參數(shù)??real?alha;?//?截距?real?bta;?//?斜率(回歸系數(shù))?real?<?lowr?=?0?>?sima;?//?誤差SD模型?beta?~?nomal(1,?0.1);?y?~?normal(apha?+?x?*?beta?,?siga);

我們將擬合該模型并將其與使用均勻先驗的均值估計進行比較。

?stan(stn_oel)plot(y?~?x)

圖 5. 北半球海冰范圍隨時間的變化（Stan?線性模型擬合）。

后驗預(yù)測發(fā)生了什么變化？模型是否更好地擬合數(shù)據(jù)？為什么模型擬合發(fā)生了變化？通過制作非常窄的先驗分布，我們的模型改變了什么？

點擊標題查閱往期內(nèi)容

視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型

左右滑動查看更多

01

02

03

04

嘗試自己將先驗更改為一些不同的數(shù)字，看看會發(fā)生什么，這是貝葉斯建模中的一個常見問題，如果您的先驗分布非常窄，但不符合您對系統(tǒng)或數(shù)據(jù)分布的理解，您可以運行無法有意義地解釋數(shù)據(jù)變化的模型。然而，這并不是說您不應(yīng)該選擇一些信息豐富的先驗，您確實希望使用先前的分析和對您的研究系統(tǒng)的理解來告知您的模型先驗和設(shè)計。您只需要仔細考慮您做出的每個建模決策！

6. 收斂診斷

在繼續(xù)之前，我們應(yīng)該再次檢查模型參數(shù)的?Rhat?值、有效樣本大小 (?n_eff) 和跟蹤圖，以確保模型已收斂且可靠。

n_f?是有效樣本大小的粗略度量。通常只需要擔心這個數(shù)字小于迭代次數(shù)的 1/100 或 1/1000。

對于跟蹤圖，我們可以直接從后驗查看它們：

plot(alpha,?type?=?"l")plot(beta,?type?=?"l")plot(sigma,?type?=?"l")

圖 6. alpha 的跡線圖，截距。

對于更簡單的模型，收斂通常不是問題，除非您的代碼中有錯誤，或者運行采樣器的迭代次數(shù)太少。

收斂性差

嘗試僅運行 50 次迭代的模型并檢查跟蹤圖。

fid?<-?stan(?iter?=?50)

這在預(yù)熱后也有一些“轉(zhuǎn)換”，表明模型指定錯誤，或者采樣器未能完全采樣后驗。

plot(alpha,?type?=?"l")plot(beta,?type?=?"l")plot(sigma,?type?=?"l")

圖 7. alpha 的錯誤軌跡圖，截距。

參數(shù)匯總

我們也可以直接通過后驗得到參數(shù)的匯總。讓我們還繪制非貝葉斯線性模型值，以確保我們的模型運行正確

par(mfrow?=?c(1,3))plot(dnsty(alpha)

圖 8.Stan?模型擬合的密度圖分布?與一般lm?擬合的估計值比較?。

從后驗我們可以直接計算任何參數(shù)超過或低于某個感興趣值的概率。

beta > 0 的概率：

sum(beta>0)/lenth(beta)#?0

beta > 0.2 的概率：

sum(bea>0.2)/legth(beta)#?0

診斷圖?rstan

雖然我們可以直接使用后驗，但?rstan?內(nèi)置了許多有用的功能。

plot(fit)

圖 9.Stan?模型不同鏈的跟蹤圖?。

我們還可以查看后驗密度和直方圖。

dens(it)hist(ft)

圖 10.Stan?模型中截距、斜率和殘差方差的后驗密度圖和直方圖?。

我們可以生成指示平均參數(shù)估計值和我們可能感興趣的任何置信區(qū)間的圖。請注意，beta?和?sigma?參數(shù)的 95% 置信區(qū)間?非常小，因此您只能看到點。根據(jù)您自己數(shù)據(jù)的差異，當您進行自己的分析時，您可能會看到更小或更大的置信區(qū)間。

plot(fit)

圖 11.Stan?模型的參數(shù)估計?。

后驗預(yù)測檢查

對于預(yù)測和作為模型診斷的另一種形式，?Stan?可以使用隨機數(shù)生成器在每次迭代中為每個數(shù)據(jù)點生成預(yù)測值。通過這種方式，我們可以生成預(yù)測，這些預(yù)測也代表了我們模型和數(shù)據(jù)生成過程中的不確定性?？捎糜讷@取我們想要的關(guān)于后驗的任何其他信息，或?qū)π聰?shù)據(jù)進行預(yù)測。

參數(shù)?real?alpha;?//?截距?real?beta;?//?斜率(回歸系數(shù))?y?~?nomal(x?*?eta?+?alpa,?sgma);產(chǎn)生的數(shù)量??for?(n?in?1:N)?{?yp[n]?= normlrg(x[n]?* bta + apha, sima)?。?}

請注意，GQ（生成量）塊不支持矢量化，因此我們必須將其放入循環(huán)中。但是由于它被編譯為?C++，循環(huán)實際上非?？?，并且 Stan 每次迭代只評估一次 GQ 塊，因此它不會為您的采樣增加太多時間。通常，數(shù)據(jù)生成函數(shù)將是您在模型塊中使用的分布，但帶有?_rng?后綴。

?stan(modl2GQ,?data?,?ier?=?1000,?hans?=?4,?cres?=?2,?tin?=?1)

y_rep?從后驗中提取?值。

處理y_rep?值有很多選擇?。

每一行都是模型的一次迭代（單一后驗估計）。

我們可以制作一些更漂亮的圖。這個包是ggplot2。

在200次后驗抽樣中，比較y的密度和y的密度。

poy(y,?yrep[1:200,?])

圖 12. 比較隨機后驗抽取的估計值。

在這里，我們看到數(shù)據(jù)（深藍色）與我們的后驗預(yù)測非常吻合。

我們還可以使用它來比較匯總統(tǒng)計的估計值。

pp(y?=?y,?yep?=?yrep,?tat?=?"mean")

圖 13. 比較匯總統(tǒng)計的估計值。

我們可以更改傳遞給?stat?函數(shù)的函數(shù)，甚至可以自己編寫！

我們可以調(diào)查每個數(shù)據(jù)點的平均后驗預(yù)測與每個數(shù)據(jù)點的觀察值（默認線為 1:1）

scttrg(y?=?y,?yrp?=?yrep)

圖 14. 每個數(shù)據(jù)點的平均后驗預(yù)測與每個數(shù)據(jù)點的觀測值。

所以現(xiàn)在您已經(jīng)學習了如何運行線性模型?Stan?并檢查模型收斂性。

如有任何問題，請聯(lián)系我們！

點擊文末?“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語言STAN貝葉斯線性回歸模型分析氣候變化影響北半球海冰范圍和可視化檢查模型收斂性》。

點擊標題查閱往期內(nèi)容

R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進球數(shù)
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預(yù)測心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計