Author: Jay Liao (ID: RE6094028)

Load in the packages

library(dplyr)
library(GGally)
library(corrplot)

Exercise 5.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.

library(ISLR)
data('Default')
summary(Default)

 default    student       balance           income     
 No :9667   No :7056   Min.   :   0.0   Min.   :  772  
 Yes: 333   Yes:2944   1st Qu.: 481.7   1st Qu.:21340  
                       Median : 823.6   Median :34553  
                       Mean   : 835.4   Mean   :33517  
                       3rd Qu.:1166.3   3rd Qu.:43808  
                       Max.   :2654.3   Max.   :73554

Exercise 5.5 - (a)

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

Correlation and scatter plots

ggpairs(Default[,3:4], mapping = aes(color = Default$default))

The correlation between income and balance is slight, which may not introduce the problem of collinearity.

Fit the logistic regression

fit_5.5 <- glm(default ~ income + balance, family = (binomial('logit')), data = Default)
summary(fit_5.5)


Call:
glm(formula = default ~ income + balance, family = (binomial("logit")), 
    data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4725  -0.1444  -0.0574  -0.0211   3.7245  

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

According the result of the logistic regression, both income and balance have significant regression coefficients.

Exercise 5.5 - (b)

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

[Step 1] Split the sample set into a training set and a validation set.
[Step 2] Fit a multiple logistic regression model using only the training observations.
[Step 3] 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.

[Step 1] Splitting ratio = 3:1.

set.seed(6094028)
idx_train <- sample(c(T, F), nrow(Default), replace = T, prob = c(.75, .25))

[Step 2]

Default$y <- as.numeric(Default$default == 'Yes')
fit_5.5_b <- glm(y ~ income + balance, family = (binomial('logit')),
                 data = Default[idx_train,])
summary(fit_5.5_b)


Call:
glm(formula = y ~ income + balance, family = (binomial("logit")), 
    data = Default[idx_train, ])

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5010  -0.1423  -0.0552  -0.0200   3.7542  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.176e+01  5.192e-01 -22.646  < 2e-16 ***
income       2.200e-05  5.796e-06   3.796 0.000147 ***
balance      5.744e-03  2.695e-04  21.314  < 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: 2195.3  on 7447  degrees of freedom
Residual deviance: 1176.4  on 7445  degrees of freedom
AIC: 1182.4

Number of Fisher Scoring iterations: 8

According the result of the logistic regression, both income and balance have significant regression coefficients.

[Step 3]

p_hat_train <- fit_5.5_b %>% 
  predict.glm(newdata = Default[idx_train,], 'response')
p_hat_val <- fit_5.5_b %>% 
  predict.glm(newdata = Default[!idx_train,], 'response')
y_hat_train <- as.numeric(p_hat_train > 0.5)
y_hat_val <- as.numeric(p_hat_val > 0.5)

Accuracy score

mean(y_hat_train == Default$y[idx_train]) # Training

[1] 0.9742213

mean(y_hat_val == Default$y[!idx_train]) # Validation

[1] 0.9721787

Exercise 5.5 - (c)

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.

Write a funtion to generalize

f_5.5_b <- function(seeds, tr_proportion) {
  # Training
  set.seed(seeds)
  idx_train <- sample(c(T, F), nrow(Default), replace = T,
                      prob = c(tr_proportion, 1 - tr_proportion))
  fit_5.5_b <- glm(y ~ income + balance, family = (binomial('logit')),
                 data = Default[idx_train,])
  
  # Predicting
  p_hat_train <- fit_5.5_b %>% 
    predict.glm(newdata = Default[idx_train,], 'response')
  p_hat_val <- fit_5.5_b %>% 
    predict.glm(newdata = Default[!idx_train,], 'response')
  y_hat_train <- as.numeric(p_hat_train > 0.5)
  y_hat_val <- as.numeric(p_hat_val > 0.5)
  
  return(list(y_hat_train, y_hat_val))
}

Apply the function to different splitting ratios (different training proportions)

To get more detailed results and pattern, I repeat the process in (b) eight times, instead of only three times.

tr_proportion_lst <- c(.60, .65, .70, .75, .80, .85, .90, .95)
mtx_retult <- sapply(tr_proportion_lst, function(p) {
  reuslt_lst <- f_5.5_b(seeds = 6094028, tr_proportion = p)
  y_hat_train <- reuslt_lst[[1]]
  y_hat_val <- reuslt_lst[[2]]
  return(c(mean(y_hat_train == Default$y[idx_train]),
           mean(y_hat_val == Default$y[!idx_train])))
})

df <- as.data.frame(t(mtx_retult))
colnames(df) <- c('Training accuracy', 'Validation accuracy')
rownames(df) <- tr_proportion_lst
df

Visualization

df_long <- df %>% t() %>% as.data.frame() %>%
  mutate(dataset = colnames(df)) %>% reshape2::melt()
df_long %>% mutate(variable = as.numeric(as.character(variable))) %>%
  qplot(x = variable, y = value, col = dataset, data = .) +
  geom_line() +
  labs(x = 'Training proportion', y = 'Accuracy', colour = '') +
  theme_bw() + theme(legend.position = 'top')

由此圖可知：

當訓練資料集比例較低時（i.e., $p_{tr} < .75$），訓練與驗證資料集得到預測結果的accuracy都相對低，training accuracy甚至還低於validation accuracy，顯示用太少資料來訓練模型，導致underfitting。

-而訓練資料集比例較高時，訓練與驗證資料集得到預測結果的accuracy也都相對低，training accuracy高於validation accuracy，顯示用太多資料來訓練模型，導致模型判別新資料的能力不佳，也就所謂的overfitting。

似乎存在一個最合適的訓練資料集比例，使得訓練與驗證資料集得到預測結果的accuracy都相對高，沒有發生嚴重的underfitting或overfitting。

Exercise 5.5 - (d)

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.

Write a funtion to generalize

f_5.5_d <- function(seeds, variables, tr_proportion, return_error = TRUE) {
  
  # Train
  set.seed(seeds)
  idx_train <- sample(c(T, F), nrow(Default), replace = T,
                      prob = c(tr_proportion, 1 - tr_proportion))
  if ('student' %in% variables) {
    Default$student <- as.numeric(Default$student == 'Yes')
  }
  glm_fit <- glm(y ~ ., family = binomial('logit'),
                 data = Default[idx_train, c('y', variables)])
  
  # Predict
  p_hat_train <- glm_fit %>% predict.glm(
    newdata = Default[idx_train, variables], 'response')
  p_hat_val <- glm_fit %>% predict.glm(
    newdata = Default[!idx_train, variables], 'response')
  y_hat_train <- as.numeric(p_hat_train > 0.5)
  y_hat_val <- as.numeric(p_hat_val > 0.5)
  
  # Return
  if (return_error) {
    error_train <- mean(y_hat_train != Default$y[idx_train])
    error_val <- mean(y_hat_val != Default$y[!idx_train])
    return(c(error_train, error_val))
  }
  else {
    return(list(y_hat_train, y_hat_val))
  }
}

Apply the funtion to different variable sets

error_no_student <- f_5.5_d(seeds = 6094028, tr_proportion = .75,
                            variables = c('income', 'balance'))
error_student <- f_5.5_d(seeds = 6094028, tr_proportion = .75,
                         variables = c('income', 'balance', 'student'))
df <- data.frame(error_no_student, error_student)
rownames(df) <- c('Training error', 'Testing error')
df

It seems that including a dummy variable for student does not lead to a reduction in the test error rate.

Exercise 5.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.

Exercise 5.6 - (a)

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.

glm_fit_5.6 <- Default %>% glm(default ~ balance + income,
                               family=binomial('logit'), data = .)
summary(glm_fit_5.6)


Call:
glm(formula = default ~ balance + income, family = binomial("logit"), 
    data = .)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4725  -0.1444  -0.0574  -0.0211   3.7245  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
---
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

coef(glm_fit_5.6) %>% round(.,6)

(Intercept)     balance      income 
 -11.540468    0.005647    0.000021

Exercise 5.6 - (b)

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(dta, index) {
  model <- glm(default ~ balance + income, data = dta,
               family = binomial('logit'), subset = index)
  return(coef(model))
}

Exercise 5.6 - (c)

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

library(boot)
set.seed(6094028)
t0 <- Sys.time()
b <- boot(Default, boot.fn, 1000); b


ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Default, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
         original        bias     std. error
t1* -1.154047e+01 -1.334946e-02 4.311382e-01
t2*  5.647103e-03  5.422527e-06 2.284816e-04
t3*  2.080898e-05 -3.188514e-08 4.623201e-06

plot(b)


Time cost of the bootstrap:  52.67 s

balance的標準誤 = 2.284816e-04
income的標準誤 = 4.623201e-06

Exercise 5.6 - (d)

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

summary(glm_fit_5.6)


Call:
glm(formula = default ~ balance + income, family = binomial("logit"), 
    data = .)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4725  -0.1444  -0.0574  -0.0211   3.7245  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
---
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


ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Default, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
         original        bias     std. error
t1* -1.154047e+01 -1.334946e-02 4.311382e-01
t2*  5.647103e-03  5.422527e-06 2.284816e-04
t3*  2.080898e-05 -3.188514e-08 4.623201e-06

用bootstrap function得到的標準誤與fit glm()得到的標準誤相當接近。

Exercise 5.9

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

library(MASS)
data('Boston')

Exercise 5.9 - (a)

Based on this data set, provide an estimate for the population mean of medv1. Call this estimate $\hat{\mu}$.

mean(Boston$medv)

[1] 22.53281

Exercise 5.9 - (b)

Provide an estimate of the standard error of $\hat{\mu}$. 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.

sqrt(var(Boston$medv) / length(Boston$medv))

[1] 0.4088611

Exercise 5.9 - (c)

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

boot.fn2 <- function(samp, n = length(samp), I) {
  return(sapply(1:I, function(i) {mean(sample(samp, n, replace = TRUE))}))
}

set.seed(6094028)
boot1 <- boot.fn2(Boston$medv, I = 10000)
mean(boot1)

[1] 22.52907

sqrt(var(boot1))

[1] 0.4127674

plot(boot1)

Exercise 5.9 - (d)

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{\mu} - 2SE(\hat{\mu})$, $\hat{\mu} + 2SE(\hat{\mu})$].

