---
title: "IS5 in R: Paired Samples and Blocks (Chapter 18)"
author: "Nicholas Horton (nhorton@amherst.edu)"
date: "December 19, 2020"
output: 
  pdf_document:
    fig_height: 3
    fig_width: 5
  html_document:
    fig_height: 3
    fig_width: 5
  word_document:
    fig_height: 4
    fig_width: 6
---


```{r, include = FALSE}
# Don't delete this chunk if you are using the mosaic package
# This loads the mosaic and dplyr packages
library(mosaic)
library(readr)
library(janitor)
library(tidyr)
```

```{r, include = FALSE}
# knitr settings to control how R chunks work.
require(knitr)
opts_chunk$set(
  tidy = FALSE, # display code as typed
  size = "small" # slightly smaller font for code
)
```

## Introduction and background 

This document is intended to help describe how to undertake analyses introduced as examples in the Fifth Edition of *Intro Stats* (2018) by De Veaux, Velleman, and Bock.
This file as well as the associated R Markdown reproducible analysis source file used to create it can be found at http://nhorton.people.amherst.edu/is5.

This work leverages initiatives undertaken by Project MOSAIC (http://www.mosaic-web.org), an NSF-funded effort to improve the teaching of statistics, calculus, science and computing in the undergraduate curriculum.
In particular, we utilize the `mosaic` package, which was written to simplify the use of R for introductory statistics courses.
A short summary of the R needed to teach introductory statistics can be found in the mosaic package vignettes (https://cran.r-project.org/web/packages/mosaic).
A paper describing the mosaic approach was published in the *R Journal*: https://journal.r-project.org/archive/2017/RJ-2017-024.

## Chapter 18: Paired Samples and Blocks

```{r message = FALSE, warning = FALSE}
library(mosaic)
library(readr)
library(janitor)
library(tidyr) # for the gather() function
Dexterity <- read_csv("http://nhorton.people.amherst.edu/is5/data/Dexterity.csv") %>%
  janitor::clean_names()
```

By default, `read_csv()` prints the variable names.
These messages have been suppressed using the `message=FALSE` code chunk option to save space and improve readability. 
Here we use the `clean_names()` function from the `janitor` package to sanitize the names of the columns (which would otherwise contain special characters or whitespace). 

```{r}
Dexterity %>%
  select(age_months, dominant_1, non_dominant_2, gender) %>%
  head(n = 7)
```

### Section 18.1: Paired Data

```{r warning = FALSE}
# Figure 18.1
Dexterity %>%
  select(dominant_1, non_dominant_2) %>%
  gather(key = hand_type, value = speed, dominant_1, non_dominant_2) %>%
  gf_boxplot(speed ~ hand_type, fill = ~hand_type) %>%
  gf_labs(x = "", y = "Speed (cyl/sec)") +
  ylim(.2, .9) +
  guides(fill = FALSE)
```

#### Example 18.1: Identifying Paired Data

We begin by creating the data set on page 586.

```{r}
WorkWeek <- rbind(
  data.frame(name = "Jeff", fiveday = 2798, fourday = 2914),
  data.frame(name = "Betty", fiveday = 7724, fourday = 6112),
  data.frame(name = "Roger", fiveday = 7505, fourday = 6177),
  data.frame(name = "Tom", fiveday = 838, fourday = 1102),
  data.frame(name = "Aimee", fiveday = 4592, fourday = 3281),
  data.frame(name = "Greg", fiveday = 8107, fourday = 4997),
  data.frame(name = "Larry G.", fiveday = 1228, fourday = 1695),
  data.frame(name = "Tad", fiveday = 8718, fourday = 6606),
  data.frame(name = "Larry M.", fiveday = 1097, fourday = 1063),
  data.frame(name = "Leslie", fiveday = 8089, fourday = 6392),
  data.frame(name = "Lee", fiveday = 3807, fourday = 3362)
)
WorkWeek
```

Looking at pairwise differences in `Dexterity`.

```{r}
Dexterity %>%
  select(dominant_1, non_dominant_2) %>%
  mutate(difference = dominant_1 - non_dominant_2) %>%
  head(n = 18)
```

### Section 18.2: The Paired *t*-Test

#### Example 18.2: Checking Assumptions and Conditions

We can display the distribution of the differences (see page 588).

```{r}
WorkWeek <- WorkWeek %>%
  mutate(difference = fiveday - fourday)
gf_histogram(~difference, data = WorkWeek, binwidth = 2000, center = 1000) %>%
  gf_labs(x = "Difference (mi)", y = "# of Workers")
```

#### Example 18.3: Doing a Paired *t*-Test

It is straightforward to carry out the paired t-test.

```{r}
t.test(~difference, data = WorkWeek)
```

or

```{r}
nwork <- nrow(WorkWeek)
nwork # number of pairs
dwork <- mean(~difference, data = WorkWeek)
dwork # mean of differences
swork <- sd(~difference, data = WorkWeek)
swork # SD of differences

sework <- swork / (nwork^.5)
sework # SE of differences

twork <- (dwork - 0) / sework
twork # t stat

2 * xpt(twork, df = nwork - 1, lower.tail = FALSE)
```

The `xpt()` function finds the p-value and plots it on a graph to visualize it. Here, the visualization shows a one-sided test, but in the book, it is two sided.

### Section 18.3: Confidence Intervals for Matched Pairs

We begin by reading the data.

```{r message = FALSE}
Couples <- read_csv("http://nhorton.people.amherst.edu/is5/data/Couples.csv") %>%
  filter(wAge != "*") %>%
  mutate(wAge = as.numeric(wAge))
# table on page 592
Couples %>%
  select(wAge, hAge) %>%
  mutate(difference = hAge - wAge) %>%
  head(n = 7)
```

#### Step-By-Step Example: A Paired *t*-Interval

We replicate the example from page 593.

```{r}
DexData <- Dexterity %>%
  select(dominant_1, non_dominant_2) %>%
  mutate(difference = dominant_1 - non_dominant_2) %>%
  filter(dominant_1 < 1)
# For some reason, the book has removed one observation where dominant_1 = 1,
# but has kept the count of children at 93 instead of 92
gf_histogram(~difference, data = DexData, binwidth = .05, center = .025) %>%
  gf_labs(x = "Dominant-Non-dominant", y = "# Subjects")
```

Here we display the calculations using the `t.test()` function and then by hand.


```{r}
df_stats(~difference, data = DexData)
t.test(~difference, data = DexData)
```

```{r}
ndex <- nrow(DexData) + 1 # the book kept n at 93 for some reason
ndex # number of pairs (children)
ddex <- mean(~difference, data = DexData)
ddex # mean difference
sdex <- sd(~difference, data = DexData)
sdex # standard deviation of the differences
sedex <- sdex / (ndex^.5)
sedex # standard error of the differences
df <- ndex - 1
df # this matches the book, but it should be 91
tstats <- qt(p = c(.025, .975), df = df)
tstats
medex <- tstats * sedex
medex # margin of error of the differences

ddex + medex

# Or, if you don't want to go through all those calculations:
t.test(~difference, data = DexData, df = df)
```

#### Effect Size  

#### Example 18.4: Looking at Effect Size with a Paired *t* Confidence Interval

We can verify the calculations from the example.

```{r}
tstats <- qt(p = c(.025, .975), df = nwork - 1)
tstats

me <- tstats * sework
me # margin of error

dwork + me # confidence interval
```

### Section 18.4: Blocking

#### What's Independent?

#### Random Matters: A Bootstrapped Paired Data Confidence Interval and Hypothesis Test

Our usual approach to bootstrapping works here.

```{r}
set.seed(2345)
numsim <- 5000

# What does do() do?
mean(~difference, data = resample(DexData)) # One mean of a random resample
mean(~difference, data = resample(DexData)) # Another mean of a random resample

do(2) * mean(~difference, data = resample(DexData)) # Calculates two means

# We need numsim means
DexBoots <- do(numsim) * mean(~difference, data = resample(DexData))
```

For more information about `resample()`, refer to the [resample vignette in mosaic](https://cran.r-project.org/web/packages/mosaic/vignettes/Resampling.pdf).  

```{r warning = FALSE}
qdata(~mean, p = c(.025, .975), data = DexBoots)
DexBoots <- DexBoots %>%
  mutate(interval = ifelse(mean > .0245 & mean < .0783, "Within 95% Confidence",
    "Outside 95% Confidence"
  ))

# Figure 18.4, page 597
gf_histogram(~mean, fill = ~interval, data = DexBoots, binwidth = .002, center = .001) %>%
  gf_vline(xintercept = .0245) %>%
  gf_vline(xintercept = .0783) %>%
  gf_labs(x = "Differences of Bootstrapped Means", y = "# of Trials") +
  guides(fill = FALSE)

# Figure 18.5
gf_histogram(~ (mean - ddex), data = DexBoots, binwidth = .002, center = .001) %>%
  gf_vline(xintercept = ddex) %>%
  gf_vline(xintercept = -ddex) %>%
  gf_labs(x = "Differences of Bootstrapped Means - Actual Difference", y = "# of Trials")
```

```{r}
df_stats(~ (mean - ddex), data = DexBoots)
```

With `df_stats()`, we can see that our minimum is within the interval, but our maximum isn't.

```{r}
DexBoots %>%
  filter((mean - ddex) > ddex)
```

Like the book, there is one instance (out of 5,000), so we estimate the P-value as 1/5,000 (the book says 50,000, which is incorrect), or .0002.
