Moses Test of Extreme Reactions

# ==============================================================================
# Moses Test: ADHD Medication Side Effect Variability
# Research Question: Does medication produce more extreme/variable reactions?
# ==============================================================================

library(tidyverse)      # Data manipulation and visualization
library(DescTools)      # MosesTest() function
library(car)            # Levene's test for comparison
library(effectsize)     # Effect size calculations

set.seed(2025)

# Simulate data:
# Placebo: Low variability (homogeneous minor side effects)
# Medication: High variability (some severe, some minimal side effects)
data <- data.frame(
  group = rep(c("Placebo", "Medication"), each=25),
  side_effect_severity = c(
    # Placebo: narrow spread around 2-3
    pmin(10, pmax(0, rnorm(25, mean=2.5, sd=1.2))),
    # Medication: wide spread (bimodal: responders + severe reactors)
    pmin(10, pmax(0, c(
      rnorm(12, mean=1.5, sd=1.5),  # low side effects
      rnorm(13, mean=6.5, sd=2.0)   # high side effects
    )))
  )
)

data$group <- factor(data$group, levels = c("Placebo", "Medication"))

# ==============================================================================
# STEP 1: DESCRIPTIVE STATISTICS
# ==============================================================================

cat("\n=== DESCRIPTIVE STATISTICS ===\n")

descriptives <- data %>%
  group_by(group) %>%
  summarise(
    n = n(),
    mean = mean(side_effect_severity),
    sd = sd(side_effect_severity),
    median = median(side_effect_severity),
    IQR = IQR(side_effect_severity),
    min = min(side_effect_severity),
    max = max(side_effect_severity),
    range = max - min,
    cv = (sd / mean) * 100  # coefficient of variation
  )
print(descriptives)

cat("\nKey observation: Check if SDs/IQRs differ(scale) while medians similar(location).\n")

# ==============================================================================
# STEP 2: ASSUMPTION CHECKING
# ==============================================================================

cat("\n=== ASSUMPTION CHECKS ===\n")

# Assumption 1: Independence
cat("1. Independence: Between-subjects RCT design confirmed(independent groups).\n")

# Assumption 2: Adequate sample size
cat("\n2. Sample Size Check:\n")
n_placebo <- sum(data$group == "Placebo")
n_medication <- sum(data$group == "Medication")
cat(sprintf("   Placebo(control): n = %d(≥10 ✓)\n", n_placebo))
cat(sprintf("   Medication(experimental): n = %d(≥10 ✓)\n", n_medication))
cat("   Adequate sample sizes for Moses test.\n")

# Assumption 3: Check if medians are similar (critical for Moses test interpretation)
cat("\n3. Location Similarity Check(for clear scale interpretation):\n")
placebo_vals <- data %>% filter(group == "Placebo") %>% pull(side_effect_severity)
medication_vals <- data %>% filter(group == "Medication") %>% pull(side_effect_severity)

mdn_placebo <- median(placebo_vals)
mdn_medication <- median(medication_vals)
mdn_diff <- abs(mdn_placebo - mdn_medication)
pooled_iqr <- (IQR(placebo_vals) + IQR(medication_vals)) / 2

cat(sprintf("   Placebo median: %.2f\n", mdn_placebo))
cat(sprintf("   Medication median: %.2f\n", mdn_medication))
cat(sprintf("   Median difference: %.2f\n", mdn_diff))
cat(sprintf("   Relative difference: %.2f (< 0.5 suggests similar location)\n", 
    mdn_diff / pooled_iqr))

# Mann-Whitney to formally test location
mw_result <- wilcox.test(medication_vals, placebo_vals)
cat(sprintf("\n   Mann-Whitney U test: W = %.1f, p = %.4f\n", 
    mw_result$statistic, mw_result$p.value))
if(mw_result$p.value > 0.05) {
  cat("   Locations are similar(p > .05) → Moses test isolates SCALE difference. ✓\n")
} else {
  cat("   Locations differ(p < .05) → Moses result may reflect both location AND scale.\n")
}

# ==============================================================================
# STEP 3: VISUALIZATIONS
# ==============================================================================

cat("\n=== GENERATING VISUALIZATIONS ===\n")

# Visualization 1: Boxplots (focus on spread)
p1 <- ggplot(data, aes(x=group, y=side_effect_severity, fill=group)) +
  geom_boxplot(alpha=0.7, outlier.shape=NA) +
  geom_jitter(width=0.2, alpha=0.5, size=2.5) +
  stat_summary(fun=median, geom="point", shape=23, size=4, fill="red") +
  labs(title="Side Effect Severity: Placebo vs. Medication",
       subtitle="Moses test: Is medication group more variable(wider box/whiskers)?",
       x="Group", y="Side Effect Severity(0-10)") +
  theme_minimal() +
  theme(legend.position="none")
