CRD Advanced: Contrasts and Multiple Comparisons

# Load radon data from GitHub
radon_url <- "https://raw.githubusercontent.com/Jaisai0611/CRD-example/main/EX%203.21.csv"
radon_data <- read.csv(radon_url) |>
    pivot_longer(
        cols = starts_with("Radon"),
        names_to = "sample",
        values_to = "radon_pct"
    ) |>
    mutate(
        orifice = factor(Orifice.Diameter),
        response = radon_pct
    ) |>
    dplyr::select(orifice, sample, response)

# Summary statistics by orifice size
radon_data |>
    group_by(orifice) |>
    summarize(
        n = n(),
        mean = mean(response),
        sd = sd(response),
        .groups = "drop"
    ) |>
    knitr::kable(digits = 1)

ggplot(radon_data, aes(x = orifice, y = response, fill = orifice)) +
    geom_boxplot(alpha = 0.7) +
    geom_jitter(width = 0.15, alpha = 0.5, size = 2.5) +
    labs(
        title = "Effect of Orifice Diameter on Radon Release",
        x = "Orifice Diameter (mm)",
        y = "Radon Released (%)",
        fill = "Diameter"
    ) +
    scale_fill_viridis_d(option = "plasma") +
    theme_minimal() +
    theme(
        text = element_text(size = 14),
        legend.position = "none"
    )

# Check ANOVA assumptions
par(mfrow = c(2, 2), mar = c(4, 4, 2, 2))
plot(model_radon)
par(mfrow = c(1, 1))

# Plot EMMs with CIs
plot(emm_radon, comparisons = TRUE) +
    labs(
        title = "Radon Release by Orifice Diameter with 95% CIs",
        x = "Radon Released (%)",
        y = "Orifice Diameter (mm)"
    ) +
    theme_minimal() +
    theme(text = element_text(size = 13))

# Create simulated but realistic NPK application data
set.seed(2026)
npk_data <- data.frame(
    treatment = factor(rep(c("Control", "Jan_Plow", "Jan_Broadcast", "Apr_Broadcast"), each = 8),
        levels = c("Control", "Jan_Plow", "Jan_Broadcast", "Apr_Broadcast")
    ),
    block = factor(rep(1:8, 4)),
    yield = c(
        rnorm(8, 45, 5), # Control: lower yield
        rnorm(8, 62, 5), # Jan plowed: best
        rnorm(8, 58, 5), # Jan broadcast: good
        rnorm(8, 55, 6) # Apr broadcast: moderate
    )
)

# Summary by treatment
npk_data |>
    group_by(treatment) |>
    summarize(
        n = n(),
        mean = mean(yield),
        sd = sd(yield),
        .groups = "drop"
    ) |>
    knitr::kable(digits = 1)

ggplot(npk_data, aes(x = treatment, y = yield, fill = treatment)) +
    geom_boxplot(alpha = 0.7) +
    geom_jitter(width = 0.15, alpha = 0.5, size = 2.5) +
    labs(
        title = "Effect of NPK Application Method on Crop Yield",
        x = "Treatment",
        y = "Yield (bushels/acre)",
        fill = "Treatment"
    ) +
    scale_fill_viridis_d(option = "plasma") +
    theme_minimal() +
    theme(
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 30, hjust = 1),
        legend.position = "none"
    )

# Fit CRD model
model_npk <- lm(yield ~ treatment, data = npk_data)

# Get estimated marginal means
emm_npk <- emmeans(model_npk, ~treatment)

# Define contrasts as a list
contrast_list <- list(
    "Control vs NPK" = c(3, -1, -1, -1) / 3, # Divide to get average
    "Jan vs Apr" = c(0, 1, 1, -2) / 2, # Compare averages
    "Plow vs Broadcast" = c(0, 1, -1, 0) # Simple comparison
)

# Test contrasts
contrast_results <- contrast(emm_npk, contrast_list)
print(contrast_results)

# With confidence intervals
confint(contrast_results)

# Plot contrasts with confidence intervals
contrast_df <- confint(contrast_results) |> as.data.frame()
if (!all(c("lower.CL", "upper.CL") %in% names(contrast_df))) {
    contrast_df <- contrast_df |>
        mutate(
            lower.CL = estimate - 1.96 * SE,
            upper.CL = estimate + 1.96 * SE
        )
}

