---
title: "Homework 06 Solutions"
output: 
  bookdown::pdf_document2:
    toc: false
    keep_tex: true
---

# Problem 1

```{r message = F, warning = F}
library(tidyverse)
library(kableExtra)
cdidata <- read.table("cdi.dat",header=T)

str(cdidata)
```

## Part a

### Make a table or tables showing appropriate summary statistics for each variable in the data set. Note that summary statistics for continuous variables will be different from the summary statistics for categorical variables.

### Indicate where (in which variables) there is missing data (NA’s), if any, how much there is (in each variable) and why it might be there.



### Make some appropriate descriptive EDA plots to illustrate any important features of the variables or possible important relationships among them.


<!---

The first two code sections produce summary information and EDA as encouraged in the assignment.

```{r results='asis', message = F, warning = F, eval = F, echo = F}

# initial approach (actually used `table1`)

summary_info <- function(col, name) {
  if (inherits(col, "character")){
    info <- table(col,useNA = "ifany")/length(col) 
    df <- as.data.frame(t(matrix(info)))
    names(df) <- names(info)
    print(kable_styling(kable(df, caption = name)))
  } else if (length(unique(col)) == 2) {
    info <- table(col,useNA = "ifany")/length(col) 
    df <- as.data.frame(t(matrix(info)))
    names(df) <- names(info)
    print(kable_styling(kable(df, caption = name)))
    } else {
  info <- summary(col)
  df <- as.data.frame(t(matrix(info)))
  names(df) <- names(info)
  print(kable_styling(kable(df, caption = name)))

  }
}

for (c_name in names(maskdata))
{
  summary_info(maskdata[,c_name], name = c_name)
  }

```

```{r message = F, warning = F, results='asis'}
library(knitr)
table_data <- maskdata %>% 
  mutate(study = factor(study),
         allo.wm = factor(study),
         mask.class = factor(mask.class),
         virus = factor(virus),
         breath.num = factor(breath.num),
         female = factor(female, 
                         labels = c("male", "female")),
         `log10(nasal)` = log10(nasal),
         `log10(throat)` = log10(throat),
         `log10(aerosols)` = log10(aerosols),
         `log10(droplets)` = log10(droplets)
         ) %>%
  select(-c(nasal, throat, aerosols, droplets))

options(knitr.kable.NA = '')
table_data %>%
  select(female, study, allo.wm, mask.class, virus, breath.num,
         cough.num, `log10(nasal)`, `log10(throat)`, `log10(aerosols)`,
         `log10(droplets)`) %>% 
  summary() %>% 
  knitr::kable(format = "latex",
               caption = "Summary of variables in mask experiment.") %>%
  kableExtra::row_spec(row = 0, bold = T) %>%
  kableExtra::kable_styling(latex_options =c("scale_down")) 
```


We observe `NA` values in each of the viral load variables, the homework and project document mention that all these variables (recorded in copies of the virus per ml) have a value of `2` if the viral load was measured as 0, and that the `NA`s are associated with measurements not being taken or being lost. Table \@ref(tab:tab-na) shows that most of these variables see very few `NA` values recorded. There also doesn't seem to be much relationship between when the `NA`s occur in these variables.

```{r tab-na}
na_table <- maskdata %>% 
  dplyr::select(nasal, throat, aerosols, droplets) %>%
  sapply(function(col) table(is.na(col))) 

any_na <- sum(apply(maskdata, 1, function(row) any(is.na(row))))
na_table <- cbind(na_table, `rows with any NA` = c(nrow(maskdata) - any_na, any_na))
rownames(na_table) <- c("not NA", "is NA")

na_table %>%
  knitr::kable(caption = "Number of NAs in each viral load column.") %>%
  kableExtra::kable_styling()
```


```{r fig.height = 12, fig.width = 9}
par(mfrow = c(4,3))
for (c_name in colnames(maskdata)){
  if(length(unique(maskdata[,c_name])) <= 3){
     c_vals <- names(table(maskdata[,c_name]))
     c_counts <- table(maskdata[,c_name])
     barplot(names.arg = c_vals,height= c_counts, main = c_name)
  }else {
    hist(maskdata[,c_name], main = c_name, xlab = c_name)
  }
}
```


It's highly encouraged to examine the data more deeply and understand how the variables relate. In the following EDA, we examine the interaction between variables, this helps us prepare for answering the questions of interesting in the later subsections.

