# Statistical Modelling, Michaelmas 2026
# Carried forward from 2025; sections are organised by topic, not lecture date.
# Run from statistical-modelling/2026, where data/yule.csv is available.
# For a downloaded copy, see the Yule data link on the course webpage.

########################################
########### Nested projection ##########
########################################

?LifeCycleSavings

model.matrix(sr ~ ddpi + pop15 + pop75, data = LifeCycleSavings)

fit1 <- lm(sr ~ ddpi + pop15 + pop75, data = LifeCycleSavings)
fit1

fit2 <- lm(LifeCycleSavings$sr ~ residuals(lm(ddpi ~ pop15 + pop75, data = LifeCycleSavings)))
fit2

fit3 <- lm(residuals(lm(sr ~ pop15 + pop75, data = LifeCycleSavings)) ~ residuals(lm(ddpi ~ pop15 + pop75, data = LifeCycleSavings)))
fit3

summary(fit1)
summary(fit3)

########################################
###### Data handling; linear models ####
########################################

library(data.table)

system.time(flights <- read.csv("https://raw.githubusercontent.com/Rdatatable/data.table/master/vignettes/flights14.csv"))

system.time(flights <- fread("https://raw.githubusercontent.com/Rdatatable/data.table/master/vignettes/flights14.csv"))

class(flights)

dim(flights)

print(flights, 2)

summary(flights)

table(flights$carrier)

## DT[i, j, by]: take rows "i", calculate "j", grouped by "by"
flights[dep_delay > 120 & month == 6, table(carrier)]

delay_summary <-
    flights[, .(total = .N,
                delayed = sum(dep_delay > 120),
                percent = mean(dep_delay > 120),
                mean = mean(dep_delay)),
            by = .(carrier, dest)]
delay_summary[total > 30][order(- percent)][1:5]

## Data visualization

## plot, hist, boxplot

x <- rnorm(100)
y <- x + rnorm(100)
plot(x, y)

hist(x)
hist(flights$air_time)
boxplot(air_time ~ carrier, flights)

## ggplot2 package (not examinable)

library(ggplot2)
ggplot(delay_summary[total > 30],
       aes(x = percent, y = mean, col = carrier, size = total)) +
    geom_point(alpha = 0.5) +
    theme_minimal(base_size = 15) +
    scale_x_continuous(labels = scales::percent) +
    xlab("Delay > 2 hours") + ylab("Average delay (in minutes)")

## Normal linear models

set.seed(42)
n <- 50
a <- rnorm(n)
b <- rnorm(n)
y <- a - b + 2 * rnorm(n)
data <- data.frame(y = y, x1 = a, x2 = b)

?lm

fit1 <- lm(y ~ x1 + x2, data)
summary(fit1)

fit2 <- lm(y ~ x1, data)
anova(fit2, fit1)

fit3 <- lm(y ~ 1, data)
anova(fit3, fit1) # What's the output based on summary(fit1)?

mylm <- function(X, Y) {

    n <- nrow(X)
    p <- ncol(X)
    stopifnot(length(Y) == n)

    beta.hat <- solve(t(X) %*% X) %*% t(X) %*% Y
    R <- Y - X %*% beta.hat
    sigma.hat <- sqrt(sum(R^2) / (n-p))
    SE <- sqrt(diag(solve(t(X) %*% X))) * sigma.hat

    print("Residuals:")
    print(quantile(R))

    print("Coefficients:")
    print(cbind(beta.hat, SE))

}

X <- model.matrix(~ x1 + x2, data)
mylm(X, y)

## plot.lm

set.seed(42)
n <- 100
a <- rnorm(n)
b <- rnorm(n)
y <- a - b + 2 * rt(n, df = 5)
data <- data.frame(y = y, x1 = a, x2 = b)

fit <- lm(y ~ x1 + x2, data)

par(mfrow = c(2, 2))
plot(fit)

############################################################
############### Interpreting linear models #################
############################################################

## Practical 4: Interpreting linear models

## Yule-Simpson paradox

library(palmerpenguins)
plot(body_mass_g ~ bill_depth_mm, penguins, col = species)

(fit <- lm(body_mass_g ~ bill_depth_mm, penguins))
abline(fit, col = "black")

fit2 <- lm(body_mass_g ~ 0 + species + bill_depth_mm, penguins)
summary(fit2)

