Executive Summary
What is associated with charges across 1,338 observations
The short answer
Smoking status is the strongest driver of insurance charges: smokers face expected charges about 23,848.535 higher than non-smokers. The model explains 75.1% of the variation in charges, with 5 of 6 predictors statistically significant.
The detail
A linear regression on 1,338 observations identifies 6 predictors of medical insurance charges. Smoking status (smoker: yes) produces the largest effect: a shift of 23,848.535 in expected charges relative to the baseline (non-smoker) category, with very strong evidence behind it. Age, as a numeric predictor, also moves charges systematically. The model accounts for 0.751 deviance explained, meaning 75.1% of the total variation in charges is captured by these 6 predictors. Five predictors reach statistical significance at p < 0.05; one does not.
What this can't tell you
The model is fit to cross-sectional data—a single snapshot across individuals—so no trajectory or trend over time can be inferred. The large effect of smoking is an association observed in this sample; generalization to other populations depends on how representative these 1,338 observations are of the target population.
Analysis Overview
Gaussian (linear) regression of charges on 6 predictors across 1,338 observations.
The short answer
A linear regression model explains what drives medical insurance charges across 1,338 individuals using six predictors: age, sex, BMI, number of children, smoking status, and region. The model accounts for 75.1% of the variation in charges, meaning the predictors capture most of what makes charges differ from person to person.
The detail
The analysis uses a gaussian (linear) model—ordinary least squares—because charges is a continuous numeric outcome. All 1,338 rows had usable charge values; no rows were dropped. The model family was chosen from the outcome itself rather than assumed, ensuring the right statistical form for the question. Six predictors were modeled with no exclusions; all preprocessing steps (median imputation for missing numeric values, grouping of rare categorical levels) were applied to ensure usable data.
What this can't tell you
The cross-sectional design captures associations at one point in time; it does not track how charges evolve or whether changes in one predictor mechanically shift charges. Individual-level prediction will remain imperfect: even with 75.1% of deviance explained, substantial row-level variation remains unaccounted for, likely reflecting factors not in this dataset.
Data Quality
Outcome validation, predictor typing, and exclusions.
The short answer
All 1,338 observations had valid charge data and were retained. Numeric predictors with missing values were filled using column medians; categorical variables were grouped into standard levels with rare categories combined into "Other" and blanks labeled "Missing." No rows were dropped for missing outcomes.
The detail
The initial dataset contained 1,338 rows; the final modeled dataset also contains 1,338 rows, with 0 rows removed. Outcome validation confirmed charges is a continuous numeric measure with fractional and negative values, supporting the gaussian family choice. Numeric predictors underwent median imputation; text predictors were treated as categories with rare levels lumped and missing values explicitly labeled. All mapped predictor columns were usable in the model.
What this can't tell you
The data does not show which numeric columns had missing values or how many were imputed per column. A column-level missing-data summary would clarify whether imputation concentrated in a few predictors or was evenly distributed. Consider requesting a preprocessing detail export showing pre- and post-imputation counts by column.
Effects
Each predictor's coefficient with a 95% confidence interval and plain-language significance.
The short answer
Smoking status dominates the effects: smokers incur 23,848.535 more in expected charges than non-smokers. Age, BMI, and number of children are also strongly associated with higher charges; sex and region show weaker or no reliable associations.
The detail
Five of eight terms are significant at p < 0.05. Ranked by effect size: smoker: yes (23,848.535; p < 0.001), children (475.501; p < 0.001), bmi (339.193; p < 0.001), age (256.856; p < 0.001), region: southeast (−1,035.022; p = 0.0308), region: southwest (−960.051; p = 0.0448). Sex: male (−131.314; p = 0.693) and region: northwest (−352.964; p = 0.459) are not significant. The 95% confidence intervals for the four largest effects exclude zero cleanly; the two regional effects exclude zero but with narrower margins.
What this can't tell you
The model assumes linear effects: a one-unit increase in age, BMI, or children produces the same charge shift regardless of baseline values. If effects vary by subgroup (e.g., age matters more for smokers), this additive form would miss those interactions. A stratified or interaction-term analysis would be needed to test whether driver strength varies across subpopulations.
Strongest Predictor Effect
Model-predicted charges across the strongest numeric predictor, all else held typical.
The short answer
Age shows a tight, linear upward relationship with predicted charges. Moving from age 18 to 64 shifts predicted charges from 2305.52 to 14120.91—a consistent rise with no curvature or outlier departures.
The detail
The strongest numeric predictor is age. The model predicts charges at each observed age value (ranging from 18 to 64) while holding all other predictors at their typical values—numeric predictors at their median, categories at their mode. The relationship is strictly linear by construction: predicted charges increase from 2305.52 at age 18 to 14120.91 at age 64. No nonlinearity or threshold effects appear; each year of age is associated with the same increment in predicted charges across the entire range.
What this can't tell you
This effect assumes all other predictors remain at their typical levels. The actual effect of age on an individual's charges varies depending on their smoking status, region, and other characteristics; this chart shows the average slope. The tight linearity reflects the model's functional form, not the precision of individual predictions—substantial row-level residual variation remains, as the predicted-vs-actual chart shows.
Predicted vs Actual
How closely the model's predictions track the observed charges.
The short answer
Predictions and actual charges correlate at 0.87, indicating a tight overall relationship but with substantial scatter around the diagonal. The model captures the broad patterns but individual predictions remain imperfect.
The detail
Predicted and actual charges show a correlation of 0.87. Points scatter around the diagonal, with some observations lying well above or below the line of perfect prediction. This pattern is consistent with the 75.1% deviance explained: the model's six predictors account for most of what varies, but row-level factors not in the dataset create residual noise. Predictions range from negative values (a statistical artifact of the linear model) to high positive values tracking the actual distribution.
What this can't tell you
The scatter does not reveal whether prediction errors are systematic by subgroup (e.g., worse for smokers or high-age individuals) or randomly distributed. Residual diagnostics (residuals versus fitted values, Q-Q plots) would show whether the linear model's assumptions hold or whether transformations or robust methods are needed. A residual-by-group breakdown would clarify whether any predictor combination yields consistently poor predictions.
Model Diagnostics
The family decision trail and fit quality for the gaussian (linear) model.
| Metric | Value | Interpretation |
|---|---|---|
| Outcome type detected | charges is a continuous numeric measure (it has fractional or negative values), so ordinary least squares applies directly | What the outcome column's own values say about which model family applies. |
| Family chosen | Gaussian (linear) | Chosen because the outcome is a continuous numeric measure, so ordinary least squares — the gaussian family — is the right fit. |
| Deviance explained | 0.751 | The predictors account for 75.1% of the deviance in charges (identical to R-squared here). |
| AIC | 27115.5 | Model comparison score (lower is better); only comparable between models fit to this same data. |
| Baseline (intercept) | -11938.539 | The expected charges when every numeric predictor is zero and categories are at their reference level. |
| Rows used / dropped | 1,338 / 0 | Rows with a valid charges value were modeled; rows with a blank or unusable outcome were dropped, never imputed. |
The short answer
The gaussian (linear) model was selected because charges is a continuous numeric measure with fractional values, making ordinary least squares the appropriate method. The model accounts for 75.1% of deviance; all 1,338 observations were retained with no rows dropped.
The detail
Outcome type detected: charges is a continuous numeric measure with fractional or negative values, so ordinary least squares applies directly. Family chosen: Gaussian (linear). Deviance explained: 0.751 (identical to R-squared for gaussian models). AIC: 27,115.5 for model comparison. Baseline (intercept): −11,938.539, the expected charges when all numeric predictors are zero and categorical predictors are at reference level. Rows used / dropped: 1,338 / 0.
What this can't tell you
The negative intercept reflects the model's linear extrapolation to zero values and does not imply negative real-world charges. The diagnostics table does not report assumption checks (normality of residuals, homoscedasticity, linearity). Consider requesting a residual diagnostic plot to confirm whether the linear model's assumptions hold or whether weighted or robust methods would be more appropriate for this outcome distribution.
GLM Explorer — The Right Regression Family, Chosen From Your Data
Models any outcome on the predictor columns the user selects, choosing the correct generalized linear model family from the outcome itself: a two-value outcome gets logistic regression (odds ratios), a non-negative whole-number outcome gets Poisson regression (rate ratios, upgraded to quasi-Poisson when overdispersed), and anything else gets ordinary least squares (plain coefficients). The decision and its evidence are reported, never hidden.
Why This Method?
Plain linear regression quietly assumes a continuous, unbounded outcome. Run it on a yes/no flag or an event count and the coefficients, intervals, and predictions are all wrong in ways that look fine. A GLM keeps the same "outcome ~ drivers" question but swaps in the distribution and link the outcome actually follows — this module inspects the outcome, picks the family, and shows its work.
What This Analysis Covers
- Automatic family selection with an explicit decision trail
- Effect sizes in the family's native language (coefficients, odds
ratios, or rate ratios) with 95% confidence intervals
- Partial effect of the strongest predictor, others held typical
- Predicted-vs-actual fit (binned calibration for logistic)
- Diagnostics: deviance explained, AIC, dispersion, rows used
Standard Library
Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {outcome, base_1..base_N}. 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))Step 1: Row accounting + semantic column discovery
initial_rows <- nrow(df)
if (!"outcome" %in% names(df)) {
stop("column_mapping must map an 'outcome' column (the value to explain)")
}
pred_cols <- grep("^base_[0-9]+$", names(df), value = TRUE)
pred_cols <- pred_cols[order(as.integer(sub("^base_", "", pred_cols)))]
if (length(pred_cols) == 0) {
stop("column_mapping must map at least one predictor column(base_1)")
}
outcome_name <- humanize_semantic("outcome", col_map)
predictor_names <- setNames(humanize_semantic(pred_cols, col_map), pred_cols)Step 2: Inspect the outcome and choose the family
The decision trail: (1) exactly two distinct non-blank values -> binomial (logistic); (2) numeric by the 95% rule AND all non-negative whole numbers -> Poisson candidate (the overdispersion check after the fit decides Poisson vs quasi-Poisson); (3) numeric otherwise -> gaussian (ordinary least squares); (4) anything else cannot be modeled and stops with a humanized message.
raw_chr <- trimws(as.character(df$outcome))
nonblank <- !is.na(df$outcome) & raw_chr != "" & tolower(raw_chr) != "na"
n_nonblank <- sum(nonblank)
if (n_nonblank < 30) {
stop(sprintf(
"Only %d rows have a value in the outcome column('%s'). At least 30 are required to fit a regression.",
n_nonblank, outcome_name))
}
vals <- raw_chr[nonblank]
lv_norm <- unique(tolower(vals))
family_key <- NULL
success_display <- NA_character_
vm_ratio <- NA_real_
outcome_evidence <- ""
if (length(lv_norm) == 1) {
stop(sprintf(
"The outcome column('%s') has only one distinct value, so there is nothing to explain.",
outcome_name))
} else if (length(lv_norm) == 2) {Two distinct values -> binomial. Recognized yes-side tokens map to 1; otherwise the alphabetically-second level is modeled as the event.
pos_tokens <- c("1", "yes", "true", "y")
neg_tokens <- c("0", "no", "false", "n")
if (sum(lv_norm %in% pos_tokens) == 1 && sum(lv_norm %in% neg_tokens) == 1) {
success_norm <- lv_norm[lv_norm %in% pos_tokens]
} else {
success_norm <- sort(lv_norm)[2]
}
success_display <- vals[match(success_norm, tolower(vals))]
y <- as.integer(tolower(raw_chr) == success_norm)
y[!nonblank] <- NA_integer_
df$outcome <- y
family_key <- "binomial"
outcome_evidence <- paste0(
outcome_name, " takes exactly 2 distinct values, so it is a yes/no ",
"outcome; the analysis models the probability that ", outcome_name,
" = '", success_display, "'")
} else {
conv <- suppressWarnings(as.numeric(raw_chr))
n_num <- sum(nonblank & !is.na(conv))
if (n_num < 0.95 * n_nonblank) {
stop(sprintf(
paste0("The outcome column('%s') has %d distinct text values. This ",
"analysis needs an outcome that is numeric(continuous or a ",
"count) or has exactly two categories(yes/no). For a ",
"multi-category outcome, use the classification drivers tool ",
"instead."),
outcome_name, length(lv_norm)))
}
conv[!nonblank] <- NA_real_
df$outcome <- conv
y_ok <- conv[!is.na(conv)]
if (isTRUE(stats::var(y_ok) == 0)) {
stop(sprintf(
"The outcome column('%s') is constant after cleaning, so there is nothing to explain.",
outcome_name))
}
is_countlike <- all(y_ok >= 0) && all(abs(y_ok - round(y_ok)) <= 1e-6)
if (is_countlike) {
family_key <- "poisson"
vm_ratio <- stats::var(y_ok) / mean(y_ok)
outcome_evidence <- paste0(
outcome_name, " contains only non-negative whole numbers ",
"(variance/mean ratio = ", round(vm_ratio, 2),
"), so it is treated as an event count")
} else {
family_key <- "gaussian"
outcome_evidence <- paste0(
outcome_name, " is a continuous numeric measure(it has ",
"fractional or negative values), so ordinary least squares ",
"applies directly")
}
}Drop rows with a missing outcome (never impute the outcome).
df <- df[!is.na(df$outcome), , drop = FALSE]
if (nrow(df) < 30) {
stop(sprintf(
"Only %d rows have a usable value in the outcome column('%s') after validation. At least 30 are required.",
nrow(df), outcome_name))
}
if (family_key == "binomial") {
class_tab <- table(df$outcome)
if (length(class_tab) < 2 || min(class_tab) < 5) {
stop(sprintf(
"The outcome column('%s') needs at least 5 rows in each of its two categories to fit a logistic regression; the rarer category has %d.",
outcome_name, if (length(class_tab) < 2) 0L else as.integer(min(class_tab))))
}
}Step 3: Type each predictor — numeric if >=95% of values convert,
else categorical. Impute numeric NA with the median; label blank categories "Missing"; drop identifier-like near-unique text columns; lump categorical levels beyond 10 into "Other".
dropped_predictors <- character(0)
for (dc in pred_cols) {
v <- df[[dc]]
if (!is.numeric(v)) {
conv_d <- suppressWarnings(as.numeric(as.character(v)))
n_orig <- sum(!is.na(v) & trimws(as.character(v)) != "")
if (n_orig > 0 && sum(!is.na(conv_d)) >= 0.95 * n_orig) {
df[[dc]] <- conv_d
}
}
v <- df[[dc]]
if (is.numeric(v)) {
med <- median(v, na.rm = TRUE)
if (is.na(med)) { dropped_predictors <- c(dropped_predictors, dc); next }
v[is.na(v)] <- med
df[[dc]] <- v
} else {
v <- as.character(v)
v[is.na(v) | trimws(v) == ""] <- "Missing"
if (length(unique(v)) > nrow(df) / 2) {Near-unique text column (an ID or free text, not a predictor) — exclude BEFORE lumping, or the lump step disguises identifiers as 10 categories.
dropped_predictors <- c(dropped_predictors, dc)
next
}
tab <- sort(table(v), decreasing = TRUE)
if (length(tab) > 10) {
keep <- names(tab)[1:10]
v[!(v %in% keep)] <- "Other"
}
df[[dc]] <- factor(v)
}
}Step 4: Drop zero-variance predictors
for (dc in setdiff(pred_cols, dropped_predictors)) {
v <- df[[dc]]
zero_var <- if (is.numeric(v)) {
isTRUE(stats::var(v, na.rm = TRUE) == 0) || is.na(stats::var(v, na.rm = TRUE))
} else {
length(unique(v)) <= 1
}
if (zero_var) dropped_predictors <- c(dropped_predictors, dc)
}
model_predictors <- setdiff(pred_cols, dropped_predictors)
if (length(model_predictors) == 0) {
stop("No usable predictor columns remained after cleaning(all were constant, empty, or identifier-like).")
}
df_clean <- df[, c(model_predictors, "outcome"), drop = FALSE]
final_rows <- nrow(df_clean)
rows_removed <- initial_rows - final_rowsStep 5: Guard — need clearly more rows than model terms
n_terms_est <- sum(sapply(model_predictors, function(dc) {
if (is.factor(df_clean[[dc]])) nlevels(df_clean[[dc]]) - 1 else 1
}))
while (n_terms_est >= final_rows - 5 && length(model_predictors) > 1) {
drop_dc <- model_predictors[length(model_predictors)]
dropped_predictors <- c(dropped_predictors, drop_dc)
model_predictors <- model_predictors[-length(model_predictors)]
df_clean <- df_clean[, c(model_predictors, "outcome"), drop = FALSE]
n_terms_est <- sum(sapply(model_predictors, function(dc) {
if (is.factor(df_clean[[dc]])) nlevels(df_clean[[dc]]) - 1 else 1
}))
}Step 6: Fit the chosen family; on the Poisson path, test
overdispersion (Pearson chi-squared / residual df, with its chi-squared p-value) and upgrade to quasi-Poisson when the counts vary more than Poisson allows (dispersion > 1.5 AND p < 0.05).
dispersion <- NA_real_
disp_p <- NA_real_
overdispersed <- FALSE
if (family_key == "binomial") {
model <- glm(outcome ~ ., family = binomial(), data = df_clean)
family_label <- "Binomial(logistic)"
effect_label <- "odds ratio"
family_reason <- paste0(
"the outcome has exactly two values, and a straight-line model on a ",
"yes/no outcome predicts impossible probabilities; logistic ",
"regression models the odds instead")
} else if (family_key == "poisson") {
fit_pois <- glm(outcome ~ ., family = poisson(), data = df_clean)
df_res <- df.residual(fit_pois)
pearson_stat <- sum(residuals(fit_pois, type = "pearson")^2)
dispersion <- if (df_res > 0) pearson_stat / df_res else NA_real_
disp_p <- if (df_res > 0) {
stats::pchisq(pearson_stat, df = df_res, lower.tail = FALSE)
} else NA_real_
overdispersed <- is.finite(dispersion) && dispersion > 1.5 &&
is.finite(disp_p) && disp_p < 0.05
if (overdispersed) {
model <- glm(outcome ~ ., family = quasipoisson(), data = df_clean)
family_key <- "quasipoisson"
family_label <- "Quasi-Poisson"
effect_label <- "rate ratio"
family_reason <- paste0(
"the outcome is a count, but its Pearson dispersion is ",
round(dispersion, 2), " (chi-squared p ",
format_p(disp_p), ") — well above the 1.0 plain Poisson assumes — ",
"so quasi-Poisson widens the standard errors to stay honest")
} else {
model <- fit_pois
family_label <- "Poisson"
effect_label <- "rate ratio"
family_reason <- paste0(
"the outcome is a non-negative whole-number count and its Pearson ",
"dispersion is ", round(dispersion, 2),
", close to the 1.0 a Poisson model assumes")
}
} else {
model <- glm(outcome ~ ., family = gaussian(), data = df_clean)
family_label <- "Gaussian(linear)"
effect_label <- "coefficient"
family_reason <- paste0(
"the outcome is a continuous numeric measure, so ordinary least ",
"squares — the gaussian family — is the right fit")
}
exp_scale <- family_key %in% c("binomial", "poisson", "quasipoisson")
dev_explained <- 1 - model$deviance / model$null.deviance
if (!is.finite(dev_explained)) dev_explained <- NA_real_
aic <- tryCatch(as.numeric(stats::AIC(model)), error = function(e) NA_real_)
if (family_key == "quasipoisson") aic <- NA_real_Step 7: Effects table — 95% CIs (Wald, tryCatch-guarded);
exponentiated to odds/rate ratios for the log/logit families.
tidy_df <- broom::tidy(model)
ci <- tryCatch(suppressMessages(stats::confint.default(model)),
error = function(e) NULL)
if (!is.null(ci) && all(tidy_df$term %in% rownames(ci))) {
ci_lo_raw <- ci[tidy_df$term, 1]
ci_hi_raw <- ci[tidy_df$term, 2]
} else {
ci_lo_raw <- tidy_df$estimate - 1.96 * tidy_df$std.error
ci_hi_raw <- tidy_df$estimate + 1.96 * tidy_df$std.error
}
scale_est <- function(x) {
out <- if (exp_scale) exp(x) else x
out[!is.finite(out)] <- NA_real_
round(out, 3)
}
effects_all <- data.frame(
term_raw = tidy_df$term,
term = humanize_term(tidy_df$term, model_predictors, col_map),
estimate = scale_est(tidy_df$estimate),
ci_low = scale_est(ci_lo_raw),
ci_high = scale_est(ci_hi_raw),
p_value = format_p(tidy_df$p.value),
significance = sig_word(tidy_df$p.value),
stringsAsFactors = FALSE
)
rownames(effects_all) <- NULLIntercept: reported on the response scale in the diagnostics table and in prose; the chart excludes it so effect bars share one scale.
int_idx <- which(tidy_df$term == "(Intercept)")
intercept_est <- if (length(int_idx) == 1) tidy_df$estimate[int_idx] else NA_real_
baseline_response <- if (!is.finite(intercept_est)) NA_real_ else {
switch(family_key,
binomial = stats::plogis(intercept_est),
poisson = exp(intercept_est),
quasipoisson = exp(intercept_est),
intercept_est)
}
coefficients_df <- effects_all[effects_all$term_raw != "(Intercept)",
c("term", "estimate", "ci_low", "ci_high",
"p_value", "significance"), drop = FALSE]
rownames(coefficients_df) <- NULLStep 8: Predictor importance — max |z| per predictor, NA-filtered
term_raw <- tidy_df$term[tidy_df$term != "(Intercept)"]
term_z <- tidy_df$statistic[tidy_df$term != "(Intercept)"]
term_est <- tidy_df$estimate[tidy_df$term != "(Intercept)"]
term_p <- tidy_df$p.value[tidy_df$term != "(Intercept)"]
imp_rows <- lapply(model_predictors, function(dc) {
idx <- which(startsWith(term_raw, dc))Aliased (perfectly collinear) terms carry NA statistics — filter NA BEFORE which.max so an all-NA predictor is excluded instead of crashing (LAT-1445 class).
idx <- idx[!is.na(term_z[idx])]
if (length(idx) == 0) return(NULL)
best <- idx[which.max(abs(term_z[idx]))]
data.frame(
semantic = dc,
driver = predictor_names[[dc]],
abs_z = abs(term_z[best]),
direction = if (is.numeric(df_clean[[dc]])) {
if (term_est[best] > 0) "positive" else "negative"
} else "categorical",
p_raw = term_p[best],
best_term = term_raw[best],
best_est = term_est[best],
stringsAsFactors = FALSE
)
})
importance_df <- do.call(rbind, Filter(Negate(is.null), imp_rows))
if (!is.null(importance_df) && nrow(importance_df) > 0) {
importance_df <- importance_df[order(-importance_df$abs_z), , drop = FALSE]
max_z <- max(importance_df$abs_z, na.rm = TRUE)
importance_df$importance_score <- if (is.finite(max_z) && max_z > 0) {
round(100 * importance_df$abs_z / max_z, 1)
} else 0
importance_df$significance <- sig_word(importance_df$p_raw)
rownames(importance_df) <- NULL
} else {
importance_df <- data.frame(
semantic = character(0), driver = character(0), abs_z = numeric(0),
direction = character(0), p_raw = numeric(0), best_term = character(0),
best_est = numeric(0), importance_score = numeric(0),
significance = character(0), stringsAsFactors = FALSE)
}
top_driver_name <- if (nrow(importance_df) > 0) importance_df$driver[1] else "Unknown"
top_term_effect <- if (nrow(importance_df) > 0) {
round(if (exp_scale) exp(importance_df$best_est[1]) else importance_df$best_est[1], 3)
} else NA_real_
top_term_p <- if (nrow(importance_df) > 0) importance_df$p_raw[1] else NA_real_
top_is_numeric <- if (nrow(importance_df) > 0) importance_df$direction[1] != "categorical" else FALSE
top_term_label <- if (nrow(importance_df) > 0) {
humanize_term(importance_df$best_term[1], model_predictors, col_map)
} else "Unknown"One computed clause describing the top effect in the family's language.
top_effect_phrase <- if (!is.finite(top_term_effect)) {
"no predictor effect was estimable"
} else if (family_key == "binomial") {
if (top_is_numeric) {
paste0("each additional unit of ", top_driver_name,
" multiplies the odds of ", outcome_name, " = '",
success_display, "' by ", top_term_effect)
} else {
paste0(top_term_label, " multiplies the odds of ", outcome_name,
" = '", success_display, "' by ", top_term_effect,
" versus the baseline category")
}
} else if (exp_scale) {
if (top_is_numeric) {
paste0("each additional unit of ", top_driver_name,
" multiplies the expected ", outcome_name, " by ", top_term_effect)
} else {
paste0(top_term_label, " multiplies the expected ", outcome_name,
" by ", top_term_effect, " versus the baseline category")
}
} else {
if (top_is_numeric) {
paste0("each additional unit of ", top_driver_name,
if (top_term_effect >= 0) " adds " else " removes ",
abs(top_term_effect),
if (top_term_effect >= 0) " to " else " from ",
"the expected ", outcome_name)
} else {
paste0(top_term_label, " shifts the expected ", outcome_name, " by ",
top_term_effect, " versus the baseline category")
}
}Step 9: Partial effect of the strongest NUMERIC predictor —
predictions across its observed range with every other predictor held typical (numeric at median, categorical at its modal level). Computed from the fitted model, never canned.
num_rank <- importance_df$semantic[
sapply(importance_df$semantic, function(s) is.numeric(df_clean[[s]]))]
partial_sem <- if (length(num_rank) > 0) num_rank[1] else NA_character_
if (!is.na(partial_sem)) {
v_top <- df_clean[[partial_sem]]
xs <- seq(min(v_top), max(v_top), length.out = 60)
nd <- data.frame(row.names = seq_len(60))
for (dc in model_predictors) {
col <- df_clean[[dc]]
if (is.numeric(col)) {
nd[[dc]] <- rep(median(col), 60)
} else {
modal <- names(sort(table(col), decreasing = TRUE))[1]
nd[[dc]] <- factor(rep(modal, 60), levels = levels(col))
}
}
nd[[partial_sem]] <- xs
preds <- tryCatch(
as.numeric(predict(model, newdata = nd, type = "response")),
error = function(e) rep(NA_real_, 60))
partial_effect_df <- data.frame(
driver_value = round(xs, 4),
predicted_outcome = round(preds, 4),
stringsAsFactors = FALSE
)
partial_effect_df <- partial_effect_df[is.finite(partial_effect_df$predicted_outcome), ,
drop = FALSE]
rownames(partial_effect_df) <- NULL
partial_driver_name <- predictor_names[[partial_sem]]
} else {
partial_effect_df <- data.frame(driver_value = numeric(0),
predicted_outcome = numeric(0),
stringsAsFactors = FALSE)
partial_driver_name <- NA_character_
}Step 10: Predicted vs actual — response-scale fit check.
For binomial, raw 0/1 points are unreadable, so the standard binned calibration is computed instead: rows are grouped into up to 10 predicted-probability bins; each point is (mean predicted, observed share). The diagonal is perfect calibration either way.
fitted_vals <- as.numeric(fitted(model))
pred_actual_cor <- suppressWarnings(
stats::cor(fitted_vals, df_clean$outcome, use = "complete.obs"))
if (family_key == "binomial") {
br <- unique(stats::quantile(fitted_vals, probs = seq(0, 1, length.out = 11),
na.rm = TRUE))
if (length(br) < 3) {
br <- unique(c(min(fitted_vals) - 1e-9, stats::median(fitted_vals),
max(fitted_vals) + 1e-9))
}
bins <- cut(fitted_vals, breaks = br, include.lowest = TRUE)
pred_actual_df <- do.call(rbind, lapply(levels(bins), function(b) {
sel <- bins == b
if (sum(sel, na.rm = TRUE) == 0) return(NULL)
data.frame(predicted = round(mean(fitted_vals[sel], na.rm = TRUE), 4),
actual = round(mean(df_clean$outcome[sel], na.rm = TRUE), 4),
stringsAsFactors = FALSE)
}))
rownames(pred_actual_df) <- NULL
pred_actual_kind <- "calibration"
} else {
set.seed(42)
sidx <- if (final_rows > 1000) sample(final_rows, 1000) else seq_len(final_rows)
pred_actual_df <- data.frame(
predicted = round(fitted_vals[sidx], 4),
actual = round(as.numeric(df_clean$outcome[sidx]), 4),
stringsAsFactors = FALSE
)
pred_actual_df <- pred_actual_df[order(pred_actual_df$predicted), , drop = FALSE]
rownames(pred_actual_df) <- NULL
pred_actual_kind <- "rows"
}Step 11: Significant predictors + KPI metrics
n_sig <- sum(importance_df$p_raw < 0.05, na.rm = TRUE)
metrics <- list(
`Observations` = final_rows,
`Model Family` = family_label,
`Deviance Explained` = round(dev_explained, 3),
`Significant Predictors` = as.integer(n_sig),
`Top Predictor` = top_driver_name
)
if (is.finite(aic)) metrics$AIC <- round(aic, 1)Step 12: Diagnostics table — the family decision trail plus fit
quality, all values as strings (no e-notation).
baseline_phrase <- if (!is.finite(baseline_response)) {
"Not estimable."
} else if (family_key == "binomial") {
paste0("With every numeric predictor at zero and categories at their ",
"reference level, the predicted probability of ", outcome_name,
" = '", success_display, "' is ",
round(100 * baseline_response, 1), "%.")
} else if (exp_scale) {
paste0("With every numeric predictor at zero and categories at their ",
"reference level, the expected ", outcome_name, " is ",
round(baseline_response, 2), ".")
} else {
paste0("The expected ", outcome_name, " when every numeric predictor ",
"is zero and categories are at their reference level.")
}
diag_metric <- c("Outcome type detected", "Family chosen")
diag_value <- c(outcome_evidence, family_label)
diag_interp <- c(
"What the outcome column's own values say about which model family applies.",
paste0("Chosen because ", family_reason, "."))
if (family_key %in% c("poisson", "quasipoisson")) {
diag_metric <- c(diag_metric, "Overdispersion(Pearson chi-squared / df)")
diag_value <- c(diag_value, paste0(
formatC(dispersion, digits = 2, format = "f"),
" (chi-squared p ", format_p(disp_p), ")"))
diag_interp <- c(diag_interp, if (overdispersed) {
"Above 1.5 with a significant chi-squared test: the counts vary more than plain Poisson allows, so Poisson p-values would be overconfident."
} else {
"Near 1: the variance tracks the mean, so the plain Poisson model is adequate."
})
}
diag_metric <- c(diag_metric, "Deviance explained")
diag_value <- c(diag_value, formatC(dev_explained, digits = 3, format = "f"))
diag_interp <- c(diag_interp, paste0(
"The predictors account for ", round(100 * dev_explained, 1),
"% of the deviance in ", outcome_name,
if (family_key == "gaussian") " (identical to R-squared here)." else "."))
diag_metric <- c(diag_metric, "AIC")
diag_value <- c(diag_value, if (is.finite(aic)) {
formatC(aic, digits = 1, format = "f")
} else {
"not defined"
})
diag_interp <- c(diag_interp, if (is.finite(aic)) {
"Model comparison score(lower is better); only comparable between models fit to this same data."
} else {
"AIC is not defined for quasi-likelihood models — quasi-Poisson has no full likelihood."
})
diag_metric <- c(diag_metric, "Baseline(intercept)")
diag_value <- c(diag_value, if (is.finite(baseline_response)) {
formatC(baseline_response, digits = 3, format = "f")
} else "not estimable")
diag_interp <- c(diag_interp, baseline_phrase)
diag_metric <- c(diag_metric, "Rows used / dropped")
diag_value <- c(diag_value, paste0(
format(final_rows, big.mark = ","), " / ",
format(rows_removed, big.mark = ",")))
diag_interp <- c(diag_interp, paste0(
"Rows with a valid ", outcome_name, " value were modeled; rows with a ",
"blank or unusable outcome were dropped, never imputed."))
diagnostics_df <- data.frame(metric = diag_metric, value = diag_value,
interpretation = diag_interp,
stringsAsFactors = FALSE)Step 13: json_output machine channel
json_output <- list(
answer = paste0(
family_label, " regression of ", outcome_name, " on ",
length(model_predictors),
if (length(model_predictors) == 1) " predictor" else " predictors",
" across ", format(final_rows, big.mark = ","), " rows. ",
"Family chosen because ", family_reason, ". ",
top_driver_name, " is the strongest predictor: ", top_effect_phrase,
" (", sig_word(top_term_p), "). ",
n_sig, " of ", nrow(importance_df),
if (nrow(importance_df) == 1) " predictor is" else " predictors are",
" significant at p < 0.05; the model explains ",
round(100 * dev_explained, 1), "% of the deviance in ",
outcome_name, "."
),
cards = lapply(
c("tldr", "overview", "preprocessing", "coefficients",
"effect_of_top_driver", "predicted_vs_actual", "diagnostics"),
function(cid) list(id = cid, metrics = metrics)
)
)Machine-readable verdict for the overdispersion diagnostic (LAT-1783) — only the Poisson path actually COMPUTES a diagnostic; gaussian and binomial fits have no genuine assumption test here, so they emit nothing rather than inventing a check. As in standard_count_regression, overdispersion is "warn" not "fail": the module SELF-CORRECTS by switching to quasi-Poisson, so downstream should surface the model change, not cap claim strength.
assumption_checks <- if (family_key %in% c("poisson", "quasipoisson")) {
list(list(
name = "overdispersion",
verdict = if (!is.finite(dispersion)) "warn"
else if (overdispersed) "warn" else "pass",
statistic = if (!is.finite(dispersion)) "dispersion not estimable"
else sprintf("Pearson dispersion = %s, chi-squared p %s(threshold 1.5)",
round(dispersion, 2), format_p(disp_p))
))
} else {
list()
}
list(
initial_rows = initial_rows,
final_rows = final_rows,
rows_removed = rows_removed,
outcome_name = outcome_name,
predictor_names = predictor_names,
model_predictors = model_predictors,
dropped_predictors = dropped_predictors,
family_key = family_key,
family_label = family_label,
family_reason = family_reason,
outcome_evidence = outcome_evidence,
effect_label = effect_label,
exp_scale = exp_scale,
success_display = success_display,
vm_ratio = vm_ratio,
dispersion = dispersion,
disp_p = disp_p,
overdispersed = overdispersed,
dev_explained = dev_explained,
aic = aic,
model = model,
df_clean = df_clean,
coefficients_df = coefficients_df,
effects_all = effects_all,
intercept_est = intercept_est,
baseline_response = baseline_response,
importance_df = importance_df,
top_driver_name = top_driver_name,
top_term_effect = top_term_effect,
top_term_p = top_term_p,
top_is_numeric = top_is_numeric,
top_term_label = top_term_label,
top_effect_phrase = top_effect_phrase,
partial_effect_df = partial_effect_df,
partial_driver_name = partial_driver_name,
pred_actual_df = pred_actual_df,
pred_actual_kind = pred_actual_kind,
pred_actual_cor = pred_actual_cor,
diagnostics_df = diagnostics_df,
n_sig = n_sig,
metrics = metrics,
json_output = json_output,
assumption_checks = assumption_checks
)
}Compute shared resources
shared <- compute_shared(df, params, col_map)