Assignment #4

Conceptual

Question 3

We now review k-fold cross-validation.

1. Explain how k-fold cross-validation is implemented.

   • k-fold cross-validation begins by diving the data into k groups or folds (k being set by the user) of approximately equal size. After dividing the data into equal groups, the model uses the first fold as the validation set. All other folds are used to train the model. The mean squared error, \(MSE_i\), is calculated on the observations in the held-out fold. This process will be repeated k times, and for each iteration a different fold will be used as the validation set. Lastly, the cross validation estimate is equal to the average of the k individual \(MSE\) values.

2. What are the advantages and disadvantages of k-fold cross-validation relative to:

   • The validation set approach?
     • Advantage of k-fold over validation set approach:
       • The k-fold method use more points for the training models compared to the validation set approace. Meaning the k-fold may provide a more accurate MSE
     • Disadvantage of k-fold over validation set approach:
       • The validation set approach runs faster then the k-fold as it does not go into as exstensive detail.
   • LOOCV?
     • Advantage of k-fold over LOOCV:
       • The k-fold method less variablility compared to LOOCV, as LOOCV is based on a single observation. Also, LOOCV may have computational problems if n is large, compared to k-fold which would run better.
     • Disadvantage of k-fold overLOOCV:
       • The LOOCV will have a lower bias then the k-fold since it uses less data to train the model.

Applied

Question 5

---

In Chapter 4, we used logistic regression to predict the probability of default using income and balance on the Default data set. We will now estimate the test error of this logistic regression model using the validation set approach. Do not forget to set a random seed before beginning your analysis.

# install.packages("ISLR2")
# install.packages("boot")
library(ISLR2)
library(boot)

default_data <- Default

1. Fit a logistic regression model that uses income and balance to predict default.

default_glm <- glm(default ~ income + balance,
                   data = default_data,
                   family = "binomial")
default_glm

## 
## Call:  glm(formula = default ~ income + balance, family = "binomial", 
##     data = default_data)
## 
## Coefficients:
## (Intercept)       income      balance  
##  -1.154e+01    2.081e-05    5.647e-03  
## 
## Degrees of Freedom: 9999 Total (i.e. Null);  9997 Residual
## Null Deviance:       2921 
## Residual Deviance: 1579  AIC: 1585

2. Using the validation set approach, estimate the test error of this model. In order to do this, you must perform the following steps:

   1. Split the sample set into a training set and a validation set.

# converting default column from "Yes" or "No" to 1 for Yes and 0 for no
default_data$default <- ifelse(default_data$default == "Yes", 1, 0)

set.seed(1)
default_train <- sample(1000, 500)
default_test <- default_data[-default_train,]

ii. Fit a multiple logistic regression model using only the training observations.

default_train_glm <- glm(default ~ income + balance,
                         data = default_data,
                         subset = default_train,
                         family = "binomial")
default_train_glm

## 
## Call:  glm(formula = default ~ income + balance, family = "binomial", 
##     data = default_data, subset = default_train)
## 
## Coefficients:
## (Intercept)       income      balance  
##  -1.177e+01    3.176e-05    5.711e-03  
## 
## Degrees of Freedom: 499 Total (i.e. Null);  497 Residual
## Null Deviance:       134.7 
## Residual Deviance: 78.33     AIC: 84.33

iii. Obtain a prediction of default status for each individual in the validation set by computing the posterior probability of default for that individual, and classifying the individual to the default category if the posterior probability is greater than 0.5.

default_preds <- predict(default_train_glm, default_data, type = "response")
default.preds <- ifelse(default_preds > 0.5, 1, 0)

iv. Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.

mean((default.preds != default_data$default)[-default_train])

## [1] 0.02652632

3. Repeat the process in (b) three times, using three different splits of the observations into a training set and a validation set. Comment on the results obtained.

# training set will only have a 1/4 of the data points

set.seed(1)
default_train_250 <- sample(1000, 250)
default_test_250 <- default_data[-default_train_250,]

default_train_250_glm <- glm(default ~ income + balance,
                         data = default_data,
                         subset = default_train_250,
                         family = "binomial")
default_train_250_glm

## 
## Call:  glm(formula = default ~ income + balance, family = "binomial", 
##     data = default_data, subset = default_train_250)
## 
## Coefficients:
## (Intercept)       income      balance  
##  -9.242e+00    9.932e-06    4.596e-03  
## 
## Degrees of Freedom: 249 Total (i.e. Null);  247 Residual
## Null Deviance:       70.81 
## Residual Deviance: 48.89     AIC: 54.89

default_250_preds <- predict(default_train_250_glm, default_data, type = "response")
default_250.preds <- ifelse(default_250_preds > 0.5, 1, 0)

mean((default_250.preds != default_data$default)[-default_train_250])

## [1] 0.02666667

# training set will contain 3/4 of the data points

set.seed(1)
default_train_750 <- sample(1000, 750)
default_test_750 <- default_data[-default_train_750,]

