survival_by_gender <- table(titanic$survived, titanic$gender)
survival_by_gender

##        
##         Female Male
##   FALSE    130 1366
##   TRUE     357  355

barplot(survival_by_gender,
        main = "Passengers on the Titanic",
        legend.text = c("Death", "Survivor"),
        args.legend = list(x = "topleft"))

men <- titanic[titanic$gender == "Male", ]
mean(!men$survived)

## [1] 0.7937246

women <- titanic[titanic$gender == "Female", ]
mean(!women$survived)

## [1] 0.2669405

# Alternative
aggregate(!survived ~ gender, data = titanic, mean)

##   gender !survived
## 1 Female 0.2669405
## 2   Male 0.7937246

79% of men died versus 27% of women.

1. Men dying with a three times higher probability than women is a meaningful difference. High male mortality is linked to the deaths of 1366 men. Therefore, the observed mortality ratio is practically significant.
2. Whether this difference is statistically significant depends on how likely such a difference is when survival does not depend on gender. We can find an answer by performing a permutation test.

1. The mortality ratio can take any nonnegative value:

   • It can be zero if all men survived and at least one woman died.
   • It can have an infinitely large positive value if at least one man died and all women survived.

   If all men and all women survived, the mortality ratio is undefined (zero divided by zero).

2. If men and women had died with the same probability, the response would have been approximately 1.

male_mortality <- mean(!men$survived)
female_mortality <- mean(!women$survived)
mortality_ratio <- male_mortality / female_mortality
mortality_ratio

## [1] 2.973414

The observed mortality ratio is approximately 2.97.

titanic$gender_p <- sample(titanic$gender)  # Permute gender
men_p <- titanic[titanic$gender_p == "Male", ]
women_p <- titanic[titanic$gender_p == "Female", ]
male_mortality_p <- mean(!men_p$survived)
female_mortality_p <- mean(!women_p$survived)
response <- male_mortality_p / female_mortality_p
response

## [1] 1.027646

Is it plausible that the response of the permuted data has a value close to 1?

Yes, the permutation broke the association between gender and survival; thus, the male and female mortalities are approximately equal.

n_trials <- 10000
responses <- numeric(n_trials)
for (trial in 1:n_trials) {
  titanic$gender_p <- sample(titanic$gender)  # Permute gender
  men_p <- titanic[titanic$gender_p == "Male", ]
  women_p <- titanic[titanic$gender_p == "Female", ]
  male_mortality_p <- mean(!men_p$survived)
  female_mortality_p <- mean(!women_p$survived)
  responses[trial] <- male_mortality_p / female_mortality_p
}

hist(responses)

The observed mortality ratio of 2.97 is statistically significant because all responses in the permutation test were smaller.

titanic$is_adult <- (titanic$age >= 18)
adults <- titanic[which(titanic$is_adult), ]  # Use which() to remove NAs
younger <- titanic[which(!titanic$is_adult), ]
adult_mortality <- mean(!adults$survived) 
younger_mortality <- mean(!younger$survived)
mortality_ratio <- adult_mortality / younger_mortality
mortality_ratio

## [1] 1.188829

It is a subjective judgement call whether to consider an adult mortality as practically significant if it is 19% higher than the mortality of younger people.

Compared to the effect of gender on survival, the effect of age is relatively small.

Still, the age-dependent difference had practical significance for those younger people whose lives were saved as a result.

n_trials <- 10000
responses <- numeric(n_trials)
for (trial in 1:n_trials) {
  titanic$is_adult_p <- sample(titanic$is_adult)
  adults_p <- titanic[which(titanic$is_adult_p), ]
  younger_p <- titanic[which(!titanic$is_adult_p), ]
  adult_mortality_p <- mean(!adults_p$survived)
  younger_mortality_p <- mean(!younger_p$survived)
  responses[trial] <- adult_mortality_p / younger_mortality_p
}

hist(responses)

# Insert vertical line at the observed value
abline(v = mortality_ratio, col = "blue", lwd = 3)

mean(responses > mortality_ratio)

## [1] 8e-04

The mortality ratio was more extreme than the observed value in only 0.08% of the permutations. We conclude that adults were statistically significantly more likely to die.

We cannot conclude that being an adult was the cause of a higher mortality. What might be a lurking variable?

slope <- coefficients(model)[2]
slope

##   hurdles 
## -0.238086

For every extra second an athlete needs in the hurdles event, the linear model predicts that the athlete throws the javelin 23.8 centimetres shorter.

summary(model)$r.squared

## [1] 0.0002317609

The model only predicts 0.02% of the variability of the javelin throw distances. However, a low \(R^2\) does not necessarily imply that the slope is insignificant (neither practically nor statistically).

Keep in mind that being better at hurdles is equivalent to having a smaller value in hurdles. However, being better in javelin throw translates into having a larger value in javelin. Thus, we might expect the slope of the linear model to be negative.

However, the observed effect is small. Even a one-second difference in hurdle time (which is a big difference in practice) only leads to a difference of a few centimetres in the predicted javelin throw distance (which has a standard deviation of 5.75 metres). Therefore, the observed slope is not practically significant.

n_trials <- 10000
responses <- numeric(n_trials)
for (trial in 1:n_trials) {
  hept$hurdles_p <- sample(hept$hurdles)  # Permute hurdles
  model_p <- lm(jt ~ hurdles_p, data = hept)
  responses[trial] <- coefficients(model_p)[2]
}

hist(responses)
abline(v = slope, col = "blue", lwd = 3)

What is the proportion of trials that produced a response that was even more strongly negative than the observed slope?

mean(responses < slope)

## [1] 0.4785

Almost half of the responses was more extreme than what we observed.

Therefore, there is no statistically significant evidence that better hurdle runners tend to be better javelin throwers.