---
title: "The Bechdel Test"
author: "Solutions"
output:
  html_document:
    highlight: pygments
    theme: cerulean
    toc: true
    toc_float: true
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(eval = TRUE)
library(readr)
library(dplyr)
library(ggplot2)
```

## The Bechdel Test

A movie passes the Bechdel test if it satisfies 3 rules:

1. it has at least two women;
2. the women talk to each other; and
3. they talk to each other about something or someone other than a man.

The Bechdel test originated in this comic by Alison Bechdel (image source http://dykestowatchoutfor.com/the-rule): 

![](https://mhc-stat140-2017.github.io/labs/20170913_bechdel_movies/The-Rule-cleaned-up.jpg)

## Data Source

The data we're going to work with today have been gathered from a variety of sources by several people.  The Bechdel test ratings themselves are from www.bechdeltest.com, where the general public can rate movies according to whether they pass or fail the Bechdel test.  Some additional information about the movies comes from www.the-numbers.com.  These data were the basis of an article on the topic at http://fivethirtyeight.com/features/the-dollar-and-cents-case-against-hollywoods-exclusion-of-women/.  The data have since been added to the `fivethirtyeight` package for R.  I took those data, scraped some additional information about the movies like the MPAA rating, run time, and ratings from IMDB users from imdb.com.  Note that this is not a random sample of movies -- which movies made it into the data set was basically determined by which movies were rated by users of www.bechdeltest.com.  That means any findings from your analysis in this lab are only tentative.

## Initial Data Exploration

The following R chunk loads the data and sets the factor levels for categorical variables:

```{r load-data}
movies <- read_csv("https://mhc-stat140-2017.github.io/data/bechdel/bechdel.csv") %>%
  mutate(
    bechdel_test = factor(
      bechdel_test,
      ordered = TRUE,
      levels = c("nowomen", "notalk", "men", "dubious", "ok")),
    bechdel_test_binary = factor(
      bechdel_test_binary,
      ordered = TRUE,
      levels = c("FAIL", "PASS")),
    mpaa_rating = factor(
      mpaa_rating,
      ordered = TRUE,
      levels = c("UNRATED", "NOT RATED", "G", "PG", "TV-PG", "PG-13", "TV-14", "R", "NC-17"))
  )

dim(movies)

names(movies)

head(movies)

str(movies)
```

I think most of the variables are self-explanatory, but a couple require explanation.  The `bechdel_test` variable has five levels:

1. "nowomen" means there are not at least two women in the movie
2. "notalk" means there are at least two women in the movie, but they don't talk to each other;
3. "men" means there are at least two women in the movie, but they only talk to each other about men;
4. "dubious" means there was some disagreement among users of bechdeltest.com about whether or not the movie passed the test;
5. "ok" means that the movie passes the test.

The `bechdel_test_binary` variable has two levels:

1. "PASS" means that the movie passed the test (i.e., its value for `bechdel_test` is "ok")
2. "FAIL" means it did not pass the test (i.e., its value for `bechdel_test` is something other than "ok")

## Univariate Exploration

#### Make an appropriate plot of the imdb_rating variable

The `imdb_rating` variable is quantitative (look at the output from `str(movies)` above).  That means either a density plot or a histogram would be appropriate.

```{r}
ggplot() +
  geom_histogram(mapping = aes(x = imdb_rating), data = movies)
```

```{r}
ggplot() +
  geom_density(mapping = aes(x = imdb_rating), data = movies)
```

Note that in the code to make those two plots, everything is the same except for the type of geometry.

#### Calculate the mean, median, standard deviation and IQR of the imdb_rating variable

You should use the `summarize()` function so that the result is stored in a data frame.

```{r}
summarize(movies,
  mean_imdb_rating = mean(imdb_rating),
  median_imdb_rating = median(imdb_rating),
  sd_imdb_rating = sd(imdb_rating),
  iqr_imdb_rating = IQR(imdb_rating)
  )
```

#### Plot of bechdel_test

Make an appropriate plot of the `bechdel_test` variable.

The `bechdel_test` variable is categorical, so we make a bar plot.

```{r}
ggplot() +
  geom_bar(aes(x = bechdel_test), data = movies)
