Page's Trend Test

# Page's Trend Test: Drug Dose-Response Example
# ================================================================
# Research Question: Do increasing analgesic doses produce 
# monotonically increasing pain relief?

# Load packages
library(crank)       # page.test()
library(tidyverse)
library(DescTools)   # KendallW()

# Create data: 15 patients × 4 doses
set.seed(123)
patient_id <- rep(1:15, each = 4)
dose <- rep(c(0, 50, 100, 200), times = 15)

# Simulate pain relief with monotonic trend + noise
base_relief <- case_when(
  dose == 0   ~ rnorm(15, mean = 2, sd = 1),
  dose == 50  ~ rnorm(15, mean = 4, sd = 1.2),
  dose == 100 ~ rnorm(15, mean = 6, sd = 1.3),
  dose == 200 ~ rnorm(15, mean = 7.5, sd = 1.5)
)

pain_data <- data.frame(
  patient = factor(patient_id),
  dose = factor(dose, levels = c("0", "50", "100", "200")),
  dose_numeric = dose,
  relief = pmax(0, pmin(10, base_relief))  # Bound 0-10
)

# -------- ASSUMPTION CHECKS --------

# 1. Check design: within-subjects, ordered conditions
cat("Design check:\n")
cat("Subjects:", length(unique(pain_data$patient)), "\n")
cat("Conditions(ordered):", nlevels(pain_data$dose), "\n")
cat("Measurements per subject:", table(table(pain_data$patient)), "\n\n")

# 2. Visualize individual trends (spaghetti plot)
ggplot(pain_data, aes(x = dose, y = relief, group = patient)) +
  geom_line(alpha = 0.3) +
  geom_point(alpha = 0.5) +
  stat_summary(aes(group = 1), fun = mean, geom = "line", 
               color = "red", size = 1.5) +
  labs(title = "Individual Dose-Response Profiles",
       subtitle = "Red line = mean trend",
       x = "Dose(mg)", y = "Pain Relief(0-10)") +
  theme_minimal()

# 3. Check for monotonic pattern in means
mean_relief <- pain_data %>%
  group_by(dose) %>%
  summarise(mean_relief = mean(relief),
            median_relief = median(relief),
            se = sd(relief)/sqrt(n()))
print(mean_relief)

# 4. Calculate within-subject ranks
pain_wide <- pain_data %>%
  pivot_wider(names_from = dose, values_from = relief)

ranks <- t(apply(pain_wide[, -1], 1, rank))
colnames(ranks) <- c("Dose_0", "Dose_50", "Dose_100", "Dose_200")
cat("\nMean ranks per dose:\n")
print(colMeans(ranks))

# -------- PAGE'S TREND TEST --------

# Specify expected order: 1, 2, 3, 4 (increasing)
expected_order <- c(1, 2, 3, 4)

# Run Page's L test
page_result <- page.test(
  x = pain_wide[, c("0", "50", "100", "200")],
  ranks = expected_order
)

cat("\n=== PAGE'S TREND TEST ===\n")
print(page_result)

# -------- EFFECT SIZE --------

# Standardized L (L ranges from minimum to maximum possible)
L_obs <- page_result$statistic
k <- 4  # conditions
n <- 15  # subjects

# L ranges from k(k+1)n/2 (no trend) to k²(k+1)n/4 (perfect trend)
L_min <- k * (k + 1) * n / 2
L_max <- k^2 * (k + 1) * n / 4
L_standardized <- (L_obs - L_min) / (L_max - L_min)

cat("\nEffect Size:\n")
cat("L statistic:", L_obs, "\n")
cat("Standardized L:", round(L_standardized, 3), "(0=no trend, 1=perfect)\n")

# Kendall's W (concordance coefficient)
friedman_data <- as.matrix(pain_wide[, -1])
kendall_w <- KendallW(friedman_data, correct = TRUE)
cat("Kendall's W:", round(kendall_w, 3), "(concordance)\n")

# -------- COMPARISON WITH FRIEDMAN --------

friedman_result <- friedman.test(relief ~ dose | patient, data = pain_data)
cat("\nFriedman test(omnibus):\n")
cat("Chi-square =", friedman_result$statistic, ", p =", 
    format.pval(friedman_result$p.value, digits = 3), "\n")

# -------- VISUALIZATION --------

# Plot mean trend with error bars
ggplot(mean_relief, aes(x = dose, y = mean_relief, group = 1)) +
  geom_line(size = 1.2, color = "steelblue") +
  geom_point(size = 3, color = "steelblue") +
  geom_errorbar(aes(ymin = mean_relief - se, ymax = mean_relief + se),
                width = 0.2) +
  labs(title = "Dose-Response Trend: Mean Pain Relief",
       x = "Dose(mg)", y = "Mean Pain Relief ± SE") +
  theme_minimal()

# -------- INTERPRETATION --------

cat("\n=== INTERPRETATION ===\n")
if (page_result$p.value < 0.05) {
  cat("Significant monotonic trend detected(p <", 
      format.pval(page_result$p.value, digits = 3), ")\n")
  cat("Conclusion: Pain relief increases monotonically with dose.\n")
  cat("Effect size(standardized L =", round(L_standardized, 3), 
      ") indicates",
      ifelse(L_standardized > 0.5, "strong", "moderate"), "trend.\n")
} else {
  cat("No significant trend(p =", 
      format.pval(page_result$p.value, digits = 3), ")\n")
}