The relationship between `allo.wm` and `mask.class`, as presented in Table \@ref(tab:tab-all-vs-mask), is very interesting, and digging into the data a bit more, it appears that even though individuals were allocated to a "wear mask" or "don't wear mask" class, but also had the opportunity to breath both with a mask and without a mask (in 2 seperate 30 minute breathing sessions). It's hard to rectify how an individual with `allo.wm = 1` (which means they were assigned to wear the mask), is listed as `mask.class = "ctrl"`. It's also unclear what the viral load `aerosol` and `droplets` data actually contains for individuals with `mask.class = "both"`. As such, we're going to use the `mask.class` feature to split of the groups (and ignore individuals that gave samples for both types as we are unclear what the data associated with that individual means). We also explored potentials for individuals with `mask.class = "both"` to have 2 rows (for both options), but that is not the case.

```{r tab-all-vs-mask, message = F}
# table(maskdata$allo.wm, maskdata$mask.class) 
maskdata %>% 
  group_by(allo.wm, mask.class) %>% 
  summarize(total = n()) %>%
  pivot_wider(names_from = mask.class, values_from = total, values_fill = 0) %>%
  knitr::kable(caption = "Split of studies (year) and types of virus") %>%
  kable_styling() %>% 
  column_spec(1,bold = T, border_right = T)
```

Table \@ref(tab:tab-study-vs-virus) helps us see that that the study year and the type of virus isn't independent, we that only 1 individual in our data set in the `flu` study year had coronavirus. 

```{r tab-study-vs-virus, message = F}
# table(maskdata$study, maskdata$virus) 
maskdata %>% 
  filter(mask.class != "both") %>% 
  group_by(study, virus) %>% 
  summarize(total = n()) %>%
  pivot_wider(names_from = virus, values_from = total, values_fill = 0) %>%
  knitr::kable(caption = "Allocated to wear a mask vs mask class recorded") %>%
  kable_styling() %>% 
  column_spec(1,bold = T, border_right = T)
```

Finally, we observed that the variable `breath.num` equals 2 if `mask.class = "both"` and otherwise it equals 1. Given this observation, we won't be using this variable due to the above discussion of `mask.class` vs `allo.wm`.

## Part b

```{r}
aerosols_base_lm <- lm(log10(aerosols-1) ~ mask.class, 
                       data = maskdata[maskdata$mask.class != "both",])
droplets_base_lm <- lm(log10(droplets-1) ~ mask.class, 
                       data = maskdata[maskdata$mask.class != "both",])
```

### aerosol model

```{r}
summary(aerosols_base_lm)
```

The above model's coefficients suggest that wearing a mask reduces the percentage of viral load in aerosol released in 30 minutes of breathing by `r round(-1*coef(aerosols_base_lm)[2]*100,2)`%. This statistically significant result. For both of the models in this section, we are basically doing a 1 way ANOVA (with a single binary covariate), as such the diagnostic plots won't make that much sense. We decided to use $log_{10}(\text{viral load}-1)$ as the viral load variables are highly skewed, and the associated paper we are comparing against also examined these loads in the log scale^[Given that the minimum value (when no viral load was actually recorded) was 2, we did a $log(x-1)$ transform so those values would be transformed to $0$.].

### droplets model

```{r}
summary(droplets_base_lm)
```

The above model's coefficients suggest that wearing a mask reduces the percentage of viral load in droplets released in 30 minutes of breathing buy `r round(-1*coef(droplets_base_lm)[2]*100,2)`%. This statistically significant result. 


## Part c

In order to only make 2 models, but still be able to answer the question of the efficiency of masks conditional on virus type we build models of the type 

$$
\log_{10}(\text{viral load}) \sim \text{virus} + \text{mask.class}:\text{virus} - 1 \;.
$$

This allows us to examine the efficiency of the masks per virus type using the interaction term.

```{r}
aerosols_interact_lm <- lm(log10(aerosols - 1) ~ virus + mask.class:virus - 1, 
                       data = maskdata[maskdata$mask.class != "both",])

droplets_interact_lm <- lm(log10(droplets - 1) ~ virus + mask.class:virus - 1, 
                       data = maskdata[maskdata$mask.class != "both",])
```

### aerosol model

```{r}
summary(aerosols_interact_lm)
```

Interestingly, when we start focusing on the efficacy of masks conditional on the virus type, the evidence is a bit more murky. The above model's coefficients suggest that wearing a mask when you have a **cold** reduces the percentage of viral load in aerosol released in 30 minutes of breathing by `r round(-1*coef(aerosols_interact_lm)[4]*100,2)`% (this is **not** statistically significant result). Similarly when you have the **flu**, wearning a mask reduces the precentage of viral load in aerosol released in 30 minutes of breathing by `r round(-1*coef(aerosols_interact_lm)[5]*100,2)`% (this is also **not** statistically significant result). If you have a coronavirus, mask wearing is associated with a reduction of your viral load in aerosol released in 30 minutes by `r round(-1*coef(aerosols_interact_lm)[6]*100,2)`% (this **is** a statistically significant result).


