# Association of athletic and IQ scores


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

We use random pairings of the athletic and IQ scores to find the
null-world distribution of the sum of the products of these scores.

We load the file containing the data:

```{r}
# Read the data file containing athletic and IQ scores.
ath_iq_df <- read.csv('data/athletic_iq.csv')
# Show the data frame.
ath_iq_df
```

```{r}
# Turn athletic and IQ scores into vectors.
ath = ath_iq_df$athletic_score
iq = ath_iq_df$iq_score
```

```{r}
# Multiply the two sets of scores together.
actual_prod <- ath * iq
# Sum the results — the "observed value."
actual_sum <- sum(actual_prod)
actual_sum
```

```{r}
# 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) {
    # Shuffle the IQ scores so we can pair them against athletic scores.
    shuffled <- sample(iq)
    # Multiply the shuffled IQ scores by the athletic scores. (Note that we
    # could shuffle the athletic scores too but it would not achieve any
    # greater randomization, because shuffling one of the set of scores
    # already gives a random pairing between the sets of scores.
    fake_prod <- ath * shuffled
    # Sum the randomized multiplications.
    fake_sum <- sum(fake_prod)
    # Subtract the sum from the sum of the "observed" multiplication.
    diff <- actual_sum - fake_sum
    # Keep track of the result in results array.
    results[i] <- diff
    # End one trial, go back and repeat until 10000 trials are complete.
}

# Obtain a histogram of the trial results.
hist(results, breaks=25,
    main='Random sums of products',
    xlab='Observed sum minus sum of random pairing')
```

We see that obtaining a chance trial result as great as that observed
was rare. R will calculate this proportion for us:

```{r}
# Determine in how many trials the random sum of products was less than
# the observed sum of products.
k <- sum(results <= 0)
# Convert to a proportion.
kk <- k / n_trials
# Print the result.
message('Proportion of random pairings giving sum <= observed: ', kk)
```
