---
title: "SDM4 in R: Inferences for Regression (Chapter 25)"
author: "Nicholas Horton (nhorton@amherst.edu) and Sarah McDonald"
date: "June 16, 2018"
output:
pdf_document:
fig_height: 2.8
fig_width: 7
html_document:
fig_height: 3
fig_width: 5
word_document:
fig_height: 4
fig_width: 6
---
```{r, include = FALSE}
# Don't delete this chunk if you are using the mosaic package
# This loads the mosaic and dplyr packages
require(mosaic)
options(digits=3)
```
```{r, include = FALSE}
# knitr settings to control how R chunks work.
require(knitr)
opts_chunk$set(
tidy = FALSE, # display code as typed
size = "small" # slightly smaller font for code
)
```
## Introduction and background
This document is intended to help describe how to undertake analyses introduced
as examples in the Fourth Edition of \emph{Stats: Data and Models} (2014) by De Veaux, Velleman, and Bock.
More information about the book can be found at http://wps.aw.com/aw_deveaux_stats_series. This
file as well as the associated R Markdown reproducible analysis source file used to create it can be found at http://nhorton.people.amherst.edu/sdm4.
This work leverages initiatives undertaken by Project MOSAIC (http://www.mosaic-web.org), an NSF-funded effort to improve the teaching of statistics, calculus, science and computing in the undergraduate curriculum. In particular, we utilize the `mosaic` package, which was written to simplify the use of R for introductory statistics courses. A short summary of the R needed to teach introductory statistics can be found in the mosaic package vignettes (http://cran.r-project.org/web/packages/mosaic).
A paper describing the mosaic approach was published in the *R Journal*: https://journal.r-project.org/archive/2017/RJ-2017-024.
## Chapter 25: Inferences for Regression
### Section 25.1: The population and the sample
```{r message = FALSE}
library(mosaic)
library(readr)
BodyFat <- read_csv("http://nhorton.people.amherst.edu/sdm4/data/Body_fat_complete.csv")
dim(BodyFat)
glimpse(BodyFat)
```
We can confirm the coefficients from the model on page 690.
```{r}
BodyFatmod <- lm(PctBF ~ waist, data = BodyFat)
coef(BodyFatmod)
```
### Section 25.2: Assumptions and conditions
We can regenerate the output and figures for the example on pages 692-696.
```{r}
msummary(BodyFatmod)
rsquared(BodyFatmod)
confint(BodyFatmod) # see page 700
```
```{r message = FALSE}
# Figure 25.4
gf_point(PctBF ~ waist, xlab = "Waist (in.)",
data = BodyFat) %>% # see smoothers on p.92-93
gf_lm() %>%
gf_smooth(se = FALSE)
# Figure 25.5
gf_point(resid(BodyFatmod) ~ waist, xlab = "Waist (in.)",
data = BodyFat) %>%
gf_lm() %>%
gf_smooth(se = FALSE)
# equiv of Figure 25.6 note that Figure 25.6 refers to the diamonds dataset
gf_point(resid(BodyFatmod) ~ fitted(BodyFatmod), xlab = "Predicted values",
ylab = "Residuals",
type = c("p", "r", "smooth"), data = BodyFat) %>%
gf_lm() %>%
gf_smooth(se = FALSE)
# Figure on bottom of page 695
gf_qq(~ resid(BodyFatmod)) %>%
gf_qqline()
```
#### Section 25.6: Confidence intervals for predicted values
We can reproduce Figure 25.12 (page 707) using layers in ggformula.
```{r}
library(broom)
gf_point(PctBF ~ waist, xlab = "Waist (in.)",
panel = panel.lmbands, lwd = 2, cex = 0.2, data = BodyFat) %>%
gf_lm(interval = "confidence", fill = "red") %>%
gf_lm(interval = "prediction", fill = "navy")
```
```{r}
Craters <- read.csv("http://nhorton.people.amherst.edu/sdm4/data/Craters.csv")
dim(Craters)
Craters <- mutate(Craters,
logDiam = log(Diam.km.),
logAge = log(age..Ma.))
Cratermod <- lm(logDiam ~ logAge, data = Craters)
favstats(~ logAge, data = Craters) # note example in book has n=39
```
```{r}
confpred <- predict(Cratermod, interval = "confidence")
intpred <- predict(Cratermod, interval = "prediction")
select(Craters, -Name) %>%
head(., 3)
head(confpred, 3)
head(intpred, 3)
```
#### Section 25.7: Logistic regression
The Pima Indian dataset example is given on pages 708-712.
```{r message = FALSE}
Pima <- read_csv("http://nhorton.people.amherst.edu/sdm4/data/Pima_Indians_Diabetes.csv")
Diabetes <- filter(Pima, BMI>0) # get rid of missing values for BMI
gf_boxplot(BMI ~ as.factor(Diabetes), data = Pima)
```
```{r fig.keep = "last"}
pimamod <- glm(Diabetes ~ BMI, family = "binomial", data = Pima)
f2 <- makeFun(pimamod)
gf_point(Diabetes ~ BMI, data = Pima) %>%
gf_function(f2, add = TRUE)
msummary(pimamod)
```