|
| 1 | +--- |
| 2 | +title: "Tidy inference for paired data" |
| 3 | +description: "Conducting tests for paired independence on tidy data with infer." |
| 4 | +output: rmarkdown::html_vignette |
| 5 | +vignette: | |
| 6 | + %\VignetteIndexEntry{Tidy inference for paired data} |
| 7 | + %\VignetteEngine{knitr::rmarkdown} |
| 8 | + \usepackage[utf8]{inputenc} |
| 9 | +--- |
| 10 | + |
| 11 | +```{r settings, include=FALSE} |
| 12 | +knitr::opts_chunk$set(fig.width = 6, fig.height = 4.5) |
| 13 | +options(digits = 4) |
| 14 | +``` |
| 15 | + |
| 16 | +```{r load-packages, echo = FALSE, message = FALSE, warning = FALSE} |
| 17 | +library(ggplot2) |
| 18 | +library(dplyr) |
| 19 | +library(infer) |
| 20 | +``` |
| 21 | + |
| 22 | +### Introduction |
| 23 | + |
| 24 | +In this vignette, we'll walk through conducting a randomization-based paired test of independence with infer. |
| 25 | + |
| 26 | +Throughout this vignette, we'll make use of the `gss` dataset supplied by `infer`, which contains a sample of data from the General Social Survey. See `?gss` for more information on the variables included and their source. Note that this data (and our examples on it) are for demonstration purposes only, and will not necessarily provide accurate estimates unless weighted properly. For these examples, let's suppose that this dataset is a representative sample of a population we want to learn about: American adults. The data looks like this: |
| 27 | + |
| 28 | +```{r glimpse-gss-actual, warning = FALSE, message = FALSE} |
| 29 | +dplyr::glimpse(gss) |
| 30 | +``` |
| 31 | + |
| 32 | +Two sets of observations are paired if each observation in one column has a special correspondence or connection with exactly one observation in the other. For the purposes of this vignette, we'll simulate an additional data variable with a natural pairing: suppose that each of these survey respondents had provided the number of `hours` worked per week when surveyed 5 years prior, encoded as `hours_previous`. |
| 33 | + |
| 34 | +```{r} |
| 35 | +set.seed(1) |
| 36 | +
|
| 37 | +gss_paired <- gss %>% |
| 38 | + mutate( |
| 39 | + hours_previous = hours + 5 - rpois(nrow(.), 4.8), |
| 40 | + diff = hours - hours_previous |
| 41 | + ) |
| 42 | +
|
| 43 | +gss_paired %>% |
| 44 | + select(hours, hours_previous, diff) |
| 45 | +``` |
| 46 | + |
| 47 | +The number of `hours` worked per week by a particular respondent has a special correspondence with the number of hours worked 5 years prior `hours_previous` by that same respondent. We'd like to test the null hypothesis that the `"mean"` hours worked per week did not change between the sampled time and five years prior. |
| 48 | + |
| 49 | +To carry out inference on paired data with infer, we pre-compute the difference between paired values at the beginning of the analysis, and use those differences as our values of interest. |
| 50 | + |
| 51 | +Here, we pre-compute the difference between paired observations as `diff`. The distribution of `diff` in the observed data looks like this: |
| 52 | + |
| 53 | +```{r plot-diff, echo = FALSE} |
| 54 | +unique_diff <- unique(gss_paired$diff) |
| 55 | +gss_paired %>% |
| 56 | + ggplot2::ggplot() + |
| 57 | + ggplot2::aes(x = diff) + |
| 58 | + ggplot2::geom_histogram(bins = diff(range(unique_diff))) + |
| 59 | + ggplot2::labs(x = "diff: Difference in Number of Hours Worked", |
| 60 | + y = "Number of Responses") + |
| 61 | + ggplot2::scale_x_continuous(breaks = c(range(unique_diff), 0)) |
| 62 | +``` |
| 63 | + |
| 64 | +From the looks of the distribution, most respondents worked a similar number of hours worked per week as they had 5 hours prior, though it seems like there may be a slight decline of number of hours worked per week in aggregate. (We know that the true effect is -.2 since we've simulated this data.) |
| 65 | + |
| 66 | +We calculate the observed statistic in the paired setting in the same way that we would outside of the paired setting. Using `specify()` and `calculate()`: |
| 67 | + |
| 68 | +```{r calc-obs-mean} |
| 69 | +# calculate the observed statistic |
| 70 | +observed_statistic <- |
| 71 | + gss_paired %>% |
| 72 | + specify(response = diff) %>% |
| 73 | + calculate(stat = "mean") |
| 74 | +``` |
| 75 | + |
| 76 | +The observed statistic is `r observed_statistic`. Now, we want to compare this statistic to a null distribution, generated under the assumption that the true difference was actually zero, to get a sense of how likely it would be for us to see this observed difference if there were truly no change in hours worked per week in the population. |
| 77 | + |
| 78 | +Tests for paired data are carried out via the `null = "paired independence"` argument to `hypothesize()`. |
| 79 | + |
| 80 | +```{r generate-null} |
| 81 | +# generate the null distribution |
| 82 | +null_dist <- |
| 83 | + gss_paired %>% |
| 84 | + specify(response = diff) %>% |
| 85 | + hypothesize(null = "paired independence") %>% |
| 86 | + generate(reps = 1000, type = "permute") %>% |
| 87 | + calculate(stat = "mean") |
| 88 | + |
| 89 | +null_dist |
| 90 | +``` |
| 91 | + |
| 92 | +For each replicate, `generate()` carries out `type = "permute"` with `null = "paired independence"` by: |
| 93 | + |
| 94 | +* Randomly sampling a vector of signs (i.e. -1 or 1), probability .5 for either, with length equal to the input data, and |
| 95 | +* Multiplying the response variable by the vector of signs, "flipping" the observed values for a random subset of value in each replicate |
| 96 | + |
| 97 | +To get a sense for what this distribution looks like, and where our observed statistic falls, we can use `visualize()`: |
| 98 | + |
| 99 | +```{r visualize} |
| 100 | +# visualize the null distribution and test statistic |
| 101 | +null_dist %>% |
| 102 | + visualize() + |
| 103 | + shade_p_value(observed_statistic, |
| 104 | + direction = "two-sided") |
| 105 | +``` |
| 106 | + |
| 107 | +It looks like our observed mean of `r observed_statistic` would be relatively unlikely if there were truly no change in mean number of hours worked per week over this time period. |
| 108 | + |
| 109 | +More exactly, we can calculate the p-value: |
| 110 | + |
| 111 | +```{r p-value} |
| 112 | +# calculate the p value from the test statistic and null distribution |
| 113 | +p_value <- null_dist %>% |
| 114 | + get_p_value(obs_stat = observed_statistic, |
| 115 | + direction = "two-sided") |
| 116 | +
|
| 117 | +p_value |
| 118 | +``` |
| 119 | + |
| 120 | +Thus, if the change in mean number of hours worked per week over this time period were truly zero, our approximation of the probability that we would see a test statistic as or more extreme than `r observed_statistic` is approximately `r p_value`. |
| 121 | + |
| 122 | +We can also generate a bootstrap confidence interval for the mean paired difference using `type = "bootstrap"` in `generate()`. As before, we use the pre-computed differences when generating bootstrap resamples: |
| 123 | + |
| 124 | +```{r generate-boot} |
| 125 | +# generate a bootstrap distribution |
| 126 | +boot_dist <- |
| 127 | + gss_paired %>% |
| 128 | + specify(response = diff) %>% |
| 129 | + hypothesize(null = "paired independence") %>% |
| 130 | + generate(reps = 1000, type = "bootstrap") %>% |
| 131 | + calculate(stat = "mean") |
| 132 | + |
| 133 | +visualize(boot_dist) |
| 134 | +``` |
| 135 | + |
| 136 | +Note that, unlike the null distribution of test statistics generated earlier with `type = "permute"`, this distribution is centered at `observed_statistic`. |
| 137 | + |
| 138 | +Calculating a confidence interval: |
| 139 | + |
| 140 | +```{r confidence-interval} |
| 141 | +# calculate the confidence from the bootstrap distribution |
| 142 | +confidence_interval <- boot_dist %>% |
| 143 | + get_confidence_interval(level = .95) |
| 144 | +
|
| 145 | +confidence_interval |
| 146 | +``` |
| 147 | + |
| 148 | +By default, `get_confidence_interval()` constructs the lower and upper bounds by taking the observations at the $(1 - .95) / 2$ and $1 - ((1-.95) / 2)$th percentiles. To instead build the confidence interval using the standard error of the bootstrap distribution, we can write: |
| 149 | + |
| 150 | +```{r} |
| 151 | +boot_dist %>% |
| 152 | + get_confidence_interval(type = "se", |
| 153 | + point_estimate = observed_statistic, |
| 154 | + level = .95) |
| 155 | +``` |
| 156 | + |
| 157 | +To learn more about randomization-based inference for paired observations, see the relevant chapter in [Introduction to Modern Statistics](https://openintro-ims.netlify.app/inference-paired-means.html). |
0 commit comments