## Chapter 13
##
## Code to support Video Vignette
##
## Videos and supporting code are not a complete portrayal
## of Chapter content.

##
## a: log and log-log transforms
##

## data on used pickups on Craigslist
pickups <- read.csv("pickups.csv")
fit <- lm(price ~ year, data=pickups)
## will use lm library routine rather than Ch7 longhand throughout

## not an ideal fit
plot(pickups$year, pickups$price, xlab="year", ylab="price")
abline(fit, lwd=2)
legend("topleft", "least squares line", lwd=2, bty="n")

## possibly easier to see as residuals versus fitted values
## (a popular diagnostic)
plot(fit$fitted, fit$resid, 
     xlab="fitted values (yhat)", ylab="residuals (ei)")
abline(h=0, col=2, lty=2, lwd=2)
legend("topleft", "zero residual", col=2, lty=2, lwd=2, bty="n")

## look how the y-values (price) are bunched up on the y-axis
## (that means more small y-values than large ones)
## taking the logarithm can undo that bunching
pickups$lprice <- log(pickups$price)     ## augment pickups df
lfit <- lm(lprice ~ year, data=pickups)
## square root also works

## much better fit
plot(pickups$year, pickups$lprice, xlab="year", ylab="log price")
abline(lfit, lwd=2)
legend("topleft", "least squares line", lwd=2, bty="n")

## again possibly a cleaner visual
plot(lfit$fitted, lfit$resid, 
     xlab="fitted values (log yhat)", ylab="residuals (ei)")
abline(h=0, col=2, lty=2, lwd=2)
legend("bottomleft", "zero residual", col=2, lty=2, lwd=2, bty="n")

## a grid of predictive locations and their predictions
xp <- seq(min(pickups$year), max(pickups$year), length=1000)
p <- predict(lfit, newdata=data.frame(year=xp), interval="prediction")
## again we will use libraries rather than Ch7

## putting everything back on the original (untransformed) scale
plot(pickups$year, pickups$price, ylim=range(exp(p)), 
     xlab="year", ylab="price")
lines(xp, exp(p[,1]), lwd=2)                ## exp undoes log
lines(xp, exp(p[,2]), col=2, lty=2, lwd=2)
lines(xp, exp(p[,3]), col=2, lty=2, lwd=2)
legend("topleft", c("mean", "95% PI"), col=1:2, lty=1:2, lwd=2, bty="n")

## data on canned food sales from Consolidated GFOods
confood <- read.csv("confood.csv")

## doesn't look like a good candidate for straight SLR
plot(confood$price, confood$sales, 
     xlab="price ($)", ylab="sales volume (units)")

## but how about if you log both variabes?
lprice <- log(confood$price)
lsales <- log(confood$sales)
plot(lprice, lsales, xlab="log price", ylab="log sales")
legend("topright", "log-log fit", lwd=2, bty="n")

## pretty good fit
lprice <- log(confood$price)
lsales <- log(confood$sales)
llfit <- lm(lsales ~ lprice)
abline(llfit, lwd=2)

## testing hypotheses
summary(llfit)
## easy reject of elasticity of zero

## how about testing for 5% elasticity
tstat <- (-5.1477 - (-5))/0.5098
2*pt(-abs(tstat), df=length(lprice)-2)
## can't reject

## when predicting, make sure to log x on the way in
xp <- seq(min(confood$price), max(confood$price), length=1000)
p <- predict(llfit, newdata=data.frame(lprice=log(xp)), 
             interval="prediction")

## and exp y on the way out
plot(confood$price, confood$sales, ylim=range(exp(p)), 
     xlab="price ($)", ylab="sales volume (units)")
lines(xp, exp(p[,1]), lwd=2)
lines(xp, exp(p[,2]), col=2, lty=2, lwd=2)
lines(xp, exp(p[,3]), col=2, lty=2, lwd=2)
legend("topright", c("mean", "95% PI"), col=1:2, lty=1:2, lwd=2, bty="n")

##
## b: polynomial regression
##

