Standard Ancova
Executive Summary

Executive Summary

Does Final Score differ by Study Program once Baseline Score is adjusted for?

Observations
240
Groups Compared
3
Adjusted Gap (best-worst)
7.32
Group Test
p < 0.001
Slopes Check
passed
Covariate Slope
1.94
After adjusting for Baseline Score, Final Score differs between the Study Program groups (F = 15.27, p < 0.001). The adjustment changes the story: on raw means Program A looks highest, but once Baseline Score is accounted for, Program B leads — the groups differ in Baseline Score and part of the raw gap reflects that, not the Study Program themselves. The adjusted gap between Program B and Program A is 7.32 (95% CI 4.64 to 10.00), with 3 of 3 pairwise adjusted difference(s) significant after Holm correction. Within groups, each 1-unit increase in Baseline Score is associated with a 1.940 change in Final Score (95% CI 1.835 to 2.046).
Suggested Interpretation

Program B produces the best adjusted final scores. After accounting for baseline differences, Final Score differs significantly across programs (F = 15.27, p < 0.001). Program B's adjusted mean is 128.795, ahead of Program C (124.458) and Program A (121.475)—a gap of 7.32 between the highest and lowest. Notably, raw means tell a different story: Program A appears strongest at 142.097, but this reflects its students' higher baseline scores. The equal-slopes check passed (p = 0.589), confirming the adjustment is valid.

Overview

Analysis Overview

Analysis of covariance: Final Score across 3 Study Program groups, adjusted for Baseline Score (240 observations).

N Observations240
N Groups3
N Pairwise3
Suggested Interpretation

ANCOVA compares groups when they start from unequal footing. Here, 240 students across three Study Programs are assessed on Final Score, but the programs differ in their students' Baseline Scores. The analysis estimates what each program's mean Final Score would be if all students began at the same baseline level (60.38), removing the covariate's confounding influence. This reveals whether program differences reflect the program itself or simply imbalanced starting points.

Data Preparation

Data Quality

Row and group cleaning applied before fitting the model.

Initial Rows240
Final Rows240
Rows Removed0
Groups Dropped0
Levels Lumped0
Suggested Interpretation

The short answer

All 240 rows of data were usable; no observations were dropped for missing outcome or covariate values, and all three program groups met the minimum size threshold.

The detail

240 rows loaded and 240 rows used. No rows were missing Final Score or Baseline Score values. All groups had at least 3 observations. The final comparison covers 3 groups: Program A, Program B, Program C. No groups were dropped and no rare levels were lumped.

What this can't tell you

Data quality checks do not address measurement validity—whether Baseline Score and Final Score were recorded accurately or whether they measure what the analysis assumes they measure.

Visualization

Raw vs Adjusted Means

Each Study Program group's raw mean of Final Score next to its Baseline Score-adjusted (least-squares) mean.

Suggested Interpretation

Program A moves most under adjustment, dropping 20.62 points from raw to adjusted mean (142.097 to 121.475), revealing that its apparent raw advantage owed heavily to higher baseline scores. Program B, raw mean 109.189, gains ground under adjustment (128.795), rising to first place. Program C's raw mean (123.442) adjusts to 124.458, barely shifting. The ranking flips completely: unadjusted, Program A leads; adjusted, Program B leads. This reversal shows covariate imbalance was the dominant driver of raw differences.

Visualization

Outcome vs Covariate by Group

Final Score against Baseline Score, colored by Study Program.

Suggested Interpretation

The scatter plot reveals why adjustment matters: Program A students average Baseline Score 71.01, while Program B students average 50.28—a substantial gap. Program C sits between them. Within each group, Final Score rises approximately 1.940 points per Baseline Score unit, and these trends are visually parallel across programs. The equal-slopes diagnostic confirmed parallelism (p = 0.589), validating the use of a single shared adjustment factor across all three programs.

Data Table

ANCOVA Table

F-tests for the covariate, the adjusted group effect, and the equal-slopes diagnostic.

TermDfSum SqMean SqF ValueP Value
Baseline Score17.558e+047.558e+043039< 0.001
Study Program (adjusted for Baseline Score)2759.7379.915.27< 0.001
Study Program x Baseline Score interaction (slopes check)226.4613.230.530.589
Residuals236586924.87—
Suggested Interpretation

The short answer

The adjusted group comparison is highly significant (p < 0.001), confirming that program differences in final scores persist after removing the effect of baseline scores. The baseline score itself is a strong predictor (p < 0.001), and the equal-slopes check passed (p = 0.589), validating the adjustment approach.

The detail

Baseline Score explains substantial variance (F = 3038.967, p < 0.001). Study Program, adjusted for Baseline Score, yields F = 15.274, p < 0.001, with partial eta-squared of 0.115, meaning the program effect explains 11.5% of remaining variation. The Study Program × Baseline Score interaction (equal-slopes diagnostic) has F = 0.53, p = 0.589, indicating no significant deviation from parallel slopes across groups.

