# In-class demos for 2019-01-22: In-sample vs. generalization error, # cross-validation # "Synthetic" (i.e., fake) data example, so we know what the real answer is, # and can tell which methods do and do not work # 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) #10k points y.new = 7*x.new^2 - 0.5*x.new + rnorm(length(x.new)) # 10k data points is very "busy", 1k gets the point across 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) # Use predict() on new data to see what the models predict 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! # ATTN: Can you add a good legend? # 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 lm 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) # ATTN: Can you 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 # ATTN: Can you 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)