# Hitters data - PCA regression # Geir's original code ------------------------------------------------------------------------ # Load the ISLR library library(ISLR) # See the variable names in the Hitters data set names(Hitters) # Get the dimensions of the Hitters data set. 263 observations with 20 columns (19 explanatory variables) dim(Hitters) # Get the number of observations with missing entries (none are missing) sum(is.na(Hitters$Salary)) # Code that removes missing entries, but nothing was missing so this has no effect. # can verify it by 'all.equal(Hitters, na.omit(Hitters))'. Hitters = na.omit(Hitters) # Get the dimensions again dim(Hitters) # See that no entries are missing (obviously) sum(is.na(Hitters)) # Extract the columns that are numeric. # ind = c(1:13,16:18) # I don't like to enter this manually. So I rather use. # Get the indices of the columns that are numeric and is not the response ind = seq(ncol(Hitters))[sapply(Hitters, is.numeric) & !names(Hitters) %in% "Salary"] ind # Get the 16 principal components for the 16 continuous variables Z = prcomp(Hitters[,ind], retx=TRUE, scale=TRUE)$x # Create a new data frame using the Principal components for the continuous variables # and including the original categorical features. Hitters2 = data.frame(Z = Z, League = Hitters$League, Division = Hitters$Division, NewLeague = Hitters$NewLeague, Salary = Hitters$Salary) # We fit a PCA model. Note that we use LM as we have allready computed the PCs. fit.pca = lm(Salary ~ ., data = Hitters2) ## Dividing into training/test set # Set seed for reproducibility set.seed(1) # Sample the indices for the training observations train = sample(1:nrow(Hitters), nrow(Hitters)/2) test = (-train) # Extract the test responses y.test = Hitters$Salary[test] # Create an empty vector which we will populate with the RSS for the different number of PCs. RSS = rep(NA, length(ind)) # Iterate over the number of PCs. for(m in 1:length(ind)) { # Create new dataset with just the first m PCs and the categorical data (and response) Hitters2 = data.frame(Z = Z[,1:m], League = Hitters$League, Division = Hitters$Division, NewLeague = Hitters$NewLeague, Salary = Hitters$Salary) # Fit the PCA model on the training data fit.pca = lm(Salary ~ ., data = Hitters2[train,]) # Predict the response of the teste data pred.pca = predict(fit.pca, Hitters2[test,]) # Store the residual sum of squares RSS[m] = sqrt(mean((y.test - pred.pca)^2)) } # Plot the RSS for the different number of PCs. # It seems like 5 or 13 is the best (for this training/test split) plot(1:m, RSS) # Additional code by Lars --------------------------------------------------------------------- ## Cross validation without packages ---------------------------------------------------------- # Set the number of CV folds num_cv_folds = nrow(Hitters) num_cv_folds = 10 # Get the indices for the different folds # Many ways to do it. This is maybe not the best cv_indices = floor(seq(1, nrow(Hitters), length = num_cv_folds + 1)) # Matrix to store the RSS results cv_model_mse_matrix = matrix(NA, nrow = num_cv_folds, ncol = length(ind)) # Iterate over the number of PCs. for(m in 1:length(ind)) { # Create new dataset with just the first m PCs and the categorical data (and response) Hitters2 = data.frame(Z = Z[,1:m], League = Hitters$League, Division = Hitters$Division, NewLeague = Hitters$NewLeague, Salary = Hitters$Salary) for (cv_fold in seq(num_cv_folds)) { # Get the indices of the test observations in this fold cv_fold_indices = seq(cv_indices[cv_fold], cv_indices[cv_fold + 1]) # Get training and evaluation data cv_test_data_train = Hitters2[cv_fold_indices,] cv_train_data_train = Hitters2[-cv_fold_indices,] # Fit the PCA model on the training data in the CV fit.pca = lm(Salary ~ ., data = cv_train_data_train) # Predict the response of the test data in the CV pred.pca = predict(fit.pca, cv_test_data_train) # Store the residual sum of squares cv_model_mse_matrix[cv_fold, m] = sqrt(mean((cv_test_data_train$Salary - pred.pca)^2)) } } # Aggregate the results RSS_CV_Version_1 = colMeans(cv_model_mse_matrix) # Plot the results { matplot(1:m, cbind(RSS, RSS_CV_Version_1), type = "l", lty = 1, xlab = "Number of Principal Components", ylab = "RSS") legend("topleft", legend = c("Train/test", bquote(.(num_cv_folds)*"-fold CV (V1)")), lty = 1, col = 1:3) abline(v = which.min(RSS), col = 1, lty = 2, lwd = 2) abline(v = which.min(RSS_CV_Version_1), col = 2, lty = 2, lwd = 2) } ## Cross validation using tidyverse ---------------------------------------------------------------------------- library(tidyverse) library(modelr) # Number of folds num_folds = 50 # Create an empty vector which we will populate with the RSS for the different number of PCs. RSS_CV_Version_2 = rep(NA, length(ind)) # Iterate over the number of PCs. for(m in 1:length(ind)) { # Set seed for reproducibility in sampling of the folds, such that they are the same for each m set.seed(1) # Create new dataset with just the first m PCs and the categorical data (and response) # Fit the PCA models on the training data trained_models = data.frame(Z = Z[,1:m], League = Hitters$League, Division = Hitters$Division, NewLeague = Hitters$NewLeague, Salary = Hitters$Salary) %>% crossv_kfold(num_folds) %>% mutate(model = purrr::map(train, ~lm(Salary ~ ., data=.))) # Store the square root of the residual sum of squares RSS_CV_Version_2[m] = mean(map2_dbl(trained_models$model, trained_models$test, rmse)) } # Plot the results { matplot(1:m, cbind(RSS, RSS_CV_Version_1, RSS_CV_Version_2), type = "l", lty = 1, xlab = "Number of Principal Components", ylab = "RSS") legend("topleft", legend = c("Train/test", bquote(.(num_cv_folds)*"-fold CV (V1)"), bquote(.(num_folds)*"-fold CV (V2)")), lty = 1, col = 1:3) abline(v = which.min(RSS), col = 1, lty = 2, lwd = 2) abline(v = which.min(RSS_CV_Version_1), col = 2, lty = 2, lwd = 2) abline(v = which.min(RSS_CV_Version_2), col = 3, lty = 3, lwd = 2) }