What this can't tell you

F-tests confirm significance but do not identify which specific programs differ most; pairwise comparisons address that. The interaction test's non-significance does not prove slopes are identical, only that observed differences are consistent with chance.

Data Table

Adjusted Means & Pairwise Differences

Least-squares means of Final Score per Study Program group, plus Holm-corrected pairwise differences.

GroupNRaw MeanAdjusted MeanSECI LowCI High
Program A80142.1121.50.796119.9123
Program B80109.2128.80.776127.3130.3
Program C80123.4124.50.558123.4125.6
ComparisonAdj DifferenceCI LowCI HighP ValueHolm PSignificant
Program A - Program B-7.32-10-4.64< 0.001< 0.001yes
Program B - Program C4.3372.4856.189< 0.001< 0.001yes
Program A - Program C-2.983-4.93-1.0350.0030.003yes
Suggested Interpretation

Program B's adjusted mean Final Score is 128.795 (95% CI 127.266 to 130.324, n = 80), the highest. Program C follows at 124.458 (95% CI 123.358 to 125.558, n = 80). Program A is lowest at 121.475 (95% CI 119.907 to 123.043, n = 80). All three pairwise differences are significant after Holm correction: Program B exceeds Program A by 7.32 (p < 0.001), Program B exceeds Program C by 4.337 (p < 0.001), and Program A trails Program C by 2.983 (p = 0.003).

Data Table

Method & Assumptions

The model, the adjustment, the correction used, and the assumptions the adjusted comparison rests on.

ItemDetail
DesignObservational analysis of covariance — adjustment for one measured covariate; other confounders remain unadjusted.
ModelLinear model: Final Score ~ Baseline Score + Study Program, fit by ordinary least squares on 240 rows.
Adjusted meansLeast-squares means: the model's prediction for each Study Program group at the grand mean of Baseline Score (60.38), with model-based standard errors.
Omnibus group testF-test of the Study Program term after Baseline Score is entered (F = 15.27 on 2 and 236 df, p < 0.001).
Pairwise correctionPairwise adjusted differences are t-tests on model contrasts with Holm-corrected p-values; confidence intervals are unadjusted.
Covariate effectWithin groups, each 1-unit increase in Baseline Score is associated with a Final Score change of 1.940 (95% CI 1.835 to 2.046, p < 0.001).
Equal-slopes checkThe equal-slopes check passed (interaction p = 0.589): the Final Score-Baseline Score relationship looks similar across Study Program groups, so a single adjusted comparison is reasonable.
Causal readingAdjusted differences are associations consistent with a Study Program effect, not proof of cause — groups were not randomized here and only Baseline Score was adjusted for.
Suggested Interpretation

This is an observational ANCOVA: groups were not randomized, and only Baseline Score was adjusted for; other potential confounders remain unadjusted. Adjusted means are least-squares predictions at Baseline Score = 60.38. Pairwise tests use Holm correction; confidence intervals are unadjusted. The covariate effect is 1.940 (95% CI 1.835 to 2.046) per unit. Adjusted differences are associations consistent with a program effect, not proof of causation.

Methodology

Methodology

Statistical methodology and diagnostics for ANCOVA — Adjusted Group Comparison

Statistical Method

ANCOVA — Adjusted Group Comparison

Standard-library analysis: do your groups differ once you adjust for a covariate? Analysis of covariance fits outcome ~ group + covariate, reports raw versus adjusted (least-squares) group means side by side, the omnibus adjusted group test, pairwise adjusted differences with multiplicity correction, the covariate's slope with a confidence interval, and a homogeneity-of-slopes diagnostic that warns you when the ANCOVA assumption itself is questionable.

Data
N = 240 observations
Assumptions
  • The outcome and covariate are numeric; the group column is categorical
  • The outcome-covariate relationship is approximately linear within each group
  • The covariate's slope is the same in every group (checked by the built-in homogeneity-of-slopes diagnostic)
  • Observations are independent; residuals are roughly normal with similar spread across groups
Limitations
  • ANCOVA adjusts for one measured covariate — other confounders remain unadjusted, so group differences are associations, not proof of cause
  • If the slopes diagnostic fails, adjusted means depend on the covariate value chosen and a single adjusted comparison can mislead
  • Groups with fewer than 3 observations are excluded; group columns with many levels are lumped beyond the 8 largest into Other
  • Adjustment assumes the covariate is measured before or independently of the group effect — adjusting for a post-treatment variable can bias the comparison
Software & Citation
MCP Analytics · mcpanalytics.ai
Code Appendix

Analysis Code

Complete R source code for this analysis

ANCOVA — Adjusted Group Comparison

