a. Fit a logistic regression model that uses
income and balance to predict
default.
library(ISLR)## Warning: package 'ISLR' was built under R version 4.4.3set.seed(1234)
data(Default)model1 <- glm(default ~ income + balance, data = Default, family = binomial)
summary(model1)##
## 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: 8b. Using the validation set approach, estimate the test error of this model. In order to do this, you must perform the following steps:
Split the sample set into a training set and a validation set.
Fit a multiple logistic regression model using only the training observations.
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.
Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
n <- nrow(Default)
# Split: 70% training, 30% validation
train_index <- sample(1:n, size = round(0.7 * n))
train_data <- Default[train_index, ]
validation_data <- Default[-train_index, ]# Fit the model on training data
model_train <- glm(default ~ income + balance, data = train_data, family = binomial)# Get predicted probabilities for validation set
pred_probs <- predict(model_train, newdata = validation_data, type = "response")# Classify based on a threshold of 0.5
pred_default <- ifelse(pred_probs > 0.5, "Yes", "No")
#pred_default# Compute the validation error rate
validation_error <- mean(pred_default != validation_data$default)
cat("Validation Error Rate (income + balance):", validation_error, "\n")## Validation Error Rate (income + balance): 0.025c. 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.
set.seed(1234) # resetting seed to ensure reproducibility of splits
errors <- numeric(3)
for (i in 1:3) {
train_index <- sample(1:n, size = round(0.7 * n))
train_data <- Default[train_index, ]
validation_data <- Default[-train_index, ]
model_train <- glm(default ~ income + balance, data = train_data, family = binomial)
pred_probs <- predict(model_train, newdata = validation_data, type = "response")
pred_default <- ifelse(pred_probs > 0.5, "Yes", "No")
errors[i] <- mean(pred_default != validation_data$default)
cat("Split", i, "Validation Error Rate (income + balance):", errors[i], "\n")
}## Split 1 Validation Error Rate (income + balance): 0.025
## Split 2 Validation Error Rate (income + balance): 0.02366667
## Split 3 Validation Error Rate (income + balance): 0.021cat("Average Validation Error Rate (income + balance):", mean(errors), "\n")## Average Validation Error Rate (income + balance): 0.02322222d. 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.
errors_student <- numeric(3)
for (i in 1:3) {
train_index <- sample(1:n, size = round(0.7 * n))
train_data <- Default[train_index, ]
validation_data <- Default[-train_index, ]
model_train_student <- glm(default ~ income + balance + student, data = train_data, family = binomial)
pred_probs_student <- predict(model_train_student, newdata = validation_data, type = "response")
pred_default_student <- ifelse(pred_probs_student > 0.5, "Yes", "No")
errors_student[i] <- mean(pred_default_student != validation_data$default)
cat("Split", i, "Validation Error Rate (income + balance + student):", errors_student[i], "\n")
}## Split 1 Validation Error Rate (income + balance + student): 0.02933333
## Split 2 Validation Error Rate (income + balance + student): 0.02666667
## Split 3 Validation Error Rate (income + balance + student): 0.02433333cat("Average Validation Error Rate (with student):", mean(errors_student), "\n")## Average Validation Error Rate (with student): 0.02677778The results indicate that when you add the student dummy variable to the logistic regression model, the average validation error rate increases slightly.
This suggests that including the student variable:
– Does not reduce the test error rate.
– Slightly increases the error rate, indicating that the extra variable
may not provide predictive benefit in this context.
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(1)
fit_glm <- glm(default ~ income + balance, data = Default, family = binomial)
summary(fit_glm)##
## 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: 8b. 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.
library(boot)boot.fn <- function(data, index) {
coef(glm(default ~ income + balance, data = data, subset = index, family = binomial))[c("income", "balance")]
}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.
boot_results <- boot(Default, boot.fn, R = 1000)
print(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.680317e-07 4.866284e-06
## t2* 5.647103e-03 1.855765e-05 2.298949e-04d. Comment on the estimated standard errors obtained using
the glm() function and using your bootstrap
function.
The estimated standard errors for income and balance obtained using the summary() function (i.e., the analytically derived values) are nearly identical to the bootstrap-based estimates.
This similarity indicates that the standard formula used by glm() is reliable and that the bootstrap procedure confirms these low standard errors.
Minor numerical differences may appear due to sampling variability in bootstrap resampling, but overall the two methods agree closely
a. Fit a logistic regression model that predicts
Direction using Lag1 and
Lag2.
set.seed(12)
data(Weekly)fit_full <- glm(Direction ~ Lag1 + Lag2, data = Weekly, family = binomial)
summary(fit_full)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22122 0.06147 3.599 0.000319 ***
## Lag1 -0.03872 0.02622 -1.477 0.139672
## Lag2 0.06025 0.02655 2.270 0.023232 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1488.2 on 1086 degrees of freedom
## AIC: 1494.2
##
## Number of Fisher Scoring iterations: 4b. Fit a logistic regression model that predicts
Direction using Lag1 and
Lag2using all but the first
observation.
fit_minus1 <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-1, ], family = binomial)
summary(fit_minus1)##
## 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: 4c. 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(P(Direction="Up" | Lag1 , Lag2)>0.5.
Was this observation correctly classified?
pred_prob_first <- predict(fit_minus1, newdata = Weekly[1, ], type = "response")
# Predict "Up" if probability > 0.5, otherwise "Down"
pred_first <- ifelse(pred_prob_first > 0.5, "Up", "Down")
cat("For the first observation:\n")## For the first observation:cat("Predicted Direction:", pred_first, "\n")## Predicted Direction: Upcat("Actual Direction: ", as.character(Weekly$Direction[1]), "\n")## Actual Direction: Down# compare whether the predicted direction matches the actual.
if(pred_first == Weekly$Direction[1]){
cat("The first observation was correctly classified.\n\n")
} else {
cat("The first observation was misclassified.\n\n")
}## The first observation was misclassified.d. Write a for loop from i=1i=1 to i=ni=n, where nn is the number of observations in the data set, that performs each of the following steps:
Fit a logistic regression model using all but
the iith observation to predict Direction
using Lag1 and Lag2 .
Compute the posterior probability of the market moving up for the iith observation.
Use the posterior probability for the iith observation in order to predict whether or not the market moves up.
Determine whether or not an error was made in predicting the direction for the iith observation. If an error was made, then indicate this as a 1, and otherwise indicate it as a 0.
n <- nrow(Weekly)
cv_errors <- rep(0, n) # to store error indicators: 1 if misclassified, 0 otherwise
for (i in 1:n) {
# Fit the logistic regression model using all but the i-th observation
fit_loocv <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-i, ], family = binomial)
# Compute predicted probability for the i-th observation
prob_i <- predict(fit_loocv, newdata = Weekly[i, ], type = "response")
# Classify as "Up" if probability > 0.5 otherwise "Down"
pred_direction <- ifelse(prob_i > 0.5, "Up", "Down")
# Check if the classification is an error (1 if error, 0 if correct)
cv_errors[i] <- ifelse(pred_direction == Weekly$Direction[i], 0, 1)
}e. Take the average of the n numbers obtained in (d) in order to obtain the LOOCV estimate for the test error. Comment on the results.
loocv_error <- mean(cv_errors)
cat("LOOCV Estimated Test Error:", loocv_error, "\n")## LOOCV Estimated Test Error: 0.4499541a. Based on this data set, provide an estimate for the
population mean of medv. Call this estimate \(\mu\)
library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.4.3##
## Attaching package: 'ISLR2'## The following objects are masked from 'package:ISLR':
##
## Auto, Creditlibrary(boot)
data(Boston)mu_hat <- mean(Boston$medv)
print(paste("Estimate for population mean of medv:", mu_hat))## [1] "Estimate for population mean of medv: 22.5328063241107"b. Provide an estimate of the standard error of \[\mu\]. Interpret this result.
se_formula <- sd(Boston$medv)/sqrt(nrow(Boston))
print(paste("Standard error of mu_hat using formula:", se_formula))## [1] "Standard error of mu_hat using formula: 0.408861147497535"c. Now estimate the standard error of \[\mu\] using the bootstrap. How does this compare to your answer from (b)?
set.seed(1)
boot_mean <- function(data, indices) {
return(mean(data[indices]))
}
boot_results <- boot(Boston$medv, boot_mean, R=1000)
se_boot <- sd(boot_results$t)
print(paste("Bootstrap estimate of standard error:", se_boot))## [1] "Bootstrap estimate of standard error: 0.410662153886544"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).
ci_lower <- mu_hat - 2*se_boot
ci_upper <- mu_hat + 2*se_boot
print(paste("95% CI using bootstrap: [", ci_lower, ",", ci_upper, "]"))## [1] "95% CI using bootstrap: [ 21.7114820163376 , 23.3541306318838 ]"# Compare with t.test
t_test_result <- t.test(Boston$medv)
print("95% CI using t.test:")## [1] "95% CI using t.test:"print(t_test_result$conf.int)## [1] 21.72953 23.33608
## attr(,"conf.level")
## [1] 0.95e. Based on this data set, provide an estimate,
\[\hat{\mu}_{med}\], for the
median value of medv in the population.
mu_med <- median(Boston$medv)
print(paste("Estimate for median of medv:", mu_med))## [1] "Estimate for median of medv: 21.2"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.
set.seed(87)
boot_median <- function(data, indices) {
return(median(data[indices]))
}
boot_median_results <- boot(Boston$medv, boot_median, R=1000)
se_median <- sd(boot_median_results$t)
print(paste("SE of median using bootstrap:", se_median))## [1] "SE of median using bootstrap: 0.399821331468316"g. Based on this data set, provide an estimate for the tenth
percentile of medv in Boston census tracts. Call this
quantity \(\hat{\mu}_{0.1}\)
.(You can use the quantile()
function.)
mu_0.1 <- quantile(Boston$medv, 0.1)
print(paste("10th percentile of medv:", mu_0.1))## [1] "10th percentile of medv: 12.75"h. Use the bootstrap to estimate the standard error of \[\hat{\mu}_{0.1}\] Comment on your findings.
set.seed(98)
boot_quantile <- function(data, indices) {
return(quantile(data[indices], 0.1))
}
boot_quantile_results <- boot(Boston$medv, boot_quantile, R=1000)
se_quantile <- sd(boot_quantile_results$t)
print(paste("SE of 10th percentile using bootstrap:", se_quantile))## [1] "SE of 10th percentile using bootstrap: 0.493784569119673"