abline(fit2$coef[1], fit2$coef[4], col = "black")
abline(fit2$coef[2], fit2$coef[4], col = "red")
abline(fit2$coef[3], fit2$coef[4], col = "green")

lm(body_mass_g ~ species * bill_depth_mm, penguins)

fit3 <- lm(body_mass_g ~ 0 + species + species:bill_depth_mm, penguins)
abline(fit3$coef[1], fit3$coef[4], col = "black")
abline(fit3$coef[2], fit3$coef[5], col = "red")
abline(fit3$coef[3], fit3$coef[6], col = "green")


## Yule on the causes of poverty

data_path <- if (file.exists("./data/yule.csv")) "./data/yule.csv" else "./yule.csv"
data <- read.table(data_path)
data

fit <- lm(paup ~ outrelief + old + pop, data)
summary(fit)

## Economic return to schooling

## Based on Card (1995)

library(ivmodel)
vars <- c("lwage", "educ", "exper", "expersq", "black", "south", "smsa", "nearc4")

summary(fit <- lm(lwage ~ educ + exper + expersq + black + south + smsa, card.data))

fit2 <- ivmodelFormula(lwage ~ educ + exper + expersq + black + south + smsa | nearc4 + exper + expersq + black + south + smsa, card.data)
fit2

card.data$educ.fitted <- fitted(lm(educ ~ nearc4 + exper + expersq + black + south + smsa, card.data))
summary(lm(lwage ~ educ.fitted + exper + expersq + black + south + smsa, card.data))

############################################################
######################## Practical 5 #######################
############################################################

n <- 10000
k <- 50
one.sim <- function() {
    pi <- rbeta(n / k, 2, 2)
    Z <- rbinom(n, 1, pi)
    sum(Z)
}
Y <- replicate(1000, one.sim())

c(mean(Y) / n, var(Y) / n)

(sigma2.hat <- var(Y) / n / 0.5^2 -1)
(sigma2.theory <- (k-1) * 0.05 / 0.5^2)

##

alpha <- 5
lambda <- 1
mu <- alpha / lambda
deviance.gamma <- function(mu1, mu2, alpha) {
    2 * alpha * (log(mu2 / mu1) + mu1 / mu2 - 1)
}
Y <- rgamma(1000000, alpha, 1)
dev <- deviance.gamma(Y, mu, alpha)
resid.dev <- sign(Y - mu) * sqrt(dev)
resid.pearson <- (Y - mu) / sqrt(alpha/lambda^2)

s <- 1:1000
pl <- qqnorm(resid.dev[s], col = scales::alpha("red", 0.4))
points(pl$x, resid.pearson[s], col = scales::alpha("blue", 0.4))
abline(0, 1, lty = "dashed")

gamma <- 2 / sqrt(alpha)
delta <- 6 / alpha
a <- 2 / sqrt(alpha) / 6
b <- (14 * gamma^2 - 9 * delta) / 72
resid.bartlett <- (resid.dev + a) / (1 + b)

points(pl$x, resid.bartlett[s], col = scales::alpha("green", 0.4))
legend("topleft", c("Pearson", "Deviance", "Bartlett"),
       fill = c("blue", "red", "green"))


gammaMLE <- function(Y, alpha, interval = FALSE, ci.method = c("asymptotic", "exact")) {

    ci.method <- match.arg(ci.method)

    lambda.hat <- alpha / mean(Y)
    if (!interval) {
        return(lambda.hat)
    } else {
        if (ci.method == "asymptotic") {
            SE <- sqrt(lambda.hat^2 / alpha / length(Y))
            return(c(lambda.hat,
                     lambda.hat - qnorm(0.975) * SE,
                     lambda.hat + qnorm(0.975) * SE))
        } else if (ci.method == "exact") {
            cdf.diff <- function(lambda, level) {pgamma(sum(Y), n * alpha, lambda) - level}
            return(c(lambda.hat,
                     uniroot(cdf.diff, c(0.1, 100), level = 0.025)$root,
                     uniroot(cdf.diff, c(0.1, 100), level = 0.975)$root))
        } else {
            stop("CI method not supported.")
        }
    }
}

one.sim <- function(n, alpha, lambda, ...) {
    Y <- rgamma(n, alpha, lambda)
    gammaMLE(Y, alpha, ...)
}

alpha <- 2
lambda <- 5

n <- 30
lambda.hat <- replicate(1000, one.sim(n, alpha, lambda))
summary(lambda.hat)

