Executive Summary

Does math score differ by test preparation course once reading score is adjusted for?

Observations

1000

Groups Compared

Adjusted Gap (best-worst)

0.67

Group Test

p = 0.259

Slopes Check

passed

Covariate Slope

0.854

After adjusting for reading score, no statistically reliable difference in math score between the test preparation course groups remains (F = 1.28, p = 0.259). The adjustment changes the story: on raw means completed looks highest, but once reading score is accounted for, none leads — the groups differ in reading score and part of the raw gap reflects that, not the test preparation course themselves. The adjusted gap between none and completed is 0.67 (95% CI -0.50 to 1.84), with 0 of 1 pairwise adjusted difference(s) significant after Holm correction. Within groups, each 1-unit increase in reading score is associated with a 0.854 change in math score (95% CI 0.816 to 0.893).

What this means— written from the computed numbers · yours to edit

The short answer

No statistically reliable difference in math scores emerges between the course groups once reading ability is adjusted for. The course did not produce a detectable improvement.

The detail

After adjustment, the test preparation course groups do not differ significantly (F = 1.28, p = 0.259). The adjusted gap between none and completed is 0.67 points (95% CI −0.50 to 1.84), with zero pairwise differences surviving Holm correction. Within groups, each 1-unit increase in reading score is associated with a 0.854-point increase in math score (95% CI 0.816 to 0.893). The equal-slopes check passed, confirming the linear adjustment is appropriate. The raw ranking (completed highest) reverses under adjustment—a signal that raw means were contaminated by reading ability imbalance.

What this can't tell you

The p-value of 0.259 means the observed difference is consistent with chance at this sample size; the pattern was not detected. The confidence interval spans zero, so the true adjusted difference could plausibly be negative, zero, or modestly positive. This analysis assumes no unmeasured confounding; if students self-selected into the course based on factors correlated with math ability, the comparison remains biased.

Overview

Analysis Overview

Analysis of covariance: math score across 2 test preparation course groups, adjusted for reading score (1,000 observations).

N Observations1000

N Groups2

N Pairwise1

What this means— written from the computed numbers · yours to edit

The test preparation course comparison is confounded by reading ability: the two groups start at different reading levels, and since reading score strongly predicts math score, the raw group means reflect both course effect and covariate imbalance. ANCOVA removes the covariate's contribution by estimating what each group's math score would be if every observation had the same reading score (the grand mean, 69.17). This adjustment reveals whether the course itself is associated with math gains after holding reading constant.

Data Preparation

Data Quality

Row and group cleaning applied before fitting the model.

Initial Rows1000

Final Rows1000

Rows Removed0

Groups Dropped0

Levels Lumped0

What this means— written from the computed numbers · yours to edit

All 1,000 rows loaded were retained; no missing values in math score or reading score, and both groups (completed, none) met the minimum 3-observation threshold. The data required no cleaning, leaving the full sample intact for the adjusted comparison.

Visualization

Raw vs Adjusted Means

Each test preparation course group's raw mean of math score next to its reading score-adjusted (least-squares) mean.

What this means— written from the computed numbers · yours to edit

The short answer

The course group that completed shows the largest swing under adjustment, dropping 4.04 points from raw to adjusted mean. This large shift exposes that the raw advantage of the completed group was largely due to its higher reading ability, not the course itself.

The detail

Raw means: completed 69.696, none 64.078. Adjusted means: completed 65.658, none 66.329. The completed group moves 4.04 points downward under adjustment—the largest shift—because it started with the highest reading score (73.89 vs. 66.53 for none). When both groups are set to the same reading level (69.17), none actually sits slightly higher. The ranking reversal is the key finding: raw data favored the completed group, but once reading ability is held constant, none edges ahead. This pattern indicates that much of the raw gap reflected reading ability imbalance, not a course effect.

What this can't tell you

The adjusted means assume a linear relationship between reading and math scores and that the relationship is the same within both groups. The equal-slopes check passed, so these assumptions hold in the data. However, adjustment does not prove the course caused the observed pattern; it only removes the confounding influence of reading score on the comparison.

Visualization

Outcome vs Covariate by Group

math score against reading score, colored by test preparation course.

What this means— written from the computed numbers · yours to edit

The short answer

The two groups occupy visibly different reading score ranges, with the completed group concentrated at higher reading levels. This imbalance is exactly what drives the raw mean difference; the scatter shows why adjustment was necessary.

The detail

Completed students average a reading score of 73.89, while none students average 66.53—a substantial gap. Points in the scatter show a consistent upward trend within each group: math score rises with reading score at a common slope of 0.854 per unit. The group trends are statistically parallel (equal-slopes assumption passed), allowing a single shared adjustment. The completed group's advantage in raw math scores is largely explained by its rightward position on the reading axis; when both groups are projected onto the same reading level, the advantage vanishes. The scatter reveals no outliers or nonlinear patterns that would invalidate the linear model.

