Lab - 6

library(ISLR2)
library(ISLR)

## 
## Attaching package: 'ISLR'

## The following objects are masked from 'package:ISLR2':
## 
##     Auto, Credit

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(boot)

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.

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

data(Default)
model_lgr <- glm(default ~ income + balance, data = Default, family = 'binomial')
summary(model_lgr)

## 
## Call:
## glm(formula = default ~ income + balance, family = "binomial", 
##     data = Default)
## 
## 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

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

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

set.seed(96)
n <- nrow(Default)
train <- sample(n, n/2)
training_set <- Default[train,]
validation_set <- Default[-train,]

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

### fitting lgr
model_vsa <- glm(default ~ income + balance, data = training_set, family = 'binomial')

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.

probs <- predict(model_vsa, validation_set, type = "response")
predict_class <- ifelse(probs > 0.5, 1, 0)
true_classes <- ifelse(validation_set$default == 'Yes', 1, 0)

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

### validation set error rate
error_rate <- mean(predict_class != true_classes)  
error_rate

## [1] 0.0272

• About 2.72 % validation data is missclassified

(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.

validation_error_rate <- function(seed = 96){
  set.seed(seed)
  n <- nrow(Default)
  train <- sample(n, n/2)
  lg <- glm(default ~ income + balance, data = Default, family = 'binomial', subset = train)
  probs <- predict(lg, Default[-train,], type = "response")  
  predict_class <- ifelse(probs > 0.5, 1, 0)
  true_classes <- ifelse(Default$default[-train] == 'Yes', 1, 0)
  error_rate <- mean(predict_class != true_classes)  # Misclassification rate
  return(error_rate)
}

seed = 96
for (i in 2:4) {
  print(paste("Validation error rate at seed:", seed * i, "is", validation_error_rate(seed * i)))
}

## [1] "Validation error rate at seed: 192 is 0.027"
## [1] "Validation error rate at seed: 288 is 0.0264"
## [1] "Validation error rate at seed: 384 is 0.026"

• the logisitic regression model appears to have a validation error rate around 0.026 to 0.0272.
• This estimate for test error obtained from validation set approach indicate model is predicting very well- 97% accuracy on unseen data.

(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.

### including student variable
validation_error_rate_with_student <- function(seed = 96){
  set.seed(seed)
  n <- nrow(Default)
  train <- sample(n, n/2)
  
  # Fit logistic regression with student variable
  lg <- glm(default ~ income + balance + student, data = Default, family = 'binomial', subset = train)  ## added student variable
  probs <- predict(lg, Default[-train,], type = "response")  # Correct subsetting for validation set
  predict_class <- ifelse(probs > 0.5, 1, 0)
  true_classes <- ifelse(Default$default[-train] == 'Yes', 1, 0)
  error_rate <- mean(predict_class != true_classes)  # Misclassification rate
  return(error_rate)
}

seed = 96
for (i in 2:4) {
  print(paste("Validation error rate (for model including dummy variable student) at seed", seed * i, "is", validation_error_rate_with_student(seed * i)))
}

## [1] "Validation error rate (for model including dummy variable student) at seed 192 is 0.0272"
## [1] "Validation error rate (for model including dummy variable student) at seed 288 is 0.0262"
## [1] "Validation error rate (for model including dummy variable student) at seed 384 is 0.0256"

• Around 2.56% to 2.72% of the validation data is misclassified by the logisitic regression when student variable is added
• This shows no improvement when compared to validation error rates of logistic regression when student variable is used.

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.

(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.

set.seed(96)

# Fit the logistic regression model
glm_fit <- glm(default ~ income + balance, data = Default, family = binomial)

# Display model summary to get standard errors
summary(glm_fit)

## 
## Call:
## glm(formula = default ~ income + balance, family = binomial, 
##     data = Default)
## 
## 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

• Std. Error of coefficient of income in model glm is 4.985e-06 and for coeff. of balance is 2.274e-04

(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(data, index) {
  # Fit logistic regression model on bootstrapped sample
  boot_model <- glm(default ~ income + balance, data = data, family = binomial, subset = index)
  
  # Return coefficients of interest (income and balance)
  return(coef(boot_model)[2:3])  
}

(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.

# Perform bootstrap with 1000 resamples
set.seed(96)  # Ensuring reproducibility
boot_results <- boot(Default, boot.fn, R = 1000)

# Display estimated standard errors from bootstrap
boot_results

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Default, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##         original        bias     std. error
## t1* 2.080898e-05 -1.770689e-07 4.838230e-06
## t2* 5.647103e-03  1.103567e-05 2.309823e-04

boot_se <- apply(boot_results$t, 2, sd)
boot_se

## [1] 4.838230e-06 2.309823e-04

• Using bootstrapping method with R = 1000, we get Std. Error of coefficient of income as 4.838230e-06 and for coeff. of balance as 2.309823e-04

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

Coefficient	SE (glm)	SE (Bootstrap)
Income	4.985e-06	4.838230e-06
Balance	2.274e-04	2.309823e-04

• Clearly, values not the exact match as we would expect but they are not to off to each other as well.
• That is a good sign, since SE estimate from summary method makes some assumptions to calculate SE, where as Bootstrapping doesn’t make such assumptions
• Since the results are consistent across both methods, we can be confident that the logistic regression model is well-specified for estimating the relationship between default probability and the predictors income and balance.

7 In Sections 5.3.2 and 5.3.3, we saw that the cv.glm() function can be used in order to compute the LOOCV test error estimate. Alternatively, one could compute those quantities using just the glm() and predict.glm() functions, and a for loop. You will now take this approach in order to compute the LOOCV error for a simple logistic regression model on the Weekly data set. Recall that in the context of classification problems, the LOOCV error is given in (5.4).

summary(Weekly)

##       Year           Lag1               Lag2               Lag3         
##  Min.   :1990   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950  
##  1st Qu.:1995   1st Qu.: -1.1540   1st Qu.: -1.1540   1st Qu.: -1.1580  
##  Median :2000   Median :  0.2410   Median :  0.2410   Median :  0.2410  
##  Mean   :2000   Mean   :  0.1506   Mean   :  0.1511   Mean   :  0.1472  
##  3rd Qu.:2005   3rd Qu.:  1.4050   3rd Qu.:  1.4090   3rd Qu.:  1.4090  
##  Max.   :2010   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260  
##       Lag4               Lag5              Volume            Today         
##  Min.   :-18.1950   Min.   :-18.1950   Min.   :0.08747   Min.   :-18.1950  
##  1st Qu.: -1.1580   1st Qu.: -1.1660   1st Qu.:0.33202   1st Qu.: -1.1540  
##  Median :  0.2380   Median :  0.2340   Median :1.00268   Median :  0.2410  
##  Mean   :  0.1458   Mean   :  0.1399   Mean   :1.57462   Mean   :  0.1499  
##  3rd Qu.:  1.4090   3rd Qu.:  1.4050   3rd Qu.:2.05373   3rd Qu.:  1.4050  
##  Max.   : 12.0260   Max.   : 12.0260   Max.   :9.32821   Max.   : 12.0260  
##  Direction 
##  Down:484  
##  Up  :605  
##            
##            
##            
## 

(a) Fit a logistic regression model that predicts Direction using Lag1and Lag2.

# set.seed(96)
lm_full <- glm(Direction ~ Lag1 + Lag2, data = Weekly, family = 'binomial')

(b) Fit a logistic regression model that predicts Direction using Lag1 and Lag2 using all but the first observation.

lm_1 <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-1,], family = 'binomial')
summary(lm_1)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = Weekly[-1, 
##     ])
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.22324    0.06150   3.630 0.000283 ***
## Lag1        -0.03843    0.02622  -1.466 0.142683    
## Lag2         0.06085    0.02656   2.291 0.021971 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1494.6  on 1087  degrees of freedom
## Residual deviance: 1486.5  on 1085  degrees of freedom
## AIC: 1492.5
## 
## Number of Fisher Scoring iterations: 4

(c) Use the model from (b) to predict the direction of the first observation. You can do this by predicting that the first observation will go up if P (Direction = “Up”|Lag1, Lag2) > 0.5. Was this observation correctly classified?

prob <- predict(lm_1, Weekly[1,], type = 'response')
pred_class <- ifelse(prob > 0.5, 'Up', 'Down')
print(paste('The Predicted Direction class of the first observation is:', pred_class, "::and the true class is:", Weekly[1, 'Direction']))

## [1] "The Predicted Direction class of the first observation is: Up ::and the true class is: Down"

• No, the model has misclassified the first observation

(d) Write a for loop from i = 1 to i= n, where n is the number ofobservations in the data set, that performs each of the following steps:

1. Fit a logistic regression model using all but the ith observation to predict Direction using Lag1 and Lag2.
2. Compute the posterior probability of the market moving up for the ith observation.
3. Use the posterior probability for the ith observation in order to predict whether or not the market moves up.
4. Determine whether or not an error was made in predicting the direction for the ith observation. If an error was made, then indicate this as a 1, and otherwise indicate it as a 0.

n <- nrow(Weekly)
errors <- rep(0, n)

for(i in 1:n){
  lm_i <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-i,], family = 'binomial')
  prob <- predict(lm_i, Weekly[i,], type = 'response')
  pred_class <- ifelse(prob > 0.5, 'Up', 'Down')
  true_class <- Weekly[i, 'Direction']
  errors[i] <- ifelse(pred_class != true_class, 1, 0)
}

print(paste('Total errors are:', sum(errors), "out of", n))

## [1] "Total errors are: 490 out of 1089"

(e) Take the average of the n numbers obtained in (d)iv in order to obtain the LOOCV estimate for the test error. Comment on the results.

validation_error_rate <- sum(errors)/n
validation_error_rate

## [1] 0.4499541

• The estimated test error rate obtained from LOOCV method is ~0.45.
• This imply that the model is performing slight better than randomly predicting the Market direction

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

(a) Based on this data set, provide an estimate for the populationmean of medv. Call this estimateˆµ.

data(Boston)
mu_hat <- mean(Boston$medv)
mu_hat

## [1] 22.53281

• The estimated mean of medv (median value of owner-occupied homes) is 22.53.
  
• This represents the average home value (in $1000s) in the Boston housing dataset.

(b) Provide an estimate of the standard error ofˆµ. 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.

n <- nrow(Boston)
mu_sd <- sd(Boston$medv)

mu_se <- mu_sd/sqrt(n)
mu_se

## [1] 0.4088611

• Since, the standard error represents the variability of the sample mean and a small SE suggests a more precise estimate of the population mean.

(c) Now estimate the standard error ofˆ µ using the bootstrap. How does this compare to your answer from (b)?

# Define bootstrap function
boot_mean <- function(data, index) {
  return(mean(data[index]))
}

# Perform bootstrap with 1000 resamples
set.seed(96)
boot_res <- boot(Boston$medv, boot_mean, R = 1000)

# Bootstrap SE
boot_se_mu <- sd(boot_res$t)
boot_se_mu

## [1] 0.3996253

• The bootstrap SE (0.3996) is very close to the analytically computed SE, confirming the reliability of both methods.
  
• This similarity indicates that the sample is large enough and follows normality assumptions well.

(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 \([ˆ µ− 2SE(ˆ µ),ˆµ + 2SE(ˆ µ)].\)

# Approximate 95% Confidence Interval
ci_boot <- c(mu_hat - 2 * boot_se_mu, mu_hat + 2 * boot_se_mu)
ci_boot

## [1] 21.73356 23.33206

t.test(Boston$medv)$conf.int

## [1] 21.72953 23.33608
## attr(,"conf.level")
## [1] 0.95

• The confidence interval obtained from t-test has more wider range than bootstrap confindence interval

(e) Based on this data set, provide an estimate,ˆµmed, for the median value of medv in the population.

mu_med <- median(Boston$medv)
mu_med

## [1] 21.2

(f) We now would like to estimate the standard error ofˆµ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 <- function(data, index) {
  return(median(data[index]))
}

set.seed(96)
boot_med_res <- boot(Boston$medv, boot_median, R = 1000)

boot_se_med <- sd(boot_med_res$t)
boot_se_med

## [1] 0.3649334

• Surprisingly, the SE for median is lower than that of SE of mean of medv.
• This suggests that the median is a more stable measure of central tendency than the mean, particularly in the presence of outliers.

(g) Based on this data set, provide an estimate for the tenth percentile of medv in Boston census tracts. Call this quantityˆµ0.1. (You can use the quantile() function.)

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

##   10% 
## 12.75

• The 10th percentile value of medv is 12.75.
  
• This means that 10% of the houses in the dataset are valued below $12,750.
  
• This helps in understanding the lower-end of housing prices.

(h) Use the bootstrap to estimate the standard error ofˆ µ0.1. Comment on your findings.

boot_percentile <- function(data, index) {
  return(quantile(data[index], 0.1))
}

set.seed(96)
boot_pctl_res <- boot(Boston$medv, boot_percentile, R = 1000)


boot_se_0.1 <- sd(boot_pctl_res$t)
boot_se_0.1

## [1] 0.4872043

• The SE for the 10th percentile is 0.4872, which is higher than that for the mean and median.
  

• This suggests that the lower tail of the distribution has more variability, meaning housing values in the bottom 10% fluctuate more than central values.

• Summary Table:

Estimate	Value	Bootstrap SE
Mean (µ)	22.53281	0.3996253
Median (µmed)	21.2	0.3649334
10th Percentile (µ0.1)	12.75	0.4872043

• The mean and median estimates are stable, with a slight skew in home values.
  
• The bootstrap SE confirms the accuracy of analytical SE for the mean.
  
• Lower-end home prices are more variable, as indicated by a higher SE for the 10th percentile.
  
• The bootstrap method provides reliable standard errors, even for statistics without simple formulas (median & percentile).