Project Presentation.Rmd

---
title: "Predictive models for cumulative COVID-19 cases in USA"
author: "Minsu Kim, Ariane Stark"
output: 
  beamer_presentation:
    theme: "Madrid"
    slide_level: 2
    toc: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
knitr::opts_chunk$set(comment = NA)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)

# Ariane packages
library(readr)
library(readxl)
library(tidyverse)
library(rpart) 
library(rpart.plot)
library(leaps)
library(glmnet)
library(cowplot)
library(kableExtra)

# Minsu packages
# install.packages("reshape2")
library(reshape2)
library(ggplot2)
library(scales)
#install.packages("patchwork")
library(patchwork)
library(class)
library(MASS)
# install.packages("gridExtra")
library(gridExtra)
library(grid)
library(tree)
```


```{r include=FALSE}
dat <- read_csv("case_and_demographics.csv") %>%
  subset(select = -c(1)) %>% 
  subset(select = c(1,3:4,2,5:6,43,56:61 )) %>% 
  mutate(Cum_Case_Median = cut_number(cum_case,2, labels= c(0,1))) %>% 
  na.omit() #removes Rio Arriba County NM which has NA in the pov. dem.
```

# Introduction 

## Background

- COVID-19 is a coronavirus identified in 2019 and has caused a pandemic of respiratory illness.
- It has been spreading rapidly since the first case of COVID-19 was reported Dec. 1, 2019.
- Vaccines have recently been available, but it is too early to be relieved due to concerns about the mutant virus.
- There is still a need to keep an eye on the trends of cases.

\ \ 

\begin{columns}[onlytextwidth,T]
\begin{column}{.45\linewidth}
    \includegraphics[width=0.8\linewidth]{CoronaVirus.jpeg}
    \end{column}
    \begin{column}{.45\linewidth}
    \includegraphics[width=0.8\linewidth]{covid19_county_view.png}
    \end{column}
 \end{columns}   
    \begingroup
  \fontsize{5}{7}\selectfont
  (Figure: https://www.denhamblythe.com/covid-19-health-statement)
  (Figure:https://covid.cdc.gov/covid-data-tracker/#datatracker-home)
  \endgroup

## Data Set Overview

The data set contains cumulative COVID-19 case counts as of April 1, 2021 for each of the 3,220 counties in the United States. Also we have racial and ethnicity demographics, poverty demographics, population demographics, and age demographics from government census data. Note Rio Arriba County, New Mexico is missing poverty demographics and has thus been excluded from analysis, leaving our data set with 9 predictors and 3219 observations.

\begingroup
\fontsize{8}{10}\selectfont
(Source:https://covid19.census.gov)
\endgroup


## Our Variables
- Cumulative COVID-19 Cases
- Total Population
- Population Density
- Median Age
- Total Number of White Residents
- Total Number of Black Residents
- Total Number of Other (Not strictly White or Black) Residents
- Total Number of Hispanic Residents
- Total Number of Households Above the Poverty Line
- Total Number of Households Below the Poverty Line

## Our Variables: Plots (Response)
```{r fig.height=4}
dat %>% ggplot() +
  geom_boxplot(aes(x=cum_case)) +
  scale_x_log10() +
  xlab("Log Base 10 Scale of Cumulative Cases")
```


```{r}
print(summary(dat$cum_case))
```


## Our Variables: Plots (Predictors)

```{r, fig.align='center', out.width = '90%'}
dat <- read_csv("case_and_demographics.csv") %>%
  subset(select = -c(1)) %>% 
  subset(select = c(1,3:4,2,5:6,43,56:61 )) %>% 
  na.omit() #row 1817 is removed 
dat1 <- dat[, -c(1:3)]

## distribution plots
d <- melt(dat1[, -c(1)])
raw.plot <- ggplot(d, aes(x = value)) +
  facet_wrap(~variable, scales = "free_x") +
  geom_histogram()
raw.plot  #all but median_age are highly right skewed 
```



## Study Objectives

- Select the demographics that yields a good predictive model for cumulative COVID-19 cases at the county level.

- Fit classification models to predict whether a given county has an instance rate higher than the median of the instance rates.


# Simple Linear Regression

```{r, fig.align='center', out.width = '90%'}
#scatter plots with simple linear regression
dat1 %>% gather(-cum_case, key = "var", value = "value") %>% 
  ggplot(aes(x = value, y = cum_case)) +
    geom_point() +
    geom_smooth(method='lm')+
    facet_wrap(~ var, scales = "free") +
    theme_bw()