## data on Soybean growing time versus harvest weight
library(nlme)    ## don't forget install.packages("nmle")
data(Soybean)
soy <- Soybean[Soybean$Time < 65 & Soybean$Variety == "P",]

## visualizing
plot(soy$Time, soy$weight, xlab="time (days)", ylab="weight (g)")

## there is clearly bunching on the y-axis, and we know
## what to do when that happens
x <- soy$Time
y <- log(soy$weight) 

## (log) linear model fit
lfit <- lm(y ~ x)
xp <- seq(min(x), max(x), length=1000)
p <- predict(lfit, newdata=data.frame(x=xp), interval="prediction")

## augment plot on original scale with fit
lines(xp, exp(p[,1]), lwd=2)
lines(xp, exp(p[,2]), col=2, lty=2, lwd=2)
lines(xp, exp(p[,3]), col=2, lty=2, lwd=2)
legend("topleft", c("mean", "95% PI"), col=1:2, lty=1:2, lwd=2, bty="n")

## looks good but a residual analysis indicates missing structure
plot(lfit$fitted, lfit$residual,
     xlab="fitted values (yhat)", ylab="residuals (e)")
abline(h=0, col=2, lty=2, lwd=2)
legend("bottom", "zero residual", col=2, lty=2, lwd=2, bty="n")

## try a degree-2 polynomial (for the log response)
x2 <- x^2
lfit2 <- lm(y ~ x + x2)
summary(lfit2)
## seems to be a good idea

## predicting
xp2 <- xp^2
p2 <- predict(lfit2, newdata=data.frame(x=xp, x2=xp2), 
              interval="prediction")

## inspecting the new fit in log space
## left panel, view in log y space
plot(x, y)
lines(xp, p2[,1], lwd=2)
lines(xp, p2[,2], col=2, lty=2, lwd=2)
lines(xp, p2[,3], col=2, lty=2, lwd=2)

## now back on the original, un-logged scale
plot(soy$Time, soy$weight,
     xlab="time (days)", ylab="weight (g)")
lines(xp, exp(p2[,1]), lwd=2)
lines(xp, exp(p2[,2]), col=2, lty=2, lwd=2)
lines(xp, exp(p2[,3]), col=2, lty=2, lwd=2)
legend("topleft", c("mean", "95% PI"), col=1:2, lty=1:2, lwd=2, bty="n")

## it could be that an even higher degree will lead to
## an even better fit
x3 <- x^3
lfit3 <- lm(y ~ x + x2 + x3)
summary(lfit3)
## Nope

##
## c: multiple linear regression
##

## continuing the polynomial regression example
X <- cbind(1, x, x2)

## calculating MLE for beta
XtX <- t(X) %*% X             ## t is transpose, %*% is matrix product
XtXi <- solve(XtX)            ## solve is inv with only one argument
Xty  <- t(X) %*% y            ## %*% also for matrix-vector product
bhat <- drop(XtXi %*% Xty)    ## drop from a 1xp matrix to p-vector
rbind(coef(lfit2), bhat)

## MLE for sigma2
p <- ncol(X)
n <- nrow(X)
s2 <- sum((y - X %*% bhat)^2)/(n - p)
c(hand=s2, lib=summary(lfit2)$sigma^2)

## completing the summary table
se <- sqrt(s2*diag(XtXi))
t <- bhat/se
pval <- 2*pt(-abs(t), n - p)
tab <- cbind(bhat, se, t, pval)
colnames(tab) <- c("est", "stderr", "t stat", "p val")
rownames(tab) <- c("(Intercept)", "x", "x2")

## compare
tab
summary(lfit2)

##
## d: dummies and interactions
##

## dummy variables

## a categorical predictor, back to pickup trucks
table(pickups$make)

## building our own dummy variables
X <- cbind(1, pickups$year, pickups$make == "Ford", pickups$make == "GMC")
y <- log(pickups$price)

## MLE
bhat <- drop(solve(t(X) %*% X) %*% t(X) %*% y)

## or R will make them for you
fit <- lm(lprice ~ year + make, data=pickups)
rbind(coef(fit), bhat)

