## ----echo=FALSE----------------------------------------------------------
library(knitr)
opts_chunk$set(size="small",background="white", highlight=FALSE, cache=TRUE, autdep=TRUE, tidy=TRUE, options(show.signif.stars=FALSE))

## ----eval=FALSE----------------------------------------------------------
## df$x.sq <- df$x^2
## lm(y~x+x.sq, data=df)

## ----eval=FALSE----------------------------------------------------------
## lm(y ~ x + I(x^2), data=df)

## ----eval=FALSE----------------------------------------------------------
## lm(y ~ poly(x,2), data=df)

## ------------------------------------------------------------------------
# Load the data
mobility <- read.csv("http://www.stat.cmu.edu/~cshalizi/mreg/15/dap/1/mobility.csv")
mob.quad <- lm(Mobility ~ Commute + poly(Latitude,2)+Longitude, data=mobility)

## ------------------------------------------------------------------------
summary(mob.quad)

## ------------------------------------------------------------------------
predict(mob.quad, newdata=data.frame(Commute=0.298,
                      Latitude=40.57, Longitude=-79.58))

## ----hypothetical-pghs, echo=FALSE---------------------------------------
hypothetical.pghs <- data.frame(Commute=0.287,
                                Latitude=seq(from=min(mobility$Latitude),
                                  to=max(mobility$Latitude), length.out=100),
                                Longitude=-79.92)
plot(hypothetical.pghs$Latitude, predict(mob.quad, newdata=hypothetical.pghs),
     xlab="Latitude", ylab="Expected mobility", type="l")

## ----eval=FALSE----------------------------------------------------------
## hypothetical.pghs <- data.frame(Commute=0.287,
##                                 Latitude=seq(from=min(mobility$Latitude),
##                                   to=max(mobility$Latitude), length.out=100),
##                                 Longitude=-79.92)
## plot(hypothetical.pghs$Latitude, predict(mob.quad, newdata=hypothetical.pghs),
##      xlab="Latitude", ylab="Expected mobility", type="l")

## ------------------------------------------------------------------------
library(MASS)
data(cats)
Hwt.lm <- lm(Hwt ~ Sex+Bwt, data=cats)
summary(Hwt.lm)

## ------------------------------------------------------------------------
# How many levels does State have?
nlevels(mobility$State)
# What are they?
levels(mobility$State)

## ------------------------------------------------------------------------
mob.state <- lm(Mobility ~ Commute + State, data=mobility)
signif(coefficients(mob.state),3)

## ------------------------------------------------------------------------
# Set up a function to make maps
  # "Terrain" color levels set based on quantiles of the variable being plotted
# Inputs: vector to be mapped over the data frame; number of levels to
  # use for colors; other plotting arguments
# Outputs: invisibly, list giving cut-points and the level each observation
  # was assigned
mapper <- function(z, levels, ...) {
    # Which quantiles do we need?
    probs <- seq(from=0, to=1, length.out=(levels+1))
    # What are those quantiles?
    z.quantiles <- quantile(z, probs)
    # Assign each observation to its quantile
    z.categories <- cut(z, z.quantiles, include.lowest=TRUE)
    # Make up a color scale
    shades <- terrain.colors(levels)
    plot(x=mobility$Longitude, y=mobility$Latitude,
         col=shades[z.categories], ...)
    invisible(list(quantiles=z.quantiles, categories=z.categories))
}

## ----basic-residuals, echo=FALSE-----------------------------------------
residuals.basic.map <- mapper(residuals(lm(Mobility~Commute, data=mobility)),
                             levels=4, pch=19,
                              cex=0.5, xlab="Longitude", ylab="Latitude",
                              main="Mobility")
legend("topright", legend=levels(residuals.basic.map$categories), pch=19,
       col=terrain.colors(4), cex=0.8)

## ----eval=FALSE----------------------------------------------------------
## 

## ----state-level-residuals, echo=FALSE-----------------------------------
residuals.state.map <- mapper(residuals(mob.state), levels=4, pch=19,
                              cex=0.5, xlab="Longitude", ylab="Latitude",
                              main="Mobility")
legend("topright", legend=levels(residuals.state.map$categories), pch=19,
       col=terrain.colors(4), cex=0.8)

## ----eval=FALSE----------------------------------------------------------
## residuals.state.map <- mapper(residuals(mob.state), levels=4, pch=19,
##                               cex=0.5, xlab="Longitude", ylab="Latitude",
##                               main="Mobility")
## legend("topright", legend=levels(residuals.state.map$categories), pch=19,
##        col=terrain.colors(4), cex=0.8)

## ------------------------------------------------------------------------
# The states of the Confederacy
Confederacy <- c("AR", "AL", "FL", "GA", "LA", "MS", "NC", "SC", "TN", "TX", "VA")
mobility$Dixie <- mobility$State %in% Confederacy

## ------------------------------------------------------------------------
mob.dixie <- lm(Mobility ~ Commute + Dixie, data=mobility)
signif(coefficients(summary(mob.dixie)),3)

## ----residuals-dixie, echo=FALSE-----------------------------------------
residuals.dixie.map <- mapper(residuals(mob.dixie), levels=4, pch=19,
                              cex=0.5, xlab="Longitude", ylab="Latitude",
                              main="Mobility")
legend("topright", legend=levels(residuals.dixie.map$categories), pch=19,
       col=terrain.colors(4), cex=0.8)

## ----eval=FALSE----------------------------------------------------------
## residuals.dixie.map <- mapper(residuals(mob.dixie), levels=4, pch=19,
##                               cex=0.5, xlab="Longitude", ylab="Latitude",
##                               main="Mobility")
## legend("topright", legend=levels(residuals.dixie.map$categories), pch=19,
##        col=terrain.colors(4), cex=0.8)