Does a numeric outcome differ between groups once a numeric covariate is adjusted for? Fits outcome ~ group + covariate, reports raw versus adjusted (least-squares) means side by side, the omnibus adjusted group F-test, Holm-corrected pairwise adjusted differences, the covariate's slope with a confidence interval, and a homogeneity-of-slopes diagnostic.

Why This Method?

Raw group means compare apples to oranges when the groups sit at different covariate levels. ANCOVA removes the covariate's contribution and compares the groups at a common covariate value — the fair comparison — while an equal-slopes check warns when that adjustment itself is on shaky ground.

What This Analysis Covers

  • Raw vs adjusted (least-squares) means per group, side by side
  • The outcome plotted against the covariate, colored by group
  • The ANCOVA F-table including the equal-slopes diagnostic
  • Pairwise adjusted differences with Holm correction

Standard Library

Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {outcome, group, covariate}. All narrative is derived from the user's own column names and computed values.

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

Core Analysis Pipeline

compute_shared <- function(df, params, col_map = list()) {
  # === SHARED EXPORTS ===
  #   initial_rows/final_rows/rows_removed  $ row accounting
  #   n_na_rows            $ rows dropped for missing outcome/covariate
  #   outcome_h / group_h / covariate_h  $ humanized user names
  #   k / group_levels     $ number of groups after cleaning + their names
  #   dropped_groups_df    $ data.frame(group, n) — groups dropped (n < 3)
  #   lumped_levels        $ character — levels folded into "Other"
  #   adjusted_means_df    $ group, n, raw_mean, adjusted_mean, se, ci_low, ci_high
  #   pairwise_df          $ comparison, adj_difference, ci_low, ci_high, p_value, holm_p, significant
  #   means_comparison_df  $ group_name, mean_type, mean_value (grouped bar)
  #   scatter_df           $ covariate_value, outcome_value, group_label (<=1500)
  #   ancova_df            $ term, df, sum_sq, mean_sq, f_value, p_value (display strings for p)
  #   methods_df           $ item, detail
  #   group_F / group_p    $ omnibus adjusted group test
  #   slope / slope_ci / slope_p  $ covariate effect
  #   slopes_F / slopes_p / slopes_ok  $ homogeneity-of-slopes diagnostic
  #   raw_best/raw_worst/adj_best/adj_worst  $ extreme groups by each mean
  #   ordering_changed     $ TRUE if raw and adjusted rankings differ
  #   adj_gap / adj_gap_ci $ adjusted difference best-vs-worst + its 95% CI
  #   metrics / json_output
  # === /SHARED EXPORTS ===

Honest p-value formatting: never e-notation, never bare stars.

p_prose <- function(p) {
    if (is.na(p)) return("p not computable")
    if (p < 0.001) return("p < 0.001")
    paste0("p = ", sprintf("%.3f", p))
  }
  p_cell <- function(p) {
    ifelse(is.na(p), "",
           ifelse(p < 0.001, "< 0.001", sprintf("%.3f", p)))
  }
  fmt_n <- function(v, d = 2) formatC(round(v, d), format = "f", digits = d, big.mark = "")

Step 1: Resolve mapped columns (humanized for all prose)

initial_rows <- nrow(df)
  outcome_h   <- humanize_semantic("outcome", col_map)
  group_h     <- humanize_semantic("group", col_map)
  covariate_h <- humanize_semantic("covariate", col_map)
  missing_keys <- setdiff(c("outcome", "group", "covariate"), names(df))
  if (length(missing_keys) > 0) {
    stop(sprintf("ANCOVA needs &#x27;%s' (numeric outcome), '%s' (the groups), and '%s' (numeric covariate) all mapped.",
                 outcome_h, group_h, covariate_h))
  }

Step 2: Coerce outcome + covariate to numeric (95% rule)

coerce_num <- function(v, label) {
    if (is.numeric(v)) return(v)
    ch <- as.character(v)
    non_blank <- !is.na(ch) & trimws(ch) != ""
    conv <- suppressWarnings(as.numeric(ch))
    if (sum(non_blank) == 0 ||
        sum(!is.na(conv[non_blank])) < 0.95 * sum(non_blank)) {
      stop(sprintf("The column &#x27;%s' does not look numeric — fewer than 95%% of its values parse as numbers. ANCOVA needs a numeric column here.",
                   label))
    }
    conv
  }
  df$outcome   <- coerce_num(df$outcome, outcome_h)
  df$covariate <- coerce_num(df$covariate, covariate_h)

  g <- as.character(df$group)
  g[is.na(g) | trimws(g) == ""] <- "Missing"

  keep <- !is.na(df$outcome) & !is.na(df$covariate)
  n_na_rows <- sum(!keep)
  df <- df[keep, , drop = FALSE]
  g  <- g[keep]
  if (nrow(df) == 0) {
    stop(sprintf("No rows with usable numeric values in both &#x27;%s' and '%s' remained after cleaning.",
                 outcome_h, covariate_h))
  }

Step 3: Clean the groups — drop n<3 (reported), lump beyond 8 levels

tab <- table(g)
  small <- names(tab)[tab < 3]
  dropped_groups_df <- data.frame(group = character(0), n = integer(0),
                                  stringsAsFactors = FALSE)
  if (length(small) > 0) {
    dropped_groups_df <- data.frame(group = small, n = as.integer(tab[small]),
                                    stringsAsFactors = FALSE)
    sel <- !(g %in% small)
    df <- df[sel, , drop = FALSE]
    g  <- g[sel]
  }

  lumped_levels <- character(0)
  tab <- sort(table(g), decreasing = TRUE)
  if (length(tab) > 8) {
    keep_lv <- names(tab)[1:8]
    lumped_levels <- setdiff(names(tab), keep_lv)
    g[g %in% lumped_levels] <- "Other"
  }

Re-check after lumping ("Other" itself could be tiny)

tab <- table(g)
  small2 <- names(tab)[tab < 3]
  if (length(small2) > 0) {
    dropped_groups_df <- rbind(dropped_groups_df,
                               data.frame(group = small2, n = as.integer(tab[small2]),
                                          stringsAsFactors = FALSE))
    sel <- !(g %in% small2)
    df <- df[sel, , drop = FALSE]
    g  <- g[sel]
  }

  gf <- factor(g)
  k  <- nlevels(gf)
  if (k < 2) {
    stop(sprintf("ANCOVA needs at least 2 groups in &#x27;%s' with 3 or more rows each; only %d usable group(s) remained after cleaning. Check that '%s' really splits the data into groups.",
                 group_h, k, group_h))
  }
  y <- df$outcome
  x <- df$covariate
  final_rows <- length(y)
  rows_removed <- initial_rows - final_rows
  if (final_rows < 10) {
    stop(sprintf("Only %d usable rows remained — at least 10 are needed for an adjusted group comparison.", final_rows))
  }
  if (isTRUE(stats::var(y) == 0)) {
    stop(sprintf("The outcome &#x27;%s' has no variation at all (every value is identical) — there is nothing to compare.", outcome_h))
  }
  if (isTRUE(stats::var(x) == 0)) {
    stop(sprintf("The covariate &#x27;%s' has no variation at all (every value is identical) — there is nothing to adjust for. Use the plain group-comparison tool instead.", covariate_h))
  }
  group_levels <- levels(gf)

Step 4: Fit the ANCOVA model — outcome ~ covariate + group

Sequential anova() with the covariate entered first makes the group row the covariate-ADJUSTED group test (SS for group given covariate).

dat <- data.frame(y = y, x = x, gf = gf)
  fit <- stats::lm(y ~ x + gf, data = dat)
  if (is.na(stats::coef(fit)["x"])) {
    stop(sprintf("The covariate &#x27;%s' is completely determined by the '%s' groups (no within-group variation), so its slope cannot be estimated. ANCOVA cannot separate the two — use the plain group-comparison tool.",
                 covariate_h, group_h))
  }
  df_res <- stats::df.residual(fit)
  if (df_res < 3) {
    stop(sprintf("Too few rows for the number of &#x27;%s' groups — only %d residual degrees of freedom remain. Add data or reduce the number of groups.",
                 group_h, df_res))
  }
  a_tab <- stats::anova(fit)
  cov_F <- a_tab["x", "F value"];  cov_p_seq <- a_tab["x", "Pr(>F)"]
  group_F <- a_tab["gf", "F value"]; group_p <- a_tab["gf", "Pr(>F)"]
  ss_g <- a_tab["gf", "Sum Sq"]; ss_r <- a_tab["Residuals", "Sum Sq"]
  partial_eta2 <- if (!is.na(ss_g) && !is.na(ss_r) && (ss_g + ss_r) > 0)
    ss_g / (ss_g + ss_r) else NA_real_

Covariate effect: model slope with 95% CI (partial, i.e. within groups)

slope <- unname(stats::coef(fit)["x"])
  slope_ci <- suppressWarnings(as.numeric(stats::confint(fit, "x", level = 0.95)))
  slope_p <- summary(fit)$coefficients["x", 4]

Step 5: Homogeneity-of-slopes diagnostic — outcome ~ group * covariate

fit_int <- stats::lm(y ~ x * gf, data = dat)
  cmp <- tryCatch(stats::anova(fit, fit_int), error = function(e) NULL)
  slopes_F <- if (!is.null(cmp)) cmp$F[2] else NA_real_
  slopes_p <- if (!is.null(cmp)) cmp$`Pr(>F)`[2] else NA_real_
  slopes_ok <- is.na(slopes_p) || slopes_p >= 0.05
  slopes_verdict <- if (is.na(slopes_p)) {
    "The equal-slopes check could not be computed; treat the adjusted means with some caution."
  } else if (slopes_ok) {
    sprintf("The equal-slopes check passed(interaction %s): the %s-%s relationship looks similar across %s groups, so a single adjusted comparison is reasonable.",
            p_prose(slopes_p), outcome_h, covariate_h, group_h)
  } else {
    sprintf("WARNING — the equal-slopes check failed(interaction %s): the %s-%s relationship appears to DIFFER between %s groups, so ANCOVA&#x27;s core assumption is questionable and the adjusted means below can mislead. Group differences depend on the %s level; read the scatter chart's group trends instead of relying on a single adjusted number.",
            p_prose(slopes_p), outcome_h, covariate_h, group_h, covariate_h)
  }

Step 6: Raw and adjusted (least-squares) means

Adjusted mean per group = model prediction at the covariate's grand mean, with model-based standard errors.

xbar <- mean(x)
  nd <- data.frame(x = rep(xbar, k), gf = factor(group_levels, levels = group_levels))
  pr <- stats::predict(fit, newdata = nd, se.fit = TRUE)
  tcrit <- stats::qt(0.975, df_res)
  raw_means <- sapply(group_levels, function(l) mean(y[gf == l]))
  grp_n <- sapply(group_levels, function(l) sum(gf == l))
  adj_means <- as.numeric(pr$fit)
  adj_se <- as.numeric(pr$se.fit)
  adjusted_means_df <- data.frame(
    group = group_levels,
    n = as.integer(grp_n),
    raw_mean = round(raw_means, 3),
    adjusted_mean = round(adj_means, 3),
    se = round(adj_se, 3),
    ci_low = round(adj_means - tcrit * adj_se, 3),
    ci_high = round(adj_means + tcrit * adj_se, 3),
    stringsAsFactors = FALSE
  )
  rownames(adjusted_means_df) <- NULL

Raw vs adjusted ordering — the tool's story (NA-safe: means over groups with n>=3 and no NA outcome cannot be NA, but guard anyway)

ok_m <- !is.na(raw_means) & !is.na(adj_means)
  raw_ord <- group_levels[ok_m][order(-raw_means[ok_m])]
  adj_ord <- group_levels[ok_m][order(-adj_means[ok_m])]
  raw_best <- raw_ord[1]; raw_worst <- raw_ord[length(raw_ord)]
  adj_best <- adj_ord[1]; adj_worst <- adj_ord[length(adj_ord)]
  ordering_changed <- !identical(raw_ord, adj_ord)
  adj_gap <- adj_means[match(adj_best, group_levels)] -
             adj_means[match(adj_worst, group_levels)]

Step 7: Pairwise adjusted differences — model contrasts, Holm-corrected

Difference of adjusted means = difference of group coefficients (the covariate term cancels), tested on the model's vcov with residual df. Correction: Holm (p.adjust) — stated explicitly in all prose; a classical Tukey HSD on covariate-adjusted means is not available without extra packages, and Holm is the conservative classical alternative.

b <- stats::coef(fit)
  V <- stats::vcov(fit)
  cvec <- function(l) {
    cv <- stats::setNames(rep(0, length(b)), names(b))
    dn <- paste0("gf", l)
    if (dn %in% names(b)) cv[dn] <- 1
    cv
  }
  pair_rows <- list()
  for (i in seq_len(k - 1)) {
    for (j in (i + 1):k) {
      cw <- cvec(group_levels[i]) - cvec(group_levels[j])
      est <- sum(cw * b)
      se <- sqrt(as.numeric(t(cw) %*% V %*% cw))
      tv <- if (se > 0) est / se else NA_real_
      pv <- if (!is.na(tv)) 2 * stats::pt(-abs(tv), df_res) else NA_real_
      pair_rows[[length(pair_rows) + 1]] <- data.frame(
        comparison = paste0(group_levels[i], " - ", group_levels[j]),
        adj_difference = round(est, 3),
        ci_low = round(est - tcrit * se, 3),
        ci_high = round(est + tcrit * se, 3),
        p_raw = pv,
        stringsAsFactors = FALSE
      )
    }
  }
  pairwise_all <- do.call(rbind, pair_rows)
  pairwise_all$holm_p_num <- stats::p.adjust(pairwise_all$p_raw, method = "holm")
  pairwise_all$p_value <- p_cell(pairwise_all$p_raw)
  pairwise_all$holm_p <- p_cell(pairwise_all$holm_p_num)
  pairwise_all$significant <- ifelse(is.na(pairwise_all$holm_p_num), "",
                              ifelse(pairwise_all$holm_p_num < 0.05, "yes", "no"))
  pairwise_all <- pairwise_all[order(-abs(pairwise_all$adj_difference)), , drop = FALSE]
  rownames(pairwise_all) <- NULL
  n_sig_pairs <- sum(pairwise_all$significant == "yes", na.rm = TRUE)
  pairwise_df <- head(pairwise_all[, c("comparison", "adj_difference", "ci_low",
                                       "ci_high", "p_value", "holm_p",
                                       "significant")], 15)

CI of the best-vs-worst adjusted gap, from its pairwise row

gap_row <- pairwise_all[
    (grepl(adj_best, pairwise_all$comparison, fixed = TRUE) &
       grepl(adj_worst, pairwise_all$comparison, fixed = TRUE)), , drop = FALSE]
  adj_gap_ci <- if (nrow(gap_row) > 0) {
    r1 <- gap_row[1, ]
    if (r1$adj_difference >= 0) c(r1$ci_low, r1$ci_high) else c(-r1$ci_high, -r1$ci_low)
  } else c(NA_real_, NA_real_)
  gap_holm_p <- if (nrow(gap_row) > 0) gap_row$holm_p_num[1] else NA_real_

Step 8: Chart datasets

Grouped bar — raw vs adjusted means, the central story

means_comparison_df <- rbind(
    data.frame(group_name = group_levels, mean_type = "Raw mean",
               mean_value = round(raw_means, 3), stringsAsFactors = FALSE),
    data.frame(group_name = group_levels, mean_type = "Adjusted mean",
               mean_value = round(adj_means, 3), stringsAsFactors = FALSE)
  )
  rownames(means_comparison_df) <- NULL

Scatter — outcome vs covariate colored by group, <= 1500 sampled rows

set.seed(42)
  sidx <- if (final_rows > 1500) sample(final_rows, 1500) else seq_len(final_rows)
  scatter_df <- data.frame(
    covariate_value = x[sidx],
    outcome_value = y[sidx],
    group_label = as.character(gf[sidx]),
    stringsAsFactors = FALSE
  )

ANCOVA table (F-tests; p shown as text so no e-notation ever renders)

ancova_df <- data.frame(
    term = c(covariate_h, paste0(group_h, " (adjusted for ", covariate_h, ")"),
             paste0(group_h, " x ", covariate_h, " interaction(slopes check)"),
             "Residuals"),
    df = c(a_tab["x", "Df"], a_tab["gf", "Df"],
           if (!is.null(cmp)) cmp$Df[2] else NA_integer_,
           a_tab["Residuals", "Df"]),
    sum_sq = round(c(a_tab["x", "Sum Sq"], a_tab["gf", "Sum Sq"],
                     if (!is.null(cmp)) cmp$`Sum of Sq`[2] else NA_real_,
                     a_tab["Residuals", "Sum Sq"]), 3),
    mean_sq = round(c(a_tab["x", "Mean Sq"], a_tab["gf", "Mean Sq"],
                      if (!is.null(cmp) && !is.na(cmp$Df[2]) && cmp$Df[2] > 0)
                        cmp$`Sum of Sq`[2] / cmp$Df[2] else NA_real_,
                      a_tab["Residuals", "Mean Sq"]), 3),
    f_value = round(c(cov_F, group_F, slopes_F, NA_real_), 3),
    p_value = c(p_cell(cov_p_seq), p_cell(group_p), p_cell(slopes_p), ""),
    stringsAsFactors = FALSE
  )

Methods table

methods_df <- data.frame(
    item = c("Design", "Model", "Adjusted means", "Omnibus group test",
             "Pairwise correction", "Covariate effect", "Equal-slopes check",
             "Causal reading"),
    detail = c(
      "Observational analysis of covariance — adjustment for one measured covariate; other confounders remain unadjusted.",
      sprintf("Linear model: %s ~ %s + %s, fit by ordinary least squares on %s rows.",
              outcome_h, covariate_h, group_h, format(final_rows, big.mark = ",")),
      sprintf("Least-squares means: the model&#x27;s prediction for each %s group at the grand mean of %s (%s), with model-based standard errors.",
              group_h, covariate_h, fmt_n(xbar)),
      sprintf("F-test of the %s term after %s is entered(F = %s on %d and %d df, %s).",
              group_h, covariate_h, fmt_n(group_F), a_tab["gf", "Df"], df_res, p_prose(group_p)),
      "Pairwise adjusted differences are t-tests on model contrasts with Holm-corrected p-values; confidence intervals are unadjusted.",
      sprintf("Within groups, each 1-unit increase in %s is associated with a %s change of %s(95%% CI %s to %s, %s).",
              covariate_h, outcome_h, fmt_n(slope, 3), fmt_n(slope_ci[1], 3),
              fmt_n(slope_ci[2], 3), p_prose(slope_p)),
      slopes_verdict,
      sprintf("Adjusted differences are associations consistent with a %s effect, not proof of cause — groups were not randomized here and only %s was adjusted for.",
              group_h, covariate_h)
    ),
    stringsAsFactors = FALSE
  )

  metrics <- list(
    `Observations`            = final_rows,
    `Groups Compared`         = k,
    `Adjusted Gap(best-worst)` = round(adj_gap, 2),
    `Group Test`              = p_prose(group_p),
    `Slopes Check`            = if (slopes_ok) "passed" else "questionable",
    `Covariate Slope`         = round(slope, 3)
  )

  flip_phrase <- if (ordering_changed) {
    paste0("Adjustment changes the ranking: on raw means ", raw_best,
           " looks best, but adjusted for ", covariate_h, ", ", adj_best,
           " leads. ")
  } else {
    paste0("The group ranking is the same before and after adjustment(",
           adj_best, " highest). ")
  }
  json_output <- list(
    answer = paste0(
      "ANCOVA of ", outcome_h, " across ", k, " ", group_h, " groups, adjusting for ",
      covariate_h, " (", format(final_rows, big.mark = ","), " rows): after adjustment the groups ",
      if (!is.na(group_p) && group_p < 0.05) "differ" else "do not differ reliably",
      " (F = ", fmt_n(group_F), ", ", p_prose(group_p), "). ", flip_phrase,
      "The adjusted gap between ", adj_best, " and ", adj_worst, " is ",
      fmt_n(adj_gap), if (!is.na(adj_gap_ci[1])) paste0(
        " (95% CI ", fmt_n(adj_gap_ci[1]), " to ", fmt_n(adj_gap_ci[2]), ")") else "", ". ",
      "Each 1-unit increase in ", covariate_h, " is associated with a ",
      fmt_n(slope, 3), " change in ", outcome_h,
      " (95% CI ", fmt_n(slope_ci[1], 3), " to ", fmt_n(slope_ci[2], 3), "). ",
      if (slopes_ok) "The equal-slopes check passed. "
      else paste0("CAUTION: the equal-slopes check failed(", p_prose(slopes_p),
                  ") — the adjusted comparison is questionable. "),
      n_sig_pairs, " of ", nrow(pairwise_all),
      " pairwise adjusted difference(s) significant after Holm correction."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "means_comparison",
        "covariate_scatter", "ancova_table", "adjusted_means", "methods"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows,
    rows_removed = rows_removed, n_na_rows = n_na_rows,
    outcome_h = outcome_h, group_h = group_h, covariate_h = covariate_h,
    k = k, group_levels = group_levels,
    dropped_groups_df = dropped_groups_df, lumped_levels = lumped_levels,
    adjusted_means_df = adjusted_means_df,
    pairwise_df = pairwise_df, n_sig_pairs = n_sig_pairs,
    n_pairs = nrow(pairwise_all),
    means_comparison_df = means_comparison_df,
    scatter_df = scatter_df,
    ancova_df = ancova_df, methods_df = methods_df,
    group_F = group_F, group_p = group_p, partial_eta2 = partial_eta2,
    slope = slope, slope_ci = slope_ci, slope_p = slope_p,
    slopes_F = slopes_F, slopes_p = slopes_p, slopes_ok = slopes_ok,
    slopes_verdict = slopes_verdict,
    raw_best = raw_best, raw_worst = raw_worst,
    adj_best = adj_best, adj_worst = adj_worst,
    ordering_changed = ordering_changed,
    adj_gap = adj_gap, adj_gap_ci = adj_gap_ci, gap_holm_p = gap_holm_p,
    xbar = xbar, p_prose = p_prose, fmt_n = fmt_n,
    metrics = metrics, json_output = json_output
  )
}

Per-group covariate means — quantify the imbalance the chart shows

sc <- shared$scatter_df
  cov_means <- tapply(sc$covariate_value, sc$group_label, mean)
  cov_means <- cov_means[!is.na(cov_means)]
  hi_g <- names(cov_means)[which.max(cov_means)]
  lo_g <- names(cov_means)[which.min(cov_means)]
  imb_note <- if (length(cov_means) >= 2 && hi_g != lo_g) {
    paste0("The groups sit at visibly different ", shared$covariate_h,
           " levels — ", hi_g, " averages ", shared$fmt_n(max(cov_means)),
           " while ", lo_g, " averages ", shared$fmt_n(min(cov_means)),
           " — which is exactly the imbalance the adjustment corrects. ")
  } else ""
  slope_note <- paste0(
    "The common within-group trend is ", shared$fmt_n(shared$slope, 3), " in ",
    shared$outcome_h, " per unit of ", shared$covariate_h, ". ")
  parallel_note <- if (shared$slopes_ok) {
    "The group trends look statistically parallel(equal-slopes check passed), which is what allows one shared adjustment."
  } else {
    paste0("CAUTION: the group trends are NOT parallel(equal-slopes check failed, ",
           shared$p_prose(shared$slopes_p), ") — which group is ahead depends on the ",
           shared$covariate_h, " level, so read the trends here rather than any ",
           "single adjusted number.")
  }
  list(
    title = "Outcome vs Covariate by Group",
    description = paste0(shared$outcome_h, " against ", shared$covariate_h,
                         ", colored by ", shared$group_h, "."),
    text = paste0(
      "Each point is one observation. ", imb_note, slope_note, parallel_note
    ),
    chart_labels = list(
      covariate_value = shared$covariate_h,
      outcome_value = shared$outcome_h,
      group_label = shared$group_h
    ),
    data = list(covariate_scatter = shared$scatter_df)
  )
}