CI_b <- round(mean(boot1) + c(-2, 2)*sqrt(var(boot1)), 4)
cat(paste0('95% CI = [', CI_b[1], ', ', CI_b[2], ']\n'))

95% CI = [21.7035, 23.3546]

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

Exercise 5.9 - (e)

Based on this dataset, provide an estimate, $\hat{\mu}_{med}$, for the median value of medv in the population.

median(Boston$medv)

[1] 21.2

Exercise 5.9 - (f)

We now would like to estimate the standard error of $\hat{\mu}_{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.fn3 <- function(samp, n = length(samp), I) {
  return(sapply(1:I, function(i) {median(sample(samp, n, replace = TRUE))}))
}

set.seed(6094028)
boot2 <- boot.fn3(Boston$medv, I = 10000)
mean(boot2)

[1] 21.18475

sqrt(var(boot2))

[1] 0.3821231

plot(boot2)

Exercise 5.9 - (g)

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

quantile(Boston$medv, 0.1)

  10% 
12.75

Exercise 5.9 - (h)

Use the bootstrap to estimate the standard error of $\hat{\mu}_{0.1}$. Comment on your findings.

boot.fn4 <- function(samp, n = length(samp), I) {
  return(sapply(1:I, function(i) {quantile(sample(samp, n, replace = TRUE),  0.1)}))
}

set.seed(6094028)
boot3 <- boot.fn4(Boston$medv, I = 10000)
mean(boot3)

[1] 12.75516

sqrt(var(boot3))

[1] 0.5046076

plot(boot3)

Machine Learning - HW6

Machine Learning - HW6

Author: Jay Liao (ID: RE6094028)

Load in the packages

Exercise 5.5

Exercise 5.5 - (a)

Correlation and scatter plots

Fit the logistic regression

Exercise 5.5 - (b)

[Step 1] Splitting ratio = 3:1.

[Step 2]

[Step 3]

Accuracy score

Exercise 5.5 - (c)

Write a funtion to generalize

Apply the function to different splitting ratios (different training proportions)

Visualization

Exercise 5.5 - (d)

Write a funtion to generalize

Apply the funtion to different variable sets

Exercise 5.6

Exercise 5.6 - (a)

Exercise 5.6 - (b)

Exercise 5.6 - (c)

Exercise 5.6 - (d)

Exercise 5.9

Exercise 5.9 - (a)

Exercise 5.9 - (b)

Exercise 5.9 - (c)

Exercise 5.9 - (d)

Exercise 5.9 - (e)

Exercise 5.9 - (f)

Exercise 5.9 - (g)

Exercise 5.9 - (h)