# Sample size for pig rations with bootstrap


Find this notebook on the web at
<a class="quarto-xref" href="https://resampling-stats.github.io/latest-r/how_big_sample.html#nte-sample_size_rations">Note <span>30.6</span></a>.

Use the bootstrap (resampling with replacement) procedure to find the
plausible differences between two rations in the benchmark (combined)
universe.

In due course, we then use these plausible values to put ourselves in
the situation where the distributions of rations are distinct and
different, and see how often we get a simulated difference drawn from
these distributions that is implausible on criterion above.

First we need to get the measured data from the data file using R:

```{r}
# Estimated weights for ration A.
a_weights <- c(36, 35, 34, 33, 33, 32, 32, 31, 31, 30, 29, 28)
# Estimated weights for ration B.
b_weights <- c(32, 31, 30, 29, 29, 28, 28, 27, 27, 26, 25, 24)

# Make a combined (benchmark) universe from the weights.
both = c(a_weights, b_weights)

# Set the number of trials
n_trials <- 10000

# An empty array to store the trial results.
results <- numeric(n_trials)

# Do 10,000 experiments.
for (i in 1:n_trials) {
    # Take a "resample" of 12 with replacement from both and put it in fake_a
    fake_a <- sample(both, size=12, replace=TRUE)
    # Likewise to make fake_b
    fake_b <- sample(both, size=12, replace=TRUE)
    # Mean the first "resample."
    fake_a_mean <- mean(fake_a)
    # Mean the second "resample."
    fake_b_mean <- mean(fake_b)
    # Calculate the difference between the two resamples.
    fake_diff <- fake_a_mean - fake_b_mean
    # Keep track of each trial result.
    results[i] <- fake_diff
    # End one experiment, go back and repeat until all trials are complete,
    # then proceed.
}
# Produce a histogram of trial results.
hist(results, breaks=25,
      xlab='Second resample mean minus first',
      main='Distribution difference in means of resamples')

# Get the 95% percentile.  Only 5% of results are above this value, by chance.
q_95 <- quantile(results, 0.95)

message('95% quantile for resampled mean difference: ', q_95)
```

```{r}
# Check the quantile.
k <- sum(results > q_95)
kk = k / n_trials

message('Proportion > 95% quantile: ', kk)
```

We find that in only of the trials the difference between the two
observed groups of 12 pigs was more than 2.083. Now we investigate how
often samples of 12 pigs, drawn from the *separate* universes, will show
a mean difference larger than 2.083 pounds. We do so by making up a deck
of 25 or 50 cards for *each* of the 12 hypothesized A’s and each of the
12 B’s, with the ration name and the weight gain written on it — that
is, a deck of, say, 300 cards for each ration. Then from each deck we
draw a set of 12 cards at random, record the group averages, and find
the difference.

Here is the same work done with more runs on the computer. In this
version we are sampling from the separate A and B universes we have
estimated.

```{r}
# A new empty vector to store the trial results.
results <- numeric(n_trials)

# Do 10,000 experiments.
for (i in 1:n_trials) {
    # Take a "resample" of 12 with replacement from A and put it in fake_a.
    # Notice we are sampling from "a_weights" this time.
    fake_a <- sample(a_weights, size=12, replace=TRUE)
    # Likewise to make fake_b
    # Notice we are sampling from "b_weights" this time.
    fake_b <- sample(b_weights, size=12, replace=TRUE)
    # Mean the first "resample."
    fake_a_mean <- mean(fake_a)
    # Mean the second "resample."
    fake_b_mean <- mean(fake_b)
    # Calculate the difference between the two resamples.
    fake_diff <- fake_a_mean - fake_b_mean
    # Keep track of each trial result.
    results[i] <- fake_diff
    # End one experiment, go back and repeat until all trials are complete,
    # then proceed.
}
# Produce a histogram of trial results.
hist(results, breaks=25,
      main='Second resample mean minus first, in separate universes',
      xlab='Distribution difference in means of resamples')

# What proportion of the differences are implausible in the combined universe?
k <- sum(results >= q_95)
kk <- k / n_trials

message('p for separate universe results >= 95% quantile: ', kk)
```

If `kk` is close to one, we know that the sample size is large enough
that samples drawn from the universes we have hypothesized will not
mislead us into thinking that they could come from the same universe.

Therefore, two samples of twelve pigs each are clearly large enough,
and, in fact, even smaller samples might be sufficient if the universes
are really like those we guessed at. If, on the other hand, the
differences in the guessed universes had been smaller, then twelve-pig
groups would have seemed too small and we would then have had to try out
larger sample sizes, say forty-eight pigs in each group and perhaps 200
pigs in each group if forty-eight were not enough. And so on until the
sample size is large enough to promise the accuracy we want. (In that
case, the decks would also have to be much larger, of course.)

If we had guessed different universes for the two rations, then the
sample sizes required would have been larger or smaller. If we had
guessed the averages for the two samples to be closer together, then we
would have needed larger samples. Also, if we had guessed the weight
gains *within* each universe to be less spread out, the samples could
have been smaller and vice versa.