```{r echo = F, eval = F}
# this is if we just did it with the interaction term.
aerosols_interact_lm_hidden <- lm(log10(aerosols) ~ mask.class:virus - 1, 
                       data = maskdata[maskdata$mask.class != "both",])

summ <- summary(aerosols_interact_lm_hidden)
coef_info_aerosols_interact <- summ$coefficients
#vcov(aerosols_interact_lm)

cold_mask_coef <- diff(coef_info_aerosols_interact[c(1,2), "Estimate"])
flu_mask_coef <- diff(coef_info_aerosols_interact[c(3,4), "Estimate"])
hcov_mask_coef <- diff(coef_info_aerosols_interact[c(5,6), "Estimate"])

cold_mask_se <- sqrt(c(1,-1,0,0,0,0) %*% 
                       vcov(aerosols_interact_lm_hidden) %*% 
                       c(1,-1,0,0,0,0))
flu_mask_se <- sqrt(c(0,0,1,-1,0,0) %*% 
                      vcov(aerosols_interact_lm_hidden) %*% 
                      c(0,0,1,-1,0,0))
hcov_mask_se <- sqrt(c(0,0,0,0,1,-1) %*% 
                       vcov(aerosols_interact_lm) %*% 
                       caerosols_interact_lm_hidden)

n_obs <- nrow(maskdata[maskdata$mask.class != "both",])

p_cold_mask_norm <- 2*(1-pnorm(abs(cold_mask_coef)/ cold_mask_se))
p_flu_mask_norm <- 2*(1-pnorm(abs(flu_mask_coef)/ flu_mask_se))
p_hcov_mask_norm <- 2*(1-pnorm(abs(hcov_mask_coef)/ hcov_mask_se))
```

If we interpret the question as *"Is the efficacy of masks different between the different virus types"*, then we'd be looking to see if $\beta_\text{mask:cold} = \beta_\text{mask:flu} = \beta_\text{mask:coronavirus}$. To do this type of test we'd actually want to look at a slightly different model, namely

$$
\log_{10}(\text{viral load}) \sim \text{virus} +\text{mask.class} +\text{mask.class}:\text{virus} - 1 \;.
$$

```{r}
aerosols_interact_lm2 <- lm(log10(aerosols - 1) ~ virus + mask.class + 
                              mask.class:virus - 1,
   data = maskdata[maskdata$mask.class != "both",])
```

We can then examine amount of variability explained by the interaction term to determine if the efficacy for masks is different across the virus types. In the following summary, we see that the associated p-value for this question is not significant (`r round(anova(aerosols_interact_lm2)[["Pr(>F)"]][3],3)`), suggesting that we cannot conclude that efficacy of masks is different across virus strands.

```{r}
anova(aerosols_interact_lm2)
```

### droplets model

```{r}
summary(droplets_interact_lm)
```

A similar translation of this summary table can be made. We leave this to the reader.

```{r echo = F, eval = F}
# this is if we just did it with the interaction term.
droplets_interact_lm_hidden <- lm(log10(droplets) ~ mask.class:virus - 1, 
                       data = maskdata[maskdata$mask.class != "both",])

summ <- summary(droplets_interact_lm_hidden)
coef_info_droplets_interact <- summ$coefficients
#vcov(droplets_interact_lm)


cold_mask_coef <- diff(coef_info_droplets_interact[c(1,2), "Estimate"])
flu_mask_coef <- diff(coef_info_droplets_interact[c(3,4), "Estimate"])
hcov_mask_coef <- diff(coef_info_droplets_interact[c(5,6), "Estimate"])

cold_mask_se <- sqrt(c(1,-1,0,0,0,0) %*% 
                       vcov(droplets_interact_lm_hidden) %*% 
                       c(1,-1,0,0,0,0))
flu_mask_se <- sqrt(c(0,0,1,-1,0,0) %*% 
                      vcov(droplets_interact_lm_hidden) %*% 
                      c(0,0,1,-1,0,0))
hcov_mask_se <- sqrt(c(0,0,0,0,1,-1) %*% 
                       vcov(droplets_interact_lm_hidden) %*% 
                       c(0,0,0,0,1,-1))

n_obs <- nrow(maskdata[maskdata$mask.class != "both",])

p_cold_mask_norm <- 2*(1-pnorm(abs(cold_mask_coef)/ cold_mask_se))
p_flu_mask_norm <- 2*(1-pnorm(abs(flu_mask_coef)/ flu_mask_se))
p_hcov_mask_norm <- 2*(1-pnorm(abs(hcov_mask_coef)/ hcov_mask_se))
```