```


# Subset Selection

## Cross Validation of Best Subset, Forward and Backward Stepwise 

10-fold cross validation of data was done to find the optimal number of parameters for the model selected using each method. This is done to balance the model fitting the data well and the model having decent predictive accuracy.


```{r}
predict.regsubsets <- function (object ,newdata ,id ,...){
  form<-as.formula(object$call [[2]])
  mat<-model.matrix(form,newdata)
  coefi<-coef(object ,id=id)
  xvars<-names(coefi)
  mat[,xvars]%*%coefi
}
```

```{r}
set.seed(100)

k<-10

folds<-sample(1:k,nrow(dat[,c(4:13)]),replace=TRUE)
cv.errors.best_subset <- matrix(NA,k,9, dimnames=list(NULL, paste(1:9)))
cv.errors.forward <- matrix(NA,k,9, dimnames=list(NULL, paste(1:9)))
cv.errors.backward <- matrix(NA,k,9, dimnames=list(NULL, paste(1:9)))


for(j in 1:k){
  # best subset select
  best_subset.fit <- regsubsets(cum_case~.,data = dat[folds!=j,c(4:13)],
                                method = "exhaustive", nvmax = 9)
  # forward stepwise select
  forward.fit <- regsubsets(cum_case~.,data = dat[folds!=j,c(4:13)],
                                method = "forward", nvmax = 9)
  # backward stepwise select
  backward.fit <- regsubsets(cum_case~.,data = dat[folds!=j,c(4:13)],
                                method = "backward", nvmax = 9)
  
  for(i in 1:9){
    # best subset error
    pred.best <- 
      predict.regsubsets(best_subset.fit, dat[folds==j,c(4:13)], id=i)
    cv.errors.best_subset[j,i] <-
      mean((dat$cum_case[folds==j]-pred.best)^2)
    
    # forward stepwise error
    pred.forward <- 
      predict.regsubsets(forward.fit, dat[folds==j,c(4:13)], id=i)
    cv.errors.forward[j,i] <-
      mean((dat$cum_case[folds==j]-pred.forward)^2)
    
    # backward stepwise error
    pred.backward <- 
      predict.regsubsets(backward.fit, dat[folds==j,c(4:13)], id=i)
    cv.errors.backward[j,i] <-
      mean((dat$cum_case[folds==j]-pred.backward)^2)
  
    
  } }
```

##  

```{r}
#best subset CV model selection
mean.cv.errors.best <- apply(cv.errors.best_subset,2,mean)

plot.best <- ggplot(mapping=aes(x=(1:length(mean.cv.errors.best)),
                           y = mean.cv.errors.best)) +
  geom_point(size=2) +
  geom_line()+
  geom_point(mapping = aes(x=which.min(mean.cv.errors.best),
       y=mean.cv.errors.best[which.min(mean.cv.errors.best)]),
       color="red", size=2.2)+
  scale_x_continuous(breaks = seq(1,9,2))+
  xlab("Number of Variables")+
  ylab("Mean C.V. Error") +
  labs(title="Best Subset Selection")

#forward stepwise CV model selection
mean.cv.errors.forward <- apply(cv.errors.forward,2,mean)

plot.forward <- ggplot(mapping=aes(x=(1:length(mean.cv.errors.forward)),
                           y = mean.cv.errors.forward)) +
  geom_point(size=2) +
  geom_line()+
  geom_point(mapping = aes(x=which.min(mean.cv.errors.forward),
       y=mean.cv.errors.forward[which.min(mean.cv.errors.forward)]),
       color="red", size=2.2)+
  scale_x_continuous(breaks = seq(1,9,2))+
  xlab("Number of Variables")+
  ylab("Mean C.V. Error") +
  labs(title="Forward Stepwise Selection")


#backward stepwise CV model selection
mean.cv.errors.backward <- apply(cv.errors.backward,2,mean)
plot.backward <- ggplot(mapping=aes(
  x=(1:length(mean.cv.errors.backward)),
                           y = mean.cv.errors.backward)) +
  geom_point(size=2) +
  geom_line()+
  geom_point(mapping = aes(x=which.min(mean.cv.errors.backward),
       y=mean.cv.errors.backward[which.min(mean.cv.errors.backward)]),
       color="red", size=2.2)+
  scale_x_continuous(breaks = seq(1,9,2))+
  xlab("Number of Variables")+
  ylab("Mean C.V. Error") +
  labs(title="Backward Stepwise Selection")


