SMHS rANOVA - Revision history

Dinov: /* References */

2025-12-09T17:35:22Z

‎References

Dinov: /* Problems and Exercises */

2025-12-09T17:34:48Z

‎Problems and Exercises

Dinov: /* Violation of Sphericity */

2025-12-09T17:33:51Z

‎Violation of Sphericity

Dinov: /* Violation of Sphericity */

2025-12-09T17:32:44Z

‎Violation of Sphericity

Dinov: /* Violation of Sphericity */

2025-12-09T17:18:45Z

‎Violation of Sphericity

Dinov: /* Software Implementation */

2025-12-09T17:03:44Z

‎Software Implementation

Dinov: /* Handling Autocorrelated Errors */

2025-12-09T16:57:15Z

‎Handling Autocorrelated Errors

Dinov: /* Handling Autocorrelated Errors */

2025-12-09T16:55:39Z

‎Handling Autocorrelated Errors

Dinov: /* Handling Autocorrelated Errors */

2025-12-09T16:54:03Z

‎Handling Autocorrelated Errors

Dinov: /* Example 2: Clinical Trial with Missing Data (Using Mixed Models) */

2025-12-09T16:42:02Z

‎Example 2: Clinical Trial with Missing Data (Using Mixed Models)

@@ Line 1,040: / Line 1,040: @@
 . Maxwell, S. E., & Delaney, H. D. (2004). *Designing Experiments and Analyzing Data: A Model Comparison Perspective* (2nd ed.). Psychology Press.
 . Pinheiro, J. C., & Bates, D. M. (2000). *Mixed-Effects Models in S and S-PLUS*. Springer.
 . Verbeke, G., & Molenberghs, G. (2000). *Linear Mixed Models for Longitudinal Data*. Springer.
 . Field, A., Miles, J., & Field, Z. (2012). *Discovering Statistics Using R*. Sage Publications.
 . Gueorguieva, R., & Krystal, J. H. (2004). Move over ANOVA: Progress in analyzing repeated-measures data and its reflection in papers published in the Archives of General Psychiatry. *Archives of General Psychiatry, 61*(3), 310-317.

@@ Line 967: / Line 967: @@
 . '''Applied Problems''':
-    a) Analyze the [SOCR LDT dataset](http://wiki.socr.umich.edu/index.php/SOCR_Data_072108_LDT) using rANOVA
+    a) Analyze the [http://wiki.socr.umich.edu/index.php/SOCR_Data_072108_LDT SOCR LDT dataset] using rANOVA
-    b) Conduct a repeated measures analysis of the [ChickWeight dataset](https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/ChickWeight.html) with diet as a between-subjects factor
+    b) Conduct a repeated measures analysis of the [https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/ChickWeight.html ChickWeight dataset] with diet as a between-subjects factor
     c) Perform power analysis for a planned longitudinal study with 5 time points, expecting a medium effect (f² = 0.15), and correlation ρ = 0.6

@@ Line 923: / Line 923: @@
    # 3. Use linear mixed models with unstructured covariance
    cat("\n3. Linear mixed model with unstructured covariance:\n")
-   library(nlme)
+   # library(nlme) # Already loaded globally
    # Convert grouping variable to factor (required for lme grouping)