default_train_750_glm <- glm(default ~ income + balance,
                         data = default_data,
                         subset = default_train_750,
                         family = "binomial")
default_train_750_glm

## 
## Call:  glm(formula = default ~ income + balance, family = "binomial", 
##     data = default_data, subset = default_train_750)
## 
## Coefficients:
## (Intercept)       income      balance  
##  -1.224e+01    4.369e-05    5.816e-03  
## 
## Degrees of Freedom: 749 Total (i.e. Null);  747 Residual
## Null Deviance:       212.4 
## Residual Deviance: 119.5     AIC: 125.5

default_750_preds <- predict(default_train_750_glm, default_data, type = "response")
default_750.preds <- ifelse(default_750_preds > 0.5, 1, 0)

mean((default_750.preds != default_data$default)[-default_train_750])

## [1] 0.0267027

# training set will contain 632 points and set seed will be at 5

set.seed(5)
default_train_632 <- sample(1000, 632)
default_test_632 <- default_data[-default_train_632,]

default_train_632_glm <- glm(default ~ income + balance,
                         data = default_data,
                         subset = default_train_632,
                         family = "binomial")
default_train_632_glm

## 
## Call:  glm(formula = default ~ income + balance, family = "binomial", 
##     data = default_data, subset = default_train_632)
## 
## Coefficients:
## (Intercept)       income      balance  
##  -1.065e+01    3.445e-05    4.966e-03  
## 
## Degrees of Freedom: 631 Total (i.e. Null);  629 Residual
## Null Deviance:       177.5 
## Residual Deviance: 112   AIC: 118

default_632_preds <- predict(default_train_632_glm, default_data, type = "response")
default_632.preds <- ifelse(default_632_preds > 0.5, 1, 0)

mean((default_632.preds != default_data$default)[-default_train_632])

## [1] 0.02690009

The main thing I noticed when changing the training set size, is that the test error did not change much. The initial test (half the data points in the training set) the error was 0.02652632 but as the size change the only by the ten-thousands place (example: 0.02690009).

4. Now consider a logistic regression model that predicts the probability of default using income, balance, and a dummy variable for student. Estimate the test error for this model using the validation set approach. Comment on whether or not including a dummy variable for student leads to a reduction in the test error rate.

# creating dummy variable for student
default_data$student <- ifelse(default_data$student == "Yes", 1, 0)

updated_default_glm <- glm(default ~ income + balance + student,
                   data = default_data,
                   family = "binomial")
updated_default_glm

## 
## Call:  glm(formula = default ~ income + balance + student, family = "binomial", 
##     data = default_data)
## 
## Coefficients:
## (Intercept)       income      balance      student  
##  -1.087e+01    3.033e-06    5.737e-03   -6.468e-01  
## 
## Degrees of Freedom: 9999 Total (i.e. Null);  9996 Residual
## Null Deviance:       2921 
## Residual Deviance: 1572  AIC: 1580

set.seed(1)
updated_default_train <- sample(1000, 500)
updated_default_test <- default_data[-updated_default_train,]

updated_default_train_glm <- glm(default ~ income + balance + student,
                         data = default_data,
                         subset = updated_default_train,
                         family = "binomial")
updated_default_train_glm

## 
## Call:  glm(formula = default ~ income + balance + student, family = "binomial", 
##     data = default_data, subset = updated_default_train)
## 
## Coefficients:
## (Intercept)       income      balance      student  
##  -1.072e+01    1.949e-06    5.931e-03   -1.151e+00  
## 
## Degrees of Freedom: 499 Total (i.e. Null);  496 Residual
## Null Deviance:       134.7 
## Residual Deviance: 77.14     AIC: 85.14

updated_default_preds <- predict(updated_default_train_glm, default_data, type = "response")
updated_default.preds <- ifelse(updated_default_preds > 0.5, 1, 0)

mean((updated_default.preds != default_data$default)[-updated_default_train])

## [1] 0.02673684

After adding creating a dummy variable for student and adding it into logistical regression, it did not have a significant impact to the test error. After using the updated model and training the set, the error lied in the same range of the initial tests in the previous questions.

Question 6

---

We continue to consider the use of a logistic regression model to predict the probability of default using income and balance on the Default data set. In particular, we will now compute estimates for the standard errors of the income and balance logistic regression coefficients in two different ways: (1) using the bootstrap, and (2) using the standard formula for computing the standard errors in the glm() function. Do not forget to set a random seed before beginning your analysis.

1. Using the summary() and glm() functions, determine the estimated standard errors for the coefficients associated with income and balance in a multiple logistic regression model that uses both predictors.

default.glm <- glm(default ~ income + balance,
                   data = default_data,
                   family = "binomial")
summary(default.glm)

## 
## Call:
## glm(formula = default ~ income + balance, family = "binomial", 
##     data = default_data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
## income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
## balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1579.0  on 9997  degrees of freedom
## AIC: 1585
## 
## Number of Fisher Scoring iterations: 8