What this can't tell you

The scatter shows association, not causation. The fact that reading and math scores move together within each group does not prove reading causes math scores; both may reflect a common underlying ability. Similarly, the completed group's higher reading level could reflect selection bias (stronger students chose the course) rather than a course effect on reading.

Data Table

ANCOVA Table

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

Term	Df	Sum Sq	Mean Sq	F Value	P Value
reading score	1	1.535e+05	1.535e+05	2013	< 0.001
test preparation course (adjusted for reading score)	1	97.31	97.31	1.276	0.259
test preparation course x reading score interaction (slopes check)	1	19.6	19.6	0.257	0.613
Residuals	997	7.606e+04	76.29	—

What this means— written from the computed numbers · yours to edit

The test preparation course row (the headline) reports F = 1.276, p = 0.259, with partial eta-squared of 0.1%—the course term explains negligible remaining variation after reading score is entered. The reading score row confirms the covariate matters (F = 2012.537, p < 0.001), validating the adjustment. The interaction row (equal-slopes check) yielded F = 0.257, p = 0.613, showing no evidence that the math-reading relationship differs between course groups, so one shared slope is appropriate.

Data Table

Adjusted Means & Pairwise Differences

Least-squares means of math score per test preparation course group, plus Holm-corrected pairwise differences.

Group	N	Raw Mean	Adjusted Mean	SE	CI Low	CI High
completed	358	69.7	65.66	0.471	64.73	66.58
none	642	64.08	66.33	0.349	65.64	67.01

Comparison	Adj Difference	CI Low	CI High	P Value	Holm P	Significant
completed - none	-0.671	-1.836	0.495	0.259	0.259	no

What this means— written from the computed numbers · yours to edit

At the grand mean reading score of 69.17, none has the highest adjusted mean (66.329, 95% CI 65.645 to 67.013, n = 642) and completed the lowest (65.658, 95% CI 64.735 to 66.582, n = 358). The adjusted difference is −0.671 (completed minus none), with 95% CI −1.836 to 0.495. The Holm-corrected p-value is 0.259, so this difference is not statistically significant—the confidence interval includes zero.

Data Table

Method & Assumptions

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

Item	Detail
Design	Observational analysis of covariance — adjustment for one measured covariate; other confounders remain unadjusted.
Model	Linear model: math score ~ reading score + test preparation course, fit by ordinary least squares on 1,000 rows.
Adjusted means	Least-squares means: the model's prediction for each test preparation course group at the grand mean of reading score (69.17), with model-based standard errors.
Omnibus group test	F-test of the test preparation course term after reading score is entered (F = 1.28 on 1 and 997 df, p = 0.259).
Pairwise correction	Pairwise adjusted differences are t-tests on model contrasts with Holm-corrected p-values; confidence intervals are unadjusted.
Covariate effect	Within groups, each 1-unit increase in reading score is associated with a math score change of 0.854 (95% CI 0.816 to 0.893, p < 0.001).
Equal-slopes check	The equal-slopes check passed (interaction p = 0.613): the math score-reading score relationship looks similar across test preparation course groups, so a single adjusted comparison is reasonable.
Causal reading	Adjusted differences are associations consistent with a test preparation course effect, not proof of cause — groups were not randomized here and only reading score was adjusted for.

What this means— written from the computed numbers · yours to edit

This is an observational ANCOVA adjusting for one measured covariate (reading score); other confounders remain unadjusted. The linear model fits math score ~ reading score + test preparation course on 1,000 rows by ordinary least squares. Adjusted means are least-squares predictions at the grand mean reading score (69.17). Pairwise differences use Holm-corrected p-values; confidence intervals are unadjusted. The equal-slopes check passed (p = 0.613), confirming the adjustment is valid. Adjusted differences are associations consistent with a course effect, not causal proof, because groups were not randomized and only reading score was adjusted for.

Rate this report Was this the answer you needed?

The exact source that produced this report — yours to keep, read, and re-run.

Download PDF

How this was computed method · R source · citation

The code that did it

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)
  )
}

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

Download PDF · Print · Share ·Refresh this analysis ·Edit this report

What's wrong with this card?

Executive Summary

The short answer

The detail

What this can't tell you

Analysis Overview

Data Quality

Raw vs Adjusted Means

The short answer

The detail

What this can't tell you

Outcome vs Covariate by Group

The short answer

The detail

What this can't tell you

ANCOVA Table

Adjusted Means & Pairwise Differences

Method & Assumptions

ANCOVA — Adjusted Group Comparison

Why This Method?

What This Analysis Covers

Standard Library

Core Analysis Pipeline

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

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

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

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

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

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

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

Step 8: Chart datasets

Cite this analysis

Report an Issue