```{r setup, include=FALSE}
library(magrittr)
library(tidyverse)
library(vroom)
library(conflicted)

filter <- dplyr::filter

knitr::opts_chunk$set(echo = TRUE)
```

## 2021-07-23 Review Day

Review a couple of concepts have generally not been understood well or been confusing.

---

## Q-Q Plots

Q-Q Plots are only for find out whether your data is normally distributed. They're not really exploratory plots because (1) you can only plot a single variable at a time so you can't use Q-Q plots to look at the relationships between the variables in your data and (2) just because variables are normally distributed that doesn't mean they're necessarily alike

We'll use some song data from Spotify to look at this (because it has a bunch of numeric variables), specifically:

- **valence**: *positivity of the song*; A measure from 0.0 to 1.0 describing the musical positiveness conveyed by a track. Tracks with high valence sound more positive (e.g. happy, cheerful, euphoric), while tracks with low valence sound more negative (e.g. sad, depressed, angry).
- **tempo**: *speed of the song*; The overall estimated tempo of a track in beats per minute (BPM). In musical terminology, tempo is the speed or pace of a given piece and derives directly from the average beat duration.

```{r}
spotify_songs <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-01-21/spotify_songs.csv')
```

---

Different variables can be fairly normally distributed, but still look different from one another and have different properties

```{r, fig.width = 1, fig.height = 1}
ggplot(spotify_songs, aes(sample = valence)) +
  geom_qq() +
  geom_qq_line()

ggplot(spotify_songs, aes(sample = tempo)) +
  geom_qq() +
  geom_qq_line()
```

But when you look at their distributions, they're quite different

```{r, fig.width = 2, fig.height = 1}
spotify_songs %>%
  pivot_longer(valence:tempo, names_to = 'variable', values_to = 'value') %>%
ggplot(aes(x = value)) +
  geom_density(aes(color = variable)) +
  facet_wrap(~ variable, scales = 'free') +
  theme_bw()
```

Also, you can't look at them together in a Q-Q plot! When you plot these two variables that are both normally distributed, looks like they have no relationship to each other

```{r, fig.width = 1, fig.height = 1}
ggplot(spotify_songs, aes(x = valence, y = tempo)) +
  geom_point() +
  geom_smooth(method = 'lm', se = F) +
  theme_bw()

lm(tempo ~ valence, data = spotify_songs) %>% broom::glance()
```

---

<br>

---

## Model Prediction

If you're going to make prediction using your model you need to either leave some data out OR have multiple datasets. We'll use the colon dataset for this demo and try to predict time to event

  + `X` (numeric) unique id for each row
  + `id` (numeric) Unique id for each patient, although each patient has two records, one for reoccurrence and one for death
  + `study` (binary) 1 for all patients in dataset as they are in the study, 0 would indicate a control but there are none in this data set
  + `rx` (categorical) treatment given, either levamisole (Lev), levamisole plus 5-FU (Lev+5FU), or just observation (Obs)
  + `sex` (binary) 1 for male, 0 for female
  + `age` (numeric) age in years of patient
  + `obstruct` (binary) indicates whether the colon was obstructed by the tumor, 0 for no, 1 for yes
  + `perfor` (binary) indicates whether the colon was perforated by the tumor, 0 for no, 1 for yes
  + `adhere` (binary) indicates whether the tumor had adhered to nearby organs, 0 for no, 1 for yes
  + `nodes` (numeric) the number of lymph nodes with detectable cancer
  + `status` (binary) censor status, patients who died without a reocurrence of colon cancer were censored to avoid biasing time to reoccurrence analyses, 0 for no censor, 1 for censor
  + `differ` (categorical) differentiation of tumor, 1 for well, 2 for good, 3 for poor
  + `extent` (categorical) extent of local spread of tumor, 1 for submucosa, 2 for muscle, 3 for serosa, 4 for contiguous structures
  + `surg` (binary) time from tumor surgery to enrollment in study, 0 for short, 1 for long
  + `node4` (binary) more than four positive lymph nodes, 0 for no, 1 for yes
  + `time` (numeric) time in days to either reoccurrence or death aka time to event
  + `etype` (binary) 1 for reoccurrence, 2 for death

### Read in Data

```{r}
vroom('R/2021-07-16_exploratory_data_analysis/colon.csv') %>%
  select(-`...1`) %>%
  na.omit() %>%
  filter(status == 0) -> colon
```

### Predict Time to Event

First, leave out some data to test your model predictions with later

```{r}
### leave some data out
set.seed(42)
colon %>%
  sample_frac(0.2) %>%
  mutate(set = 'test') -> colon_test

### now remove the hold out data from the main table
colon %>%
  anti_join(colon_test) %>%
  mutate(set = 'train') -> colon_train
```

Run the model on the training data

```{r}
lm(time ~ age + obstruct + perfor + adhere + nodes + 
     status + differ + extent, data = colon_train) -> model
```

Check accuracy

```{r}
### can run prediction on original data
colon_train %>%
  mutate(prediction = predict(model, newdata = .)) %>%
  mutate(diff = prediction - time) %>%
  summarize(average_error = mean(diff),
            mean_absolute_error = mean(abs(diff))) 

### need to have held out data as well because prediction always does better
### can run prediction on original data
colon_test %>%
  mutate(prediction = predict(model, newdata = .)) %>%
  mutate(diff = prediction - time) %>%
  summarize(average_error = mean(diff),
            mean_absolute_error = mean(abs(diff))) 
```

And compare accuracy between training and testing data visually