## interactions

## data on grades in an MBA program
grades <- read.csv("grades.csv")

## MBA GPA depends on Bachelors GPA
fit1 <- lm(MBAGPA ~ BachGPA, data=grades)
summary(fit1)
## yes

## what about also including age?
fit2 <- lm(MBAGPA ~ BachGPA + Age, data=grades)
summary(fit2)
## yes, but does it make sense that older people get lower grades?

## how about an interaction term
fit3 <- lm(MBAGPA ~ BachGPA*Age, data=grades)
## or lm(MBAGPA ~ BachGPA + Age + BachGPA:Age, data=grades)
summary(fit3)
## this doesn't look like a better fit?

## lets get rid of the age main effect, since it has the 
## highest p-value
fit4 <- lm(MBAGPA ~ BachGPA*Age - Age, data=grades)
## same as lm(MBAGPA ~ BachGPA + BachGPA:Age)
summary(fit4)
## everything is useful
## apparently, the older you are the less your Bachelors GPA matters

##
## e: model selection
##

## read in data about worker's satisfaction with boss
boss <- read.csv("boss.csv")

## predictors (like fair raises, opportunity to learn new things, 
## special privlidges) are highly "multicollinear" (more in book)
cor(boss[,3:5])

## full MLR fit
full <- lm(Y ~ ., data=boss)  ## shorthand . for X1 + ... + x6
summary(full)
## none of the slopes look statistically significant

## but if you drop several collinear predictors ...
base <- lm(Y ~ X1 + X2, data=boss)
summary(base)
## ... they look useful

## not much difference in R-squared
r2f <- summary(full)$r.squared   ## same as cor(boss$y, full$fitted)
r2b <- summary(base)$r.squared   ## same as cor(boss$y, base$fitted)
r2diff <- r2f - r2b
r2diff

## is a difference of 2.5% worth an extra ...
pf <- length(coef(full))
pb <- length(coef(base))
pdiff <- pf - pb
pdiff
## ... four (seemingly useless) predictors?

## some setup to simulate from H0
X <- as.matrix(boss[,-1])  ## discard the y column
n <- nrow(X)
Dfs <- boss

## more setup
beta <- c(as.numeric(coef(base)), rep(0, 4))  ## coefficients under H0
mu <- beta[1] + X %*% beta[-1]                ## mean/prediction under H0
s2 <- summary(base)$sigma^2                   ## est. uncertainty under H0

## simulate
N <- 100000
R2fs <- R2bs <- R2diffs <- rep(NA, N)
for(i in 1:N) {
  
  ## simulate data under H0 (via residuals)
  sigma2 <- s2  ## or sigma2 <- (n - pb)*s2/rchisq(1, n - pb)
  Ys <- rnorm(n, mu, sqrt(sigma2))
  
  ## put simulated Ys along with the data Xs
  Dfs$Y <- Ys
  
  ## and then fit both models to these simulated data
  Fs <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6, data=Dfs)     # H1
  Bs <- lm(Y ~ X1 + X2, data=Dfs)                         # H0
  
  ## extract both R2 values
  R2fs[i] <- summary(Fs)$r.squared
  R2bs[i] <- summary(Bs)$r.squared
  
  ## save their difference
  R2diffs[i] <- R2fs[i] - R2bs[i]
}

## visualize (simulation isn't speedy)
hist(R2diffs, main="", xlab="sampled R2diff values")
abline(v=r2diff, col=2, lwd=2)
legend("right", "observed", col=2, lwd=2, bty="n")

## p-value
mean(R2diffs > r2diff)
## not enough evidence to reject H0: simpler model is good enough

## the math way
f <- ((r2f - r2b)/pdiff)/((1 - r2f)/(n - pf))
f

## p-value (always right-tailed)
pf(f, pdiff, n - pf, lower.tail=FALSE)
## still fail to reject

## you won't get the same p-value because the stats are different
## but you can always modify your simulation to use f instead
Fs <- (R2diffs/pdiff)/((1 - R2fs)/(n - pf))
mean(Fs > f)
## now they match