Nothing
## ----setup, include=FALSE-----------------------------------------------------
# require(knitr)
require(MASS)
require(dplyr)
## ----install, eval=FALSE------------------------------------------------------
# install.packages("BAS")
## ----devtools, eval=FALSE-----------------------------------------------------
# devtools::install_github("merliseclyde/BAS")
## ----data---------------------------------------------------------------------
data(UScrime, package = "MASS")
## ----transform----------------------------------------------------------------
UScrime[, -2] <- log(UScrime[, -2])
## ----bas----------------------------------------------------------------------
library(BAS)
crime.ZS <- bas.lm(y ~ .,
data = UScrime,
prior = "ZS-null",
modelprior = uniform(), initprobs = "eplogp",
force.heredity = FALSE, pivot = TRUE
)
## ---- fig.show='hold'---------------------------------------------------------
plot(crime.ZS, ask = F)
## ----pip, fig.width=5, fig.height=5-------------------------------------------
plot(crime.ZS, which = 4, ask = FALSE, caption = "", sub.caption = "")
## ----print--------------------------------------------------------------------
crime.ZS
## ----summary------------------------------------------------------------------
options(width = 80)
summary(crime.ZS)
## ----image, fig.width=5, fig.height=5-----------------------------------------
image(crime.ZS, rotate = F)
## ----coef---------------------------------------------------------------------
coef.ZS <- coef(crime.ZS)
## ----plot---------------------------------------------------------------------
plot(coef.ZS, subset = c(5:6), ask = F)
## ----confint-coef-------------------------------------------------------------
confint(coef.ZS)
## ----plot-confint, fig.width=7------------------------------------------------
plot(confint(coef.ZS, parm = 2:16))
## ---- warning=FALSE, fig.width=7---------------------------------------------
plot(confint(coef(crime.ZS, estimator = "HPM")))
## ---- warning=FALSE, fig.width=7, eval=FALSE---------------------------------
# plot(confint(coef(crime.ZS, estimator = "MPM")))
## ----choice of estimator------------------------------------------------------
muhat.BMA <- fitted(crime.ZS, estimator = "BMA")
BMA <- predict(crime.ZS, estimator = "BMA")
# predict has additional slots for fitted values under BMA, predictions under each model
names(BMA)
## ---- fig.width=5, fig.height=5-----------------------------------------------
par(mar = c(9, 9, 3, 3))
plot(muhat.BMA, BMA$fit,
pch = 16,
xlab = expression(hat(mu[i])), ylab = expression(hat(Y[i]))
)
abline(0, 1)
## ----HPM----------------------------------------------------------------------
HPM <- predict(crime.ZS, estimator = "HPM")
# show the indices of variables in the best model where 0 is the intercept
HPM$bestmodel
## -----------------------------------------------------------------------------
variable.names(HPM)
## ----MPM----------------------------------------------------------------------
MPM <- predict(crime.ZS, estimator = "MPM")
variable.names(MPM)
## ----BPM----------------------------------------------------------------------
BPM <- predict(crime.ZS, estimator = "BPM")
variable.names(BPM)
## ---- fig.width=6, fig.height=6-----------------------------------------------
GGally::ggpairs(data.frame(
HPM = as.vector(HPM$fit), # this used predict so we need to extract fitted values
MPM = as.vector(MPM$fit), # this used fitted
BPM = as.vector(BPM$fit), # this used fitted
BMA = as.vector(BMA$fit)
)) # this used predict
## ----se, fig.width=7----------------------------------------------------------
BPM <- predict(crime.ZS, estimator = "BPM", se.fit = TRUE)
crime.conf.fit <- confint(BPM, parm = "mean")
crime.conf.pred <- confint(BPM, parm = "pred")
plot(crime.conf.fit)
plot(crime.conf.pred)
## ----pred---------------------------------------------------------------------
new.pred <- predict(crime.ZS, newdata = UScrime, estimator = "MPM")
## -----------------------------------------------------------------------------
system.time(
for (i in 1:10) {
crime.ZS <- bas.lm(y ~ .,
data = UScrime,
prior = "ZS-null", method = "BAS",
modelprior = uniform(), initprobs = "eplogp"
)
}
)
system.time(
for (i in 1:10) {
crime.ZS <- bas.lm(y ~ .,
data = UScrime,
prior = "ZS-null", method = "deterministic",
modelprior = uniform(), initprobs = "eplogp"
)
}
)
## ----MCMC---------------------------------------------------------------------
crime.ZS <- bas.lm(y ~ .,
data = UScrime,
prior = "ZS-null",
modelprior = uniform(),
method = "MCMC"
)
## ----diagnostics--------------------------------------------------------------
diagnostics(crime.ZS, type = "pip", pch = 16)
diagnostics(crime.ZS, type = "model", pch = 16)
## ----biggerMCMC, eval=FALSE---------------------------------------------------
# crime.ZS <- bas.lm(y ~ .,
# data = UScrime,
# prior = "ZS-null",
# modelprior = uniform(),
# method = "MCMC", MCMC.iterations = 10^6
# )
#
# diagnostics(crime.ZS, type="model", pch=16)
## ----add-out------------------------------------------------------------------
data("stackloss")
stackloss <- cbind(stackloss, diag(nrow(stackloss)))
stack.bas <- bas.lm(stack.loss ~ .,
data = stackloss,
method = "MCMC", initprobs = "marg-eplogp",
prior = "ZS-null",
modelprior = tr.poisson(4, 10),
MCMC.iterations = 200000
)
## -----------------------------------------------------------------------------
knitr::kable(as.data.frame(summary(stack.bas)))
## -----------------------------------------------------------------------------
data(ToothGrowth)
ToothGrowth$dose <- factor(ToothGrowth$dose)
levels(ToothGrowth$dose) <- c("Low", "Medium", "High")
## ---- fig.width=7-------------------------------------------------------------
TG.bas <- bas.lm(len ~ supp*dose,
data = ToothGrowth,
modelprior = uniform(), method = "BAS"
)
image(TG.bas)
## ---- fig.width=7-------------------------------------------------------------
TG.bas <- bas.lm(len ~ supp * dose,
data = ToothGrowth,
modelprior = uniform(), method = "BAS", force.heredity = TRUE
)
image(TG.bas)
## ----force-herid--------------------------------------------------------------
TG.bas <- bas.lm(len ~ supp * dose,
data = ToothGrowth,
modelprior = uniform(), method = "BAS", force.heredity = FALSE
)
TG.herid.bas <- force.heredity.bas(TG.bas)
## ----climate------------------------------------------------------------------
data(climate, package="BAS")
str(climate)
summary(climate)
## ----read-climate-------------------------------------------------------------
library(dplyr)
climate <- filter(climate, proxy != 6) %>%
mutate(proxy = factor(proxy))
## ----wtreg--------------------------------------------------------------------
climate.bas <- bas.lm(deltaT ~ proxy * poly(latitude, 2),
data = climate,
weights = 1 / sdev^2,
prior = "hyper-g-n", alpha = 3.0,
n.models = 2^20,
force.heredity=TRUE,
modelprior = uniform()
)
## ----climate-image, fig.width=7-----------------------------------------------
image(climate.bas, rotate = F)
## ----wtreg-wo-contraint, fig.width=7------------------------------------------
# May take a while to enumerate all 2^20 models
climate.bas <- bas.lm(deltaT ~ proxy * poly(latitude, 2),
data = climate,
weights = 1 / sdev^2,
prior = "hyper-g-n", alpha = 3.0,
n.models = 2^20,
modelprior = uniform(),
force.heredity = FALSE
)
image(climate.bas)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.