plot_grid(plot.best,plot.forward,plot.backward)

```



## Best Subset Overall

Recall 10 fold cross validation of best subset selection selected the model with 2 parameters as having the lowest average cross validation predictive error.
```{r}
best_subset.fit <- regsubsets(cum_case~.,data = dat[,c(4:13)],
                              method = "exhaustive", nvmax = 9)
best_subset.summary <- summary(best_subset.fit)
best_subset.summary.frame <- 
  data_frame("Parameters" = seq(1:9),
             "R^2" = round(best_subset.summary$rsq,4),
             "AdjR^2" = round(best_subset.summary$adjr2,4),
             "CP" = round(best_subset.summary$cp,4),
             "BIC" = round(best_subset.summary$bic,4))%>% 
  kable() %>% 
  kable_classic()

best_subset.summary.frame

```



## Forward Stepwise Overall

Recall 10 fold cross validation of forward stepwise selection selected the model with 2 parameters as having the lowest average cross validation predictive error.

```{r}
forward.fit <- regsubsets(cum_case~.,data = dat[,c(4:13)],
                              method = "forward", nvmax = 9)
forward.summary <- summary(forward.fit)
forward.summary.frame <- 
  data_frame("Parameters" = seq(1:9),
             "R^2" = round(forward.summary$rsq,4),
             "AdjR^2" = round(forward.summary$adjr2,4),
             "CP" = round(forward.summary$cp,4), 
             "BIC" = round(forward.summary$bic,4)) %>% 
  kable() %>% 
  kable_classic()

forward.summary.frame


```



## Backward Stepwise Overall

Recall 10 fold cross validation of backward stepwise selection selected the model with 6 parameters as having the lowest average cross validation predictive error.
```{r}
backward.fit <- regsubsets(cum_case~.,data = dat[,c(4:13)],
                              method = "backward", nvmax = 9)
backward.summary <- summary(backward.fit)
backward.summary.frame <- 
  data_frame("Parameters" = seq(1:9),
             "R^2" = round(backward.summary$rsq,4),
             "AdjR^2" = round(backward.summary$adjr2,4),
             "CP" = round(backward.summary$cp,4),
             "BIC" = round(backward.summary$bic,4))%>% 
  kable() %>% 
  kable_classic()