contrast_df |>
    ggplot(aes(x = contrast, y = estimate)) +
    geom_point(size = 4, color = "steelblue") +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
        width = 0.2, linewidth = 1.2, color = "steelblue"
    ) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 0.8) +
    coord_flip() +
    labs(
        title = "Planned Contrasts: NPK Application Methods",
        x = "Contrast",
        y = "Effect Size (bushels/acre)",
        caption = "Error bars show 95% confidence intervals. Red line at zero = no effect."
    ) +
    theme_minimal() +
    theme(text = element_text(size = 13))

# Dunnett's test: all treatments vs Control
# Control must be the reference level (first level of factor)
dunnett_results <- contrast(emm_npk, method = "trt.vs.ctrl", adjust = "dunnett")
dunnett_results |>
    as.data.frame() |>
    knitr::kable(digits = 2)

# With confidence intervals
confint(dunnett_results) |>
    as.data.frame() |>
    knitr::kable(digits = 2)

# Complete workflow: emmeans_test() - pipe-friendly!
rstatix_results <- npk_data |>
    emmeans_test(yield ~ treatment, p.adjust.method = "tukey")

# Display key columns
rstatix_results |>
    select(term, .y., group1, group2, statistic, p, p.adj, p.adj.signif) |>
    knitr::kable(digits = 3)

# rstatix provides a simple tukey_hsd() function
npk_tukey <- npk_data |>
    tukey_hsd(yield ~ treatment)

# Display pairwise comparison results
npk_tukey |>
    select(group1, group2, estimate, conf.low, conf.high, p.adj, p.adj.signif) |>
    knitr::kable(digits = 2)

# Step 1: Run statistical test using rstatix
pwc <- npk_data |>
    emmeans_test(yield ~ treatment, p.adjust.method = "tukey")

# Step 2: Prepare positions for brackets
pwc <- pwc |> add_xy_position(x = "treatment")

# Step 3: Create plot with automated significance brackets
ggplot(npk_data, aes(x = treatment, y = yield, fill = treatment)) +
    geom_boxplot(alpha = 0.7, outlier.shape = NA) +
    geom_jitter(width = 0.15, alpha = 0.5, size = 2) +
    # Step 4: Add significance brackets automatically
    stat_pvalue_manual(pwc,
        hide.ns = TRUE, label = "p.adj.signif",
        tip.length = 0.01, step.increase = 0.05
    ) +
    labs(
        title = "NPK Treatment Effects with Automated Significance Testing",
        subtitle = "Brackets show significant pairwise differences (Tukey HSD, α = 0.05)",
        x = "Treatment",
        y = "Yield (bushels/acre)",
        fill = "Treatment"
    ) +
    scale_fill_viridis_d(option = "plasma") +
    theme_minimal() +
    theme(
        text = element_text(size = 13),
        axis.text.x = element_text(angle = 30, hjust = 1),
        legend.position = "none"
    )

# Compare different adjustment methods
methods_comparison <- list(
    Tukey = as.data.frame(pairs(emm_npk, adjust = "tukey")),
    Bonferroni = as.data.frame(pairs(emm_npk, adjust = "bonferroni")),
    None = as.data.frame(pairs(emm_npk, adjust = "none")) # Fisher's LSD
) |>
    bind_rows(.id = "method")

# Plot p-values
methods_comparison |>
    as.data.frame() |>
    ggplot(aes(x = contrast, y = p.value, color = method, group = method)) +
    geom_point(size = 3, position = position_dodge(width = 0.5)) +
    geom_hline(yintercept = 0.05, linetype = "dashed", color = "red") +
    coord_flip() +
    labs(
        title = "P-values Across Multiple Comparison Methods",
        subtitle = "Red line shows α = 0.05 threshold",
        x = "Pairwise Comparison",
        y = "P-value",
        color = "Adjustment"
    ) +
    scale_color_viridis_d(option = "plasma", end = 0.9) +
    theme_minimal() +
    theme(text = element_text(size = 12))

# Install easystats if needed: install.packages("easystats")
library(performance)
library(parameters)
library(effectsize)
library(report)

# Check model assumptions with performance::check_model()
check_model(model_npk, check = c("normality", "homogeneity", "qq", "outliers"))

# Extract model parameters with confidence intervals and standardization
model_parameters(model_npk, effects = "fixed", ci = 0.95)

# Calculate effect sizes - proportion of variance explained by treatment
# η² (eta-squared): Proportion of total variance explained (biased, overestimates)
# ω² (omega-squared): Less biased estimate, adjusts for sample size
# ε² (epsilon-squared): Alternative unbiased estimate
# All range from 0 (no effect) to 1 (treatment explains all variance)
eta_squared(model_npk, partial = TRUE, ci = 0.95)

omega_squared(model_npk, partial = TRUE, ci = 0.95)