R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)(SAT)建立分層模型
原文鏈接:http://tecdat.cn/?p=11974
R2WinBUGS軟件包提供了從R調(diào)用WinBUGS的便捷功能。它自動(dòng)以WinBUGS可讀的格式寫入數(shù)據(jù)和腳本,以進(jìn)行批處理(自1.4版開始)。WinBUGS流程完成后,可以通過(guò)程序包本身將結(jié)果數(shù)據(jù)讀取到R中(這提供了推斷和收斂診斷的緊湊圖形摘要),也可以使用coda程序包的功能對(duì)輸出進(jìn)行進(jìn)一步分析。?
?WinBUGS軟件可從http://www.mrc-bsu.cam.ac.uk/bugs/免費(fèi)獲得。??
R是一種“用于數(shù)據(jù)分析和圖形處理的語(yǔ)言”,是一種實(shí)現(xiàn)該語(yǔ)言的開放源代碼和免費(fèi)提供的統(tǒng)計(jì)軟件包,請(qǐng)參見http://www.R-project.org/。? R和R2WinBUGS可從CRAN 獲得,即http://CRAN.R-Project.org或其鏡像之一。R2WinBUGS可以移植到商業(yè)實(shí)施S-Plus。 如果可以使用Internet連接,則可以在R命令提示符下鍵入install.packages(“ R2WinBUGS”)來(lái)安裝R2WinBUGS。別忘了用library(R2WinBUGS)
?
例子
?
學(xué)校數(shù)據(jù)
學(xué)術(shù)能力測(cè)驗(yàn)(SAT)衡量高中生的能力,以幫助大學(xué)做出入學(xué)決定。 我們的數(shù)據(jù)來(lái)自1970年代后期進(jìn)行的一項(xiàng)實(shí)驗(yàn),來(lái)自八所高中的SAT-V(學(xué)業(yè)能力測(cè)試語(yǔ)言)。SAT-V是由教育測(cè)試服務(wù)局管理的標(biāo)準(zhǔn)多項(xiàng)選擇測(cè)試。該服務(wù)對(duì)所選學(xué)校中每所學(xué)校的教練計(jì)劃的效果很感興趣。
?
實(shí)現(xiàn)
R2WinBUGS軟件包的實(shí)現(xiàn)非常簡(jiǎn)單。\ main“函數(shù)bugs()?旨在由用戶調(diào)用。原則上,它是對(duì)?其中逐步調(diào)用的其他幾個(gè)函數(shù)的包裝,如下:
1. bugs.data.inits()寫入數(shù)據(jù)文件' data.txt”和“ inits1.txt”,“ inits2.txt” ...進(jìn)入 工作目錄。?
2. bugs.script()寫入WinBUGS用于批處理的文件“ script.txt”。
3. bugs.run()更新WinBUGS注冊(cè)表中自適應(yīng)階段的長(zhǎng)度 ,調(diào)用WinBUGS,并使用 'script.txt' 以批處理模式運(yùn)行它。
R2WinBUGS軟件包的實(shí)現(xiàn)非常簡(jiǎn)單。\ main“函數(shù)bugs()
?原則上,它是對(duì) 其中逐步調(diào)用的其他幾個(gè)函數(shù)的包裝,如下:
1. bugs.data.inits()寫入數(shù)據(jù)文件' data.txt”和“ inits1.txt”,“ inits2.txt” ...進(jìn)入 工作目錄,這些文件將由WinBUGS在批處理過(guò)程中使用。?
2 . bugs.script()寫入文件“ script”。 txt”,由WinBUGS用于批處理 。
3. bugs.run()更新WinBUGS注冊(cè)表中自適應(yīng)階段的長(zhǎng)度(使用bugs.update.settings()函數(shù)),調(diào)用WinBUGS并批量運(yùn)行“script.txt”。
4. bugs.sims()如果參數(shù)codaPkg已設(shè)置為false(默認(rèn)值)才調(diào)用。
否則,bugs()返回存儲(chǔ)數(shù)據(jù)的文件名。例如,這些可以通過(guò)打包的coda 導(dǎo)入,該軟件包提供了收斂診斷,蒙特卡洛估計(jì)的計(jì)算,跡線圖等功能。
bugs.sims()函數(shù)將WinBUGS中的模擬讀取到R中(不一定 由bugs()本身調(diào)用),將其格式化,監(jiān)視收斂,執(zhí)行收斂檢查并
計(jì)算中位數(shù)和分位數(shù)。它還為bugs()本身準(zhǔn)備輸出。
這些功能不旨在由用戶直接調(diào)用。 如果合適,參數(shù)將從bugs()傳遞給其他函數(shù)。一個(gè)簡(jiǎn)短的bugs()幫助文件,其中列出了?
軟件包R2WinBUGS 之后,請(qǐng)?jiān)赗中進(jìn)行調(diào)試?
?
例子?
我們將 R2WinBUGS提供的功能應(yīng)用于示例數(shù)據(jù)并分析輸出。
學(xué)校數(shù)據(jù)
示例數(shù)據(jù) :
?
> schools

為了對(duì)這些數(shù)據(jù)進(jìn)行建模,我們使用了Gelman等人提出的分層模型。 我們假設(shè)每所學(xué)校的觀測(cè)估計(jì)值具有正態(tài)分布,且均值theta 和反方差tau.y。逆方差為1 =σ.y2,其先驗(yàn)分布在(0,1000)上是均勻的。對(duì)于均值theta,我們采用另一個(gè)正態(tài)分布 平均mu.theta和反方差tau.theta。有關(guān)其先前的分布,請(qǐng)參見以下
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)
}
此模型必須存儲(chǔ)在單獨(dú)的文件中,例如'schools.bug'2,在適當(dāng)?shù)哪夸浿校鏲:/ schools /。在R中,用戶必須準(zhǔn)備bugs()函數(shù)所需的數(shù)據(jù)輸入。這可以是包含每個(gè)數(shù)據(jù)向量名稱的列表,例如
> J <- nrow(schools)
?使用這些數(shù)據(jù)和模型文件,我們可以運(yùn)行MCMC模擬以獲取theta, mu.theta和sigma.theta的估計(jì)值。在運(yùn)行之前,用戶必須確定要運(yùn)行多少個(gè)鏈 (n.chain = 3)和迭代次數(shù)(n.iter = 1000)。另外,用戶必須指定鏈的初始值,例如通過(guò)編寫函數(shù)。
> inits <- function(){
+ list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
+ sigma.theta = runif(1, 0, 100))
+ }
可以開始MCMC模擬,R中的參數(shù)bugs.directory必須指向WinBUGS的安裝目錄。可以通過(guò)print(schools.sim)方便地打印school.sim對(duì)象中的結(jié)果。 泛型函數(shù)print()調(diào)用由R2WinBUGS提供的類錯(cuò)誤對(duì)象的print方法。
對(duì)于此示例,將獲得類似的結(jié)果
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).
此外,用戶可以通過(guò)鍵入plot(schools.sim)生成結(jié)果圖。 結(jié)果圖如圖4所示。在該圖中,左列顯示了以下內(nèi)容的快速摘要:
推論和收斂(所有參數(shù)的Rb都接近1.0,表明三個(gè)鏈的良好混合,因此近似收斂);右列顯示每組參數(shù)的推論。從右欄中可以看到,R2WinBUGS使用 WinBUGS中的參數(shù)名稱將輸出構(gòu)造為標(biāo)量,向量和參數(shù)數(shù)組。
