rationale.R
In MASSExtra: Some 'MASS' Enhancements

## ----presetup, include = FALSE------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE,
                      comment = "",
                      message = FALSE,
                      warning = FALSE)
old_plot_new <- getHook("plot.new")
old_before_plot_new <- getHook("before.plot.new")
setHook("plot.new", list(), "replace")
setHook("before.plot.new", rev(old_before_plot_new)[1], "replace")
old_par <- par(no.readonly = TRUE)
old_options <- options(show.signif.stars = FALSE)

library(ggplot2)
old_theme <- theme_get()

## ----packages_etc-------------------------------------------------------------
suppressPackageStartupMessages({
  library(visreg)
  library(knitr)
  library(dplyr)
  library(ggplot2)
  library(patchwork)
  library(MASSExtra)
})
options(knitr.kable.NA = "")
theme_set(theme_bw() + theme(plot.title = element_text(hjust = 0.5)))

## ----setup, fig.height = 5, fig.width = 10, out.width="100%", fig.cap="Box-cox, old and new displays"----
par(mfrow = c(1, 2))
mod0 <- lm(MPG.city ~ Weight, Cars93)
boxcox(mod0)  ## MASS
box_cox(mod0) ## MASSExtra tweak

## ---- fig.height=5, fig.width=10, out.width="100%", fig.cap="The Box-Cox transformation effect"----
p0 <- ggplot(Cars93) + aes(x = Weight) + geom_point(colour = "#2297E6")  + xlab("Weight (lbs)") +
  geom_smooth(se = FALSE, method = "loess", formula = y ~ x, size=0.7, colour = "black")
p1 <- p0 + aes(y = MPG.city) + ylab("Miles per gallon (MPG)") + ggtitle("Untransformed response")
p2 <- p0 + aes(y = bc(MPG.city, lambda(mod0))) + ggtitle("Transformed response") + 
  ylab(bquote(bc(MPG, .(round(lambda(mod0), 2)))))
p1 + p2

## ---- fig.width=8, fig.height=6, fig.align="center", out.width="75%"----------
big_model <- lm(medv ~ . + (rm + tax + lstat + dis)^2 + poly(dis, 2) + poly(rm, 2) +
                   poly(tax, 2) + poly(lstat, 2), Boston)
big_model %>% drop_term(k = "GIC") %>% plot() %>% kable(booktabs=TRUE, digits=3)

## ---- fig.width=8, fig.height=6, fig.align="center", out.width="75%"----------
base_model <- lm(medv ~ ., Boston)
gic_model <- step_GIC(base_model, scope = list(lower = ~1, upper = formula(big_model)))
drop_term(gic_model) %>% plot() %>% kable(booktabs = TRUE, digits = 3)

## ---- fig.width=14, fig.height=6, fig.align="center", out.width="100%", fig.cap="Two component visualisations for the Boston house price model"----
capture.output(suppressWarnings({
   g1 <- visreg(gic_model, "dis", plot = FALSE)
   g2 <- visreg(gic_model, "lstat", plot = FALSE)
   plot(g1, gg = TRUE) + plot(g2, gg = TRUE) 
})) -> void

## ---- fig.width=8, fig.height=6, fig.align="center", out.width="75%"----------
add_term(base_model, scope = formula(big_model), k = "gic") %>% 
   plot() %>% kable(booktabs = TRUE, digits = 3)

## -----------------------------------------------------------------------------
quine_full <- glm.nb(Days ~ Age*Eth*Sex*Lrn, data = quine)
drop_term(quine_full) %>% kable(booktabs = TRUE, digits = 4)
quine_gic <- step_GIC(quine_full)
drop_term(quine_gic) %>% kable(booktabs = TRUE, digits = 4)

## -----------------------------------------------------------------------------
mpg0 <- glm(MPG.city ~ Weight + Cylinders + EngineSize + Origin, 
            family = Gamma(link = "inverse"), data = Cars93)
drop_term(mpg0) %>% kable(booktabs = TRUE, digits = 3)
mpg_gic <- step_down(mpg0, k = "gic")      ## simple backward elimination, mainly used for GLMMs
drop_term(mpg_gic) %>% kable(booktabs = TRUE, digits = 3)

## ---- fig.width=10, fig.height=7, fig.align="center", out.width="80%", fig.cap="MPG model using a Gamma family with inverse link"----
ggplot(Cars93) + aes(x = Weight, y = MPG.city, colour = Cylinders) + geom_point() +
   geom_smooth(method = "glm", method.args = list(family = Gamma), formula = y ~ x) +
   ylab("Miles per gallon (city driving)") + scale_colour_brewer(palette = "Dark2")

## ---- fig.width=10, fig.height=7, fig.align="center", out.width="80%", fig.cap="A one-dimensional KDE of a predictor variable from the Boston data set"----

Criminality <- with(Boston, log(crim))
kcrim <- kde_1d(Criminality, n = 1024, kernel = demoKde::kernelBiweight)
kcrim
plot(kcrim)

## ---- fig.width=10, fig.height=7, fig.align="center", out.width="80%", fig.cap="A two-dimensional KDE using two predictors from the Boston data set"----
Spaciousness <- with(Boston, sqrt(rm))
kcrimrm <- kde_2d(Criminality, Spaciousness, n = 512, kernel = "opt")
kcrimrm
plot(kcrimrm, ## col = hcl.colors(25, rev = TRUE),
     xlab = expression(italic(Criminality)),
     ylab = expression(italic(Spaciousness)))
contour(kcrimrm, col = "dark green", add = TRUE)

## ---- fig.width=10, fig.height=7, fig.align="center", out.width="80%", fig.cap="A two-dimensional KDE presented as a perspective plot"----
with(kcrimrm, persp(x, 10*y, 3*z, border="transparent", col = "powder blue",
                    theta = 30, phi = 15, r = 100, scale = FALSE, shade = TRUE, 
                    xlab = "Criminality", ylab = "Spaciousness", zlab = "kde"))

## ---- echo=FALSE--------------------------------------------------------------
setdiff(getNamespaceExports("MASS"), c("lmwork", getNamespaceExports("MASSExtra"))) %>% 
  sort() %>%  noquote()

## ---- echo=FALSE--------------------------------------------------------------
intersect(getNamespaceExports("MASS"), getNamespaceExports("MASSExtra")) %>% sort() %>%  noquote()

## ---- echo=FALSE--------------------------------------------------------------
setdiff(getNamespaceExports("MASSExtra"), getNamespaceExports("MASS")) %>% 
   setdiff(getNamespaceExports("splines")) %>% 
   grep("^[.]__", ., invert = TRUE, value = TRUE) %>% ## exclude S4 class objects
   sort() %>% noquote()

## ---- echo=FALSE--------------------------------------------------------------
intersect(getNamespaceExports("MASSExtra"), getNamespaceExports("splines")) %>% sort() %>%  noquote()

## ----cleanup, include=FALSE---------------------------------------------------
setHook("plot.new", old_plot_new, "replace")
setHook("before.plot.new", old_before_plot_new, "replace")
par(old_par)
options(old_options)
theme_set(old_theme)