p7.R
In vibass: Materials for Introductory Course on Bayesian Learning

## -----------------------------------------------------------------------------
library("MASS")
data("nlschools")
summary(nlschools)

## ----message = FALSE, warning = FALSE-----------------------------------------
library("R2BayesX")
m1 <- bayesx(lang ~ IQ + GS +  SES + COMB, data = nlschools)

summary(m1)

## ----fig = TRUE, fig.width = 15, fig.height = 5-------------------------------
boxplot(lang ~ class, data = nlschools, las = 2)

## -----------------------------------------------------------------------------
m2 <- bayesx(
  lang ~ IQ + GS +  SES + COMB + sx(class, bs = "re"),
  data = nlschools
)

summary(m2)

## -----------------------------------------------------------------------------
library(lme4)
m2lmer <- lmer(
  lang ~ IQ + GS +  SES + COMB + (1 | class), # Fit a separate intercept for each level of class
  data = nlschools
)

summary(m2lmer)

## ----message = FALSE, warning = FALSE-----------------------------------------
library(spData)
data(nc.sids)
summary(nc.sids)

## -----------------------------------------------------------------------------
# Overall mortality rate
r74 <- sum(nc.sids$SID74) / sum(nc.sids$BIR74)
# Expected cases
nc.sids$EXP74 <- r74 * nc.sids$BIR74

## -----------------------------------------------------------------------------
nc.sids$SMR74 <- nc.sids$SID74 / nc.sids$EXP74

## ----fig = TRUE---------------------------------------------------------------
hist(nc.sids$SMR, xlab = "SMR")

## -----------------------------------------------------------------------------
nc.sids$NWPROP74 <- nc.sids$NWBIR74/ nc.sids$BIR74

## ----fig = TRUE---------------------------------------------------------------
plot(nc.sids$NWPROP74, nc.sids$SMR74)

# Correlation
cor(nc.sids$NWPROP74, nc.sids$SMR74)

## -----------------------------------------------------------------------------
m1nc <- bayesx(
  SID74 ~ 1 + NWPROP74,
  family = "poisson",
  offset = log(nc.sids$EXP74),
  data = nc.sids
)
summary(m1nc)

## -----------------------------------------------------------------------------
# Index for random effects
nc.sids$ID <- seq_len(nrow(nc.sids))

# Model WITH covariate
m2nc <- bayesx(
  SID74 ~  1 + NWPROP74 + sx(ID, bs = "re"),
  family = "poisson",
  offset = log(nc.sids$EXP74),
  data = as.data.frame(nc.sids)
)

summary(m2nc)

## ----fig = TRUE---------------------------------------------------------------
x.predict <- seq(0,1,length=1000)
y.predict <- exp(coef(m2nc)["(Intercept)","Mean"]+coef(m2nc)["NWPROP74","Mean"]*x.predict)
oldpar <- par(mfrow = c(1, 1))
plot(nc.sids$NWPROP74, nc.sids$SMR74)
lines(x.predict,y.predict)
par(oldpar)

## -----------------------------------------------------------------------------
# Model WITHOUT covariate
m3nc <- bayesx(
  SID74 ~  1 + sx(ID, bs = "re"),
  family = "poisson",
  offset = log(nc.sids$EXP74),
  data = as.data.frame(nc.sids)
)

summary(m3nc)

## ----fig = TRUE---------------------------------------------------------------
oldpar <- par(mfrow = c(1, 2))
boxplot(m2nc$effects$`sx(ID):re`$Mean, ylim = c(-1, 1), main = "With NWPROP74")
boxplot(m3nc$effects$`sx(ID):re`$Mean, ylim = c(-1, 1), main = "Without NWPROP74")
par(oldpar)