library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.5.3library(mlbench) # For the PimaIndiansDiabetes dataset## Warning: package 'mlbench' was built under R version 4.5.3library(tidyverse) # For data manipulation and ggplot2## Warning: package 'ggplot2' was built under R version 4.5.2## Warning: package 'dplyr' was built under R version 4.5.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(modelr) # For elegant cross-validation infrastructure
library(MASS) # For lda() and qda() functions##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## The following object is masked from 'package:ISLR2':
##
## Bostonlibrary(class) # For knn() functions
library(gtsummary) # For statistical summary tables## Warning: package 'gtsummary' was built under R version 4.5.3##
## Attaching package: 'gtsummary'
##
## The following object is masked from 'package:MASS':
##
## selectlibrary(pROC) # For validation diagnostics## Warning: package 'pROC' was built under R version 4.5.3## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, varlibrary(DT) # For cool interactive tables## Warning: package 'DT' was built under R version 4.5.3library(boot)#Chapter 5 question #3 One fold of the data is held as the test set and the remaining is set as the training set. This process is repeated K times, with each fold serving as the validation set exactly once. The prediction error is calculated for each iteration and the average of the K errors is given as the estimated test error.
#validation approach Advantages of K fold is that we can help produce a more stable and reliable estimate of test error, has lower variablity, more efficient use of limited data. Disadvantages of K fold can be more demanding computationally since we have to fit the model K times, a little more complex to setup and implement into data.
#LOOVC Advatages of K fold CV runs alot faster since the we are using the model fit only 5 or 10 times instead of n times. Along with much lower variance than LOOVC. Disadvantages of K fold CV we can come out with a slightly higher bias than LOOCV since we have less observations.
#chapter 5 question 5
data(Default)
# Fit logistic regression model
fit_log <- glm(
default ~ income + balance,
data = Default,
family = binomial
)
summary(fit_log)##
## 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#set random seed first
set.seed(101)
train_index <- sample(
1:nrow(Default),
nrow(Default) / 2
)
train <- Default[train_index, ]
test <- Default[-train_index, ]
fit_log <- glm(
default ~ income + balance,
data = train,
family = binomial
)
# Predicted probabilities
prob_default <- predict(
fit_log,
newdata = test,
type = "response"
)
# Convert probabilities into class labels
pred_default <- ifelse(prob_default > 0.5,
"Yes",
"No")
validation_error <- mean(pred_default != test$default)
validation_error## [1] 0.025#split 2
set.seed(202)
train_index <- sample(
1:nrow(Default),
nrow(Default)/2
)
train <- Default[train_index, ]
test <- Default[-train_index, ]
fit_log <- glm(default ~ income + balance,
data=train,
family=binomial)
prob <- predict(fit_log,
newdata=test,
type="response")
pred <- ifelse(prob > 0.5,"Yes","No")
mean(pred != test$default)## [1] 0.0276#split 3
set.seed(303)
train_index <- sample(
1:nrow(Default),
nrow(Default)/2
)
train <- Default[train_index, ]
test <- Default[-train_index, ]
fit_log <- glm(default ~ income + balance,
data=train,
family=binomial)
prob <- predict(fit_log,
newdata=test,
type="response")
pred <- ifelse(prob > 0.5,"Yes","No")
mean(pred != test$default)## [1] 0.028The results are very similar but not exactly identical due to the random split creating a different training and validation set each seed.
set.seed(101)
train_index <- sample(
1:nrow(Default),
nrow(Default)/2
)
train <- Default[train_index, ]
test <- Default[-train_index, ]
fit_log2 <- glm(
default ~ income + balance + student,
data = train,
family = binomial
)
prob2 <- predict(
fit_log2,
newdata = test,
type = "response"
)
pred2 <- ifelse(prob2 > 0.5,
"Yes",
"No")
validation_error2 <- mean(pred2 != test$default)
validation_error2## [1] 0.0242We a got a very similar validation error rate for the model just slighlty lower, so even adding student variable does not a very meaningful validation error to where we can say student variable has an impact.
#chapter 5 question 6
data(Default)
set.seed(101)
fit_log <- glm(
default ~ income + balance,
data = Default,
family = binomial
)
summary(fit_log)##
## 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: 8boot.fn <- function(data, index){
fit <- glm(
default ~ income + balance,
data = data,
subset = index,
family = binomial
)
return(coef(fit)[2:3])
}
set.seed(101)
boot_results <- boot(
data = Default,
statistic = boot.fn,
R = 1000
)
boot_results##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Default, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 2.080898e-05 -5.207761e-08 4.854883e-06
## t2* 5.647103e-03 1.862213e-05 2.241446e-04apply(boot_results$t, 2, sd)## [1] 4.854883e-06 2.241446e-04#chpater 5 question 9
data(Boston)
set.seed(101)
#Mean
mu_hat <- mean(Boston$medv)
mu_hat## [1] 22.53281#SD
n <- nrow(Boston)
s <- sd(Boston$medv)
se_mu <- s / sqrt(n)
se_mu## [1] 0.4088611boot.mean <- function(data, index){
mean(data[index])
}
set.seed(101)
boot_mean <- boot(
data = Boston$medv,
statistic = boot.mean,
R = 1000
)
boot_mean##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.mean, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.01271779 0.4073505boot_se <- sd(boot_mean$t)
boot_se## [1] 0.4073505#theoretical formula accurately estimates the sampling variability of the mean. when we compare the results from part b to the bootstrap theyre very close to the analytical estimate
lower <- mu_hat - 2 * boot_se
upper <- mu_hat + 2 * boot_se
c(lower, upper)## [1] 21.71811 23.34751t.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#this also confirms that the bootstrap interval is similar to ttest meaning it provides reliable estimate of the uncertainty associated with the sample mean.
median_hat <- median(Boston$medv)
median_hat## [1] 21.2boot.median <- function(data, index){
median(data[index])
}
set.seed(101)
boot_median <- boot(
data = Boston$medv,
statistic = boot.median,
R = 1000
)
boot_median##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.median, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.0084 0.3795438median_se <- sd(boot_median$t)
median_se## [1] 0.3795438