# ========================================================
# ORDINAL LOGISTIC REGRESSION IN R
# PHAA outcomes data (the Lesson 4 dataset), 795 adults
# ========================================================
# Question: are age and smoking associated with how people
# rate their own health (Poor, Fair, Good, Very good,
# Excellent)?
# ========================================================


# ========================================================
# SECTION 1: PACKAGES AND DATA
# ========================================================

# install.packages("brant")   # run once, if brant is not installed
library(MASS)    # polr() fits the ordinal logistic model
library(brant)   # brant() checks the proportional-odds assumption

# The data file is stored on the course website
phaa_url <- "https://www.sfu-epi.ca/r-activities/data/phaa_outcomes.csv"
phaa <- read.csv(phaa_url)
dim(phaa)

# Keep the people with no missing values
phaa <- na.omit(phaa)
nrow(phaa)


# ========================================================
# SECTION 2: PREPARE THE OUTCOME AND THE PREDICTORS
# ========================================================

# --- 2.1 The outcome: self-rated health
# read.csv() reads categories as text, so table() sorts them
# alphabetically
table(phaa$self_rated_health)

# Set the categories in their natural order, lowest first
phaa$self_rated_health <- factor(phaa$self_rated_health,
  levels = c("Poor", "Fair", "Good", "Very good", "Excellent"),
  ordered = TRUE)
table(phaa$self_rated_health)

# --- 2.2 The predictors
summary(phaa$age)
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
table(phaa$smoker)


# ========================================================
# SECTION 3: LOOK AT THE DATA BEFORE MODELLING
# ========================================================

# --- 3.1 Counts in every category, by smoking status
table(phaa$smoker, phaa$self_rated_health)

# --- 3.2 Row percentages: each smoking group adds to 100
srh_pct <- round(100 * prop.table(table(phaa$smoker,
                                        phaa$self_rated_health),
                                  margin = 1), 1)
srh_pct

barplot(srh_pct, beside = TRUE,
        main = "Self-rated health by smoking status",
        xlab = "Self-rated health",
        ylab = "Percent of smoking group",
        col = c("grey80", "grey30"), ylim = c(0, 40),
        legend.text = c("Non-smoker", "Smoker"))

# --- 3.3 Age in each category of self-rated health
boxplot(age ~ self_rated_health, data = phaa,
        main = "Age by self-rated health",
        xlab = "Self-rated health",
        ylab = "Age (years)",
        col = "grey80")


# ========================================================
# SECTION 4: FIT THE ORDINAL LOGISTIC MODEL
# ========================================================

ord <- polr(self_rated_health ~ age + smoker, data = phaa, Hess = TRUE)
summary(ord)


# ========================================================
# SECTION 5: ODDS RATIOS
# ========================================================

# Odds ratios with 95% confidence intervals
exp(cbind(OR = coef(ord), confint(ord)))

# Odds ratio for a 10-year difference in age
exp(10 * coef(ord)["age"])


# ========================================================
# SECTION 6: PREDICTED PROBABILITIES
# ========================================================

# Two 40-year-olds who differ only in smoking status
new_people <- data.frame(age = 40, smoker = c("No", "Yes"))
probs_40 <- predict(ord, newdata = new_people, type = "probs")
round(probs_40, 2)
rowSums(probs_40)    # each person's probabilities add to 1

barplot(probs_40, beside = TRUE,
        main = "Predicted self-rated health at age 40",
        xlab = "Self-rated health",
        ylab = "Predicted probability",
        col = c("grey80", "grey30"), ylim = c(0, 0.4),
        legend.text = c("Non-smoker", "Smoker"))


# ========================================================
# SECTION 7: CHECK THE ASSUMPTIONS
# ========================================================

# --- 7.1 Independence: one row per person
sum(duplicated(phaa$id))

# --- 7.2 Enough people in each category: see the table in 3.1

# --- 7.3 Multicollinearity: are the predictors strongly related?
tapply(phaa$age, phaa$smoker, mean)

# --- 7.4 Proportional odds: the Brant test
brant(ord)


# ========================================================
# SECTION 8: WHAT THE BRANT TEST COMPARES
# ========================================================

# One binary outcome for each split: 1 = at or above, 0 = below
phaa$fair_up <- 0
phaa$fair_up[phaa$self_rated_health >= "Fair"] <- 1
phaa$good_up <- 0
phaa$good_up[phaa$self_rated_health >= "Good"] <- 1
phaa$vgood_up <- 0
phaa$vgood_up[phaa$self_rated_health >= "Very good"] <- 1
phaa$exc_up <- 0
phaa$exc_up[phaa$self_rated_health >= "Excellent"] <- 1

# One binary logistic regression for each split
split1 <- glm(fair_up ~ age + smoker, data = phaa, family = binomial)
split2 <- glm(good_up ~ age + smoker, data = phaa, family = binomial)
split3 <- glm(vgood_up ~ age + smoker, data = phaa, family = binomial)
split4 <- glm(exc_up ~ age + smoker, data = phaa, family = binomial)

# The odds ratio for each predictor at each split
split_or <- exp(rbind(coef(split1), coef(split2),
                      coef(split3), coef(split4)))
round(split_or[, c("age", "smokerYes")], 3)

barplot(split_or[, "smokerYes"],
        names.arg = c("Fair+", "Good+", "Very good+", "Excellent"),
        main = "Odds ratio for smokers at each split",
        xlab = "Split (at or above versus below)",
        ylab = "Odds ratio",
        col = "grey80", ylim = c(0, 1))
abline(h = exp(coef(ord)["smokerYes"]), lty = 2, lwd = 2)