# Card: ancova_table (table)
card_ancova_table <- function(shared, df, params) {
  headline <- paste0(
    "The ", shared$group_h, " row is the headline: it tests whether the groups ",
    "differ AFTER ", shared$covariate_h, " is accounted for — here ",
    shared$p_prose(shared$group_p),
    if (!is.na(shared$partial_eta2)) paste0(
      ", with the ", shared$group_h, " term explaining ",
      shared$fmt_n(100 * shared$partial_eta2, 1),
      "% of the remaining variation(partial eta-squared)") else "", ". ")
  cov_row_note <- paste0(
    "The ", shared$covariate_h, " row confirms the covariate matters(",
    shared$p_prose(if (is.na(shared$slope_p)) NA_real_ else shared$slope_p),
    " for its slope) — if it did not, plain group comparison would suffice. ")
  slopes_note <- paste0(
    "The interaction row is the equal-slopes diagnostic: ", shared$slopes_verdict)
  list(
    title = "ANCOVA Table",
    description = "F-tests for the covariate, the adjusted group effect, and the equal-slopes diagnostic.",
    text = paste0(headline, cov_row_note, slopes_note),
    data = list(ancova_table = shared$ancova_df)
  )
}

# Card: adjusted_means (table)
card_adjusted_means <- function(shared, df, params) {
  am <- shared$adjusted_means_df
  best_row <- am[am$group == shared$adj_best, , drop = FALSE]
  text <- paste0(
    "Adjusted(least-squares) means: each ", shared$group_h,
    " group&#x27;s model-estimated mean ", shared$outcome_h, " at the grand mean of ",
    shared$covariate_h, " (", shared$fmt_n(shared$xbar), "), with model-based ",
    "standard errors and 95% confidence intervals. ",
    shared$adj_best, " is highest after adjustment(",
    shared$fmt_n(best_row$adjusted_mean[1]), ", 95% CI ",
    shared$fmt_n(best_row$ci_low[1]), " to ", shared$fmt_n(best_row$ci_high[1]),
    "); ", shared$adj_worst, " is lowest. ",
    "The pairwise table below tests each adjusted difference with Holm-corrected ",
    "p-values(confidence intervals are unadjusted): ", shared$n_sig_pairs, " of ",
    shared$n_pairs, " pair(s) remain significant. Rows marked &#x27;yes' have a ",
    "difference whose Holm-corrected p-value is below 0.05; where the raw and ",
    "adjusted means disagree, trust the adjusted column — that is the ",
    "comparison at equal ", shared$covariate_h, "."
  )
  list(
    title = "Adjusted Means & Pairwise Differences",
    description = paste0("Least-squares means of ", shared$outcome_h, " per ",
                         shared$group_h, " group, plus Holm-corrected pairwise differences."),
    text = text,
    data = list(adjusted_means = shared$adjusted_means_df,
                pairwise_adjusted = shared$pairwise_df)
  )
}

# Card: methods (table)
card_methods <- function(shared, df, params) {
  text <- paste0(
    "ANCOVA is an observational adjustment, not a randomized comparison: it ",
    "earns its adjusted reading only through the assumptions listed here. ",
    "Adjusted means are model predictions at a common ", shared$covariate_h,
    " value; pairwise tests use Holm correction. ",
    shared$slopes_verdict
  )
  list(
    title = "Method & Assumptions",
    description = "The model, the adjustment, the correction used, and the assumptions the adjusted comparison rests on.",
    text = text,
    data = list(methods = shared$methods_df)
  )
}

This is the exact R code executed to produce the analysis above. Run it with the same dataset to reproduce these results.

Your data has more stories to tell. Run any analysis on your own data — validated R modules, interactive reports, AI insights, and PDF export. 500 free credits on signup.
Try Free — No Signup Sign Up Free
Something wrong? Report it and we'll fix it — or get your credits back
Download PDF · Print · Share ·Refresh this analysis ·Edit this report

Report an Issue

Tell us what's wrong. You'll get a free re-run of this analysis so you can try again with different parameters. If the re-run still doesn't meet your expectations, we'll refund your credits.

Want to run this analysis on your own data? Upload CSV — Free Analysis See Pricing