-
Notifications
You must be signed in to change notification settings - Fork 14
Expand file tree
/
Copy path06-STAT_220_RandomizationPreLab-tidy.Rmd
More file actions
160 lines (98 loc) · 5.91 KB
/
Copy path06-STAT_220_RandomizationPreLab-tidy.Rmd
File metadata and controls
160 lines (98 loc) · 5.91 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
---
title: 'Randomization hypothesis testing pre-lab: tidy'
author: "Professor McNamara"
output:
html_document:
toc: true
toc_depth: 3
toc_float: true
---
We use three packages in this course: `Lock5Data`, `tidyverse` and `infer`. To load a package, you use the `library()` function, wrapped around the name of a package. I've put the code to load one package into the chunk below. Add the other two you need.
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, error = TRUE)
library(Lock5Data)
# put in the other two packages you need here
```
## Loading in data
As usual, we'll load the example data, `GSS_clean.csv`. It should be inside the `data` folder in your RStudio Cloud. We'll use the `read_csv()` function to read in the data.
```{r data-load}
GSS <- read_csv("data/GSS_clean.csv")
```
To make our lives easier, let's just drop NA values for the variables we're talking about today,
```{r}
GSS <- GSS |>
drop_na(
number_of_hours_worked_last_week,
self_emp_or_works_for_somebody
)
```
## Hypothesizing
Let's consider the relationship between `self_emp_or_works_for_somebody` and `number_of_hours_worked_last_week`. We might want to know if there is a significance difference in working hours between self employed people and and those employed by someone else.
We can compute the observed sample statistics,
```{r}
```
Those means **look** different to me, but we need to use statistics to be sure. So, let's perform a hypothesis test. Step one is hypothesizing.
We can write out hypotheses in RMarkdown using either words or notation.
We could write this in notation using LaTeX:
$$
H_0: \mu_{\text{self employed}} = \mu_{\text{someone else}} \\
H_A:\mu_{\text{self employed}} \neq \mu_{\text{someone else}}
$$
We could also have written this out in words:
Null hypothesis: there is no difference in average number of hours worked by people who are self employed and those who are employed by someone else.
Alternative hypothesis: there is a difference in average number of hours worked by people who are self employed and those who are employed by someone else.
## Randomization
The process of randomization is to create what the world would look like if the null hypothesis were true. In this case, that means we want to ensure there is no relationship between whether someone is self employed and the number of hours they worked. We can do this by breaking the relationship between the two variables. (See [Professor McNamara's video](https://youtu.be/12pQmmLMtTs) about making the connection between bootstrap and randomization distributions on YouTube for more.)
We can use the same sort of data pipeline we did for the bootstrap to help us with this. We need to specify our response and explanatory variables, set our hypothesis, and the generate permutation samples. "Permute" means to mix up the values of one variable, while while leaving the other one alone. Let's start by trying this one time.
```{r}
one_randomization <- GSS |>
specify(response = number_of_hours_worked_last_week,
explanatory = self_emp_or_works_for_somebody) |>
hypothesize(null = "independence") |>
generate(reps = 1, type = "permute")
```
How many observations does this new data set have? How many variables?
We can once again compute the grouped means,
```{r}
one_randomization |>
group_by(self_emp_or_works_for_somebody) |>
summarize(mean = mean(number_of_hours_worked_last_week))
```
Those means are slightly different from one another, but they look like they could just be the result of chance. (And in fact, they are! We ensured that, by permuting.)
## Setting a seed
Because these permuations are random, they look different every time we do them, just like bootstrap samples did. As before, this is good, because it will help us build up a null distribution, but it can be frustrating when you "knit" your document and the results change. To fix that, we are going to "set the seed" for the randomness. We use the `set.seed()` function, and then put our favorite number inside. I often use 42, or the day's date. It doesn't matter what you pick.
```{r}
set.seed(42)
```
Setting the seed means that the psuedo-random-number-generator in R will always start off at the same point, and so your "randomness" will be predictable.
## Randomization distribution
Just one randomization sample isn't that interesting. In order to make a distribution, we need to repeat this process many times, where many usually means at least 1,000.
Let's do that here. Once again, to avoid saving many datasets, we will add another verb to our pipeline, `calculate()`.
```{r}
randomization <- GSS |>
specify(response = number_of_hours_worked_last_week,
explanatory = self_emp_or_works_for_somebody) |>
hypothesize(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "diff in means",
order = c("Self-employed", "Someone else"))
```
How many observations does that dataset have? How many variables?
Let's look at the distribution of that statistic.
```{r}
```
How is this distribution shaped?
What is the mean of that distribution?
```{r}
```
Does that make sense?
## Evidence from randomization distributions
We could use the 95\% rule to eyeball a confidence interval that contained the middle 95\% of the distribution. My guess would be around (-2.5, 2.5). That would mean that any observed statistic that was larger than 2.5 or smaller than -2.5 would be pretty weird to see if the null hypothesis were true.
What was our observed statistic, again?
```{r}
```
Okay, so we think it is pretty weird, just by eyeballing. But, we'd rather have a specific p-value for how much of the distribution is as extreme or more extreme than that observed value. We can use the `get_p_value()` function to tell us how much of the distribution is to the left of the value,
```{r}
```
What do we conclude, at a $\alpha=0.05$ level?
How would this change if we were doing a one-sided test?