Statistical Inference and Data Exploration for Archaeologists

# Statistical Inference and Data Exploration for Archaeologists
## A Guide to using R for Beginners
### Ben Marwick, June 2020

---

.pull-right[
# Find me at...
.left[
[<svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 512 512"><path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"/></svg> @benmarwick](http://twitter.com/benmarwick)  
[<svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 496 512"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg> @benmarwick](http://github.com/benmarwick)  
[<svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 512 512"><path d="M326.612 185.391c59.747 59.809 58.927 155.698.36 214.59-.11.12-.24.25-.36.37l-67.2 67.2c-59.27 59.27-155.699 59.262-214.96 0-59.27-59.26-59.27-155.7 0-214.96l37.106-37.106c9.84-9.84 26.786-3.3 27.294 10.606.648 17.722 3.826 35.527 9.69 52.721 1.986 5.822.567 12.262-3.783 16.612l-13.087 13.087c-28.026 28.026-28.905 73.66-1.155 101.96 28.024 28.579 74.086 28.749 102.325.51l67.2-67.19c28.191-28.191 28.073-73.757 0-101.83-3.701-3.694-7.429-6.564-10.341-8.569a16.037 16.037 0 0 1-6.947-12.606c-.396-10.567 3.348-21.456 11.698-29.806l21.054-21.055c5.521-5.521 14.182-6.199 20.584-1.731a152.482 152.482 0 0 1 20.522 17.197zM467.547 44.449c-59.261-59.262-155.69-59.27-214.96 0l-67.2 67.2c-.12.12-.25.25-.36.37-58.566 58.892-59.387 154.781.36 214.59a152.454 152.454 0 0 0 20.521 17.196c6.402 4.468 15.064 3.789 20.584-1.731l21.054-21.055c8.35-8.35 12.094-19.239 11.698-29.806a16.037 16.037 0 0 0-6.947-12.606c-2.912-2.005-6.64-4.875-10.341-8.569-28.073-28.073-28.191-73.639 0-101.83l67.2-67.19c28.239-28.239 74.3-28.069 102.325.51 27.75 28.3 26.872 73.934-1.155 101.96l-13.087 13.087c-4.35 4.35-5.769 10.79-3.783 16.612 5.864 17.194 9.042 34.999 9.69 52.721.509 13.906 17.454 20.446 27.294 10.606l37.106-37.106c59.271-59.259 59.271-155.699.001-214.959z"/></svg> faculty.washington.edu/bmarwick](http://faculty.washington.edu/bmarwick)  
[<svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 512 512"><path d="M476 3.2L12.5 270.6c-18.1 10.4-15.8 35.6 2.2 43.2L121 358.4l287.3-253.2c5.5-4.9 13.3 2.6 8.6 8.3L176 407v80.5c0 23.6 28.5 32.9 42.5 15.8L282 426l124.6 52.2c14.2 6 30.4-2.9 33-18.2l72-432C515 7.8 493.3-6.8 476 3.2z"/></svg> bmarwick@uw.edu](mailto:bmarwick@uw.edu)
]
]

---

# Welcome

- Statistical inference is the process of making conclusions from a sample of data

- One kind of statistical inference is hypothesis testing, where we use a sample of data to determine the strength of evidence for or against our theory

- Data exploration is searching for patterns, groups and characteristics in our sample

- We are going to do statistical inference and data exploration using the tidyverse (to catch up, see [tidyverse for archaeologists](https://github.com/benmarwick/tidyverse-for-archaeology))

- When we write R code to analyse archaeological data, our analysis is **transparent and reproducible**. These are two vital characteristics of good scientific work.

- We should **share our code and data** with our research publications for others to study and learn

---

# What are we going to do today?

I will demonstrate and you will practice...

- Hypothesis testing: chi-square, t-test, ANOVA

- Exploration: Principal Components Analysis

- Exploration: k-means clustering

---

# Let's import some archaeological data 💽

---

---

```r
*library(rio)
```
]
 
.right-output-import-auto[

]

---
class: split-40
count: false

```r
library(rio)

# read in the data
*j_data <- import("https://bit.ly/j_data_xlsx", setclass = "tbl_df")
```
]
 
.right-output-import-auto[

]

---
class: split-40
count: false

```r
library(rio)

# read in the data
j_data <- import("https://bit.ly/j_data_xlsx", setclass = "tbl_df")

# take a look
*j_data
```
]
 
.right-output-import-auto[

```
# A tibble: 9,752 x 48
   Site  Square  Spit Group Artno Material Colour Weight Length Artclas Cortex Cortype Initiat Breaks Noseg
   <chr> <chr>  <dbl> <dbl> <dbl> <chr>    <chr>   <dbl>  <dbl> <chr>    <dbl> <chr>   <chr>   <chr>  <dbl>
 1 J     B          1     1     1 Chert    Brown   1.83   22.0  Flake        0 <NA>    Hertz   <NA>      NA
 2 J     B          1     1     2 Chert    Yellow  0.61   14.6  RetF         0 <NA>    Hertz   <NA>       3
 3 J     B          1     1     3 Chert    Dk Gr…  3.81   29.6  Core         0 <NA>    <NA>    <NA>      NA
 4 J     B          1     1     4 Chert    Grey    0.13    5.83 RetFrag      0 <NA>    <NA>    Both       1
 5 J     B          1     1     5 Chert    Grey    0.05    7.07 Flake        0 <NA>    Hertz   <NA>      NA
 6 J     B          1     1     6 Chert    Lt Br…  0.580  11.1  Core        NA <NA>    <NA>    <NA>      NA
 7 J     B          1     1     7 Chert    Dk Br…  0.03    5.72 Flake        0 <NA>    Hertz   <NA>      NA
 8 J     B          1     1     8 Chert    Cream   0.02    8.37 FFrag        0 <NA>    <NA>    <NA>      NA
 9 J     B          1     1     9 Chert    Grey    0.21   10.7  Flake-b      0 <NA>    Hertz   Long      NA
10 J     B          1     1    10 Chert    Grey    0.25   14.3  Flake        0 <NA>    Hertz   <NA>      NA
# … with 9,742 more rows, and 33 more variables: Platwid <dbl>, Plat <chr>, Focal <lgl>, Overhang <chr>,
#   NoDS <dbl>, NoPS <dbl>, Rtch <chr>, Plat_2 <chr>, Focal_2 <lgl>, Overhang_2 <chr>, NoDS_2 <dbl>,
#   NoPS_2 <dbl>, Rtch_2 <chr>, RetOri <chr>, Retype <chr>, Retloc <chr>, Retlen <dbl>, Retdep <dbl>,
#   Portion <chr>, Heat <dbl>, EdgeDam <dbl>, Weathering <dbl>, Corerot <dbl>, RetScNo <dbl>, Term <chr>,
#   Janus <lgl>, Recycled <lgl>, Redirect <lgl>, Width <dbl>, Thick <dbl>, Platthic <dbl>, area <dbl>,
#   elongation <dbl>
```
]

---

background-image: url(https://images.unsplash.com/photo-1515490480959-ce9152f7ea2b?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=2850&q=80)
background-size: cover

---

# Your turn: import the data

```r
library(rio)

# read in the data
j_data <- import("https://bit.ly/j_data_xlsx", setclass = "tbl_df")

# take a look
j_data
```

<div class="countdown" id="timer_5edc2aa6" style="right:0;bottom:0;" data-warnwhen="0">
<code class="countdown-time"><span class="countdown-digits minutes">02</span><span class="countdown-digits colon">:</span><span class="countdown-digits seconds">00</span></code>
</div>
---

# Let's do some hypothesis testing 🧪 We will use the `infer` package and the [Modern Dive](https://moderndive.com/) book by Chester Ismay and Albert Y. Kim

---

# This really is only one test

[![Only One Test](figures/downey.png)](http://allendowney.blogspot.com/2016/06/there-is-still-only-one-test.html)

---

# Here are the basic steps 👣

[![](figures/infer-diagram.png)](https://infer.netlify.app/)

---

Start with our data

---

Specify the variables in our data frame that we want to investigate

---

---

Declare the null hypothesis, the observed effect was simply due to random chance

---

---

Generate many datasets according to the null hypothesis

---

---

Calculate the distribution of the statistic from the generated datasets to create the null distribution.

---

Visualize the statistic relative to the null distribution of our generated data

---

# Let's do a chi-square test for the independence of raw material type and artefact type of our stone artefacts.

# Is the association between artefact type and raw material type random or independent?

---

class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
*library(tidyverse)
```
]
 
.right-output-test-chi-square-1-auto[

]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
*j_data
```
]
 
.right-output-test-chi-square-1-auto[

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
* filter(Material %in% c("Volcanic", "Quartzite"))
```
]
 
.right-output-test-chi-square-1-auto[

```
# A tibble: 475 x 48
   Site  Square  Spit Group Artno Material Colour Weight Length Artclas Cortex Cortype Initiat Breaks Noseg
   <chr> <chr>  <dbl> <dbl> <dbl> <chr>    <chr>   <dbl>  <dbl> <chr>    <dbl> <chr>   <chr>   <chr>  <dbl>
 1 J     B          1     1    15 Volcanic Grey     1.07   6.36 Flake-b     20 Angul   Hertz   Long      NA
 2 J     B          1     1    33 Volcanic Brown    0.23   8.5  Flake        0 <NA>    Bending <NA>      NA
 3 J     B          1     1    38 Volcanic Dk Gr…   0.06   4.24 Flake       60 Angul   Hertz   <NA>      NA
 4 J     B          1     1    59 Volcanic Dk Gr…   0.43  12.9  Shatt        0 <NA>    <NA>    <NA>      NA
 5 J     B          1     1    70 Volcanic Dk Br…   0.31   7.43 Hshat       NA <NA>    <NA>    <NA>      NA
 6 J     B          1     1    71 Volcanic Dk Br…   0.26   6.4  Hshat       NA <NA>    <NA>    <NA>      NA
 7 J     B          1     1    72 Volcanic Dk Br…   0.06   7.3  Hshat       NA <NA>    <NA>    <NA>      NA
 8 J     B          1     1    73 Volcanic Dk Br…   0.16   7.84 Hshat       NA <NA>    <NA>    <NA>      NA
 9 J     B          1     1    78 Volcanic Lt Br…   0.02   6.47 Hshat       NA <NA>    <NA>    <NA>      NA
10 J     B          2     1    88 Volcanic Grey     0.38  13.3  Flake       40 Angul   Hertz   <NA>      NA
# … with 465 more rows, and 33 more variables: Platwid <dbl>, Plat <chr>, Focal <lgl>, Overhang <chr>,
#   NoDS <dbl>, NoPS <dbl>, Rtch <chr>, Plat_2 <chr>, Focal_2 <lgl>, Overhang_2 <chr>, NoDS_2 <dbl>,
#   NoPS_2 <dbl>, Rtch_2 <chr>, RetOri <chr>, Retype <chr>, Retloc <chr>, Retlen <dbl>, Retdep <dbl>,
#   Portion <chr>, Heat <dbl>, EdgeDam <dbl>, Weathering <dbl>, Corerot <dbl>, RetScNo <dbl>, Term <chr>,
#   Janus <lgl>, Recycled <lgl>, Redirect <lgl>, Width <dbl>, Thick <dbl>, Platthic <dbl>, area <dbl>,
#   elongation <dbl>
```
]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
  filter(Material %in% c("Volcanic", "Quartzite")) %>%  
* filter(Artclas %in% c("Flake", "Hshat", "FFrag"))
```
]
 
.right-output-test-chi-square-1-auto[

```
# A tibble: 297 x 48
   Site  Square  Spit Group Artno Material Colour Weight Length Artclas Cortex Cortype Initiat Breaks Noseg
   <chr> <chr>  <dbl> <dbl> <dbl> <chr>    <chr>   <dbl>  <dbl> <chr>    <dbl> <chr>   <chr>   <chr>  <dbl>
 1 J     B          1     1    33 Volcanic Brown    0.23   8.5  Flake        0 <NA>    Bending <NA>      NA
 2 J     B          1     1    38 Volcanic Dk Gr…   0.06   4.24 Flake       60 Angul   Hertz   <NA>      NA
 3 J     B          1     1    70 Volcanic Dk Br…   0.31   7.43 Hshat       NA <NA>    <NA>    <NA>      NA
 4 J     B          1     1    71 Volcanic Dk Br…   0.26   6.4  Hshat       NA <NA>    <NA>    <NA>      NA
 5 J     B          1     1    72 Volcanic Dk Br…   0.06   7.3  Hshat       NA <NA>    <NA>    <NA>      NA
 6 J     B          1     1    73 Volcanic Dk Br…   0.16   7.84 Hshat       NA <NA>    <NA>    <NA>      NA
 7 J     B          1     1    78 Volcanic Lt Br…   0.02   6.47 Hshat       NA <NA>    <NA>    <NA>      NA
 8 J     B          2     1    88 Volcanic Grey     0.38  13.3  Flake       40 Angul   Hertz   <NA>      NA
 9 J     B          3     1   130 Volcanic Grey    17.8   38.4  Flake        0 <NA>    Hertz   <NA>      NA
10 J     B          3     1   131 Volcanic Grey    32.1   57.4  Flake       40 Round   <NA>    <NA>      NA
# … with 287 more rows, and 33 more variables: Platwid <dbl>, Plat <chr>, Focal <lgl>, Overhang <chr>,
#   NoDS <dbl>, NoPS <dbl>, Rtch <chr>, Plat_2 <chr>, Focal_2 <lgl>, Overhang_2 <chr>, NoDS_2 <dbl>,
#   NoPS_2 <dbl>, Rtch_2 <chr>, RetOri <chr>, Retype <chr>, Retloc <chr>, Retlen <dbl>, Retdep <dbl>,
#   Portion <chr>, Heat <dbl>, EdgeDam <dbl>, Weathering <dbl>, Corerot <dbl>, RetScNo <dbl>, Term <chr>,
#   Janus <lgl>, Recycled <lgl>, Redirect <lgl>, Width <dbl>, Thick <dbl>, Platthic <dbl>, area <dbl>,
#   elongation <dbl>
```
]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
  filter(Material %in% c("Volcanic", "Quartzite")) %>%  
  filter(Artclas %in% c("Flake", "Hshat", "FFrag"))

# plot the data
*ggplot(j_data1)
```
]
 
.right-output-test-chi-square-1-auto[

<img src="stat-inference-and-exploration-for-archaeologists_files/figure-html/test-chi-square-1_auto_5_output-1.png" width="432" />
]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
  filter(Material %in% c("Volcanic", "Quartzite")) %>%  
  filter(Artclas %in% c("Flake", "Hshat", "FFrag"))

# plot the data
ggplot(j_data1) +  
*      aes(Material,
*          fill = Artclas)
```
]
 
.right-output-test-chi-square-1-auto[

<img src="stat-inference-and-exploration-for-archaeologists_files/figure-html/test-chi-square-1_auto_6_output-1.png" width="432" />
]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
  filter(Material %in% c("Volcanic", "Quartzite")) %>%  
  filter(Artclas %in% c("Flake", "Hshat", "FFrag"))

# plot the data
ggplot(j_data1) +  
       aes(Material,  
           fill = Artclas) +  
* geom_bar()
```
]
 
.right-output-test-chi-square-1-auto[

<img src="stat-inference-and-exploration-for-archaeologists_files/figure-html/test-chi-square-1_auto_7_output-1.png" width="432" />
]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
  filter(Material %in% c("Volcanic", "Quartzite")) %>%  
  filter(Artclas %in% c("Flake", "Hshat", "FFrag"))

# plot the data
ggplot(j_data1) +  
       aes(Material,  
           fill = Artclas) +  
  geom_bar() +  
* theme_bw(base_size = 20)
```
]
 
.right-output-test-chi-square-1-auto[

<img src="stat-inference-and-exploration-for-archaeologists_files/figure-html/test-chi-square-1_auto_8_output-1.png" width="432" />
]

---
class: split-40
count: false
 
Prepare our data for a chi-square test: Is flake type independant from raw material type?
.left-code-test-chi-square-1-auto[

```r
library(tidyverse)

# filter the data to
# focus on our question
j_data %>%  
  filter(Material %in% c("Volcanic", "Quartzite")) %>%  
  filter(Artclas %in% c("Flake", "Hshat", "FFrag"))

# plot the data
ggplot(j_data1) +  
       aes(Material,  
           fill = Artclas) +  
  geom_bar() +  
  theme_bw(base_size = 20)  
```
]
 
.right-output-test-chi-square-1-auto[

<img src="stat-inference-and-exploration-for-archaeologists_files/figure-html/test-chi-square-1_auto_9_output-1.png" width="432" />
]

---

class: split-40
count: false
 
Compute the observed statistic, and use our real data to generate many fake datasets based on our null hypothesis
.left-code-test-chi-square-2-auto[

```r
*library(infer)
```
]
 
.right-output-test-chi-square-2-auto[

]