par(mfrow = c(1, 2))
hist(sqrt(n) * (lambda.hat - lambda), main = "n=30")
var(sqrt(n) * (lambda.hat - lambda))

n <- 10000
lambda.hat <- replicate(1000, one.sim(n, alpha, lambda))
hist(sqrt(n) * (lambda.hat - lambda), main = "n=10000")

## Monte-Carlo estimate of limiting variance
var(sqrt(n) * (lambda.hat - lambda))

## Theoretical variance
lambda^2 / alpha

## Asymptotic confidence interval
n <- 30
lambda.hat <- replicate(10000, one.sim(n, alpha, lambda, interval = TRUE))

mean((lambda.hat[2, ] < lambda) & (lambda.hat[3, ] > lambda))

## Exact distribution
par(mfrow = c(1, 1))
qqplot(lambda.hat[1, ], n * alpha / rgamma(1000, n * alpha, lambda))
abline(0, 1, col = "red")

## Exact confidence interval
n <- 30
lambda.hat <- replicate(10000, one.sim(n, alpha, lambda, interval = TRUE, ci.method = "exact"))

mean((lambda.hat[2, ] < lambda) & (lambda.hat[3, ] > lambda))

############################################################
## Practical 6: Binomial GLMs

path <- "http://www.statslab.cam.ac.uk/~rds37/teaching/statistical_modelling/"
Myopia <- read.csv(file.path(path, "Myopia.csv"))
head(Myopia)

fit1 <- glm(myopic ~ ., data = Myopia, family = binomial)
summary(fit1)

## Test whether the intercept only model is correct

pchisq(480.08 - 439.60, df = 617 - 609, lower.tail = FALSE)

fit2 <- glm(myopic ~ 1, data = Myopia, family = binomial)

anova(fit2, fit1)

## Smoing

Smoking <- read.csv(paste0(path, "Smoking.csv"))

Smoking$Total <- Smoking$Survived + Smoking$Died
Smoking$propDied <- Smoking$Died / Smoking$Total

plot(propDied ~ Age.group, data = Smoking, subset = (Smoker == "Yes"))
points(propDied ~ Age.group, data = Smoking, subset = (Smoker == "No"), pch = 4)

logit <- function(p) {log(p / (1-p))} ## Or use qlogis

plot(logit(propDied) ~ Age.group, data = Smoking, subset = (Smoker == "Yes"))
points(logit(propDied) ~ Age.group, data = Smoking, subset = (Smoker == "No"), pch = 4)

plot(qlogis(propDied) ~ Age.group, data = Smoking, subset = (Smoker == "Yes"))
points(qlogis(propDied) ~ Age.group, data = Smoking, subset = (Smoker == "No"), pch = 4)

fit3 <- glm(propDied ~ Age.group + Smoker, data = Smoking,
            family = binomial, weights = Total)
summary(fit3)

summary(fit4 <- glm(propDied ~ Age.group + I(Age.group^2) + Smoker, data = Smoking,
            family = binomial, weights = Total))

summary(fit5 <- glm(propDied ~ factor(Age.group) + Smoker, data = Smoking,
            family = binomial, weights = Total))

## Exercises

## Exercise 4.5

fit6 <- glm(myopic ~ . - TVHR - compHR, data = Myopia, family = binomial)
anova(fit6, fit1)

## Exercise 4.6

Myopia$parentMyopic <- (Myopia$dadMyopic == "Yes") + (Myopia$mumMyopic == "Yes")

fit7 <- glm(myopic ~ . - dadMyopic - mumMyopic + parentMyopic, data = Myopia, family = binomial)
anova(fit7, fit1)

summary(fit8 <- glm(myopic ~ . - dadMyopic, data = Myopia, family = binomial))

summary(fit9 <- glm(myopic ~ . - dadMyopic - mumMyopic + I((mumMyopic == "Yes") - (dadMyopic == "Yes")), data = Myopia, family = binomial))

## Exercise 4.7
summary(fit10 <- glm(myopic ~ . - dadMyopic - mumMyopic + I((mumMyopic == "Yes") * (dadMyopic == "Yes")), data = Myopia, family = binomial))


############################################################

## Practical 7

path <- "http://www.statslab.cam.ac.uk/~rds37/teaching/statistical_modelling/"
football2014 <- read.csv(file.path(path, "football2014.csv"), stringsAsFactors = TRUE)

