# In-class demos for 2015-01-20


### Illustration of over-fitting and model selection

# 20 standard-Gaussian X's
x = rnorm(20)
# Quadratic Y's
y = 7*x^2 - 0.5*x + rnorm(20)

# Initial plot of training data plus true regression curve
plot(x,y)
curve(7*x^2-0.5*x,col="grey",add=TRUE)

# Fit polynomials and add them to the plot
# Make a list of model formulae for polynomials of degree 0 to 9
poly.formulae <- c("y~1", paste("y ~ poly(x,", 1:9, ")", sep=""))
  # because "y ~ poly(x,0)" doesn't work
poly.formulae <- sapply(poly.formulae, as.formula)
# Get evenly spaced points for pretty plotting models
df.plot <- data.frame(x=seq(min(x),max(x),length.out=200))
# Fit polynomials of order 0 to 9
  # Store them in a list
fitted.models <- list(length=length(poly.formulae))
for (model_index in 1:length(poly.formulae)) {
  fm <- lm(formula=poly.formulae[[model_index]])
  lines(df.plot$x, predict(fm,newdata=df.plot),lty=model_index)
  fitted.models[[model_index]] <- fm
}
legend("top",legend=paste("degree",0:9),lty=1:10, ncol=2, cex=0.5)

# Calculate and plot in-sample errors
mse.q <- sapply(fitted.models, function(mdl) { mean(residuals(mdl)^2) })
plot(0:9,mse.q,type="b",xlab="polynomial degree",ylab="mean squared error",
     log="y")


# Plot the old curves with testing data
x.new = rnorm(2e4)
y.new = 7*x.new^2 - 0.5*x.new + rnorm(length(x.new))
plot(x.new[1:2e3],y.new[1:2e3],xlab="x",ylab="y",pch=24,cex=0.1,col="blue",
     xlim=2*range(x.new),ylim=2*range(y.new))
curve(7*x^2-0.5*x,col="grey",add=TRUE)
newdf <- data.frame(x=seq(from=2*min(x.new),to=2*max(x.new),length.out=200))
for (model_index in 1:length(fitted.models)) {
  fm <- fitted.models[[model_index]]
  lines(newdf$x, predict(fm, newdf), lty=model_index)
}
points(x,y)
legend("top",legend=paste("degree",0:9),lty=1:10, ncol=2, cex=0.5)

# Calculate and plot the out-of-sample errors

# Calculate an already-fitted model's MSE on new data
  # Inputs: the model object
  # Outputs: the MSE
  # ATTN: Hard-wired to use "x.new" and "y.new" --- can you improve this?
gmse <- function(mdl) {
  predictions <- predict(mdl, data.frame(x=x.new))
  residuals <- y.new - predictions
  return(mean(residuals^2))
}
gmse.q <- sapply(fitted.models, gmse)
plot(0:9,mse.q,type="b",xlab="polynomial degree",
     ylab="mean squared error",log="y",ylim=c(min(mse.q),max(gmse.q)))
lines(0:9,gmse.q,lty=2,col="blue")
points(0:9,gmse.q,pch=24,col="blue")
# Note: logarithmic scale of vertical axis!


# Demonstration: R^2 is an EPIC FAIL as a goodness-of-fit statistic
# So is adjusted R^2

# Extract R^2 and adjusted R^2 from a fitted model
  # Inputs: the model object
  # Output: vector of length 2, R^2 first
extract.rsqd <- function(mdl) {
  c( summary(mdl)$r.squared, summary(mdl)$adj.r.squared)
}
rsqd.q <- sapply(fitted.models, extract.rsqd)
plot(0:9,rsqd.q[1,],type="b",xlab="polynomial degree",ylab=expression(R^2),
  ylim=c(0,1))
lines(0:9,rsqd.q[2,],type="b",lty="dashed")
legend("bottomright",legend=c(expression(R^2),expression(R[adj]^2)),
  lty=c("solid","dashed"))



# General function to do k-fold CV for a bunch of linear models
  # Inputs: dataframe to fit all models on, list or vector of model formulae,
    # number of folds of cross-validation
  # Output: vector of cross-validated MSEs for the models
cv.lm <- function(data, formulae, nfolds=5) {
  # Strip data of NA rows
    # ATTN: Better to check whether NAs are in variables used by the models
  data <- na.omit(data)
  # Make sure the formulae have type "formula"
  formulae <- sapply(formulae, as.formula)
  # Extract the name of the response variable from each formula
    # ATTN: CV doesn't make a lot of sense unless these are all the same!
  responses <- sapply(formulae, response.name)
  names(responses) <- as.character(formulae)
  n <- nrow(data)
  # Assign each data point to a fold, at random
    # see ?sample for the effect of sample(x) on a vector x
  fold.labels <- sample(rep(1:nfolds, length.out=n))
  mses <- matrix(NA, nrow=nfolds, ncol=length(formulae))
  colnames <- as.character(formulae)
  # EXERCISE: Replace the double for() loop below by defining a new
  # function and then calling outer()
  for (fold in 1:nfolds) {
    test.rows <- which(fold.labels == fold)
    train <- data[-test.rows,]
    test <- data[test.rows,]
    for (form in 1:length(formulae)) {
       # Fit the model on the training data
       current.model <- lm(formula=formulae[[form]], data=train)
       # Generate predictions on the testing data
       predictions <- predict(current.model, newdata=test)
       # Get the responses on the testing data
       test.responses <- test[,responses[form]]
       test.errors <- test.responses - predictions
       mses[fold, form] <- mean(test.errors^2)
    }
  }
  return(colMeans(mses))
}

# Extract the name of the response variable from a regression formula
  # Presumes response is the left-most variable
  # EXERCISE: Write a more robust version using terms()
  # Inputs: regression formula
  # Outputs: name of the response variable
response.name <- function(formula) {
  var.names <- all.vars(formula)
  return(var.names[1])
}


# How well does cross-validation work?
  # Remember that our original data had 20 points, and we're seeing how
  # we generalize to 20,000 points from the same distribution
# Make a little data frame out of our little data
little.df <- data.frame(x=x, y=y)
# CV for the polynomials (defaults to five-fold)
cv.q <- cv.lm(little.df, poly.formulae)
plot(0:9,mse.q,type="b",xlab="polynomial degree",
     ylab="mean squared error",log="y",ylim=c(min(mse.q),max(gmse.q)))
lines(0:9,gmse.q,lty=2,col="blue",type="b",pch=2)
lines(0:9,cv.q,lty=3,col="red",type="b",pch=3)
legend("topleft",legend=c("In-sample","Generalization","CV"),
       col=c("black","blue","red"),lty=1:3,pch=1:3)