---
title: "ADA2: Homework Ch 07b, ANCOVA 2"
author: "Name Here"
date: "mm/dd/yyyy"
output:
pdf_document:
number_sections: yes
toc: yes
html_document:
toc: true
number_sections: true
code_folding: show
---
```{R, echo=FALSE}
# I set some GLOBAL R chunk options here.
# (to hide this message add "echo=FALSE" to the code chunk options)
knitr::opts_chunk$set(comment = NA, message = FALSE, warning = FALSE, width = 100)
knitr::opts_chunk$set(fig.align = "center", fig.height = 4, fig.width = 6)
#knitr::opts_chunk$set(cache = TRUE, autodep=TRUE)
knitr::opts_chunk$set(cache = FALSE, autodep=TRUE)
```
This is a challenging dataset, in part because it's real and messy.
I will guide you through a simplified sensible analysis,
but other models are possible.
_Note that I needed to set `cache=FALSE` to assure all output was updated._
# ANCOVA model: Albuquerque NM 87108, House and Apartment listing prices
Prof Erhardt constructed a dataset of listing prices for dwellings (homes and apartments) for sale from
[Zillow.com](http://www.zillow.com/homes/for_sale/Albuquerque-NM-87108/95303_rid/any_days/35.095087-106.52167835.035021-106.633258_rect/13_zm/0_mmm/)
on Feb 26, 2016 at 1 PM for Albuquerque NM 87108.
In this assignment we'll develop a model to help understand which qualities that
contribute to a __typical dwelling's listing price__.
We will then also predict the listing prices of new listings
posted on the following day, Feb 27, 2016 by 2 PM.
Because we want to model a _typical dwelling_,
it is completely reasonable to remove "unusual" dwellings from the dataset.
Dwellings are a distribution that has a [long tail](https://en.wikipedia.org/wiki/Long_tail)!
## Unusual assignment, not top-down, but up-down-up-down
This is an unusual assignment because the workflow of this assignment isn't top-down;
instead, you'll be scrolling up and down as you make decisions about the data and model you're fitting.
Yes, I have much of the code worked out for you.
However, there are data decisions to make early in the code
(such as excluding observations, transforming variables, etc.)
that depend on the analysis (model checking) later.
Think of it as a "choose your own adventure" that I've chosen for you.
### Keep a record of your decisions
It is always desirable to make your work reproducible,
either by someone else or by your future self.
For each step you take, keep a diary of
(a) what the next minor goal is,
(b) what evidence/information you have,
(c) what decision you make, and
(d) what the outcome was.
For example, here's the first couple entries of your diary:
1. Include only "typical dwellings". Based on scatterplot, remove extreme observations. Keep only HOUSE and APARTMENT.
2. Exclude a few variables to reduce multicollinearity between predictor variables. Exclude `Baths` and `LotSize`.
## __(1 p)__ (Step 1) Restrict data to "typical" dwellings
__Step 1:__
After looking at the scatterplot below, identify what you consider to be a "typical dwelling"
and exclude observations far from that range.
For example, there are only a couple `TypeSale` that are common enough to model;
remember to run `factor()` again to remove factor levels that no longer appear.
```{R}
# read the data, skip the first two comment lines of the data file
dat.abq <- read.csv("http://statacumen.com/teach/ADA2/homework/ADA2_HW_07_HomePricesZillow_Abq87108.csv"
, skip=2, stringsAsFactors = FALSE)[,-1]
dat.abq$TypeSale <- factor(dat.abq$TypeSale)
# To help scale the intercept to a more reasonable value
# Scaling the x-variables are sometimes done to the mean of each x.
# center year at 1900
dat.abq$YearBuilt_1900 <- dat.abq$YearBuilt - 1900
dat.abq$YearBuilt <- NULL
head(dat.abq)
## RETURN HERE TO SUBSET THE DATA
# dat.abq <- subset(dat.abq, !(X > Z)) # exclude observations of this type
# dat.abq <- dat.abq[!(dat.abq$X > Z),] # exclude observations of this type
# SOLUTION
# these deletions are based only on the scatter plot in order to have
# "typical" dwellings
```
## __(1 p)__ (Step 3) Transform response, if necessary.
__Step 3:__
Does the response variable require a transformation?
If so, what transformation is recommended from the model diagnostic plots (Box-Cox)?
### Solution
[answer]
```{R}
# Price in units of $1000
dat.abq$PriceListK <- dat.abq$PriceList / 1000
dat.abq$PriceList <- NULL
# SOLUTION
```
## __(1 p)__ (Step 4) Remove extremely influential observations.
__Step 4:__
The goal is to develop a model that will work well for the typical dwellings.
If an observation is highly influential, then it's unusual.
```{R}
## Remove influential observation
#dat.abq <- dat.abq[-which(row.names(dat.abq) %in% c(...)),]
# SOLUTION
```
## Subset data for model building and prediction
Create a subset of the data for building the model,
and another subset for prediction later on.
```{R}
# remove observations with NAs
dat.abq <- na.omit(dat.abq)
# the data subset we will use to build our model
dat.sub <- subset(dat.abq, DaysListed > 0)
# the data subset we will predict from our model
dat.pred <- subset(dat.abq, DaysListed == 0)
# the prices we hope to predict well from our model
dat.pred$PriceListK_true <- dat.pred$PriceListK
# set them to NA to predict them later
dat.pred$PriceListK <- NA
```
Scatterplot of the model-building subset.
```{R, fig.height = 8, fig.width = 8}
# NOTE, this plot takes a long time if you're repeadly recompiling the document.
# comment the "print(p)" line so save some time when you're not evaluating this plot.
library(GGally)
library(ggplot2)
p <- ggpairs(dat.sub
, mapping = ggplot2::aes(colour = TypeSale, alpha = 0.5)
, lower = list(continuous = "points")
, upper = list(continuous = "cor")
)
print(p)
```
There are clearly some unusual observations.
Go back to the first code chunk and remove some observations that don't represent a "typical" dwelling.
For example, remove these dwellings (in code above):
* Super expensive dwelling
* Dwellings with huge lots
* Dwellings that were listed for years
* Because most dwellings were APARTMENTs and HOUSEs, remove the others (there are only 1 or so of each).
Discuss the observed correlations or other outstanding features in the data.
### Solution
[answer]
## __(1 p)__ (Step 2) Fit full two-way interaction model.
_You'll revisit this section after each modification of the data above._
__Step 2:__
Let's fit the full two-way interaction model and assess the assumptions.
However, some of the predictor variables are highly correlated.
Recall that the interpretation of a beta coefficient is
"the expected increase in the response for a 1-unit increase in $x$
with all other predictors held constant".
It's hard to hold one variable constant if it's correlated with another variable you're increasing.
Therefore, we'll make a decision to retain some variables but not others
depending on their correlation values.
(In the PCA chapter, we'll see another strategy.)
Somewhat arbitrarily, let's exclude `Baths` (since highly correlated with `Beds` and `Size.sqft`).
Let's also exclude `LotSize` (since highly correlated with `Size.sqft`).
Modify the code below.
Notice that because APARTMENTs don't have more than 1 Beds or Baths,
those interaction terms need to be excluded from the model;
I show you how to do this manually using the `update()` function.
Note that the formula below `y ~ (x1 + x2 + x3)^2` expands into all main effects and two-way interactions.
```{R}
## SOLUTION
lm.full <- lm(PriceListK ~ (TypeSale + Beds + Baths + Size.sqft + LotSize + DaysListed + YearBuilt_1900)^2, data = dat.sub)
#lm.full <- lm(PriceListK ~ (Beds + Baths + Size.sqft + LotSize + DaysListed + YearBuilt_1900)^2, data = dat.sub)
library(car)
try(Anova(lm.full, type=3))
## Note that this doesn't work because APARTMENTs only have 1 bed and 1 bath.
## There isn't a second level of bed or bath to estimate the interaction.
## Therefore, remove those two terms
lm.full <- update(lm.full, . ~ . - TypeSale:Beds - TypeSale:Baths)
library(car)
try(Anova(lm.full, type=3))
```
Note the diagnostic plots below are printed automatically with a new function below (`echo=FALSE`, so look in Rmd file).
Thanks to Eric Kruger for the idea of a diagnostic plotting function.
```{R, echo=FALSE}
lm.diag.plots <- function(fit, rc.mfrow=NA, which.plot=c(1,4,6), outliers=3, sw.order.of.data = FALSE) {
## lm.diag.plots() is a general function for plotting residual diagnostics for an lm() object
## Arguments:
# fit linear model object returned by lm()
# rc.mfrow number of rows and columns for the graphic plot, e.g., c(2,3)
# which.plot default plot numbers for lm()
# outliers number to identify in plots from lm() and qqPlot()
# sw.order.of.data T/F for whether to show residuals by order of data
# variable names
var.names <- names(fit$model)[-1]
# display settings
if (is.na(rc.mfrow[1])) {
rc.mfrow <- c(ceiling((length(var.names) + 3 + length(which.plot) + sw.order.of.data) / 3), 3)
}
op <- par(no.readonly = TRUE) # the whole list of settable par
par(mfrow = rc.mfrow)
# default: Fitted, Cook's distance (with cutoff), and Leverage (with cutoffs)
for(i.plot in which.plot) {
plot(fit, which = i.plot, id.n = outliers)
if (i.plot == 4) {
Di.large <- 4 / (dim(fit$model)[1] - dim(fit$model)[2] - 1)
abline(h = Di.large, col = "blue", lty = 3) # horizontal line
}
if (i.plot == 6) {
lev.large <- c(2, 3) * dim(fit$model)[2] / dim(fit$model)[1]
abline(v = lev.large[1], col = "blue", lty = 3) # horizontal line
abline(v = lev.large[2], col = "blue", lty = 2) # horizontal line
}
}
# Evaluate homoscedasticity
library(car)
# non-constant error variance test
print(ncvTest(fit))
# plot studentized residuals vs. fitted values
try(spreadLevelPlot(fit, sub = "(Homoscedasticity)"))
# Evaluate Collinearity
library(car)
vif.val <- vif(fit) # variance inflation factors
dotchart(vif.val, main = "Collinearity", xlab = "Variance Inflation Factor (VIF)", sub = "Not as useful with interactions")
abline(v = 2^2, col = "blue", lty = 2) # vertical line
# Normal quantile plot (QQ-plot)
library(car)
qqPlot(fit$residuals, las = 1, id.n = outliers, main="QQ Plot", ylab = "Residuals")
# Box-Cox transformation suggestion
# only if all values are positive
if(min(fit$model[,1] > 0)){
library(car)
boxCox(lm.full, lambda = seq(-3,3,length=101), main = "Box-Cox power transformation")
}
# residuals vs order of data
if(sw.order.of.data) {
# order of data (not always interesting)
plot(fit$residuals, main="Residuals vs Order of data", ylab = "Residuals")
abline(h = 0, col = "gray75", lty = 3) # horizontal line at zero
}
# residuals plotted vs each main effect
for(i.plot in 1:length(var.names)) {
m.lab <- paste("Residuals vs.", var.names[i.plot])
plot(fit$model[,var.names[i.plot]], fit$residuals, main=m.lab, ylab = "Residuals", xlab = var.names[i.plot])
abline(h = 0, col = "gray75", lty = 3) # horizontal line at zero
}
par(op) # reset plotting options
## Useful list of diags: http://www.statmethods.net/stats/rdiagnostics.html
} # end of reg.diag.plots()
```
```{R, fig.height = 3, fig.width = 8}
## Uncomment this line when you're ready to assess the model assumptions
#lm.diag.plots(lm.full, rc.mfrow = c(1,3))
```
After Step 2, interpret the residual plots. What are the primary issues?
### Solution
[answer]
## __(1 p)__ (Step 5) Model selection, check model assumptions.
Using `step(..., direction="both")` with the BIC criterion,
perform model selection.
### Solution
```{R, fig.height = 3, fig.width = 8}
## BIC
# option: test="F" includes additional information
# for parameter estimate tests that we're familiar with
# option: for BIC, include k=log(nrow( [data.frame name] ))
lm.red.BIC <- step(lm.full, direction="both", test="F", trace = 0
, k=log(nrow(dat.sub)))
lm.final <- lm.red.BIC
## Uncomment this line when you're ready to assess the model assumptions
#lm.diag.plots(lm.final, rc.mfrow = c(1,3))
```
Model assumptions appear to be reasonably met.
A few influential observations exist.
## __(2 p)__ (Step 6) Plot final model, interpret coefficients.
If you arrived at the same model I did,
then the code below will plot it.
Eventually (after Step 7), the fitted model equations will describe the
each dwelling `TypeSale`
and interpret the coefficients.
```{R, fig.height = 5, fig.width = 8, echo=FALSE}
library(ggplot2)
p <- ggplot(dat.sub, aes(x = Size.sqft, y = PriceListK, colour = TypeSale, shape = TypeSale))
p <- p + geom_point(size = 2, alpha=0.5)
#p <- p + expand_limits(x = 0, y = 8.5)
p <- p + geom_smooth(method = lm, se = FALSE) # , alpha=0.15, fullrange = TRUE)
p <- p + labs(title="Log Listing Price", x = "log10(Size.sqft)")
print(p)
```
```{R}
library(car)
#Anova(lm.final, type=3)
summary(lm.final)
```
Fitted model equation is
$$
\widehat{\log_{10}(\textrm{PriceList})}
=
`r signif(coef(lm.final)[1], 3)`
+ `r signif(coef(lm.final)[2], 3)` I(\textrm{TypeSale} = \textrm{HOUSE})
+ `r signif(coef(lm.final)[3], 3)` \log_{10}(\textrm{Size.sqft})
$$
### Solution
After Step 7, return and intepret the model coefficients above.
[answer]
## __(1 p)__ (Step 7) Transform predictors.
We now have enough information to see that a transformation of a predictor can be useful.
See the curvature with `Size.sqft`?
This is one of the headaches of regression modelling,
_everything depends on everything else_
and you learn as you go.
Return to the top and transform `Size.sqft` and `LotSize`.
A nice feature of this transformation is that the model interaction goes away.
Our interpretation is now on the log scale, but it's a simpler model.
## __(2 p)__ (Step 8) Predict new observations, interpret model's predictive ability.
Using the `predict()` function, we'll input the data we held out to predict earlier,
and use our final model to predict the `PriceListK` response.
Note that `10^lm.pred` is the table of values on the scale of "thousands of dollars".
Interpret the predictions below the output.
How well do you expect this model to predict? Justify your answer.
```{R}
# predict new observations, convert to data frame
lm.pred <- as.data.frame(predict(lm.final, newdata = dat.pred, interval = "prediction"))
# add column of actual list prices
lm.pred$PriceListK <- dat.pred$PriceListK_true
lm.pred
# on "thousands of dollars" scale
10^lm.pred
# attributes of the three predicted observations
dat.pred
```
### Solution
[answer]