# Homework Assignment 2: The Advantages of Backwardness # 36-402, Data Analysis, Spring 2011 # Solutions by Z. Kurtz # Setup install.packages("np") require(np) data(oecdpanel) # Take a quick look at the variables names(oecdpanel) plot(oecdpanel) ############################################################### # PROBLEM 1: # Fit a linear model of growth on initgdp lm.1 = lm(growth ~ initgdp, data = oecdpanel) lm.1\$coefficients ############################################################### # PROBLEM 2: # Use provided code to explore dependence of MSE on bandwidth # for a kernel regression of growth on initgdp # run the supplied code source("hw-02.R") # modify the inner loop of supplied code for the "whole data" case (non-CV) whole.MSE = rep(NA, length(bandwidths)) for (bw in bandwidths) { # Fit to the training set = whole data current.npr = npreg(growth ~ initgdp, data=oecdpanel,bws=bw) # Predict on the test set = whole data, so predictions are just fitted values predictions = fitted(current.npr) # What's the mean-squared error? whole.MSE[which(bandwidths == bw)] = mean((oecdpanel\$growth - predictions)^2) # Note the alternative to the paste() trick for accessing bandwidths # Also: one part of the objets returned by npreg is an MSE attribute, so # we could have done # whole.MSE(which(bandwidths == bw)] = current.npr\$MSE # but the more explicit approach will also work if we decide to start using # a different smoothing package } # display the results (bandwidths.cv.mses and whole.MSE) pdf("graphics/p2.pdf") plot(bandwidths, bandwidths.cv.mses, type = "l", lwd = 2, ylim = c(min(c(bandwidths.cv.mses, whole.MSE)), max(c(bandwidths.cv.mses, whole.MSE))), main = "MSE vs. bandwidth \n5-fold CV (black, solid) and in-sample (green, dashed)", ylab = "MSE", xlab = "Kernel bandwidth") lines(bandwidths, whole.MSE, lwd = 2, col = "green", lty=2) dev.off() ############################################################### # PROBLEM 3: pdf("graphics/p3.pdf") # Make scatterplot of initgdp versus growth: plot(oecdpanel\$initgdp, oecdpanel\$growth, lwd = 2, cex = 0.3, xlab = "initgdp", ylab = "growth", main = "growth vs. initgdp") # Add the line for the linear model: abline(lm.1, lwd = 2, col = "red") # Fit the kernel regression with the [approximately] optimal bandwidth 0.3: opt.reg = npreg(growth ~ initgdp, data=oecdpanel,bws=0.3) fitted.values = fitted(opt.reg) #equivelently: fitted.values = predict(opt.reg, newdata = oecdpanel) points(oecdpanel\$initgdp, fitted.values, col = "green", cex = 0.5, lwd = 1.8) dev.off() ############################################################### # PROBLEM 4: # Fit a linear model of growth on initgdp, popgro, and inv lm.4 = lm(growth ~ initgdp+popgro+inv, data = oecdpanel) signif(lm.4\$coefficients,3) ############################################################### # PROBLEM 5: oecd.npr <- npreg(growth ~ initgdp + popgro + inv, data=oecdpanel, tol=0.1, ftol=0.1) summary(oecd.npr) ############################################################### # PROBLEM 6: # Find and report the medians: med.popgro = median(oecdpanel\$popgro) signif(med.popgro,3) med.inv = median(oecdpanel\$inv) signif(med.inv,3) # Make a new data set med.oecdpanel that contains growth, initgdp, popgro, # and inv, but replace popgro and inv and with their respective medians # Use the actual values of initgdp, but sort them for easier plotting later initgdp.sort = sort(oecdpanel\$initgdp) med.oecdpanel = data.frame(initgdp=initgdp.sort, popgro=med.popgro,inv=med.inv) # Find the predicted growth rates under the model from problem 4 mod4.med.preds = predict(lm.4, newdata = med.oecdpanel) # Find the predicted growth rates under the model from problem 5 mod5.med.preds = predict(oecd.npr, newdata = med.oecdpanel) # Plot growth v.s. initgdp for med.data and add the fitted values for both models: pdf("graphics/p6.pdf") # Make scatterplot of initgdp versus growth: plot(oecdpanel\$initgdp, oecdpanel\$growth, lwd = 2, cex = 0.3, xlab = "initgdp", ylab = "growth", main = "Data with popgro, inv set to their medians: growth vs. initgdp") lines(initgdp.sort, mod4.med.preds, col = "red", cex = 0.5, lwd = 2.8) lines(initgdp.sort, mod5.med.preds, col = "gray", cex = 0.5, lwd = 2.8) mtext("Red straight = linear model fit; Gray curve = kernel regression fit") dev.off() ############################################################### # PROBLEM 7: # Code from problem 2 with some modifications for doing CV with both the linear # model and the kernel regression (with automatic bandwidth selection) nfolds = 5 case.folds = rep(1:nfolds,length.out=nrow(oecdpanel)) case.folds = sample(case.folds) linear.fold.mses = vector(length=nfolds) kernel.fold.mses = vector(length=nfolds) for (fold in 1:nfolds) { train = oecdpanel[case.folds!=fold,] test = oecdpanel[case.folds==fold,] # Fit the models linear.model = lm(growth ~ initgdp+popgro+inv, data=train) kernel.model = npreg(growth ~ initgdp + popgro + inv, data=train, tol=0.1, ftol=0.1) # Predict on the test set linear.predictions = predict(linear.model, newdata=test) kernel.predictions = predict(kernel.model, newdata=test) linear.fold.mses[fold] = mean((test\$growth - linear.predictions)^2) kernel.fold.mses[fold] = mean((test\$growth - kernel.predictions)^2) } cv.linear.mod = mean(linear.fold.mses) cv.kernel.mod = mean(kernel.fold.mses)