Project Code.Rmd

---
title: "Project Code"
author: "Minsu Kim and Ariane Stark"
date: "5/12/2021"
output: pdf_document
---

```{r setup, include=FALSE}
library(readr)
library(readxl)
library(tidyverse)
library(rpart) 
library(rpart.plot)
library(leaps)
library(glmnet)
library(kableExtra)
library(reshape2)
library(ggplot2)
library(scales)
library(patchwork)
library(class)
library(MASS)
library(gridExtra)
library(grid)
library(glmnet)
library(kableExtra)
library(tree)
```

```{r message=FALSE, warning=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(Prop_Pov = Total_Households_Below_Poverty/(
    Total_Households_Below_Poverty + Total_Households_Above_Poverty)) %>%
  na.omit()  #removes Rio Arriba County NM which has NA in the pov. dem.
```


# Cross Validation of Best Subset, Forward and Backward Stepwise

```{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:10,14)]),replace=TRUE)
cv.errors.best_subset <- matrix(NA,k,7, dimnames=list(NULL, paste(1:7)))
cv.errors.forward <- matrix(NA,k,7, dimnames=list(NULL, paste(1:7)))
cv.errors.backward <- matrix(NA,k,7, dimnames=list(NULL, paste(1:7)))


for(j in 1:k){
  # best subset select
  best_subset.fit <- regsubsets(cum_case~.,data = dat[folds!=j,c(4:10,14)],
                                method = "exhaustive", nvmax = 7)
  # forward stepwise select
  forward.fit <- regsubsets(cum_case~.,data = dat[folds!=j,c(4:10,14)],
                                method = "forward", nvmax = 7)
  # backward stepwise select
  backward.fit <- regsubsets(cum_case~.,data = dat[folds!=j,c(4:10,14)],
                                method = "backward", nvmax = 7)
  
  for(i in 1:7){
    # best subset error
    pred.best <- 
      predict.regsubsets(best_subset.fit, dat[folds==j,c(4:10,14)], 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:10,14)], 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:10,14)], 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)
mean.cv.errors.best
plot(mean.cv.errors.best,type="b",
     xlab="Number of Variables", 
     ylab="Best Subset Mean C.V. Error")
points(which.min(mean.cv.errors.best),
       mean.cv.errors.best[which.min(mean.cv.errors.best)],
       col="red", cex=1.5,pch=20)
```


```{r}
#forward stepwise CV model selection
mean.cv.errors.forward <- apply(cv.errors.forward,2,mean)
mean.cv.errors.forward
plot(mean.cv.errors.forward,type="b",
     xlab="Number of Variables", 
     ylab="Forward Stepwise Mean C.V. Error")
points(which.min(mean.cv.errors.forward),
       mean.cv.errors.forward[which.min(mean.cv.errors.forward)],
       col="red", cex=1.5,pch=20)
```


```{r}
#backward stepwise CV model selection
mean.cv.errors.backward <- apply(cv.errors.backward,2,mean)
mean.cv.errors.backward
plot(mean.cv.errors.backward,type="b",
     xlab="Number of Variables", 
     ylab="Backward Stepwise Mean C.V. Error")
points(which.min(mean.cv.errors.backward),
       mean.cv.errors.backward[which.min(mean.cv.errors.backward)],
       col="red", cex=1.5,pch=20)
```




# Best Subset Overall
```{r}
best_subset.fit <- regsubsets(cum_case~.,data = dat[,c(4:10,14)],
                              method = "exhaustive", nvmax = 7)
best_subset.summary <- summary(best_subset.fit)
best_subset.summary.frame <- data_frame("Parameters" = seq(1:7),
                                        "R^2"=best_subset.summary$rsq,
                                      "AdjR^2"=best_subset.summary$adjr2,
                                      "CP"=best_subset.summary$cp,
                                      "BIC"=best_subset.summary$bic)
best_subset.summary.frame %>% kable()

```
Based on Cross Validation we pick the model with 2 parameters
```{r}
coef(best_subset.fit,2)
```


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

```
Based on Cross Validation we pick the model with 2 parameters
```{r}
coef(forward.fit,2)
```


# Backward Stepwise Overall
```{r}
backward.fit <- regsubsets(cum_case~.,data = dat[,c(4:10,14)],
                              method = "backward", nvmax = 7)
backward.summary <- summary(backward.fit)
backward.summary.frame <- data_frame("Parameters" = seq(1:7),
                                        "R^2"=backward.summary$rsq,
                                      "AdjR^2"=backward.summary$adjr2,
                                      "CP"=backward.summary$cp,
                                      "BIC"=backward.summary$bic)
backward.summary.frame %>% kable() 

```
Based on Cross Validation we pick the model with 6 parameters
```{r}
coef(backward.fit,6)
```


# Ridge 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:10,14)]),replace=TRUE)
for (i in 1:k){
  train <- dat[folds!=i,c(4:10,14)] #Set the training set
  test <- dat[folds==i,c(4:10,14)] #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)

}
```

```{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))
  
```

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


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

Ridge Regression performs better when the response is a function of many predictors of roughly equal size

# Ridge With Best CV Error Lambda
```{r}
x<-model.matrix(cum_case~., data=dat[,c(4:10,14)])[,-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:7,]
ridge.coef[ridge.coef!=0]
```


# Lasso With Best CV Error Lambda
```{r}
lasso.fit<-glmnet(x,y,alpha=1)
lasso.coef<-predict(lasso.fit,type="coefficients",
                    s=cv.lam.lasso[which.min(cv.error.lasso)])[1:7,]
lasso.coef[lasso.coef!=0]
```






## Regression Trees

```{r}
tree_dat <- dat %>% 
  subset(select = c(4:10,14))

tree <- rpart(cum_case ~.,
              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
```


# Minsu


## Descriptive data
```{r}
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 
dat <- dat[, -c(1:3)]

dim(dat) #obs: 3219 counties, 9 predictors
names(dat)

# data summary 
summary(dat)

## distribution plots
d <- melt(dat[, -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 

#scatter plots with simple linear regression
dat %>% 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() +
    labs(x = "variable", y = "cum_case", 
       title = "Simple linear regressions on scatter plots")
    
```


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

```{r}
cls.dat <- dat %>% 
  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()+
    labs(title = "Distribution of each variable")
raw.plot
```

### Data transformation
```{r}
#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() +
    labs(title = "Distribution of each variable after log-transformation")
log.dat.plot 
```

### create binary response: high/low instance rate 
```{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(summary(glm.fit.all)$coefficient) #Total_Hispanic is not statistically significant, but I will keep it. 
```


```{r}
log.dat2 <- log.dat2 %>% subset(select = -c(5))  #remove Total_Hispanic
```

```{r}
#glm fit with all predictors without Total_Hispanic
glm.fit.all <- glm(rate01 ~ . , data = log.dat2, family = binomial)
summary(glm.fit.all)$coefficient
```


```{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)))
```


```{r}
#KNN with different k
set.seed(1)
knn.err <- 0
for(k in 1:70){
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))
```


```{r}
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)
```


```{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))
colnames(df.err) <- c("GLM", "LDA", "QDA", "KNN(K=50)")
kable(df.err, caption = "Misclassification Error Rate")
```