R語言用WinBUGS 軟件對學術能力測驗（SAT）創建層次（分層）貝葉斯模型

原文連接：http://tecdat.cn/?p=11974

---

R2WinBUGS軟件包提供了從R調用WinBUGS的便捷功能。它自動以WinBUGS可讀的格式寫入數據和腳本，以進行批處理（自1.4版開始）。WinBUGS流程完成後，能夠經過程序包自己將結果數據讀取到R中（這提供了推斷和收斂診斷的緊湊圖形摘要），也可使用coda程序包的功能對輸出進行進一步分析。 

 WinBUGS軟件可從http://www.mrc-bsu.cam.ac.uk/...。  

R是一種「用於數據分析和圖形處理的語言」，是一種實現該語言的開放源代碼和免費提供的統計軟件包，請參見http://www.R-project.org/。  R和R2WinBUGS可從CRAN 得到，即http://CRAN.R-Project.org或其...。 若是可使用Internet鏈接，則能夠在R命令提示符下鍵入install.packages（「 R2WinBUGS」）來安裝R2WinBUGS。別忘了用library(R2WinBUGS)

例子

學校數據

學術能力測驗（SAT）衡量高中生的能力，來幫助大學作出入學決定。 咱們的數據來自1970年代後期進行的一項實驗，來自八所高中的SAT-V（學業能力測試語言）。SAT-V是由教育測試服務局管理的標準多項選擇測試。該服務對所選學校中每所學校的教練計劃的效果很感興趣。

實現

R2WinBUGS軟件包的實現很是簡單。 main「函數bugs（） 由用戶調用。原則上，它是對 其中逐步調用的其餘幾個函數的包，以下：

1. bugs.data.inits（）寫入數據文件' data.txt」和「 inits1.txt」，「 inits2.txt」 ...進入 工做目錄。 
2. bugs.script（）寫入WinBUGS用於批處理的文件「 script.txt」。
3. bugs.run（）更新WinBUGS註冊表 ，調用WinBUGS，並使用 'script.txt' 以批處理模式運行它。
4. bugs.sims（）若是參數codaPkg已設置爲false（默認值）才調用。
   不然，bugs（）返回存儲數據的文件名。例如，這些能夠經過打包的coda 導入，該軟件包提供了收斂診斷，蒙特卡洛估計的計算，跡線圖等功能。
   bugs.sims（）函數將WinBUGS中的模擬讀取到R中，將其格式化，監視收斂，執行收斂檢查並計算中位數和分位數。它還爲bugs（）自己準備輸出。
   這些功能不禁用戶直接調用。參數將從bugs（）傳遞給其餘函數。

例子 

咱們將 R2WinBUGS提供的功能應用於示例數據並分析輸出。

學校數據

示例數據 ：

> schools

​
爲了對這些數據進行建模，咱們使用了Gelman等人提出的分層模型。 咱們假設每所學校的觀測估計值具備正態分佈，且均值theta 和方差tau.y，逆方差爲1 =σ.y2，其先驗分佈在（0,1000）上是均勻的。對於均值theta，咱們採用另外一個正態分佈 平均爲mu.theta和逆方差爲tau.theta。有關其先驗分佈，請參見如下WinBUGS代碼：

model {
for (j in 1:J)
{
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}

此模型必須存儲在單獨的文件中，例如'schools.bug'2，在適當的目錄中，例如c：/ schools /。在R中，用戶必須準備bugs（）函數所需的數據輸入。這能夠是包含每一個數據向量名稱的列表，例如

> J <- nrow(schools)

 使用這些數據和模型文件，咱們能夠運行MCMC模擬以獲取theta， mu.theta和sigma.theta的估計值。在運行以前，用戶必須肯定要運行多少個鏈 （n.chain = 3）和迭代次數（n.iter = 1000）。另外，用戶必須指定鏈的初始值，例如經過編寫函數：

> inits <- function(){
+ list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
+ sigma.theta = runif(1, 0, 100))
+ }

能夠開始MCMC模擬，R中的參數bugs.directory必須指向WinBUGS的安裝目錄。能夠經過print（schools.sim）方便地輸出school.sim對象中的結果。 
對於此示例，將得到相似的結果

Inference for Bugs model at "c:/schools/schools.bug"
3 chains, each with 1000 iterations (first 500 discarded)
n.sims = 1500 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
theta[1] 11.1 9.1 -3.0 5.0 10.0 16.0 31.8 1.1 39
theta[2] 7.6 6.6 -4.7 3.3 7.8 11.6 21.1 1.1 42
theta[3] 5.7 8.4 -12.5 0.6 6.1 10.8 21.8 1.0 150
theta[4] 7.1 7.0 -6.6 2.7 7.2 11.5 21.0 1.1 42
theta[5] 5.1 6.8 -9.5 0.7 5.2 9.7 18.1 1.0 83
theta[6] 5.7 7.3 -9.7 1.0 6.2 10.2 20.0 1.0 56
theta[7] 10.4 7.3 -2.1 5.3 9.8 15.3 25.5 1.1 27
theta[8] 8.3 8.4 -6.6 2.8 8.1 12.7 26.2 1.0 64
mu.theta 7.6 5.9 -3.0 3.7 8.0 11.0 19.5 1.1 35
sigma.theta 6.7 5.6 0.3 2.8 5.1 9.2 21.2 1.1 46
deviance 60.8 2.5 57.0 59.1 60.2 62.1 66.6 1.0 170
pD = 3 and DIC = 63.8 (using the rule, pD = var(deviance)/2)
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC is an estimate of expected predictive error (lower deviance is better).

此外，用戶能夠經過輸入plot（schools.sim）生成結果圖。 結果圖如圖所示。在該圖中，左列顯示瞭如下內容的快速摘要：
推論和收斂（全部參數的Rb都接近1.0，代表三個鏈的良好混合，所以近似收斂）；右列顯示每組參數的推論。從右欄中能夠看到，R2WinBUGS使用 WinBUGS中的參數名稱將輸出構造爲標量，向量和參數數組。
​