backward.summary.frame
```

## Model Summary
**Best Subset:**
The model has parameters: Total Population, and Total Hispanic.

Cumulative Cases = 315.11931521 + 0.07300349(Total Population) + 0.08304726(Total Hispanic) 
```{r include=FALSE}
coef(best_subset.fit,2)
```

**Forward Stepwise:**
The model has parameters: Total Population, Total White, and Total Hispanic.

Cumulative Cases = 315.11931521 + 0.07300349(Total Population) + 0.08304726(Total Hispanic) 
```{r include=FALSE}
coef(forward.fit,2)
```


## 

**Backward Stepwise:**
The model has parameters: Total Population, Total White, Total Black, Total Hispanic, Total Other, and Total Number of Households Above Poverty.

Cumulative Cases = -338.60412493 -0.38868473(Total Population) +0.58210948(Total White) +0.55858828(Total Black) +0.05514001(Total Hispanic) +0.51284419(Total Other) -0.29137162(Total Number of Households Above Poverty)
```{r include=FALSE}
coef(backward.fit,6)
```


# Shrinkage Methods

## Ridge Regression and Lasso

```{r}
set.seed(100)
k <- 10
cv.error.lasso <- rep(NA,k)
cv.error.ridge <- rep(NA,k)
cv.lam.lasso <- rep(NA,k)
cv.lam.ridge <- rep(NA,k)
set.seed(1)
folds<-sample(1:k,nrow(dat[,c(4:13)]),replace=TRUE)
for (i in 1:k){
  train <- dat[folds!=i,c(4:13)] #Set the training set
  test <- dat[folds==i,c(4:13)] #Set testing set
  x.train <- model.matrix(cum_case~., data=train)[,-1]
  y.train <- train$cum_case
  x.test <- model.matrix(cum_case~., data=test)[,-1]
  y.test <- test$cum_case
  
  #Lasso
  lasso.mod <- glmnet(x.train,y.train,alpha=1)
  cv.out <- cv.glmnet(x.train,y.train,alpha=1)
  bestlam <- cv.out$lambda.min
  cv.lam.lasso[i] <- bestlam
  lasso.pred <- predict(lasso.mod,s=bestlam ,newx=x.test)
  cv.error.lasso[i] <- mean((lasso.pred-y.test)^2)
  
  
  #Ridge
  ridge.mod <- glmnet(x.train,y.train,alpha=0)
  cv.out <- cv.glmnet(x.train,y.train,alpha=0)
  bestlam <- cv.out$lambda.min
  cv.lam.ridge[i] <- bestlam
  ridge.pred <- predict(ridge.mod,s=bestlam ,newx=x.test)
  cv.error.ridge[i] <- mean((ridge.pred-y.test)^2)

}
```

Recall that Ridge Regression performs better when the response is a function of many predictors of roughly equal size and Lasso performs better when the response is a function of a small subset of predictors.

10-fold Cross Validation was performed to find the optimal tuning parameter $\lambda$ for each method. 


## 

```{r}
data.frame("Fold" = rep(1:10,2),
           "Error" = c(cv.error.ridge,cv.error.lasso),
           "Method" = c(rep("Ridge Regresion",10),rep("Lasso",10))) %>%
  ggplot(aes(x=Fold,y=Error,color=Method))+
  geom_point()+
  geom_line() + 
  geom_point(aes(x=which.min(cv.error.ridge),
                 y=cv.error.ridge[which.min(cv.error.ridge)]),
             color="blue")+ 
  geom_point(aes(x=which.min(cv.error.lasso),
                 y=cv.error.lasso[which.min(cv.error.lasso)]),
             color="red") +
  scale_x_continuous(breaks = seq(1,10,2))
  
```


## Ridge With Lambda that Produces Smallest C.V. Error

```{r include=FALSE}
# Ridge CV Errors and Mean
cv.error.ridge[which.min(cv.error.ridge)]
mean(cv.error.ridge)
```

The smallest Cross Validation Error for Ridge Regression is 24901393 with and average of 58650401.

The model selected is:

Cumulative Cases = 2659.26735391 
-1.46826014(Population Density)
+0.01675252(Total Population)
-60.01296350(Median Age)
+0.03389108(Total White) 
+0.04046944(Total Black) 
+0.06843143 (Total Hispanic) 
+0.03043294 (Total Other) 
+0.29137162(Total Number of Households Above Poverty)
+0.27874317(Total Number of Households Below Poverty)

```{r include=FALSE}
x<-model.matrix(cum_case~., data=dat[,c(4:13)])[,-1]
y<-dat$cum_case

ridge.fit<-glmnet(x,y,alpha=0)
ridge.coef<-predict(ridge.fit,type="coefficients",
                    s=cv.lam.ridge[which.min(cv.error.ridge)])[1:9,]
ridge.coef[ridge.coef!=0]
```



## Lasso With Lambda that Produces Smallest C.V. Error

```{r include=FALSE}
# Lasso CV Errors and Mean
cv.error.lasso[which.min(cv.error.lasso)]
mean(cv.error.lasso)

```

The smallest Cross Validation Error for Lasso is 21711998 with and average of 64144413.

The model selected is:

Cumulative Cases = 286.95400883 
+0.04444873(Total Population)
+0.02943025(Total White) 
+0.01685243(Total Black) 
+0.08547591(Total Hispanic) 
+0.09854572(Total Number of Households Below Poverty)
```{r include=FALSE}
lasso.fit<-glmnet(x,y,alpha=1)
lasso.coef<-predict(lasso.fit,type="coefficients",
                    s=cv.lam.lasso[which.min(cv.error.lasso)])[1:9,]
lasso.coef[lasso.coef!=0]
```

# Tree Based Methods

## Regression Trees

```{r}
tree_dat <- dat %>% 
  subset(select = c(4:13)) %>% 
  rename("Cum. Cases" = 1) %>% 
  rename("Pop. Dens." = 2) %>% 
  rename("Total Pop." = 3) %>% 
  rename("Med. Age" = 4) %>% 
  rename("Total White" = 5) %>% 
  rename("Total Black" = 6) %>% 
  rename("Total Hisp." = 7) %>% 
  rename("Total Other" = 8) %>% 
  rename("House. Bel. Pov." = 9) %>% 
  rename("House. Ab. Pov." = 10) 
  