nrow(football2014) / 2

football2014$By <- relevel(football2014$By, "Leicester")
football2014$Against <- relevel(football2014$Against, "Leicester")

summary(fit.football <- glm(GoalsScored ~ HomeAway + By + Against, data = football2014, family = poisson))

team <- names(coef(fit.football))[3:21]
team <- substr(team, 3, nchar(team))
team_strength <- data.frame(team = team,
attack_rating = coef(fit.football)[3:21],
defense_rating = - coef(fit.football)[22:40])
team_strength$team <- with(team_strength,
reorder(team, attack_rating + defense_rating))
library(ggplot2)
ggplot(reshape2::melt(team_strength), aes(x = team, fill = variable)) +
geom_bar(aes(y = value), stat = "identity", position = "dodge") +
    theme_minimal(base_size = 25) + coord_flip()


football2014_rem <- read.csv(file.path(path, "fixtures_remaining.csv"))

pred <- predict(fit.football, newdata = football2014_rem, type = "response")

head(cbind(GoalsPred = pred, football2014_rem))

standing <- function(data) {
library(data.table)
data <- data.table(data)
out <- data[1:(nrow(data)/2), c("By", "Against", "GoalsScored")]
names(out) <- c("HomeTeam", "AwayTeam", "GoalsHome")
out$GoalsAway <- data[- (1:(nrow(data)/2)), ]$GoalsScored
out[, HomeScore := 3 * (GoalsHome > GoalsAway) + (GoalsHome == GoalsAway)]
out[, AwayScore := 3 * (GoalsAway > GoalsHome) + (GoalsHome == GoalsAway)]
out1 <- out[, .(HomeScore = sum(HomeScore)), by = HomeTeam]
out2 <- out[, .(AwayScore = sum(AwayScore)), by = AwayTeam]
out <- merge(out1, out2, by.x = "HomeTeam", by.y = "AwayTeam")
names(out)[1] <- "Team"
out[, Score := HomeScore + AwayScore]
out[order(- Score)]
}

one.sim <- function(sim) {
    football2014_rem$GoalsScored <- rpois(nrow(football2014_rem), pred)
    out <- merge(standing(football2014), standing(football2014_rem), by ="Team")
    out <- out[, .(Score = Score.x + Score.y), by = "Team"]
    out <- out[order(-Score)]
    out$Ranking <- 1:nrow(out)
    out
}

out <- do.call(rbind ,lapply(1:1000, one.sim))

out[, .(Score.med = median(Score),
        Ranking.med = median(Ranking),
        Champion.prob = mean(Ranking == 1),
        Top4.prob = mean(Ranking <= 4),
        Bottom3.prob = mean(Ranking >= 18)),
    by = Team][order(- Score.med)]

##

SD_data <- read.csv(file.path(path, "SD_match.csv"), stringsAsFactors=T)

head(SD_data)

(SD_table <- table(SD_data[, c("match", "subject_m", "subject_f")]))

head(as.data.frame(SD_table))

fit1 <- glm(Freq ~ subject_m + subject_f + match,
data = SD_table, family = poisson)
fit2 <- glm(Freq ~ (subject_m + subject_f + match)^2,
            data = SD_table, family = poisson)

anova(fit1, fit2, test = "LRT")

summary(fit2)

fit3 <- glm(Freq ~ subject_m + subject_f + match,
data = SD_table, family = quasipoisson)
fit4 <- glm(Freq ~ (subject_m + subject_f + match)^2,
data = SD_table, family = quasipoisson)

## Exercise 4.14

fit.football2 <- glm(GoalsScored ~ HomeAway*By + Against, data = football2014, family = poisson)
anova(fit.football, fit.football2, test = "LRT")

fit.football3 <- glm(GoalsScored ~ HomeAway*By + HomeAway*Against, data = football2014, family = poisson)
anova(fit.football, fit.football3, test = "LRT")

## Exercise 4.15
anova(fit1, fit2, test = "LRT")
anova(fit3, fit4, test = "LRT")

pchisq(144.43 / summary(fit4)$dispersion, df = 15, lower.tail = FALSE)

## Exercise 4.16
SD_table <- table(SD_data[, c("match", "subject_m", "subject_f", "goal_m", "goal_f")])

library(MASS)
stepAIC(glm(Freq ~ .^3, data = SD_table, family = poisson))