```{r, fig.width = 1.5, fig.height = 1}
rbind(colon_train, colon_test) %>%
  mutate(prediction = predict(model, newdata = .)) %>%
ggplot(aes(x = time, y = prediction, color = set)) +
  geom_abline(slope = 1, intercept = 0, linetype = 'dashed', color = 'gray60') +
  geom_point(alpha = 0.7) +
  coord_cartesian(xlim = c(0, 3500), ylim = c(0, 3500)) +
  labs(x = 'Actual Time to Event', y = 'Predicted Time to Event', color = '') +
  theme_bw()
```

---

<br>

---

## `readr`/`tidyr` Assignment Review: Combining Skills: Answer Questions About Data

Read in the `measles` data set at the link and answer the following questions.

- **state**: School's state	
county	School's county	y
district	School's district	y
name	School name	
type	Whether a school is public, private, charter	y
enroll	Enrollment: Depending on the state, enrollment is for kindergarten only or may extend to include other grades.
mmr	School's Measles, Mumps, and Rubella (MMR) vaccination rate	y
overall	School's overall vaccination rate	y
xmed	Percentage of students exempted from vaccination for medical reasons	y
xper	Percentage of students exempted from vaccination for personal reasons	y
xrel	Percentage of students exempted from vaccination for religious reasons	y
lat	School latitude	(only in individual state files)
lng	School longitude

```{r}
### read the data in here
read_csv("https://raw.githubusercontent.com/WSJ/measles-data/master/all-measles-rates.csv") %>%
  select(-index) -> measles
```

---

1. Which state has the lowest average percentage of MMR vaccinations?

```{r}
measles %>%
  group_by(state) %>%
  summarize(avg_mmr = mean(mmr)) %>%
  ungroup() %>% 
  arrange(avg_mmr)

### try to find mistake
measles %>%
ggplot(aes(x = mmr)) +
  geom_density()

measles %>%
  distinct(mmr) %>%
  arrange(mmr)

measles %>% filter(mmr == -1)

### without
measles %>%
  filter(mmr >= 0) %>%
  group_by(state) %>%
  summarize(avg_mmr = mean(mmr)) %>%
  ungroup() %>% 
  arrange(avg_mmr)
```

```{r}

### using numbers
measles %>% 
  filter(mmr >= 0) %>%
  group_by(state) %>%
  summarize(avg_mmr = mean(mmr)) %>%
  ungroup() %>%
  arrange(avg_mmr)

### visually
measles %>% 
  filter(mmr >= 0) %>%
  group_by(state) %>%
  summarize(avg_mmr = mean(mmr)) %>%
  ungroup() %>%
  arrange(desc(avg_mmr)) %>%
  mutate(state = factor(state, state)) %>%

ggplot(aes(x = state, y = avg_mmr)) +
  geom_col(aes(fill = avg_mmr), color = 'gray20') +
  scale_fill_gradient(high = 'dodgerblue4', low = 'firebrick3') +
  labs(x = 'US State', 
       y = 'Average Measles, Mumps, and Rubella (MMR)\nVaccination Rate') +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30))
```

2. What's the average percentage of students exempted from vaccination for each of California's counties?

```{r}
measles %>%
  distinct(overall) %>%
  arrange(overall)

### overall + exemption columns don't add up to 100, must use exemption columns to answer the question
measles %>%
  filter(overall >= 0)

### check exemption columns
# no religious exemptions because data is bad
measles %>%
  distinct(xrel) %>%
  arrange(xrel)
# only one problematic value xper == 169, remove before calculating
measles %>%
  distinct(xper) %>%
  arrange(xper)
measles %$% range(na.omit(xper))
measles %>% filter(xper >= 100)
#  medical looks good except for NAs
measles %>%
  distinct(xmed) %>%
  arrange(xmed)
measles %$% range(na.omit(xmed))

measles %>%
  filter(state == 'California') %>% 
  mutate(exempt = case_when(is.na(xmed) & !is.na(xper) ~ xper,
                            !is.na(xmed) & is.na(xper) ~ xmed,
                            TRUE ~ xmed + xper)) %>%
  select(year, name, county, exempt) %>%
  na.omit() %>%
  group_by(county) %>%
  summarize(avg_exempt = mean(exempt)) %>%
  ungroup() %>%
  arrange(desc(avg_exempt))
```

3. What's the year to year change in average vaccination rates?

```{r}
# using numbers
measles %>% 
  mutate(year = ifelse(year == '2017', '2017-18', year)) %>%
  filter(year != 'null', overall >= 0) %>% 
  group_by(year) %>%
  summarize(avg_vac = mean(overall)) %>%
  ungroup()

# visually
measles %>%
  mutate(year = ifelse(year == '2017', '2017-18', year)) %>%
  filter(year != 'null', overall >= 0) %>% 
  # group_by(year) %>%
  # summarize(avg_vac = mean(overall)) %>%

ggplot(aes(x = year, y = overall)) +
  geom_boxplot() +
  labs(x = 'School Year', y = 'Percent of Students Vaccinated') +
  theme_classic()
```

4. What's the difference in distribution between the three reasons for exemption from vaccination?

```{r}
measles %>%
  pivot_longer(xmed:xper, 
               names_to = 'exempt_reason', values_to = 'exempt_percent') %>%

ggplot(aes(x = exempt_percent, fill = exempt_reason)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(labels = c('medical', 'personal'),
                    values = c('deepskyblue3', 'forestgreen'),
                    name = 'Reason for Exemption') +
  labs(x = 'Percentage of Students Exempted from Vaccination') +
  coord_cartesian(xlim = c(0, 20)) +
  theme_classic()
```

5. What's the breakdown of average student enrollment by type of school per state?

```{r}
measles %>% #filter(enroll > 100)
  # mutate(enroll = ifelse(enroll >= 100, 100, enroll)) %>%
  select(state, type, enroll) %>% 
  na.omit() %>%
  group_by(state, type) %>%
  summarize(avg_enroll = mean(enroll)) %>%
  ungroup()
```



<br><br>



