# ========================================================
# MULTINOMIAL LOGISTIC REGRESSION IN R
# PHAA outcomes data (the Lesson 4 dataset), 795 adults
# ========================================================
# Question: are age and smoking associated with where people
# usually go for health care (family doctor, walk-in clinic,
# emergency department, or no usual place)?
# ========================================================


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

library(nnet)    # multinom() fits the multinomial logistic model

# 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: usual place of care
table(phaa$usual_care)

# Make it a factor and choose the reference category
phaa$usual_care <- relevel(factor(phaa$usual_care), ref = "Family doctor")
levels(phaa$usual_care)

# --- 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$usual_care)

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

barplot(care_pct, beside = TRUE,
        names.arg = c("Family doctor", "Emergency", "No usual place", "Walk-in"),
        main = "Usual place of care by smoking status",
        xlab = "Usual place of care",
        ylab = "Percent of smoking group",
        col = c("grey80", "grey30"), ylim = c(0, 70),
        legend.text = c("Non-smoker", "Smoker"))

# --- 3.3 Age in each category of usual care
boxplot(age ~ usual_care, data = phaa,
        names = c("Family doctor", "Emergency", "No usual place", "Walk-in"),
        cex.axis = 0.85,
        main = "Age by usual place of care",
        xlab = "Usual place of care",
        ylab = "Age (years)",
        col = "grey80")


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

mn <- multinom(usual_care ~ age + smoker, data = phaa)
summary(mn)


# ========================================================
# SECTION 5: RELATIVE RISK RATIOS
# ========================================================

# Relative risk ratios (RRR): exponentiated coefficients
round(exp(coef(mn)), 2)

# 95% confidence intervals, one block per comparison
round(exp(confint(mn)), 2)

# RRR for a 10-year difference in age
round(exp(10 * coef(mn)[, "age"]), 2)


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

# Four people: aged 30 or 60, non-smoker or smoker
new_people <- data.frame(age = c(30, 30, 60, 60),
                         smoker = c("No", "Yes", "No", "Yes"))
probs <- predict(mn, newdata = new_people, type = "probs")
round(probs, 2)
rowSums(probs)    # each person's probabilities add to 1

barplot(probs, beside = TRUE,
        names.arg = c("Family doctor", "Emergency", "No usual place", "Walk-in"),
        main = "Predicted usual place of care",
        xlab = "Usual place of care",
        ylab = "Predicted probability",
        col = c("white", "grey75", "grey50", "grey20"), ylim = c(0, 0.8),
        legend.text = c("30, non-smoker", "30, smoker",
                        "60, non-smoker", "60, smoker"))


# ========================================================
# SECTION 7: CHANGING THE REFERENCE CATEGORY
# ========================================================

# The same model with Walk-in clinic as the reference
phaa$care_walkin <- relevel(phaa$usual_care, ref = "Walk-in clinic")
mn_walkin <- multinom(care_walkin ~ age + smoker, data = phaa,
                      trace = FALSE)
round(exp(coef(mn_walkin)), 2)

# The fit and the predicted probabilities do not change
AIC(mn, mn_walkin)
round(predict(mn_walkin, newdata = new_people, type = "probs"), 2)


# ========================================================
# SECTION 8: CHECK THE ASSUMPTIONS
# ========================================================

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

# --- 8.2 Enough people in every category
table(phaa$usual_care)
table(phaa$smoker, phaa$usual_care)

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

# --- 8.4 Independence of irrelevant alternatives (IIA)
# This assumption is judged by reasoning: no two categories
# should be close substitutes for each other.