# Simulation of offender strategies using squared differences


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

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

The first step is to examine our data and calculate this measure: The
mean is 11, the differences (call *deviations from the mean*) are 6, 1,
1, and 4, the squared deviations are 36, 1, 1, and 16, and their sum is
54. In R:

```{r}
# The actual scores for each treatment.
treat_scores <- c(17, 10, 10, 7)
# The mean.
m <- mean(treat_scores)
# The four deviations from the mean
deviations <- treat_scores - m
# Squared deviations.
sq_deviations <- deviations ** 2
# Sum of squared deviations.
actual_sum_sq_deviations <- sum(sq_deviations)
# Show the result.
actual_sum_sq_deviations
```

Our experiment will be, as before, to constitute four groups of 20 at
random from a universe with a 55 percent rehabilitation rate. We then
calculate this same measure for the random groups. If it is frequently
larger than 54, then we conclude that a uniform cure rate of 55 percent
could easily have produced the observed results. The program that
follows also generates the four treatments by using a `for` loop, rather
than spelling out the `sample` command 4 times as above. In R:

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

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

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

# Do 10000 trials.
for (i in 1:n_trials) {
    # Repeat the following steps 4 times to constitute 4 groups of 20 and
    # count their success rates.
    group_scores <- numeric(4)
    for (group_no in 1:4) {
        # A treatment group
        group <- sample(c('success', 'failure'),
                        size=20,
                        prob=c(0.55, 0.45),
                        replace=TRUE)
        # Count successes.
        group_count <- sum(group == 'success')
        # Store result.
        group_scores[group_no] <- group_count
        # End the procedure for one group of 20, go back and repeat until all 4
        # are done.
    }
    # Calculate the mean
    m <- mean(group_scores)
    # Calculate the deviations of the scores from the mean of the scores.
    deviations <- group_scores - m
    # Square them, making them them positive.
    sq_deviations <- deviations ** 2
    # Sum up the squared deviations.
    sum_sq_deviations <- sum(sq_deviations)
    # Keep track of the result for each trial.
    scores[i] <- sum_sq_deviations

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

# Produce a histogram of the trial results.
hist(scores, breaks=50,
     main='Distribution of sum of squared deviations, in null universe')
```

From this histogram, we see that in only about 1 percent of the cases
did our trial sum of squared differences equal or exceed 54, confirming
our conclusion that this is an unusual result. We can have R calculate
this proportion:

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

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

The conventional way to approach this problem would be with what is
known as a [chi-squared
test](https://en.wikipedia.org/wiki/Chi-squared_test).