print(p1)

# Visualization 2: Violin plots (show full distribution shape)
p2 <- ggplot(data, aes(x=group, y=side_effect_severity, fill=group)) +
  geom_violin(alpha=0.6, trim=FALSE) +
  geom_boxplot(width=0.1, fill="white", alpha=0.8) +
  labs(title="Distribution Shapes: Medication Shows Greater Spread",
       subtitle="Wider violin = more variability(extreme reactions)",
       x="Group", y="Side Effect Severity") +
  theme_minimal() +
  theme(legend.position="none")
print(p2)

# Visualization 3: Q-Q plot (compare quantiles)
placebo_quantiles <- quantile(placebo_vals, probs = seq(0, 1, 0.1))
medication_quantiles <- quantile(medication_vals, probs = seq(0, 1, 0.1))

qq_data <- data.frame(
  placebo_q = placebo_quantiles,
  medication_q = medication_quantiles
)

p3 <- ggplot(qq_data, aes(x=placebo_q, y=medication_q)) +
  geom_point(size=3, color="blue") +
  geom_abline(slope=1, intercept=0, linetype="dashed", color="red") +
  labs(title="Q-Q Plot: Medication vs. Placebo",
       subtitle="Points above line = medication more variable at that quantile",
       x="Placebo Quantiles", y="Medication Quantiles") +
  theme_minimal()
print(p3)

# Visualization 4: Rank plot showing control span (Moses test logic)
data_ranked <- data %>%
  arrange(side_effect_severity) %>%
  mutate(rank = row_number())

# Identify control group span
control_ranks <- data_ranked %>% 
  filter(group == "Placebo") %>% 
  pull(rank)
control_span <- max(control_ranks) - min(control_ranks) + 1

p4 <- ggplot(data_ranked, aes(x=rank, y=side_effect_severity, 
                               color=group, shape=group)) +
  geom_point(size=3, alpha=0.7) +
  geom_segment(aes(x=min(control_ranks), xend=max(control_ranks),
                   y=-0.5, yend=-0.5), color="darkgreen", size=2) +
  annotate("text", x=mean(control_ranks), y=-1, 
           label=sprintf("Control span = %d", control_span), 
           color="darkgreen", fontface="bold") +
  labs(title="Moses Test Logic: Control Group Span in Combined Ranking",
       subtitle="Short span = control clustered(less variable); Medication more spread out",
       x="Rank(Combined Data)", y="Side Effect Severity") +
  theme_minimal()
print(p4)

# ==============================================================================
# STEP 4: MOSES TEST OF EXTREME REACTIONS
# ==============================================================================

cat("\n=== MOSES TEST OF EXTREME REACTIONS ===\n")

# Moses test: tests if experimental group (medication) has greater variability
# Uses DescTools package
moses_result <- MosesTest(
  x = medication_vals,  # experimental group
  y = placebo_vals,     # control group
  alternative = "greater"  # one-sided: experimental > control variability
)

print(moses_result)

cat("\n--- Test Summary ---\n")
cat(sprintf("Control group(Placebo): n = %d\n", n_placebo))
cat(sprintf("Experimental group(Medication): n = %d\n", n_medication))
cat(sprintf("Control span: %d\n", moses_result$statistic))
cat(sprintf("p-value: %.4f\n", moses_result$p.value))
cat(sprintf("Result: %s at α = .05\n",
    ifelse(moses_result$p.value < 0.05, "SIGNIFICANT", "Not significant")))

if(moses_result$p.value < 0.05) {
  cat("\nInterpretation: Medication group has significantly greater variability\n")
  cat("(more extreme reactions) than placebo. This indicates responder heterogeneity:\n")
  cat("some patients experience severe side effects, others minimal.\n")
}

# ==============================================================================
# STEP 5: EFFECT SIZE
# ==============================================================================

cat("\n=== EFFECT SIZE ===\n")

# Variance ratio (primary effect size for scale tests)
var_medication <- var(medication_vals)
var_placebo <- var(placebo_vals)
variance_ratio <- var_medication / var_placebo

cat(sprintf("Variance ratio(Medication / Placebo): %.3f\n", variance_ratio))
cat("Interpretation: ratio = 1.0 (equal variability), >1.5 (substantially more variable)\n")

# IQR ratio
iqr_ratio <- IQR(medication_vals) / IQR(placebo_vals)
cat(sprintf("\nIQR ratio(Medication / Placebo): %.3f\n", iqr_ratio))

