############################################################# ## A BMA example: Attitude data ############################################################# data("attitude", package = "datasets") library("BMS") att <- bms(attitude, mprior = "uniform", g = "UIP", user.int = FALSE) coef(att) coef(att, std.coefs = TRUE, order.by.pip = FALSE, include.constant = TRUE) summary(att) topmodels(att)[, 1:3] image(att) ############################################################# ## Model size and model priors ############################################################# sum(coef(att)[, 1]) plotModelsize(att) att_fixed <- bms(attitude, mprior = "fixed", mprior.size = 2, user.int = TRUE) plotModelsize(att_fixed) att_pip <- bms(attitude, mprior = "pip", mprior.size = c(0.01, 0.5, 0.5, 0.5, 0.5, 0.5), user.int = FALSE) summary(att_pip) plotModelsize(att_pip) att_random <- bms(attitude, mprior = "random", mprior.size = 3, user.int = FALSE) plotModelsize(att_random) plotComp(Uniform = att, fixed = att_fixed, PIP = att_pip, Random = att_random, cex = 2) ############################################################# ## MCMC samplers and more variables ############################################################# data("datafls", package = "BMS") fls1 <- bms(datafls, burn = 50000, iter = 1e+05, g = "BRIC", mprior = "uniform", nmodel = 2000, mcmc = "bd", user.int = FALSE) summary(fls1) plotConv(fls1) plotConv(fls1[1:100]) pmp(fls1)[1:5, ] colSums(pmp(fls1)) round(colSums(pmp(fls1))[[2]], 2) * 100 coef(fls1)[1:5, ] coef(fls1, exact = TRUE)[1:5, ] fls2 <- bms(datafls, burn = 20000, iter = 50000, g = "BRIC", mprior = "uniform", mcmc = "rev.jump", start.value = 0, user.int = FALSE) summary(fls2) round(cor(pmp(fls2))[2, 1], 2) round(cor(pmp(fls1))[2, 1], 2) fls_combi <- c(fls1, fls2) summary(fls_combi) ############################################################# ## Alternative formulations for Zellner's g prior ############################################################# fls_g5 <- bms(datafls, burn = 20000, iter = 50000, g = 5, mprior = "uniform", user.int = FALSE) coef(fls_g5)[1:5, ] summary(fls_g5) fls_ebl <- bms(datafls, burn = 20000, iter = 50000, g = "EBL", mprior = "uniform", nmodel = 1000, user.int = FALSE) summary(fls_ebl) fls_hyper <- bms(datafls, burn = 20000, iter = 50000, g = "hyper=UIP", mprior = "random", mprior.size = 7, nmodel = 1000, user.int = FALSE) summary(fls_hyper) gdensity(fls_hyper) round(fls_hyper$gprior.info$shrinkage.moments[[1]] - as.numeric(strsplit(summary(fls_hyper)[13], "=")[[1]][[3]]), 2) image(fls_hyper) density(fls_combi, reg = "Muslim") round(estimates(fls_combi, exact = TRUE)["Muslim", 1], 3) coef(fls_combi, exact = TRUE, condi.coef = TRUE)["Muslim", ] dmuslim <- density(fls_hyper, reg = "Muslim", addons = "Eebl") quantile(dmuslim, c(0.025, 0.975)) ############################################################# ## Predictive densities ############################################################# fcstbma <- bms(datafls[1:70, ], mprior = "uniform", burn = 20000, iter = 50000, user.int = FALSE) pdens <- pred.density(fcstbma, newdata = datafls[71:72, ]) plot(pdens, 2) quantile(pdens, c(0.05, 0.95)) pdens$dyf(datafls[71:72, 1]) plot(pdens, "ZM", realized.y = datafls["ZM", 1]) lps(pdens, datafls[71:72, 1]) ############################################################# ## Appendix ############################################################# data("attitude", package = "datasets") att_full <- zlm(attitude, g = "UIP") summary(att_full) att_best <- as.zlm(att, model = 1) summary(att_best) att_bestlm <- lm(model.frame(as.zlm(att))) summary(att_bestlm) att_learn <- bms(attitude, mprior = "uniform", fixed.reg = c("complaints", "learning")) fls_culture <- bms(datafls, fixed.reg = c(1, 8:16, 24, 26:41), mprior = "random", mprior.size = 28, mcmc = "enumeration", user.int = FALSE) coef(fls_culture)[28:41, ] plotModelsize(fls_culture, ksubset = 27:41)