@@ Line 935: / Line 935: @@
    # Fit model with unstructured covariance (dynamic formula)
    cor_form <- as.formula(paste("~", time_var, "|", subject_var))
    model_unstructured <- lme(as.formula(paste(response_var, "~", time_var)),
-                             random = ~ 1 | subject_var,
+                             random = random_form, # Use the dynamically created formula
                              correlation = corSymm(form = cor_form),
                              data = data,
@@ Line 944: / Line 948: @@
    # Compare with compound symmetry model
    model_compound <- lme(as.formula(paste(response_var, "~", time_var)),
-                         random = ~ 1 | subject_var,
+                         random = random_form, # Also use the dynamically created formula
-                         correlation = corCompSymm(form = ~ 1 | subject_var),
+                         correlation = corCompSymm(form = as.formula(paste("~ 1 |", subject_var))),
                          data = data)
    cat("\nModel comparison (unstructured vs. compound symmetry):\n")

@@ Line 894: / Line 894: @@
+)
-# Define the function (with fixes for MANOVA and correlation formula)
+# Define the function (with fixes for grouping variable and time conversion)
 handle_sphericity_violation <- function(data, subject_var, time_var, response_var) {
    cat("=== Handling Sphericity Violation ===\n")
@@ Line 924: / Line 924: @@
    cat("\n3. Linear mixed model with unstructured covariance:\n")
    library(nlme)
    # Fit model with unstructured covariance (dynamic formula)
-   cor_form <- as.formula(paste("~ as.numeric(", time_var, ") | ", subject_var))
+   cor_form <- as.formula(paste("~", time_var, "|", subject_var))
    model_unstructured <- lme(as.formula(paste(response_var, "~", time_var)),
                              random = ~ 1 | subject_var,

@@ Line 824: / Line 824: @@
 ====Violation of Sphericity====
 <pre>
 handle_sphericity_violation <- function(data, subject_var, time_var, response_var) {
@@ Line 874: / Line 877: @@
    print(anova(model_unstructured, model_compound))
+}
 </pre>

@@ Line 639: / Line 639: @@
 <pre>
 # Comprehensive rANOVA analysis function
-run_ranova_analysis <- function(data, subject_var, time_var, response_var,
+run_ranova_analysis <- function(data, subject_var, time_var, response_var, between_vars = NULL, method = "mixed") {
    cat("=== REPEATED MEASURES ANOVA ANALYSIS ===\n\n")
    cat("Subject variable:", subject_var, "\n")
@@ Line 653: / Line 651: @@
    # Check for missing data
    missing_count <- sum(is.na(data[[response_var]]))
-   cat("Missing observations:", missing_count,
+   cat("Missing observations:", missing_count, sprintf("(%.1f%%)\n", missing_count/nrow(data)*100))
    if (method == "mixed") {
@@ Line 668: / Line 665: @@
        # Add interactions with time
        interaction_str <- paste(paste(time_var, between_vars, sep = ":"), collapse = " + ")
-       formula_str <- paste(response_var, "~", time_var, "*", between_str,
+       formula_str <- paste(response_var, "~", time_var, "*", between_str, "+ (1|", subject_var, ")")
+     }
      cat("Fitting model:", formula_str, "\n\n")
      model <- lmer(as.formula(formula_str), data = data)
@@ Line 695: / Line 690: @@
      # Traditional rANOVA using aov()
      if (!is.null(between_vars)) {
-       formula_str <- paste(response_var, "~",
+       formula_str <- paste(response_var, "~", paste(between_vars, collapse = " + "), "+ Error(", subject_var, "/", time_var, ")")
      } else {
-       formula_str <- paste(response_var, "~", time_var,
+       formula_str <- paste(response_var, "~", time_var, "+ Error(", subject_var, "/", time_var, ")")
+     }
      cat("Fitting model:", formula_str, "\n\n")
      model <- aov(as.formula(formula_str), data = data)
@@ Line 713: / Line 704: @@
    cat("\n--- Post-hoc Comparisons ---\n")
    if (require(emmeans)) {
      if (method == "mixed") {
-       emm_time <- emmeans(model, specs = pairwise ~ time_var)
+       anova_res <- anova(model)
      } else {
-       # For traditional ANOVA, fit a linear model for post-hoc tests
+       anova_res <- summary(model)
-      if (is.null(between_vars)) {
+      time_p <- anova_res[[2]][[1]]$"Pr(>F)"[1]  # Assuming time is the first term in Error strata
-         lm_formula <- as.formula(paste(response_var, "~", subject_var, "+", time_var))
        } else {
-         lm_formula <- as.formula(paste(response_var, "~", subject_var, "+", time_var,
+         # For traditional ANOVA, fit a linear model for post-hoc tests
-                                      "+", paste(between_vars, collapse = " + ")))
+        if (is.null(between_vars)) {
+       }
-       lm_model <- lm(lm_formula, data = data)
+       print(emm_time$contrasts)
-       emm_time <- emmeans(lm_model, specs = pairwise ~ time_var, adjust = "tukey")
+    } else {
+     }
+   }

@@ Line 627: / Line 627: @@
                  random = ~ 1 | subject,
                  correlation = corCompSymm(form = ~ 1 | subject),
-                 data = ar_data)
+                 data = ar_data,
 # Model comparison

@@ Line 571: / Line 571: @@
 ====Handling Autocorrelated Errors====
 <pre>
 # When measurements are equally spaced in time, errors may be autocorrelated
@@ Line 614: / Line 620: @@
 cat("=== Refit Model with AR(1) Correlation Structure and Maximum Likelihood  Estimation ===\n")
 print(summary(model_ar1_ml))
-cat("\nCorrelation parameter:", intervals(model_ar1)$corStruct[2], "\n")
 # Compare with compound symmetry model
@@ Line 624: / Line 631: @@
 # Model comparison
 cat("\n=== Model Comparison ===\n")
-print(anova(model_ar1, model_cs))
+print(anova(model_ar1_ml, model_cs))
 </pre>

@@ Line 603: / Line 603: @@
                   data = ar_data)
-cat("=== Model with AR(1) Correlation Structure ===\n")
+# Refit the model using maximum likelihood (method = "ML") instead of REML.
 print(summary(model_ar1))
 cat("\nCorrelation parameter:", intervals(model_ar1)$corStruct[2], "\n")

@@ Line 434: / Line 434: @@
    # Check if confidence interval columns exist
-   if (!all(c("lower.CL", "upper.CL") %in% names(plot_data))) {
+   if (!all(c("asymp.LCL", "asymp.LCL") %in% names(plot_data))) {
      # If not, get them using summary
      plot_data <- as.data.frame(summary(emm_interaction))
@@ Line 442: / Line 442: @@
    ggplot(plot_data, aes(x = time, y = emmean, group = treatment, color = treatment)) +
      geom_line(size = 1.5) +
-     geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1) +
+     geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1) +
      geom_point(size = 3) +
      labs(title = "Treatment × Time Interaction",