```

## Multivariate Exploration

#### Make an appropriate plot using the imdb_rating and bechdel_test variables

The `imdb_rating` variables is quantitative, and the `bechdel_test` variable is categorical.  For that combination, we could make box plots, or use a density curve with the different levels of `bechdel_test` associated with different colors.

```{r}
ggplot() +
  geom_boxplot(mapping = aes(x = bechdel_test, y = imdb_rating), data = movies)
```

```{r}
ggplot() +
  geom_density(mapping = aes(x = imdb_rating, color = bechdel_test),
    data = movies)
```

#### Calculate the mean and standard deviation of the imdb_rating variable, but do it separately for each level of the bechdel_test variable

You should use the `group_by` function and the `summarize()` function so that the result is stored in a data frame.

```{r}
movies %>%
  group_by(bechdel_test) %>%
  summarize(
    mean_imdb_rating = mean(imdb_rating),
    sd_imdb_rating = sd(imdb_rating))
```


#### How would you make the plot below?

![](https://mhc-stat140-2017.github.io/labs/20170913_bechdel_movies/bechdel_plot_to_make.png)

This is a bar plot based on two categorical variables: `bechdel_test_binary` on the horizontal axis and `bechdel_test` used for the fill color.  You can tell how those two variables were used in making the plot from the labels on the horizontal axis and on the legend for the colors.

```{r}
ggplot() +
  geom_bar(mapping = aes(x = bechdel_test_binary, fill = bechdel_test),
    data = movies)
```

#### Make an appropriate plot using the budget_2013 and imdb_rating variables

Both of these variables are quantitative, so we make a scatter plot (with geometry of type `point`)

```{r}
ggplot() +
  geom_point(mapping = aes(x = budget_2013, y = imdb_rating), data = movies)
```

## Additional exploration

Make a few more plots of your choosing.  Try different plot types and see what relationships you can find among the variables in the data set.  Add new R chunks as needed.

Here are a few more plots:

#### Is there a connection between MPAA rating ("G", "R", etc.) and Bechdel test results?

Both of these variables are categorical, so let's use a bar plot.

```{r}
ggplot() +
  geom_bar(mapping = aes(x = mpaa_rating, fill = bechdel_test), data = movies)
```

The relative proportion of different Bechdel test results seems fairly consistent across the different MPAA ratings categories.

Here's another variation on this plot, where each bar is scaled to 1, and the plot shows the proportion of movies with each result for the Bechdel test **within each MPAA ratings category**.

```{r}
ggplot() +
  geom_bar(mapping = aes(x = mpaa_rating, fill = bechdel_test),
    position = "fill",
    data = movies)
```

It does look like the distribution of Bechdel test results is similar across different MPAA rating categories, at least for the categories where we had a decent sample size (G, PG, PG-13, R).

#### Is there a connection between Bechdel test results and international gross earnings?

Let's use the variable with international gross inflation-adjusted to 2013 dollars.

```{r}
ggplot() +
  geom_boxplot(mapping = aes(x = bechdel_test, y = intgross_2013), data = movies)
```

```{r}
ggplot() +
  geom_density(mapping = aes(x = intgross_2013, color = bechdel_test), data = movies)
```

If you look reaaaaaallly closely at those plots, it seems like movies that pass the Bechdel test might tend to earn less than movies that fail the Bechdel test, on average.  But from the first box plot, the highest-earning movie in our data set did pass the Bechdel test.  What was that movie?

```{r}
arrange(movies, desc(intgross_2013))
```

The Exorcist!  Surprising to me...

Anyways, does this become clearer if we look at just the binary pass/fail measure?

```{r}
ggplot() +
  geom_boxplot(mapping = aes(x = bechdel_test_binary, y = intgross_2013), data = movies)
```

```{r}
ggplot() +
  geom_density(mapping = aes(x = intgross_2013, color = bechdel_test_binary), data = movies)
```

Yeah, it does look like movies that pass the test earn slightly less than movies that fail the test, in our data set.  The difference is not huge though.

#### Bechdel test results over time

```{r}
ggplot() +
  geom_histogram(mapping = aes(x = year, fill = bechdel_test_binary),
    data = movies)
```

```{r}
ggplot() +
  geom_histogram(mapping = aes(x = year, fill = bechdel_test_binary),
    position = "fill",
    data = movies)
```

Similar to the bar charts above, `position = "fill"` scales the values to show the proportion of movies in each category of `bechdel_test_binary` within each bin of years on the horizontal axis.

It does seem like more movies are passing the test now than used to, generally.  But we're still at less than 50% passing the test in recent years.
