# Simulation of offender strategies


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

This notebook uses simulation to test the null hypothesis that the four
treatments in question were each equally effective.

```{r}
# Set the number of trials
n_trials <- 10000

# Set the sample size for each trial
sample_size <- 20

# An empty array to store the trials
scores <- numeric(n_trials)

# Do 10000 trials
for (i in 1:n_trials) {
    # The first treatment group
    a <- sample(c('success', 'failure'),
                size=sample_size,
                prob=c(0.55, 0.45),
                replace=TRUE)
    # Count successes.
    a_count <- sum(a == 'success')
    # Second, third and fourth treatment groups, and successes.
    b <- sample(c('success', 'failure'),
                size=sample_size,
                prob=c(0.55, 0.45),
                replace=TRUE)
    b_count <- sum(b == 'success')
    c <- sample(c('success', 'failure'),
                size=sample_size,
                prob=c(0.55, 0.45),
                replace=TRUE)
    c_count <- sum(c == 'success')
    d <- sample(c('success', 'failure'),
                size=sample_size,
                prob=c(0.55, 0.45),
                replace=TRUE)
    d_count <- sum(d == 'success')

    # Now find all the pairwise differences in successes among the groups.
    a_minus_b <- a_count - b_count
    a_minus_c <- a_count - c_count
    a_minus_d <- a_count - d_count
    b_minus_c <- b_count - c_count
    b_minus_d <- b_count - d_count
    c_minus_d <- c_count - d_count

    # Concatenate, or join, all the differences in a single array "diffs".
    diffs <- c(a_minus_b, a_minus_c, a_minus_d,
               b_minus_c, b_minus_d, c_minus_d)
    # Since we are interested only in the magnitude of the difference, not its
    # direction, we take the absolute value of all the differences (we remove
    # any minus signs, making all values positive).
    abs_diffs <- abs(diffs)
    # Find the largest of all the differences
    max_abs_diff <- max(abs_diffs)
    # Keep score of the largest
    scores[i] <- max_abs_diff

    # End a trial, go back and repeat until all 10000 are complete.
}

# How many of the trials yielded a maximum difference greater than the
# observed maximum difference?
k <- sum(scores >= 10)
# Convert to a proportion
kk <- k / n_trials

message('Proportion >= 10 was: ', kk)
```
