library(BMS)
data(attitude)
att = bms(attitude, mprior = "uniform", g = "UIP", user.int = F)
coef(att)
> coef(att)
PIP Post Mean Post SD Cond.Pos.Sign Idx
complaints 0.9996351 0.684449094 0.12848607 1.00000000 2
learning 0.4056392 0.096481513 0.15035025 1.00000000 4
advance 0.2129325 -0.026686161 0.09044093 0.00000107 7
privileges 0.1737658 -0.011854183 0.06055425 0.00046267 3
raises 0.1665853 0.010567022 0.08231172 0.73338938 5
critical 0.1535886 0.001034563 0.05368470 0.89769774 6
#We see that with 99:9%, virtually all of posterior model
#mass rests on models that include "complaints".
topmodels.bma(att)[, 1:5]
> topmodels.bma(att)[, 1:5]
20 28 29 24 30
complaints 1.0000000 1.0000000 1.00000000 1.00000000 1.00000000
privileges 0.0000000 0.0000000 0.00000000 0.00000000 1.00000000
learning 0.0000000 1.0000000 1.00000000 0.00000000 0.00000000
raises 0.0000000 0.0000000 0.00000000 1.00000000 0.00000000
critical 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000
advance 0.0000000 0.0000000 1.00000000 0.00000000 0.00000000
PMP (Exact) 0.2947721 0.1661679 0.06678871 0.05891111 0.05690948
PMP (MCMC) 0.2947721 0.1661679 0.06678871 0.05891111 0.05690948
#We see that the best model, with 29% posterior model probability,
#is the one that only includes complaints.
#预测
predict(att) #fitted values based on MCM frequencies
predict(att, exact=TRUE) #fitted values based on best models
#预设是 100 models
#采(20,28,29,24,30) 5 topmodels
#当然你可以试更多cases,比较预测结果
predict(att, topmodels=1:5)
进一步请参考package "BMS",
Bayesian Model Averaging with BMS.pdf
|