---
jupyter:
  jupytext:
    notebook_metadata_filter: all,-language_info
    split_at_heading: true
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.16.0
  kernelspec:
    display_name: Python 3 (ipykernel)
    language: python
    name: python3
---

# Permutations and pairs

Sometimes the difference we are interested in is a difference between pairs of values.

In the [idea of permutation](permutation_idea) page, we introduced data from an
experiment on effect of drinking beer on mosquitoes.   These data are a good example of a situation where we want to look at differences between pairs of values.

You may remember that the experiment involved two groups of participants.
A group of 24 participants drank beer for the experiment, and another group of
18 subjects drank water.   The experimenters took a measure of how attractive
each person was to mosquitoes.   Specifically they put each person in a tent,
from which there was an air tube leading to a closed box of 50 mosquitoes.
The experimenters then opened the box and counted how many mosquitoes flew down
the tube towards the tent containing the person.  This count is the "activated"
column in the dataset.

In fact, they did this procedure *twice* for each volunteer, once *before* they
drank their allocated drink, and once *after*.  The difference between before
and after is a measure of the difference to the mosquitoes after the subject
had had their drink.

Without further ado, let us load the data.

If you want to run this notebook on your own computer, Download the data from
{download}`mosquito_beer.csv <data/mosquito_beer.csv>`.

See [this
page](https://github.com/matthew-brett/datasets/tree/master/mosquito_beer) for
more details on the dataset, and [the data license page](data/LICENSE).


```{python}
# Import Numpy library, rename as "np"
import numpy as np
# Make random number generator.
rng = np.random.default_rng()
# Import Pandas library, rename as "pd"
import pandas as pd
# Safe setting for Pandas.
pd.set_option('mode.copy_on_write', True)

# Set up plotting
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
```

```{python}
mosquitoes = pd.read_csv('data/mosquito_beer.csv')
mosquitoes.head()
```

As you saw above, the experimental procedure left us with two mosquito counts
for each volunteer, one count taken *before* they had had their drink, and
another count *after*.

Here we collect those *before* and *after* values for each beer-drinking
volunteer.  Please ignore the code below, we will cover this kind of data
selection and organization later in the course.

```{python}
# Make new DataFrame with before and after for each volunteer.
# Run this cell for now.  We will cover this code later in the course.
# Just the beer drinkers.
beer = mosquitoes[mosquitoes['group'] == 'beer']
before = beer[beer['test'] == 'before']
after = beer[beer['test'] == 'after']
# Merge before and after rows for matching volunteers.
both = before.merge(after, on=['volunteer', 'group'],
                    suffixes=['_before', '_after'])
# Select the columns we're interested in.
before_after = both[['group', 'activated_before', 'activated_after']]
before_after
```

Here is our planned subtraction of the control `before` numbers from the
experimental `after` numbers.

```{python}
# Transfer to arrays for simplicity.
befores = np.array(before_after['activated_before'])
afters = np.array(before_after['activated_after'])
```

```{python}
actual_diffs = afters - befores
actual_diffs
```

Here we show the result using Pandas, of which more soon in the course:

```{python}
before_after['after_minus_before'] = actual_diffs
before_after
```

If our hypothesis is correct, we expect this difference (after minus before
counts) to be positive, on average.  Let's see what this average difference was
for our sample:

```{python}
actual_mean_diff = np.mean(actual_diffs)
actual_mean_diff
```

## Using permutation for pairs

We find that the difference is positive *for our sample*.   Our question of
course is whether this positive mean difference is compatible with sampling
variation — the differences we will expect to see given we have taken a sample
of beer-drinking people.

We now have to think about what our null world would be for such a mean
difference.

In the null world, there is 0 (not-any) average difference between the control
`before` scores and the corresponding `after` scores.  That is, the average
difference between these two scores will be 0.

How can we simulate such a world, where we expect the average difference
between this *pair* of scores to be 0?

If the null world it true, and the average difference is 0, then we can just do
a random swap of the before and after scores in the pair, and we'll still have
an observation that is valid in the null world.

That is, to make a new dataset that could occur in the null world, we could go
through row by row and, at random, swap the `before` and `after` scores.  Then
we would recalculate the mean difference, and this mean difference would be
a mean difference we might see in the null world, where there is no difference
on average between the two values in the pair.  Then we would do this same
procedure thousands of times to build up the *sampling distribution* of the
mean difference, and then we would compare our observed mean difference the
sampling distribution, to see if it was rare in the null world.

We could do this operation, of going through each row, and randomly flipping
the `before` and `after` values, but we can also simplify our task with a tiny
bit of algebra.

Let's say we have the subtraction between any two values $x$ and $y$: $x - y$, and we want the subtraction the other way round: $y - x$.  We can get the value for $y - x$ by multiplying $x - y$ by -1.

$$
-1 * (x - y) = -x + y = y - x
$$

We were thinking to randomly swap the two elements of the pair, and then
subtract the results, but we can get the same result by taking the differences
between the original pairs, and randomly choosing whether to multiply each
difference by -1.

Here we choose 1 or -1 at random for each row in our data frame.

```{python}
n = len(before_after)
# Choose 1 or -1 at random, n times.
rand_signs = rng.choice([-1, 1], size=n)
rand_signs
```

The values of -1 represent rows for which we are flipping the pairs, and values
of 1 correspond to pairs we have left in the original order.

We recalculate the differences, as we did above:

```{python}
actual_diffs = afters - befores
actual_diffs
```

This is the mean of the actual differences that we saw above.

```{python}
actual_mean_diff = np.mean(actual_diffs)
actual_mean_diff
```

Make a new set of difference corresponding to the flips encoded in the `signs` array:

```{python}
# Generate result of flipping pairs randomly
fake_diffs = rand_signs * actual_diffs
fake_diffs
```

Now we have a mean difference in the null world:

```{python}
fake_mean_diff = np.mean(fake_diffs)
fake_mean_diff
```

Here then is the procedure for one trial in the null world:

```{python}
# One trial in the null world.
rand_signs = rng.choice([-1, 1], size=n)
fake_diffs = rand_signs * actual_diffs
fake_mean_diff = np.mean(fake_diffs)
fake_mean_diff
```

Run this procedure a few times to get a feel for how much these numbers vary from trial to trial.

We repeat this procedure thousands of times to build up the sampling
distribution in the null world.

```{python}
results = np.zeros(10000)
for i in np.arange(10000):
    # Do the single trial procedure.
    rand_signs = rng.choice([-1, 1], size=n)
    fake_diffs = rand_signs * actual_diffs
    fake_mean_diff = np.mean(fake_diffs)
    # Store the result
    results[i] = fake_mean_diff

# Show the first 10 values
results[:10]
```

```{python}
plt.hist(results, bins=50)
plt.title('Sampling distribution of mean of differences')
# Show the position of the actual value on the x-axis.
plt.axvline(actual_mean_diff, color='red', label='Actual value')
# Label the actual value line.
plt.legend();
```

Finally we ask how unusual the actual value is in the sampling distribution from the null world:

```{python}
p = np.count_nonzero(results >= actual_mean_diff) / 10000
p
```

We have found that there is less than a 0.1% chance we would see the actual value, or greater, in the null world.  The actual value is surprising in the null world, and we have reason to continue to investigate causes of this value, including the presumed cause, that mosquitoes are, in fact, attracted to people who drink beer.
