431-class03-all-code.Rmd

---
title: "Analyzing our Quick Survey"
author: "YOUR NAME SHOULD BE HERE"
date: "`r Sys.Date()`"
output:
  html_document: 
    theme: paper
    highlight: textmate
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
---

```{r setup, message=FALSE}
knitr::opts_chunk$set(comment=NA)
options(width = 70)
```

# Load Packages and Set Theme {-}

```{r library_load, message = FALSE}

library(janitor)
library(patchwork)
library(tidyverse)

theme_set(theme_bw())
```

# Load Data {-}

Our data come from the Quick 15-item Survey we did in Class 02 ([pdf in Class 02 README](https://github.com/THOMASELOVE/431-classes-2022/tree/main/class02)), which we've done (in various forms) since 2014. 

```{r load_data}
## if you want to load in a data set called namebeta.csv
## and then create a tibble from it called namealpha
## then uncomment the next line by removing the #

quicksur_raw <- 
  read_csv("quick_survey_2022.csv", show_col_types = FALSE) |>
  clean_names()
```

Our data come from the Quick 15-item Survey we did in Class 02 ([pdf in Class 02 README](https://github.com/THOMASELOVE/431-classes-2022/tree/main/class02)), which we've done (in various forms) since 2014. 

# Our Questions of Interest

We'll tackle several exploratory questions of interest...

1. What is the distribution of pulse rates among students in 431 since 2014?
2. Does the distribution of student heights change materially over time?
3. Is a Normal distribution is a good model for our data?
4. Do taller people appear to have paid less for their most recent haircut?
5. Do students have a more substantial tobacco history if they prefer to speak English or a language other than English?

# Initial Cleaning

## What does the `quicksur_raw` tibble include?

```{r}
quicksur_raw

glimpse(quicksur_raw)
```

## Counting Categories

```{r}
quicksur_raw |> count(glasses)

quicksur_raw |> count(glasses, english)
```

## Favorite Color in 2022?

```{r}
quicksur_raw |>
    filter(year == "2022") |>
    tabyl(favcolor) |>
    adorn_pct_formatting()
```

## Using `summary()` on Quantities

```{r}
quicksur_raw |> 
  select(love_htcm, haircut, height_in, lastsleep) |>
  summary()
```

- Numerical summaries (five quantiles, plus the mean) for:
  - your guess of my height (in cm), the price of your last haircut, your height (in inches), and the hours of sleep you got last night
- How many observations are available for these measures?

# Manage the data into `qsdat`

## Variables we'll look at closely today

To address our Questions of Interest, we need these seven variables in our analytic data frame (tibble.)

-   `student`: student identification (numerical code)
-   `year`: indicates year when survey was taken (August)
-   `english`: y = prefers to speak English, else n
-   `smoke`: 1 = never smoker, 2 = quit, 3 = current
-   `pulse`: pulse rate (beats per minute)
-   `height_in`: student's height (in inches)
-   `haircut`: price of student's last haircut (in \$)

## Select our variables

```{r}
qsdat <- quicksur_raw |>
    select(student, year, english, smoke, 
           pulse, height_in, haircut)

qsdat
```

## Initial Numeric Summaries

- Is everything the "type" of variable it should be? 
- Are we getting the summaries we want?

```{r}
summary(qsdat)
```

## What should we be seeing?

- Categorical variables should list the categories, with associated counts. 
  - To accomplish this, the variable needs to be represented in R with a `factor`, rather than as a `character` or `numeric` variable.
- Quantitative variables should show the minimum, median, mean, maximum, etc.

```{r}
names(qsdat)
```

## Change categorical variables to factors

We want the `year` and `smoke` information treated as categorical, rather than as quantitative, and the `english` information as a factor, too. Also, do we want to summarize the student ID codes?

- We use the `mutate()` function to help with this.

```{r}
qsdat <- qsdat |>
    mutate(year = as_factor(year),
           smoke = as_factor(smoke),
           english = as_factor(english),
           student = as.character(student))
```

- Note that it's `as_factor()` but `as.character()`. Sigh.

## Next step: Recheck the summaries and do range checks

- Do these summaries make sense?
- Are the minimum and maximum values appropriate?
- How much missingness are we to deal with?

```{r}
summary(qsdat)
```

## What to look for...

- Are we getting counts for all variables that are categorical?
    - Do the category levels make sense?
- Are we getting means and medians for all variables that are quantities?
    - Do the minimum and maximum values make sense for each of these quantities?
- Which variables have missing data, as indicated by `NA's`?

## The summary for `year` is an issue

- Just to fill in the gap left by the `summary()` result, how many students responded each year?

```{r}
qsdat |> tabyl(year) |> adorn_totals() |> adorn_pct_formatting()
```

# Question 1 

1. What is the distribution of pulse rates among students in 431 since 2014?

## Histogram, first try

-   What is the distribution of student `pulse` rates?

```{r}
ggplot(data = qsdat, aes(x = pulse)) +
    geom_histogram(bins = 30, fill = "royalblue", col = "seagreen1")
```

## Describing the Pulse Rates

How might we describe this distribution?

- What is the center?
- How much of a range around that center do we see? How spread out are the data?
- What is the shape of this distribution?
    - Is it symmetric, or is it skewed to the left or to the right? 

## Fundamental Numerical Summaries

```{r}
qsdat |> select(pulse) |> summary()
```

- How do the summary statistics help us describe the data?
- Do the values make sense to you?

```{r}
mosaic::favstats(~ pulse, data = qsdat)
```

## Histogram, version 2

```{r}
dat1 <- qsdat |>
  filter(complete.cases(pulse))

ggplot(data = dat1, aes(x = pulse)) +
    geom_histogram(fill = "seagreen", col = "white", bins = 20) +
    labs(title = "Pulse Rates of Dr. Love's students",
         subtitle = "2014 - 2022",
         y = "Number of Students",
         x = "Pulse Rate (beats per minute)")
```

- How did we deal with missing data?
- How did we add axis labels and titles to the plot?
- What is the distinction between `fill` and `col`?
- How many bins should we use?

# Question 2 

2. Does the distribution of student heights change materially over time?

## Yearly Five-Number Summaries

```{r}
qsdat |>
    filter(complete.cases(height_in)) |>
    group_by(year) |>
    summarize(n = n(), min = min(height_in), q25 = quantile(height_in, 0.25),
              median = median(height_in), q75 = quantile(height_in, 0.75),
              max = max(height_in))
```

- Does the distribution of heights change materially in 2014-2022?

## Five-Number Summary

- Key summaries based on percentiles / quantiles
    - minimum = 0th, maximum = 100th, median = 50th
    - quartiles (25th, 50th and 75th percentiles)
    - Range is maximum - minimum
    - IQR (inter-quartile range) is 75th - 25th percentile
- These summaries are generally more resistant to outliers than mean, standard deviation
- Form the elements of a boxplot (box-and-whisker plot)

## Comparison Boxplot of Heights by Year

```{r}
dat2 <- qsdat |>
    filter(complete.cases(height_in)) 

ggplot(data = dat2, aes(x = year, y = height_in)) +
    geom_boxplot() +
    labs(title = "Heights of Dr. Love's students, by year",
         subtitle = "2014 - 2022", x = "Year", y = "Height (in inches)")
```

## Thinking about the Boxplot

- Box covers the middle half of the data (25th and 75th percentiles), and the solid line indicates the median
- Whiskers extend from the quartiles to the most extreme values that are not judged by **Tukey's** "fences" method to be candidate outliers
    - Fences are drawn at 25th percentile - 1.5 IQR and 75th percentile + 1.5 IQR
- Are any values candidate outliers by this method? For which years?
- Was it important to change `year` to a factor earlier?

## Adding a Violin to the Boxplot

- When we'd like to better understand the shape of a distribution, we can amplify the boxplot.

```{r}
dat2 <- qsdat |>
    filter(complete.cases(height_in))

ggplot(data = dat2, aes(x = year, y = height_in)) +
    geom_violin() +
    geom_boxplot(aes(fill = year), width = 0.3) +
    guides(fill = "none") +
    scale_fill_viridis_d() +
    labs(title = "Heights of Dr. Love's students, by year",
         subtitle = "2014 - 2022", x = "Year", y = "Height (in inches)")
```

- How did we change the boxplot when we added the violin?
- What would happen if we added the boxplot first and the violin second?
- What does `guides(fill = "none")` do?
- What does `scale_fill_viridis_d()` do?

## Table of Means and Standard Deviations

```{r}
qsdat |>
    filter(complete.cases(height_in)) |>
    group_by(year) |>
    summarize(n = n(), mean = mean(height_in), sd = sd(height_in))
```

# Question 3 

3. Is a Normal distribution is a good model for our data?

- A Normal distribution is completely specified by its mean and standard deviation. The "bell shape" doesn't change.

## Summarizing Quantitative Data

If the data followed a Normal model, 

- we would be justified in using the sample **mean** to describe the center, and
- in using the sample **standard deviation** to describe the spread (variation.)

But it is often the case that these measures aren't robust enough, because the data show meaningful skew (asymmetry), or the data have lighter or heavier tails than a Normal model would predict.


## The Empirical Rule for Approximately Normal Distributions 

If the data followed a Normal distribution,

- approximately 68% of the data would be within 1 SD of the mean, 
- approximately 95% of the data would be within 2 SD of the mean, while 
- essentially all (99.7%) of the data would be within 3 SD of the mean.

## The Empirical Rule and 2022 Student Heights

In 2022, we had 54 students whose `height_in` was available, with mean 68.4 inches (173.7 cm) and standard deviation 3.7 inches (9.4 cm).

What do the histogram and boxplot (seen earlier) suggest about whether a Normal model with this mean and standard deviation would hold well for these 54 student heights?

## Histogram of 2022 Student Heights

```{r}
dat3 <- qsdat |>
  filter(complete.cases(height_in)) |>
  filter(year == "2022")

ggplot(data = dat3, aes(x = height_in)) +
    geom_histogram(fill = "salmon", col = "white", binwidth = 1) +
    labs(title = "Heights of Dr. Love's students",
         subtitle = "2022 (n = 54 students with height data)",
         y = "Number of Students", x = "Height (inches)")
```

- How did we use the two `filter()` statements?
- Why might I have changed from specifying `bins` to `binwidth` here?

## Checking the 1-SD Empirical Rule

- Of the 54 students in 2022 with heights, how many were within 1 SD of the mean?
  - Mean = 68.4, SD = 3.7.
  - 68.4 - 3.7 = 64.7 inches and 68.4 + 3.7 = 72.1 inches

```{r}
qsdat |> filter(complete.cases(height_in)) |>
    filter(year == "2022") |>
    count(height_in >= 64.7 & height_in <= 72.1)
```

## 2-SD Empirical Rule

- How many of the 54 `height_in` values gathered in 2022 were between 68.4 - 2(3.7) = 61.0 and 68.4 + 2(3.7) = 75.8 inches?


```{r}
qsdat |> filter(complete.cases(height_in)) |>
    filter(year == "2022") |>
    count(height_in >= 61.0 & height_in <= 75.8)
```

## 3-SD Empirical Rule

- How many of the 54 `height_in` values gathered in 2022 were between 68.4 - 3(3.7) = 57.3 and 68.4 + 3(3.7) = 79.5 inches?

```{r}
qsdat |> filter(complete.cases(height_in)) |>
    filter(year == "2022") |>
    count(height_in >= 57.3 & height_in <= 79.5)
```

## Empirical Rule Table for 2022 data

- $\bar{x}$ = sample mean, $s$ = sample SD
- For `height_in`: $n$ = 54 with data, $\bar{x} = 68.4, s = 3.7$
- For `pulse`: $n$ = 52 with data, $\bar{x} = 75.4, s = 11.2$

Range | "Normal" | `height_in` | `pulse`
:----: | :---: | :-----: | :-----:
$\bar{x} \pm s$ | ~68% | $\frac{39}{54}$ = 72.2% | $\frac{43}{52}$ = 82.7% 
$\bar{x} \pm 2\times s$ | ~95% | $\frac{52}{54}$ = 96.3% | $\frac{51}{52}$ = 98.1%
$\bar{x} \pm 3\times s$ | ~99.7% | $\frac{54}{54}$ = 100% | $\frac{52}{52}$ = 100%

## Boxplots of Height and of Pulse Rate

```{r}
dat4 <- qsdat |> filter(complete.cases(height_in), year == "2022")

p4 <- ggplot(data = dat4, aes(x = "height (inches)", y = height_in)) +
  geom_violin() + geom_boxplot(width = 0.3, fill = "tomato") +
  labs(title = "Boxplot of 2022 Student Heights", x = "")

dat5 <- qsdat |> filter(complete.cases(pulse), year == "2022")

p5 <- ggplot(data = dat5, aes(x = "pulse rate (beats/minute)", y = pulse)) +
  geom_violin() + geom_boxplot(width = 0.3, fill = "dodgerblue") +
  labs(title = "Boxplot of 2022 Pulse Rates", x = "")

p4 + p5 + 
  plot_annotation(title = "2022 Quick Survey Data")
```

## Normality and Mean/SD as summaries

If the data are approximately Normally distributed (like `height_in` and `pulse`) we can safely use the sample mean and standard deviation as summaries. If not "Normal", then ... 

- The median is a more robust summary of the center.
- For spread, we often use the 25th and 75th percentiles.

```{r}
dat3 <- qsdat |> filter(year == "2022")
mosaic::favstats(~ height_in, data = dat3)
mosaic::favstats(~ pulse, data = dat3)
```

## A new quantitative variable

Let's look at haircut prices, across all years.

```{r}
mosaic::favstats(~ haircut, data = qsdat)
```

Does it seem like the Normal model will be a good fit for these prices?

- Why or why not?
- What more information do you need to make a decision?

## 2022 Haircut Prices {.tabset}

### Unsorted

```{r}
qsdat |> filter(year == "2022") |> 
  select(haircut) |> 
  as.vector() ## just to print it here horizontally
```

### Sorted

```{r}
qsdat |> filter(year == "2022") |> 
  select(haircut) |> arrange(haircut) |> 
  as.vector() ## just to print it here horizontally
```

### Counts

```{r}
qsdat |> filter(year == "2022") |> 
  count(haircut) 
```

## 2022 Haircut Prices, tabulated

```{r}
qsdat |> filter(year == "2022") |> tabyl(haircut) |> adorn_pct_formatting()
```

## Normality of Haircut prices?

```{r}
dat6 <- qsdat |> filter(complete.cases(haircut))

p6a <- ggplot(data = dat6, aes(x = haircut)) +
  geom_histogram(binwidth = 5, fill = "purple", col = "white") +
  labs(x = "Haircut Price (in $)")

p6b <- ggplot(data = dat6, aes(x = haircut, y = "Price")) +
  geom_violin(fill = "plum") + geom_boxplot(width = 0.3) +
  labs(y = "", x = "Haircut Prices in $")

p6a + p6b +
  plot_annotation(
    title = "Histogram and Boxplot of Haircut Prices",
    subtitle = "2014-2022 Students of Dr. Love in 431")
```

- Do you think that the distribution of these prices follows a Normal model?

## Stem-and-Leaf of Haircut Prices

```{r}
stem(qsdat$haircut, scale = 2) # scale makes plot twice as long as default
```

## Empirical Rule Table for Haircut Prices

Let's look across all years, as well as just in 2022

```{r}
mosaic::favstats(~ haircut, data = qsdat)
mosaic::favstats(~ haircut, data = qsdat |> filter(year == "2022"))
```

Range | "Normal" | 2014-2022 | 2022
----: | :---: | :-----: | :-----:
$\bar{x} \pm s$ | ~68% | $\frac{438}{485}$ = 90.3% | $\frac{50}{54}$ = 92.6% 
$\bar{x} \pm 2\times s$ | ~95% | $\frac{465}{485}$ = 95.6% | $\frac{52}{54}$ = 96.3%
$\bar{x} \pm 3\times s$ | ~99.7% | $\frac{478}{485}$ = 98.6% | $\frac{52}{54}$ = 96.3%

## How did I calculate those fractions?

```{r}
# haircut price mean = 30.17 and sd = 31.41 across 2014-2022

qsdat |> count(haircut >= 30.17 - 31.41 & haircut <= 30.17 + 31.41)
qsdat |> count(haircut >= 30.17 - 2*31.41 & haircut <= 30.17 + 2*31.41)
qsdat |> count(haircut >= 30.17 - 3*31.41 & haircut <= 30.17 + 3*31.41)

# haircut price mean = 36.26 and sd = 41.82 in 2022 alone

qsdat |> filter(year == "2022") |> 
  count(haircut >= 36.26 - 41.82 & haircut <= 36.26 + 41.82)
qsdat |> filter(year == "2022") |> 
  count(haircut >= 36.26 - 2*41.82 & haircut <= 36.26 + 2*41.82)
qsdat |> filter(year == "2022") |> 
  count(haircut >= 36.26 - 3*41.82 & haircut <= 36.26 + 3*41.82)

```

# Question 4 

4. Do taller people appear to have paid less for their most recent haircut?

## Do tall people pay less for haircuts?

Why might we think that they do, before we see the data?

- Convert our student heights from inches to centimeters...

```{r}
qsdat <- qsdat |> mutate(height_cm = height_in * 2.54)

qsdat |> select(student, height_in, height_cm) |> head()
```

## Initial Numerical Summaries

```{r}
qsdat |> filter(complete.cases(haircut, height_cm)) |>
  summarize(n = n(), median(haircut), median(height_cm), median(height_in))
```

## A First Scatterplot

- We'll include the straight line from a linear model, in red.

```{r}
dat7 <- qsdat |> filter(complete.cases(height_cm, haircut)) 

ggplot(dat7, aes(x = height_cm, y = haircut)) +
    geom_point(alpha = 0.3) + 
    geom_smooth(method = "lm", col = "red",
                formula = y ~ x, se = TRUE) +
    labs(x = "Height (in cm)",
         y = "Price of last haircut (in $)",
         title = "Do taller people pay less for haircuts?")
```

## What is the (Pearson) correlation of height and haircut price?

```{r}
dat7 <- qsdat |> filter(complete.cases(height_cm, haircut)) 

dat7 |> 
    select(height_in, height_cm, haircut) |>
    cor() 
```

## What is the straight line regression model?

```{r}
dat7 <- qsdat |> filter(complete.cases(height_cm, haircut)) 

mod1 <- lm(haircut ~ height_cm, data = dat7)

mod1
```

## Summarizing our model `mod1`

```{r}
summary(mod1)
```

## Compare `lm` fit to `loess` smooth curve?

```{r}
dat7 <- qsdat |> filter(complete.cases(height_cm, haircut)) 

ggplot(dat7, aes(x = height_cm, y = haircut)) +
    geom_point(alpha = 0.5) + 
    geom_smooth(method = "lm", col = "red",
                formula = y ~ x, se = FALSE) +
    geom_smooth(method = "loess", col = "blue",
                formula = y ~ x, se = FALSE) +
    labs(x = "Height (in cm)",
         y = "Price of last haircut (in $)",
         title = "Do taller people pay less for haircuts?")
```

- Does a linear model appear to fit these data well?
- Do taller people pay less for their haircuts?

# Question 5 

5. Do students have a more substantial tobacco history if they prefer to speak English or a language other than English?

## Restrict ourselves to 2022 data

- Do students in the 2022 class have a more substantial history of tobacco use if they prefer to speak a language other than English?

```{r}
dat9 <- qsdat |> 
    filter(year == "2022") |>
    select(student, year, english, smoke)
```

```{r}
summary(dat9)
```

No missing data.

## Tabulating the categorical variables individually

```{r}
dat9 |> tabyl(english)

dat9 |> tabyl(smoke) |> adorn_pct_formatting()
```

## Cross-Classification (2 rows $\times$ 3 columns)

```{r}
dat9 |> tabyl(english, smoke)
```

## Recode the `smoke` levels to more meaningful names in `tobacco`

```{r}
dat9 <- dat9 |> 
    mutate(tobacco = fct_recode(smoke, 
            "Never" = "1", "Quit" = "2", "Current" = "3"))
```

### Check our work?

```{r}
dat9 |> count(smoke, tobacco)
```

- Everyone with `smoke` = 1 has `tobacco` as Never, etc.

## Restate the cross-tabulation 

Now we'll use this new variable, and this time, add row and column totals.

```{r}
dat9 |> tabyl(english, tobacco) |> 
    adorn_totals(where = c("row", "col"))
```

- What can we conclude about this association?

## How about in 2014-2022?

```{r}
dat8 <- qsdat |> 
  filter(complete.cases(english, smoke)) |>
  mutate(tobacco = fct_recode(smoke, 
            "Never" = "1", "Quit" = "2", "Current" = "3"))

dat8 |> 
  tabyl(english, tobacco) |> 
  adorn_totals(where = c("row", "col"))
```

- Now, what is your conclusion?

# Cleaning up the temporary objects

```{r}
rm(mod1,
   p4, p5, p6a, p6b,
   dat1, dat2, dat3, dat4, dat5, dat6, dat7, dat8, dat9
   )

## this just leaves
## qsdat and quicksur_raw in my Global Environment
```

# Session Information

Don't forget to close your file with the session information.

```{r}
sessionInfo()
```