2. Write a function, boot.fn(), that takes as input the Default data set as well as an index of the observations, and that outputs the coefficient estimates for income and balance in the multiple logistic regression model.

boot.fn <- function(data, index) {
  coef(glm(default ~ income + balance,
                  data = data,
                  family = "binomial",
                  subset = index))
}

boot.fn(default_data, 1:1000)

##   (Intercept)        income       balance 
## -1.155359e+01  3.074546e-05  5.609161e-03

3. Use the boot() function together with your boot.fn() function to estimate the standard errors of the logistic regression coefficients for income and balance.

set.seed(1)
boot(default_data, boot.fn, R = 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = default_data, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##          original        bias     std. error
## t1* -1.154047e+01 -3.945460e-02 4.344722e-01
## t2*  2.080898e-05  1.680317e-07 4.866284e-06
## t3*  5.647103e-03  1.855765e-05 2.298949e-04

4. Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.

When comparing the estimated standard errors between glm() and bootstrap function, it can be noticed that the standard errors are very close to one another. The estimated standard error for the income intercept are within .120 of each other and the estimated standard error for the balance intercept are within .02 of each other.

Question 9

---

We will now consider the Boston housing data set, from the ISLR2 library.

1. Based on this data set, provide an estimate for the population mean of medv. Call this estimate \(\hat{µ}\)

# u hat
mean_of_medv <- mean(Boston$medv)
mean_of_medv

## [1] 22.53281

2. Provide an estimate of the standard error of \(\hat{µ}\). Interpret this result.
   • Hint: We can compute the standard error of the sample mean by dividing the sample standard deviation by the square root of the number of observations.

# 
num <- sd(Boston$medv)
denom <- sqrt(nrow(Boston))

se_u_hat <- num / denom
se_u_hat

## [1] 0.4088611

3. Now estimate the standard error of \(\hat{µ}\) using the bootstrap. How does this compare to your answer from (b)?

boston.boot.fn <- function(data, index) {
  mean(data$medv[index])
}

boston.boot.fn(Boston, 1:506)

## [1] 22.53281

set.seed(1)
boot(Boston, boston.boot.fn, R = 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston, statistic = boston.boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 22.53281 0.007650791   0.4106622

Using bootstap, the standard error was only .01 greater than the standard error calculated in part b.

4. Based on your bootstrap estimate from (c), provide a 95 % confidence interval for the mean of medv. Compare it to the results obtained using t.test(Boston$medv).
   • Hint: You can approximate a 95 % confidence interval using the formula \([\hat{µ} - 2SE(\hat{µ}), \hat{µ} + 2SE(\hat{µ})]\)

set.seed(1)
boot_medv <- boot(Boston, boston.boot.fn, R = 1000)

confidence_lower <- mean_of_medv - (2*(sd(boot_medv$t)))
confidence_upper <- mean_of_medv + (2*(sd(boot_medv$t)))

c(confidence_lower, confidence_upper)

## [1] 21.71148 23.35413

t.test(Boston$medv)

## 
##  One Sample t-test
## 
## data:  Boston$medv
## t = 55.111, df = 505, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  21.72953 23.33608
## sample estimates:
## mean of x 
##  22.53281

The confidence interval provided by the T test and the confidence interval calculated are nearly identical. The lower limit of the confidence interval have a difference of .05 and the upper limit of the confidence interval have a difference of .04.

5. Based on this data set, provide an estimate,\(\hat{µ}_(med)\), for the median value of medv in the population.

median_medv <- median(Boston$medv)
median_medv

## [1] 21.2

6. We now would like to estimate the standard error of \(\hat{µ}_(med)\). Unfortunately, there is no simple formula for computing the standard error of the median. Instead, estimate the standard error of the median using the bootstrap. Comment on your findings.

boot_median.fn <- function(data, index) {
  median(data$medv[index])
}

set.seed(1)
boot(Boston, boot_median.fn, R = 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston, statistic = boot_median.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     21.2 0.02295   0.3778075

When using the median, bootstrap calculated a lower standard error than when using the mean of medv. The standard error using median medv is lower by .04.

7. Based on this data set, provide an estimate for the tenth percentile of medv in Boston census tracts. Call this quantity \(\hat{µ}_(0.1)\). (You can use the quantile() function.)

u_hat_0.1 <- quantile(Boston$medv, 0.1)
u_hat_0.1

##   10% 
## 12.75

8. Use the bootstrap to estimate the standard error of \(\hat{µ}_(0.1)\). Comment on your findings.

boot_10_percent.fn <- function(data, index) {
  quantile(data$medv[index], 0.1)
}

set.seed(1)
boot(Boston, boot_10_percent.fn, R = 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston, statistic = boot_10_percent.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    12.75  0.0339   0.4767526

After using bootstrap to find the standard error of \(\hat{µ}\), the tenth percentile provides the largest standard error. This standard error is .47 and is higher than the inital test on the mean by .06.