```{r}
droplets_interact_lm2 <- lm(log10(droplets - 1) ~ virus + mask.class + 
                              mask.class:virus - 1,
   data = maskdata[maskdata$mask.class != "both",])
anova(droplets_interact_lm2)
```

## Part d

As mentioned in the EDA in part a, we won't include # of breaths as it is only 2 when with `mask.class == "both"` (else it is 1), and we're not looking at the data where `mask.class == "both"`. 

Our model focuses on inference (focusing on potential collaborators' interest in the best estimate of the impact of masks on reducing the viral load across these populations). We include `age` and `female` features to account for demographic covariates that might impact an individual's production viral load production (which interestingly seems to have significant coefficients). We include also include nasal and throat recorded viral load, as it is natural to assume that an individual's viral load produced by breathing is associated with the amount of viral load found in the nose and throat. And finally, we include `study` to account for any measurement / experiment differences across the studies.

```{r}
aerosols_interact_lm <- lm(log10(aerosols) ~ virus + mask.class:virus + age + 
                             female + log10(nasal) + log10(throat) + study - 1, 
                       data = maskdata[maskdata$mask.class != "both",])

summary(aerosols_interact_lm)
```
Under this model, it appears that we can reject hypothesis that masks don't reduce the amount of viral load for those with a cold or a flu (and the estimated impact is `r -round(coef(aerosols_interact_lm)[9]*100,2)` and `r -round(coef(aerosols_interact_lm)[11]*100,2)`% respectively - aka pretty high.

Our data has a lot of left censured variables (basically all the viral load recordings), and the diagnostic plots in Figures \@ref(fig:fig-left-censured-aerosol) and \@ref(fig:fig-diagostics-aerosol) still highlight such things. We also examined added variable plots to understanding if additional patterns between variables and residuals where present but didn't find any. VIF analysis (when accounting for subspace structure in the data for the categorical variables) the generalized VIF values are not too worrisome related to the inflation of standard errors for the $\beta$ values. Additionally, given the small size of the data set, and the nature of the variables there wasn't too much imagined could be done / variable that could be removed in the situation.  




```{r fig-left-censured-aerosol, fig.cap = "True vs Estimated log10(aerosol viral load)."}
plot(x = aerosols_interact_lm$fitted.values,
     y = (aerosols_interact_lm$residuals + aerosols_interact_lm$fitted.values),
     ylab = "log10(aerosol viral load)",
     xlab = "estimated log10(aerosol viral load)")
```

```{r fig-diagostics-aerosol, fig.cap = "Diagnostics for final aerosol viral load model."}
par(mfrow = c(2,2))
plot(aerosols_interact_lm)
```



```{r fig.height = 8, message = F, warning = F, eval = F, echo = F}
library(car)
avPlots(aerosols_interact_lm, layout = c(4,3))

vif(aerosols_interact_lm)
```

-->

# Problem 2

```{r message = F, warning = F}
beauty <- read_csv("../ProfEvaltnsBeautyPublic.csv") %>%
  dplyr::select(-profevaluation, -profnumber, -multipleclass, 
                -starts_with("class"))
```

```{r warning =F , message = F , echo = F}
beauty2 <- read_csv("../ProfEvaltnsBeautyPublic.csv") %>%
  dplyr::select( -profnumber, -multipleclass, 
                -starts_with("class"))
```


## Part a

```{r}
library(leaps)
all_subsets <- regsubsets(courseevaluation ~.,
                          data = beauty,
                          nvmax = 10, intercept = T,
                          method = "exhaustive")

summ <- summary(all_subsets)

best_model_index <- which.min(summ$bic)

dropped <- summ$outmat %>% apply(2, function(col) {mean(col == " ") == 1})
all_subset_info <- cbind(summ$outmat[,!dropped], BIC = summ$bic) 
rownames(all_subset_info) <- 1:nrow(all_subset_info)
all_subset_info %>% data.frame %>%
  rownames_to_column(var = "num variables") %>%
  knitr::kable() %>% 
  kableExtra::kable_styling(latex_options = "scale_down") %>%
  row_spec(best_model_index, background = "gray")
```

```{r}
var_names <- colnames(all_subset_info)[all_subset_info[best_model_index,] == "*"]

all_subsets_model <- lm(as.formula(paste("courseevaluation ~",
                                         paste0(var_names, collapse = " + "))),
                        data = beauty)

summary(all_subsets_model)
```

## Part b

We'll show forward regression.

```{r}
forward_step <- step(lm(courseevaluation ~1, data = beauty),
  scope  = as.formula(paste("courseevaluation ~",
                      paste0(names(beauty)[names(beauty) != "courseevaluation"], 
                             collapse = " + "))),
                          direction = "forward", trace = 0)

summary(forward_step)
```

## Part c

```{r message = F, warning = F}
library(glmnet)
cv_lasso_selection <- cv.glmnet(x = as.matrix(beauty[,names(beauty) != "courseevaluation"]),
          y = beauty$courseevaluation, alpha = 1)

```
```{r}
lasso <- glmnet(x = as.matrix(beauty[,names(beauty) != "courseevaluation"]),
                y = beauty$courseevaluation,
                alpha = 1, lambda = cv_lasso_selection$lambda.1se)


coef_lasso <- as.matrix(coef(lasso))
coef_lasso_df <- data.frame(columns = names(coef_lasso[coef_lasso[,1] != 0,]),
         beta = coef_lasso[coef_lasso[,1] != 0,]) 
rownames(coef_lasso_df) <- NULL

coef_lasso_df
```

## Part d

```{r echo = T}
#basically coding up a document term matrix function for some reason

dtm_coef <- function(names_list){
  all_words <- unique(unlist(names_list))
  matrix_info <- matrix(nrow = length(names_list), ncol = length(all_words))
  
  r_idx <- 1
  for (col_names in names_list){
    matrix_info[r_idx,] <- all_words %in% col_names
    r_idx <- r_idx + 1
  }
  
  colnames(matrix_info) <- all_words
  rownames(matrix_info) <- names(names_list)
  
  return(matrix_info)
}

variables_in_model <- list(
  "best subset" = names(coef(all_subsets_model)),
  "lasso" = coef_lasso_df$columns,
  "forward stepwise" = names(forward_step$coefficients),
  "homework 05, 2c" = c(names(beauty[,names(beauty) != "courseevaluation"]), 
                        "profevaluation", "(Intercept)")
)
```

```{r tab-feature, echo = T}
t(dtm_coef(variables_in_model)) %>%
  data.frame(check.names = F) %>% 
  rownames_to_column(var = "features") %>%
  mutate(`best subset` = ifelse(`best subset`,  
                                cell_spec(`best subset`, "latex", color = "black"),
                                cell_spec(`best subset`, "latex", color = "red")),
         `lasso` = ifelse(`lasso`,  
                                cell_spec(`lasso`, "latex", color = "black"),
                                cell_spec(`lasso`, "latex", color = "red")),
         `forward stepwise` = ifelse(`forward stepwise`,  
                                cell_spec(`forward stepwise`, "latex", color = "black"),
                                cell_spec(`forward stepwise`, "latex", color = "red")),
         `homework 05, 2c` = ifelse(`homework 05, 2c`,  
                                cell_spec(`homework 05, 2c`, "latex", color = "black"),
                                cell_spec(`homework 05, 2c`, "latex", color = "red"))) %>%
  knitr::kable(format = "latex", caption = "Features include in each model.",
               escape = F) %>%
  kableExtra::kable_styling(latex_options = c("scale_down"),)
```

The number of features (excluding the intercept) in our models increase from 6 in the best subset model, lasso with 10, forward stepwise regression selecting 14 and our full model with 30 features in homework 05. If we take out the lasso model these models are nested inside each other. **It's hard to evaluate prediction models as we left nothing to evaluate them on.** We observe that the lasso had actually higher training mse than the best subset model - but that's not to crazy because the penalty means that we're get more training error allowed. 

```{r echo = F}
full_model <- lm(courseevaluation~., data = beauty2)
training_mse <- rep(NA, 4)
m_idx <- 1
for (model in list(all_subsets_model, lasso, forward_step, full_model)){
  if (inherits(lasso, "glmnet")){
    newx <- as.matrix(beauty[,names(beauty) != "courseevaluation"])
    yhat <- predict(object = model, newx = newx)
  } else {
    yhat <- predict(object = model)
  }
  y <- beauty$courseevaluation
  
  training_mse[m_idx] <- mean((y-yhat)^2)
  m_idx <- m_idx + 1
}

```


