---
title: "Logistic wrap-up + project time (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')

# Use a color blind friendly color palette throughout doc
library(tidyverse)
cb_palette <- c("black", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
scale_colour_discrete <- function(...) scale_colour_manual(values = cb_palette, ...)
scale_fill_discrete   <- function(...) scale_fill_manual(values = cb_palette, ...)
theme_set(theme_bw())
```


\
\


# Warm-up {-}


Fake news is...a problem.
Using data obtained from [https://www.kaggle.com/mdepak/fakenewsnet](https://www.kaggle.com/mdepak/fakenewsnet), and cleaned using code provided at [https://www.kaggle.com/kumudchauhan/fake-news-analysis-and-classification](https://www.kaggle.com/kumudchauhan/fake-news-analysis-and-classification), we'll try to predict if an article is fake news based on various features.
NOTE: These data are from real articles posted online.

```{r}
# Load packages and data
library(tidyverse)
fake_news <- read_csv("https://ajohns24.github.io/data/fake_news.csv") %>% 
  mutate(title_has_excl = title_excl > 0,
         is_fake = (type == "fake")) %>% 
  select(is_fake, title_has_excl, negative)
```




\
\
\
\




**EXAMPLE 1: Intuition**

We'll explore the relationship of whether or not an article `is_fake` (TRUE or FALSE) with `title_has_excl` (whether the title has an exclamation point) and `negative` (the percent of words that have negative sentiment):

```{r}
head(fake_news)
```

The model formula for this relationship is:

log(odds the article is fake) = $\beta_0$ + $\beta_1$ `negative` + $\beta_2$ `title_has_exclTRUE`

Intuitively:

- Do you think $\beta_1$, the `negative` coefficient, will be positive or negative?

- Do you think $\beta_2$, the `title_has_exclTRUE` coefficient, will be positive or negative?




\
\
\
\




**PAUSE: Build some plots and models of the relationships of interest**

```{r}
# Plot relationship of is_fake with title_has_excl
fake_news %>% 
  ggplot(aes(x = title_has_excl, fill = is_fake)) + 
  geom_bar(position = "fill")
```


```{r}
# Plot relationship of is_fake with negative sentiment
fake_news %>% 
  mutate(negative_brackets = cut(negative, 15)) %>% 
  group_by(negative_brackets) %>% 
  summarize(fake_rate = mean(is_fake)) %>% 
  ggplot(aes(x = negative_brackets, y = fake_rate)) + 
  geom_point() + 
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
```


```{r}
# Model the relationship of is_fake with negative sentiment & title_has_excl
news_model <- glm(is_fake ~ title_has_excl + negative, fake_news, family = "binomial")
coef(summary(news_model))
```

```{r}
# Plot the model on the probability scale
fake_news %>% 
  ggplot(aes(y = as.numeric(is_fake), x = negative, color = title_has_excl)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial")) +
  labs(y = "probability of being fake") + 
  lims(y = c(0,1))
```




\
\
\
\



**EXAMPLE 2: Interpret the coefficients of the estimated logistic regression model**

log(odds the article is fake) = -1.8791 + 2.5265 `title_has_exclTRUE` + 0.3829 `negative` 

```{r}
# Exponentiated coefficients
exp(-1.8791)
exp(2.5265)
exp(0.3829)
```


Intercept:


`title_has_exclTRUE` coefficient:


`negative` coefficient (on your own):








\
\
\
\



**EXAMPLE 3: How good is the model at distinguishing between real and fake news?**

We can use our model to calculate the probability that an article is fake based on its negative sentiment score and whether or not its title includes an exclamation mark.
Below are the probability calculations for the articles in our sample, split up by those that were & weren't actually fake.
Discuss.

```{r}
fake_news %>% 
  mutate(.fitted = predict(news_model, newdata = ., type = "response")) %>% 
  ggplot(aes(x = is_fake, y = .fitted)) +
  geom_boxplot()
```




\
\
\
\



**EXAMPLE 4: Evaluating the model (on your own)**

Let's explore using a probability threshold of 0.5 to make a binary prediction of an article's status:

- If the article's probability of being fake is greater than or equal to 50%, we'll predict that it's fake.
- Otherwise, we'll predict that it's real.

We can visualize this threshold on our predicted probability boxplots:

```{r}
fake_news %>% 
  mutate(.fitted = predict(news_model, newdata = ., type = "response")) %>% 
  ggplot(aes(x = is_fake, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.5, color = "red")
```

The resulting confusion matrix is below:

```{r}
# Rows = actual / observed to be fake (FALSE or TRUE)
# Columns = predicted to be fake (FALSE or TRUE)
library(janitor)
threshold <- 0.5
fake_news %>% 
  mutate(.fitted = predict(news_model, newdata = ., type = "response")) %>% 
  mutate(predictFake = .fitted >= threshold) %>% 
  tabyl(is_fake, predictFake) %>%
  adorn_totals(c("row", "col"))
```

```{r}
# Calculate the overall accuracy

```


```{r}
# Calculate the sensitivity (true positive rate)

```


```{r}
# Calculate the specificity (true negative rate)

```


**QUESTION**

In general, sensitivity & specificity are inversely related -- increasing one leads to a decrease in the other.
In this context, would you rather have high sensitivity or high specificity?
Why?



\
\
\
\



**EXAMPLE 5: Calculations (on your own)**

Suppose you come across a news article that has an exclamation point in the title, and 4% of its words have a negative sentiment.
Calculate the *probability* that this article is fake.
You should make this calculation by hand, but can check it using `predict()`:

```{r}
predict(news_model, newdata = data.frame(title_has_excl = TRUE, negative = 4), type = "response")
```





\
\
\
\



**EXAMPLE 6: Interaction (on your own)**

```{r}
data(penguins)
penguins_sub <- penguins %>% 
  filter(species %in% c("Adelie", "Gentoo"))

penguins_sub %>% 
  ggplot(aes(y = body_mass, x = flipper_len, color = species)) + 
  geom_smooth(method = "lm", se = FALSE)
```

a. In their relationship with `body_mass`, `flipper_len` and `species` *interact*. What does that *mean*?



b. The model formula can be expressed as:   

    E[`body_mass` | `flipper_len`, `species`] = $\beta_0$ + $\beta_1$ `flipper_len` + $\beta_2$ `speciesGentoo` + $\beta_3$ `flipper_len * speciesGentoo`
    
    Indicate whether the coefficients below are positive or negative and what your reasoning is:
    
    - $\beta_1$
    
    - $\beta_2$
    
    - $\beta_3$
    

\
\
\
\



**EXAMPLE 7: More interaction (on your own)**

```{r}
summary(lm(body_mass ~ flipper_len*species, penguins_sub))
```

a. Write out *separate sub-model formulas* of `body_mass` vs `flipper_len` for Adelies and Gentoos.

b. Interpret the `flipper_len` and `flipper_len:speciesGentoo` coefficients.

c. Predict the body mass of a Gentoo penguin with 220mm long flippers.







\
\
\
\





# Project intro {-}

- [Project info & description](https://docs.google.com/document/d/1p4uHC-kORKvikBnIPSk3gpStzp4pg1lLfOmcC6x0iPk/edit?usp=drive_link)

- [Project group preference form](https://moodle.macalester.edu/mod/assign/view.php?id=236820)
    - You cannot fill this out until after you've reviewed the data options!
    - Due Monday by 11:59pm. Late forms will not be accepted for credit. I need this info in order to form suggested groups for Tuesday's class.



\
\
\
\



# Choose your own adventure {-}

- Finish the above examples. Let me know if you want to sit and work on them together!!

- Finish PP 5 (due today)

- Review the data options for the project & fill out the Project group preference form.

- Start reviewing for Quiz 2 / putting together your note card.




## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your 'Activities' subfolder within your 'STAT155' folder and locate the HTML file. Open and inspect this in a web browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' tab in the lower left pane.
- Check the wrap-up announcements and solutions in the online course manual, at the bottom of the corresponding chapter.
- Work on homework!