tree <- rpart(`Cum. Cases` ~.,
              data=tree_dat,control=rpart.control(cp=.0001))

#printcp(tree)

best <- tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]
pruned_tree <- prune(tree, cp=best)
prp(pruned_tree,
    faclen=0, #use full names for factor labels
    extra=1, #display number of obs. for each terminal node
    roundint=F, #don't round to integers in output
    branch = 1,
    varlen = 0,
    digits=5,) #display 5 decimal places in output
```


# Classification Analysis

**Question**: Does a given county have an instance rate higher than the median of the instance rates?

\ \ 

- Binary qualitative response: 
\begin{equation*}
\text{Instance rate} =
\begin{cases}
1,\quad \text{above the median of the instance rate}\\
0,\quad \text{otherwise}
\end{cases}
\end{equation*}

- Consider '$\color{blue}{\text{Poverty Rate}}$' instead of 'Total_Households_Above_Poverty' and 'Total_Households_Below_Poverty'

## Revisit Plot of Predictors 

```{r, fig.align='center', out.width = '90%'}
cls.dat <- dat1 %>% 
  mutate(case_rate = cum_case/Total_Pop,
         Total_Households = Total_Households_Above_Poverty + Total_Households_Below_Poverty) %>% 
  mutate(poverty_rate = Total_Households_Below_Poverty/Total_Households) %>% 
  subset(select = c(11, 2,5:8,4, 13 )) # 7 predictors: POP_DENSITY,Median_Age,Total_White,Total_Black,Total_Hispanic,Total_Other,poverty_rate  

## distribution plots
d <- melt(cls.dat[, -c(1)])
raw.plot <- ggplot(d, aes(x = value)) +
  facet_wrap(~variable, scales = "free_x") +
  geom_histogram()
raw.plot
```


## Data Transformation

```{r, fig.align='center', out.width = '90%'}
#data log transformation
log.dat <- cls.dat
log.dat[, -c(1, 7)] <- log10(cls.dat[, -c(1, 7)]+1) #log transformation except cum_case and Median Age

#distribution plots of predictors after log transformation and un-transformed Median_age
d.log <- melt(log.dat[, -c(1)])
log.dat.plot <- ggplot(d.log, aes(x = value)) +
  facet_wrap(~variable, scales = "free_x") +
  geom_histogram() 
log.dat.plot 
```

## Multiple Logistic Regression 

- Estimated coefficients of the logistic regression model with all predictors
```{r}
rate01 <- rep(0, length(log.dat$case_rate)) 
rate01[log.dat$case_rate > median(log.dat$case_rate)] <- 1 #high instance rate==1, low ins.rate==0
log.dat2 <- data.frame(rate01, log.dat[,-c(1)])

#glm fit with all predictors
glm.fit.all <- glm(rate01 ~ . , data = log.dat2, family = binomial)
kable(round(summary(glm.fit.all)$coefficient, 4)) 
```

- 'Total_Hispanic' is not statistically significant. 

## Misclassification Error Rate 

- Logistic Regression, LDA, QDA, KNN, and Decision Trees are performed.

- The 5-fold cross-validation is used for each classifier.

\ \ 

```{r}
#Randomly split the data set into 5 subsets with approximately equal size for 5-fold cross-validation
n <- nrow(log.dat2)
set.seed(1)
sample.id <- sample(rep(1:5, times = ceiling(n/5))[1:n])
training.list <- testing.list <- list()
for (i in 1:5) {
training.list[[i]] = log.dat2[sample.id != i, ]
testing.list[[i]] = log.dat2[sample.id == i, ]
}

