### ---------------------------------------------------------------------------- ### STK2100, spring 2026 ### Solution to modified Extra exercise 10 ### By Vera Kvisgaard 06.05.25 ### Last modified by Berit Omli Øksnes 10.04.26 ### ---------------------------------------------------------------------------- # Clear workspace rm(list=ls()) ### a) Read data and remove missing values ------------------------------------- path <- "https://www.uio.no/studier/emner/matnat/math/STK2100/data/heart.csv" df <- read.csv(path, stringsAsFactors=TRUE) sapply(df, class) head(df) # Exclude rows with missing values df <- df[!(rowSums(is.na(df)) > 0), ] N <- nrow(df) p <- ncol(df) - 1 # Draw train-test split set.seed(1234) train <- sample(1:N, round(N*0.7)) N_train <- length(train) N_test <- N - length(train) df_train <- df[train, ] df_test <- df[-train, ] # Convert to matrices for methods that prefer matrices to data frames... x.tr <- vapply(df_train[, -ncol(df_train)], as.numeric, numeric(N_train)) x.test <- vapply(df_test[, -ncol(df_train)], as.numeric, numeric(N_test)) y.tr <- df_train$AHD y.test <- df_test$AHD ### b) reproduce Figure 8.6 in the ISLR book ----------------------------------- library(tree) # set.seed(007) par(mfrow=c(1,3)) # Define hyperparameters and fit full tree control <- tree.control(nobs=N_train, mincut=1, mindev=0.000001) fit <- tree(AHD~., data=df_train, control=control) ## Reproduce top figure: unpruned tree plot(fit) text(fit, digits=2, pretty=3, cex=0.5) # Define function that returns a sequence of pruned trees prune_tree_sequence <- function(tree) { trees <- list(tree) ntrees <- 1 for (size in seq(summary(tree)$size, 2)) { # prune tree tree <- prune.tree(tree, best=size, # number of terminal nodes method="misclass") # Check if tree is changed # - prune.tree() returns the larger tree if no subtree of size `size` is found if (identical(trees[[ntrees]], tree)) next ntrees <- ntrees + 1 trees[[ntrees]] <- tree } return(trees) } # Define function that computes misclassification error of tree compute_errorrate <- function(fit, df) { # Compute predicted class probabilities probs <- predict(fit, df) # Define predictions and true values y_hat <- (probs[, 2] > 1/2) y <- (df$AHD == "Yes") # Return misclassification error return(mean(y_hat != y)) } maxsize <- 20 err <- matrix(NA, nrow=maxsize, ncol=3) colnames(err) <- c("Train", "Test", "CV") ## Test and train error for tree-sequence fits <- prune_tree_sequence(fit) for (fit in fits) { size <- summary(fit)$size err[size, "Train"] <- compute_errorrate(fit, df_train) err[size, "Test"] <- compute_errorrate(fit, df_test) } ## CV estimate of test error for tree-sequence: # Partition the data into K folds # For each fold, iterate over the tree-sizes, prune tree and compute test error # Draw partition of training data into K folds K <- 10 fold <- sample(1:K, N_train, replace=TRUE) err_cv <- matrix(nrow=maxsize, ncol=K) for (k in 1:K) { # Define training data for current fold train_k <- which(fold != k) df_train_k <- df_train[train_k,] df_test_k <- df_test[-train_k,] # Initiate a deep tree on current training data control <- tree.control(nobs=length(train_k), mincut=1, mindev=0.000001) fit_k <- tree(AHD~., data=df_train_k, control=control) # compute the sequence of pruned trees fits_k <- prune_tree_sequence(fit_k) # estimate test error on hold-out sample for each tree for (fit in fits_k) { size <- summary(fit)$size err_cv[size, k] <- compute_errorrate(fit, df_test_k) } } # Average test-error over all folds err[, "CV"] <- rowMeans(err_cv, na.rm = TRUE)[1:nrow(err)] ## Reproduce bottom left figure: # CV error, training, and test error, for different sizes of the pruned tree. sizes <- 1:nrow(err) j <- 1 indx <- !is.na(err[, j]) plot(sizes[indx], err[indx, j], col=j, pch=j, ylim=c(0, 0.6), xlab="Tree Size", ylab="Error") lines(sizes[indx], err[indx, j], col=j, pch=j) for (j in 2:ncol(err)) { indx <- !is.na(err[, j]) points(sizes[indx], err[indx, j], col=j, pch=j) lines(sizes[indx], err[indx, j], col=j, pch=j) } legend("topright", legend=colnames(err), col=1:ncol(err), lty=1, pch=1:3) ## Reproduce bottom right figure: # The pruned tree corresponding to the minimal cross-validation error. best_CV_size <- which.min(err[,"CV"]) fit_CV <- fits[[which(sapply(fits, function(x) summary(x)$size)==best_CV_size)]] plot(fit_CV) text(fit_CV, digits=2, pretty=3, cex=0.5) par(mfrow=c(1,1)) ### c)-d) Bagging -------------------------------------------------------------- library(randomForest) B <- 300 fit <- randomForest(x=x.tr, y=df_train$AHD, xtest=x.test, ytest=df_test$AHD, ntree=B, mtry=p) err <- list(Test=fit$test$err.rate[, 1], # error on test data OOB=fit$err.rate[, 1]) # error computed on OOB samples plot(err[[1]], col="black", type="l", ylim=c(0.1, 0.3), xlab="Number of trees", ylab="Error") lines(err[[2]], col="green", type="l", ylim=c(0.1, 0.3)) legend("topright", legend=names(err), col=c("black", "green"), lty = 1) ### e) Boosting (xgboost) ------------------------------------------------------ ## Fit boosting trees to the data, first with stumps (J=2) as the weak learner library(xgboost) # Use 10-fold cross-validation to select M and the learning rate is set to 0.1 fit.boost.1.cv <- xgb.cv(data=xgb.DMatrix(data=x.tr,label=(y.tr=="Yes")), nrounds=1000, eta=0.1, nfold=10, metrics=list("error"), max_depth=1, objective="binary:logistic") best.boost.1 <- which.min(fit.boost.1.cv$evaluation_log$test_error_mean) plot(fit.boost.1.cv$evaluation_log$train_error_mean, type="l", lwd=2, xlab="Stage", ylab="Misclassification error") lines(fit.boost.1.cv$evaluation_log$test_error_mean, col=2, lwd=2) abline(v=best.boost.1, lty=2, lwd=2) legend("topright", c("training","test"), col=c(1,2), lty=rep(1,2))