原文鏈接：http://tecdat.cn/?p=2857
本文是通過對area,perimeter,campactness幾個變量的貝葉斯建模,來查看他們對groovelength這個變量的影響.
並且對比rjagsR2jags和內置貝葉斯預測函數的結果。
讀取數據
read data seed=read.csv("seeds_dataset.csv")seed=seed[,1:7] 查看數據的結構str(seed) 'data.frame': 209 obs. of 7 variables: $ area : num 14.9 14.3 13.8 16.1 14.4 ... $ perimeter : num 14.6 14.1 13.9 15 14.2 ... $ campactness : num 0.881 0.905 0.895 0.903 0.895 ... $ length : num 5.55 5.29 5.32 5.66 5.39 ... $ width : num 3.33 3.34 3.38 3.56 3.31 ... $ asymmetry : num 1.02 2.7 2.26 1.35 2.46 ... $ groovelength: num 4.96 4.83 4.8 5.17 4.96 ...
建立回歸模型
Do a linear model lm(formula = groovelength ~ area + perimeter + campactness, data = seed) Residuals: Min 1Q Median 3Q Max -0.66375 -0.10094 0.00175 0.11081 0.45132 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 19.46173 2.45031 7.943 1.29e-13 *** area 0.49724 0.08721 5.701 4.10e-08 *** perimeter -0.63162 0.18179 -3.474 0.000624 *** campactness -14.05218 1.34325 -10.461 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1608 on 205 degrees of freedom Multiple R-squared: 0.895, Adjusted R-squared: 0.8934 F-statistic: 582.4 on 3 and 205 DF, p-value: < 2.2e-16
從回歸模型的結果來看，三的自變量對因變量都有顯著的意義。其中，area有正向的意義。而其他兩個變量是負向的影響。從r方的結果來看，達到了0.895，模型具有較好的解釋度。
使用arm軟件建立貝葉斯回歸模型
Bayesian analysis With bayesglm bayesglm(formula = groovelength ~ area + perimeter + campactness, data = seed) Deviance Residuals: Min 1Q Median 3Q Max -0.66331 -0.09974 -0.00002 0.11110 0.44841 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 18.90538 2.41549 7.827 2.63e-13 *** area 0.47826 0.08604 5.559 8.40e-08 *** perimeter -0.59252 0.17937 -3.303 0.00113 ** campactness -13.74353 1.32463 -10.375 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for gaussian family taken to be 0.02584982) Null deviance: 50.4491 on 208 degrees of freedom Residual deviance: 5.2992 on 205 degrees of freedom AIC: -164.91 Number of Fisher Scoring iterations: 6
從內置貝葉斯模型的結果來看，3個變量同樣是非常顯著，因此模型的結果和回歸模型類似。然后我們使用BUGS/JAGS軟件包來建立貝葉斯模型
使用 BUGS/JAGS軟件包來建立貝葉斯模型library(R2jags)
library(coda) 建立貝葉斯模型 jags(model.file='bayes.bug',parameters=c("area","perimeter","campactness","int"),data =list('a' =seed$area, 'b' =seed$perimeter, 'c' =seed$campactness, 'N'=N, 'y'=groovelength),n.chains =4,inits=NULL)
查看模型結果：
module glm loaded Compiling model graph Resolving undeclared variables Allocating nodes Graph information: Observed stochastic nodes: 209 Unobserved stochastic nodes: 5 Total graph size: 1608 Initializing model library('R2jags')bb <-jags1$BUGSoutput extract the "BUGS output" componentmm <-as.mcmc.bugs(bb) convert it to an "mcmc" object that coda can handleplot(jags1) large-format graph
從上面的圖中，我們可以看到哥哥自變量的中位數和置信區間。從置信區間來看，各個變量的取值和內置貝葉斯模型的結果類似。內置貝葉斯結果的值全部落入在了置信區間內。
然后繪制每次迭代中各個變量參數的軌跡圖
trace + density plots, same as above prettier trace plot
可以看到每個變量的參數都在一定區間內波動。同時可以看到誤差在一定的迭代次數之后趨於收斂。
然后繪制每個變量參數的密度圖
prettier density plot
可以看到每個變量的參數，的密度分布近似於正態分布。同時我們可以看到分布的均值和內置貝葉斯模型，得到的結果類似。
然后繪制每個變量參數的置信區間
estimate + credible interval plot
從結果來看，可以看到各個變量參數的置信區間，campatness和int的置信區間較大，而其他兩個變量的置信區間較小。
從上面的實驗結果對比，我們可以看到，三個自變量對因變量均有重要的影響。area,perimeter,campactness幾個變量他們對groovelength這個變量均有重要的影響.同時我們可以認為回歸模型的結果和內置貝葉斯模型的結果相似。然后我們使用rjags&R2jags軟件包來對數據進行貝葉斯型的建立，從結果來看，同樣和之前得到的模型結果相差不大。並且我們通過模型的迭代，可以得到每個參數的置信區間。