glm.miscls <- 0
lda.miscls <- 0
qda.miscls <- 0
for (i in 1:5){
  #glm.log fit 
glm.log <- glm(rate01 ~ . , data = training.list[[i]], family = binomial)
glm.probs <- predict(glm.log, testing.list[[i]], type = "response")
glm.pred <- rep(0, length(glm.probs))
glm.pred[glm.probs > 0.5] <- 1
glm.miscls[i] <- mean(glm.pred != testing.list[[i]]$rate01)

#LDA
lda.fit <- lda(rate01 ~  . , data = training.list[[i]])
lda.pred <- predict(lda.fit, testing.list[[i]])
lda.class <- lda.pred$class
lda.miscls[i] <- mean(lda.class != testing.list[[i]]$rate01)

#QDA
qda.fit <- qda(rate01 ~  . , data = training.list[[i]])
qda.pred <- predict(qda.fit, testing.list[[i]])
qda.class <- qda.pred$class
qda.miscls[i]  <- mean(qda.class != testing.list[[i]]$rate01)

}
glm.err <- mean(glm.miscls)
lda.err <- mean(lda.miscls)
qda.err <- mean(qda.miscls)

cat(paste("    glm.misclassification =", round(glm.err, 4)))
cat(paste("    LDA.misclassification =", round(lda.err, 4)))
cat(paste("    QDA.misclassification =", round(qda.err, 4)))
```
\ \ 

- These all three classifiers have about 63% of overall correct prediction rates. 

## K-Nearest Neighbors (KNN)

```{r, fig.align='center', out.width = '90%'}
#KNN with different k
knn.err <- 0
for(k in 1:50){
knn.miscls <- 0
for (i in 1:5){
knn.pred <- knn(train = as.matrix(training.list[[i]][, -c(1)]),
    test = as.matrix(testing.list[[i]][, -c(1)]),
    cl = training.list[[i]][,1], k = k)
    knn.miscls[i] <- mean(knn.pred != testing.list[[i]][,1])
}
knn.err[k] <- mean(knn.miscls)
}
k.best <- which(knn.err == min(knn.err))

df <- data.frame(x=1:length(knn.err), y=knn.err)
df %>% ggplot(mapping= aes(x, y)) + 
  geom_point(shape=21) + 
  geom_line() +
  geom_point(aes(x=k.best, y=knn.err[k.best]) , color="red", size=2)+
  labs(x="K", y="Error Rate", title = "KNN Misclassification Error Rate")+
  geom_label(data=df %>% filter(y==min(y)), aes(label = paste0("(K=",k.best,", ", round(y,3),")")), hjust = -0.1, vjust= .75)
```

## Tree-based Classification 

```{r, fig.align='center', out.width = '90%'}
log.dat2$high <- ifelse(log.dat2$rate01 ==1, "Higher", "Lower")
tree.rate <- tree(as.factor(high) ~ . - rate01, log.dat2)

plot(tree.rate) 
text(tree.rate, pretty = 0)
```

## Tree Misclassification Error Rate 

```{r, fig.align='center', out.width = '90%'}
set.seed(1)
cv.tree.rate <- cv.tree(tree.rate, FUN = prune.misclass, K = 5)
df.tree <- data.frame(size = cv.tree.rate$size, err = cv.tree.rate$dev / 3219)
tree.err <- min(df.tree$err)
df.tree %>% ggplot(aes(x = size, y = err)) +
  geom_point(shape = 21) +
  geom_line() +
  labs(x = "number of terminal nodes", y = "Error Rate", 
       title = "Tree Misclassification Error Rate") +
  geom_point(aes(x = 4, y = min(err)), color = "red", size = 2) +
  geom_label(data = df.tree %>% filter(size == 4), 
             aes(label = sprintf("%.4f", err)), hjust = 0.5, vjust = -0.5)
```


## Comparison of Error Rates

```{r}
df.err <- data.frame(glm = round(glm.err, 4), LDA = round(lda.err,4), QDA=round(qda.err,4), 
                     KNN=round(knn.err[k.best],4), tree=round(tree.err, 4))
colnames(df.err) <- c("GLM", "LDA", "QDA", "KNN(K=40)", "Tree")
kable(df.err, caption = "Misclassification Error Rate")
```

- All of the classifiers produce similar error rates.


# Conclusion

## Results

- Selection methods without cross-validation each selected the same model at each parameter level, however their cross validation had Backward Stepwise selecting a model with more parameters.

- Ridge Regression performed better on average but Lasso had a slightly smaller minimum error during cross validation.

- In prediction of whether a given county is going to have a higher instance rate than the its median, all the classifiers employed in this project (Logistic regression, LDA, QDA, KNN, and Decision Tree) result in the similar accuracy of prediction; about 60% - 63% of overall correct prediction rates.