# Coefficient of variation ratio
cv_placebo <- (sd(placebo_vals) / mean(placebo_vals)) * 100
cv_medication <- (sd(medication_vals) / mean(medication_vals)) * 100
cv_ratio <- cv_medication / cv_placebo

cat(sprintf("\nCoefficient of Variation:\n"))
cat(sprintf("   Placebo: %.1f%%\n", cv_placebo))
cat(sprintf("   Medication: %.1f%%\n", cv_medication))
cat(sprintf("   CV ratio: %.2f\n", cv_ratio))

# ==============================================================================
# STEP 6: COMPARE WITH LEVENE'S TEST
# ==============================================================================

cat("\n=== COMPARISON: LEVENE'S TEST FOR VARIANCE EQUALITY ===\n")
cat("Levene's test(parametric) vs. Moses test(nonparametric)\n\n")

levene_result <- leveneTest(side_effect_severity ~ group, data=data)
print(levene_result)

cat("\n--- Comparison Summary ---\n")
cat(sprintf("Moses test: p = %.4f (%s)\n", 
    moses_result$p.value,
    ifelse(moses_result$p.value < 0.05, "significant", "not significant")))
cat(sprintf("Levene's test: p = %.4f (%s)\n",
    levene_result$`Pr(>F)`[1],
    ifelse(levene_result$`Pr(>F)`[1] < 0.05, "significant", "not significant")))

if(moses_result$p.value < 0.05 && levene_result$`Pr(>F)`[1] < 0.05) {
  cat("\nBoth tests agree: variability differs significantly.\n")
  cat("Convergent evidence strengthens conclusion.\n")
}

# ==============================================================================
# STEP 7: INTERPRETATION
# ==============================================================================

cat("\n=== INTERPRETATION ===\n")

if(moses_result$p.value < 0.05) {
  cat("\nThe Moses test showed significantly greater variability in the medication\n")
  cat(sprintf("group compared to placebo, p = %.4f. The medication group(SD=%.2f,\n",
      moses_result$p.value, sd(medication_vals)))
  cat(sprintf("IQR=%.2f, range=%.1f-%.1f) exhibited more extreme reactions than placebo\n",
      IQR(medication_vals), min(medication_vals), max(medication_vals)))
  cat(sprintf("(SD=%.2f, IQR=%.2f, range=%.1f-%.1f).\n",
      sd(placebo_vals), IQR(placebo_vals), min(placebo_vals), max(placebo_vals)))
  cat(sprintf("\nVariance ratio: %.2f (medication %.0fx more variable than placebo).\n",
      variance_ratio, variance_ratio))
  cat("\nClinical significance: Responder heterogeneity suggests some patients are\n")
  cat("highly susceptible to side effects while others tolerate medication well.\n")
  cat("Consider investigating predictors of adverse reactions for personalized dosing.\n")
} else {
  cat("\nThe Moses test showed no significant difference in variability between\n")
  cat("groups. Both medication and placebo produced similar variability in side effects.\n")
}

# ==============================================================================
# APA REPORTING TEMPLATE
# ==============================================================================

cat("\n=== APA-STYLE REPORTING ===\n")
cat("A Moses test of extreme reactions was conducted to compare side effect\n")
cat("variability between placebo(control, n=25) and ADHD medication(experimental,\n")
cat("n=25). Moses test evaluates whether the experimental group exhibits greater\n")
cat("variability(more extreme scores) than control. Groups had similar medians\n")
cat(sprintf("(Placebo: Mdn=%.1f, Medication: Mdn=%.1f, Mann-Whitney p=%.3f), isolating\n",
    mdn_placebo, mdn_medication, mw_result$p.value))
cat(sprintf("scale from location effects. Results showed %s, p = %.4f.\n",
    ifelse(moses_result$p.value < 0.05, 
           "significantly greater variability in medication group",
           "no significant variability difference"),
    moses_result$p.value))

if(moses_result$p.value < 0.05) {
  cat(sprintf("\nThe medication group(SD=%.2f, IQR=%.2f) was %.1fx more variable than\n",
      sd(medication_vals), IQR(medication_vals), variance_ratio))
  cat(sprintf("placebo(SD=%.2f, IQR=%.2f), as confirmed by Levene's test(F=%.2f,\n",
      sd(placebo_vals), IQR(placebo_vals), levene_result$`F value`[1]))
  cat(sprintf("p<.001). This responder heterogeneity—with some patients experiencing\n"))
  cat("severe side effects and others minimal reactions—suggests individual\n")
  cat("differences in drug metabolism or sensitivity. Clinical implication: identify\n")
  cat("predictive biomarkers for personalized treatment planning.\n")
}