Multiple linear regression, and when to switch

Regression analysis relates one outcome to several predictors at once, and the first decision is not which software menu to open but what kind of outcome you measured. A continuous measurement points to multiple linear regression, a yes or no outcome to logistic regression, a count to Poisson or negative binomial regression, and ordered categories to ordinal logistic regression.

This guide gives a decision table from outcome type to model, explains what each model estimates and which diagnostics belong to it, and works one example in full: a multiple linear regression on R’s built-in trees data, run through a deterministic engine. The other models are shown as runnable R code without quoted results, and each has its own guide with a worked example.

What a regression model does

Every regression model in this guide has the same three parts.

An outcome. The one variable you want to explain or predict: a blood pressure reading, a pass or fail, a number of events, a rating from “poor” to “excellent”.

Predictors. The variables you think are related to it, which can be continuous (age, dose, temperature) or categorical (group, site, treatment arm). Predictors are sometimes called covariates, explanatory variables or independent variables. “Independent” here only means “on the right-hand side of the model”; predictors are usually correlated with each other, and handling that correlation is the reason to fit them together.

A link between them. A rule that turns a weighted sum of the predictors into a statement about the outcome. In linear regression the weighted sum is the mean of the outcome directly. In the other models it is the mean on a transformed scale: the log odds of an event, the log of an expected count, or the log odds of being at or below a category.

The weights are the coefficients, and they are the main result. A coefficient answers the question “if this predictor were one unit higher and every other predictor in the model stayed the same, how would the outcome differ on average?” That “holding the others fixed” clause is what separates a coefficient from a simple correlation, and it is the reason to use regression rather than a series of two-variable tests.

Linear, logistic, Poisson, negative binomial and ordinal regression are close relatives. Nelder and Wedderburn showed in 1972 that models for normal, binomial and Poisson outcomes are special cases of one family, generalised linear models, each defined by a distribution for the outcome and a link function; logistic regression is the binomial case with a log odds link. McCullagh’s 1980 paper extended the same machinery to ordered categories. That shared structure is why the workflow below is the same for every model: choose the model from the outcome, fit it, read the coefficients with their intervals, then check the assumptions that belong to that model.

Choose the model from the outcome

Before looking at any results, answer two questions about how the data were collected.

1. What kind of outcome did you measure? This picks the model family.

Outcome	Example	Model	Coefficient on the reported scale
Continuous measurement	Weight, blood pressure, test score, timber volume	Multiple linear regression	Difference in the mean outcome per unit of the predictor
Binary (two categories)	Passed or failed, event or no event	Logistic regression	Odds ratio
Count of events	Number of seizures, visits, defects	Poisson regression	Rate ratio
Count whose variance exceeds its mean	Counts with a long right tail	Negative binomial regression	Rate ratio
Three or more ordered categories	Poor, fair, good; disease stage	Ordinal (proportional odds) logistic regression	Odds ratio for being in a higher category
Three or more unordered categories	Choice of transport, blood group	Multinomial logistic regression	Odds ratio per category against a reference

2. Does the design add structure the outcome type does not show? Some designs change the model whatever the outcome.

• Repeated measurements or clusters (several measurements per person, pupils within schools, samples within batches): the observations are not independent. Use a mixed model with a random effect for the cluster, which exists for continuous, binary and count outcomes. The repeated measures ANOVA or mixed model guide explains the choice.
• Time until an event, with some subjects not yet having had it: survival methods such as Cox regression, because a logistic model on “event yes or no” throws away the timing and mishandles censored follow-up. See the Kaplan-Meier and Cox regression guide.
• Values below a detection limit: a censored regression rather than substituting half the limit. See the guide to below-detection-limit data.
• A comparison of group means adjusted for a baseline value: this is a linear regression, usually presented as ANCOVA. The ANCOVA guide covers adjusted means and the parallel-slopes assumption.

The engine that ran the worked example below has a reference page for each model in the table: linear regression, logistic regression, Poisson regression, negative binomial regression and ordinal logistic regression. If your question is a comparison between groups rather than a model with several predictors, the guide to choosing a statistical test starts from the design instead.

Two shortcuts are worth resisting. The first is turning a continuous outcome into a binary one (for example “above or below the median”) so that logistic regression can be used. That throws away information and usually needs a much larger sample for the same precision. The second is treating every numeric outcome as continuous. A count of rare events, with most values at zero or one, is not well described by a straight line with normally distributed errors, and a linear model can predict negative counts.

Multiple linear regression

Multiple linear regression is the default when the outcome is a continuous measurement. It models the mean of the outcome as an intercept plus one coefficient times each predictor:

fit <- lm(outcome ~ predictor1 + predictor2, data = d)
summary(fit)   # coefficients, standard errors, t tests, R-squared
confint(fit)   # 95% confidence intervals for the coefficients

The intercept is the predicted mean outcome when every predictor equals zero. Often that combination is impossible or far outside the data, in which case the intercept is a mathematical anchor and should not be interpreted on its own. Centring a predictor (subtracting its mean) moves the intercept to a meaningful point without changing any slope.

Each slope is the difference in the mean outcome for a one-unit difference in that predictor, holding the other predictors fixed. For a categorical predictor, R creates one dummy-coded coefficient per non-reference level, each a difference from the reference level. If one predictor’s slope may depend on another, add an interaction (predictor1 * group); each main-effect coefficient then describes that predictor’s association only where the other predictor is zero or at its reference level, which is another reason to centre.

The overall F-test, on the last line of summary(fit), asks whether all the slopes are zero together. It cannot say which predictor matters; the coefficient tests and intervals do that.

R-squared is the proportion of the variation in the outcome that the model accounts for in this sample. It describes fit, not correctness: a model can have a high R-squared and still have the wrong shape, and a low R-squared is normal for outcomes with a lot of individual variation even when the coefficients are estimated precisely. R-squared also never decreases when a predictor is added, so it rewards overfitting; in Austin and Steyerberg’s simulations the sample R-squared needed many more subjects per predictor than the coefficients did to be estimated without bias, while adjusted R-squared behaved well. The residual standard error describes the same fit in the outcome’s own units: the typical size of a miss.

The assumptions, in order of importance

Textbooks often list the assumptions with normality first. It belongs last. In order of how much damage a violation does:

1. The model has the right form. The mean outcome changes linearly with each predictor, on the scale you entered it, and the effect of one predictor does not depend on another unless you included an interaction. A missing curve or interaction biases the coefficients themselves, and no amount of data fixes it.
2. The observations are independent. This comes from the design, not from a test. Repeated measurements, clustering or matching make the standard errors wrong, usually too small.
3. The spread of the residuals is roughly constant across the range of fitted values. When it is not, the coefficients stay unbiased but the standard errors and p-values are off. Heteroscedasticity-robust standard errors or a transformation are the usual remedies.
4. The residuals are roughly normal. This matters for exact small-sample inference. Lumley and colleagues showed in 2002 that in large samples linear regression is valid for any distribution of the outcome, because the estimates’ sampling distribution becomes close to normal. With small samples and heavy tails or strong skew, it matters more.

Two further checks are not assumptions but belong in every analysis. Influential observations: a few rows can move the coefficients a long way, and Cook’s 1977 distance measures how much each row shifts the whole set of estimates. Collinearity: when predictors are strongly correlated with each other, their individual coefficients become imprecise, even if the model as a whole predicts well. A variance inflation factor shows how much each coefficient’s variance is inflated by that correlation, and a correlation matrix of the predictors shows which pairs overlap before you fit.

How many observations?

Rules of thumb for linear regression vary widely, and figures such as 10 or 20 observations per predictor are often quoted. Austin and Steyerberg’s 2015 simulations found that about two subjects per predictor was enough for unbiased coefficients, accurate standard errors and confidence intervals with close to their stated coverage, under the conditions they simulated. That is a floor for estimation, not a target: it says nothing about having the power to detect an effect of a given size, and models built by selecting predictors on the same data need far more. Plan a sample size from the smallest effect that would matter in your field, not from an effect estimated in a pilot.

Worked example: timber volume from diameter and height

R ships a dataset called trees with measurements of 31 felled black cherry trees, from a forestry text by Meyer, reproduced in the Minitab Student Handbook. For each tree it records the trunk diameter in inches (labelled Girth in the data, although its documentation says it is a diameter), the height in feet and the volume of usable timber in cubic feet. Each row is a different tree, so the rows are independent units, and the outcome is a continuous measurement: the design multiple linear regression is for.

The research question is how timber volume relates to diameter and height, each adjusted for the other. That is a practical question in forestry, because diameter and height are easy to measure on a standing tree and volume is not.

Seeded and run in R; the method, linear regression with volume as the outcome and diameter and height as predictors, was chosen before the run. The full output and the R script that reproduces it appear below the article.

The result

Term	Coefficient	95% confidence interval	t on 28 df	p-value
Intercept	-57.99	-75.68 to -40.29
Diameter (Girth), per inch	4.708	4.167 to 5.249	17.82	< 0.001
Height, per foot	0.3393	0.07265 to 0.6059	2.607	0.0145

All 31 trees were analysed, with no missing values. The model’s R-squared was 0.948.

Reading it

Diameter. Among trees of the same height, each additional inch of diameter was associated with 4.708 more cubic feet of timber on average, with a 95% confidence interval from 4.167 to 5.249. The standard error was 0.2643, small relative to the estimate, so the interval is narrow.

Height. Among trees of the same diameter, each additional foot of height was associated with 0.3393 more cubic feet of timber on average, with a 95% confidence interval from 0.07265 to 0.6059 and p = 0.0145. The interval excludes zero but is wide relative to the estimate: the upper end is nearly twice the estimate and the lower end about a fifth of it. The data say height adds volume once diameter is accounted for, but they pin down how much only loosely.

Note the phrase “among trees of the same height”. The diameter coefficient is not the difference in volume between thick and thin trees in general. To the extent that thicker trees also tend to be taller, a model with diameter alone would credit diameter with some of height’s contribution. The two-predictor model separates the two, and that separation is the point of multiple regression.

The coefficients are associations, not effects of intervention. This is observational data on felled trees. The model describes how volume varies with diameter and height across these 31 trees. It does not say what would happen to a tree’s volume if something made it taller.

The intercept. The intercept of -57.99 is the predicted volume of a tree with zero diameter and zero height, which does not exist. It is far outside the data, and its negative value is a sign that the straight-line form is a local approximation over the range of trees measured, not a description of how trees grow. Do not interpret it, and do not use the model to predict volumes for trees much smaller or larger than these.

R-squared. Diameter and height together account for 94.8% of the variation in volume across these trees. That is very high, but it does not show that the linear form is right, as the next section explains.

Reading the diagnostics

The engine ran one formal assumption check and raised two data flags. Reading them in context is more useful than treating any one as a pass or fail gate.

Normality of the residuals (Shapiro-Wilk): not rejected. W was 0.974 with p = 0.644, so the test did not detect a departure from normality. With 31 observations such a test has limited power, so this is “not detected”, not evidence that the residuals are normal. As explained above, it is also the least important assumption.

Influential observations: four flagged. The engine reported that 4 of the 31 trees had a Cook’s distance above the conventional 4/n screening threshold, the largest being 0.605. The 4/n cut-off picks out rows for a closer look; it is not a test. Cook built the distance from confidence ellipsoids for the coefficients, so it measures how far deleting one row would move the whole set of estimates relative to their uncertainty. The right response is to identify those trees, check them for recording errors, and refit without them as a sensitivity analysis to see whether the coefficients move materially. Deleting them because they are inconvenient is not a valid response.

Extreme values: one flagged. One volume lies beyond Tukey’s fences, which sit one and a half interquartile ranges outside the quartiles. That is descriptive: nothing was removed, and an extreme outcome value is neither a problem if the model predicts it well nor an error unless a check shows it is one.

What the engine did not check, and you should. The engine’s assumption output for linear regression covers residual normality, influence and extreme values. It does not test the two assumptions that matter most here: the form of the model and constant spread. Both are judged from a plot of the residuals against the fitted values, which takes two lines:

fit <- lm(Volume ~ Girth + Height, data = trees)
plot(fit, which = 1)   # residuals vs fitted: look for a curve or a funnel
plot(fit, which = 5)   # residuals vs leverage, with Cook's distance contours
car::vif(fit)          # variance inflation factors for the two predictors

There is a reason to suspect curvature in this particular model before looking. A cylinder of diameter D and height H has volume pi times D squared times H divided by 4, and a cone a third of that, so if a trunk were either shape its volume would be proportional to the square of its diameter times its height. That relationship is multiplicative, not additive, and a straight-line model of volume on diameter and height can only approximate it. The usual alternative is to model the logarithms, which is also what the examples on R’s own help page for the trees data do:

fit_log <- lm(log(Volume) ~ log(Girth) + log(Height), data = trees)
summary(fit_log)
plot(fit_log, which = 1)

On the log scale, the geometric argument predicts coefficients near 2 for diameter and near 1 for height, and each coefficient is read as a percentage change in volume per percentage change in the predictor. Whether the log model fits these trees better is an empirical question the residual plots answer. This guide’s engine run fitted only the untransformed model, so it reports no numbers for the log version; run the code above and compare the two residual plots yourself.

The residuals-versus-fitted plot is the single most informative diagnostic for a linear regression, and subject-matter knowledge, here the geometry of a trunk, is often what tells you which shape to try.

Sensitivity: what this sample could detect

The engine also reports a sensitivity figure for the design. With 31 trees and two predictors, a test of the two predictors together at the 0.05 level had 80% power only for effects of Cohen’s f2 = 0.3466 or larger. Detecting a conventional medium effect of f2 = 0.15 with 80% power would need 68 trees in total. Here f2 measures the share of variance the predictors explain relative to what they leave unexplained, for the model as a whole, not for one coefficient.

This is a sensitivity analysis of the design, computed from the sample size, the number of predictors and the significance level. It is not “observed power” computed from the result, which adds nothing to the p-value.

Logistic regression: binary outcomes

When the outcome has two categories, such as an event that happened or did not, a linear model is the wrong tool: it can predict probabilities below zero or above one, and its constant-variance assumption cannot hold for a binary outcome. Logistic regression models the log odds of the event instead:

fit <- glm(event ~ age + group, data = d, family = binomial())
exp(cbind(OR = coef(fit), confint(fit)))   # odds ratios with profile-likelihood intervals

The coefficients are differences in log odds, and exponentiating them gives odds ratios: how many times higher the odds of the event are per unit of the predictor, holding the others fixed. An odds ratio is not a risk ratio, and the two diverge when the event is common; the odds ratio or relative risk guide explains when that matters.

The diagnostics that belong to logistic regression are different from those for linear regression:

• Events per variable. What limits a logistic model is the number of events in the smaller outcome category, not the total sample size. Too few events per estimated coefficient give biased, unstable estimates.
• Separation. When a predictor, or a combination of predictors, perfectly divides events from non-events, the maximum likelihood estimate does not exist: the coefficient runs off towards infinity and its standard error explodes. The output can then show a huge odds ratio with a large p-value, which looks like “no effect” but is the opposite. Firth’s penalised likelihood gives finite estimates in that situation, as Heinze and Schemper showed in 2002.
• Linearity on the log odds scale for each continuous predictor, checked by adding a spline or a squared term and comparing.
• Calibration: whether predicted probabilities match observed proportions across the range.

There is no residual normality or constant-variance assumption to check. The logistic regression guide works an example in full, including events per variable and what to do about separation.

Poisson and negative binomial regression: counts

When the outcome is a count of events (zero, one, two and so on, with no upper limit), Poisson regression models the log of the expected count:

fit_pois <- glm(events ~ group + age, data = d, family = poisson())
exp(cbind(RR = coef(fit_pois), confint(fit_pois)))   # rate ratios

Exponentiated coefficients are rate ratios: how many times higher the expected count is per unit of the predictor. When units were observed for different lengths of time or over different areas, add the log of the exposure as an offset, + offset(log(exposure)), so the model compares rates rather than raw counts.

The Poisson model makes one strong assumption: the variance of the counts equals their mean, given the predictors. Real counts often vary more than that, which is called overdispersion. It typically arises because units differ in ways the predictors do not capture, or because events cluster. When the mean is correctly specified, overdispersion does not stop the Poisson coefficients from converging on the right values as the sample grows: Gourieroux, Monfort and Trognon showed in 1984 that the Poisson estimates stay consistent even when the counts do not follow a Poisson distribution. The standard errors are another matter. They are computed as if the variance equalled the mean, so when the variance is larger they are too small, the confidence intervals too narrow and the p-values too small. Gardner, Mulvey and Shaw’s 1995 review makes the same point: an overdispersed Poisson model is likely to mislead unless its inferential statistics are corrected.

Checking for it is quick. A Pearson residual is the observed count minus the fitted count, divided by the square root of the fitted count. If the Poisson model holds, each of these has a variance of about 1, so the sum of their squares divided by the residual degrees of freedom should be near 1; values well above 1 indicate overdispersion. This is a rough guide rather than a test, and it is least reliable when the expected counts are very small. Cameron and Trivedi’s 1990 paper describes regression-based tests that formalise the check.

sum(residuals(fit_pois, type = "pearson")^2) / df.residual(fit_pois)   # near 1 if Poisson fits
fit_nb <- MASS::glm.nb(events ~ group + age, data = d)                 # negative binomial

Negative binomial regression adds a dispersion parameter that lets the variance grow faster than the mean, and Hilbe’s book treats it in depth. Its coefficients are read as rate ratios in the same way, and its standard errors account for the extra variation. When the counts are overdispersed, it is usually the better default. Two cautions apply before reaching for it automatically:

• Look for the cause. Overdispersion from clustering (several counts per person or site) is better handled by a mixed model with a random effect for the cluster, which the negative binomial does not replace.
• Excess zeros. If there are far more zeros than either model expects, for example because some units cannot have the event at all, a zero-inflated or hurdle model may describe the process better. Compare the observed number of zeros with the number the fitted model predicts.

A quasi-Poisson model, which scales the Poisson standard errors by the estimated dispersion, is a third option with the same coefficients as Poisson. The Poisson or negative binomial guide works a full example with a dispersion check.

Ordinal logistic regression: ordered categories

When the outcome is a set of ordered categories, such as a rating from “strongly disagree” to “strongly agree” or a stage from I to IV, the order carries information, and the spacing between categories is unknown. Two common approaches each discard something: treating the categories as numbers 1 to 5 in a linear model assumes the gaps are equal, and a multinomial model ignores the order.

Ordinal logistic regression in its proportional odds form, described by McCullagh in 1980, uses the order without assuming equal spacing. It models the cumulative odds of being at or below each category, with a separate threshold for each cut-point and one set of coefficients shared across all of them:

d$rating <- factor(d$rating, levels = c("poor", "fair", "good", "excellent"), ordered = TRUE)
fit <- MASS::polr(rating ~ group + age, data = d, Hess = TRUE)
summary(fit)
exp(cbind(OR = coef(fit), confint(fit)))   # odds ratios for a higher category

Note the sign. R’s polr writes the model as logit P(Y <= k) = zeta_k - eta, with a minus sign in front of the linear predictor eta, so a positive coefficient means higher categories are more likely. Its documentation notes that other software commonly uses the opposite sign. Check the convention before writing the direction of an effect.

The shared coefficients are the model’s key assumption, the proportional odds assumption: a predictor shifts the odds of being above any cut-point by the same factor, whether the cut is between “poor” and “fair” or between “good” and “excellent”. If a predictor mainly moves people from “good” to “excellent” and does nothing at the bottom of the scale, the assumption fails, and the single odds ratio averages over effects that differ.

How to check it:

• Fit the binary logistic models one cut-point at a time and compare the coefficients across cuts. Brant’s 1990 test formalises that comparison. Similar coefficients support the assumption; a coefficient that changes sign across cuts clearly breaks it.
• Treat the formal test with care. In large samples it rejects for small, unimportant differences; in small samples it misses large ones. The size and direction of the differences between cut-points matter more than the p-value.

When the assumption fails for one predictor, a partial proportional odds model relaxes it for that predictor only and keeps the order. Moving to a multinomial model gives up the order entirely and estimates many more coefficients, so it is the last resort for an ordinal outcome rather than the first. The ordinal or multinomial logistic regression guide works an example and covers the reference category in the multinomial model.

How to report a regression

A complete report names the model and the outcome, lists every predictor in the model, gives each coefficient of interest with its confidence interval and p-value, says how many observations were analysed and how missing data were handled, and describes the checks you ran. For the worked example, a results sentence would read:

In a multiple linear regression of timber volume on diameter and height in 31 black cherry trees, volume was higher by 4.71 cubic feet per inch of diameter (95% CI 4.17 to 5.25, p < 0.001) among trees of the same height, and by 0.34 cubic feet per foot of height (95% CI 0.07 to 0.61, t(28) = 2.61, p = 0.014) among trees of the same diameter. The model explained 95% of the variation in volume (R-squared 0.95).

Points to get right:

• Say “adjusted for” or “among units with the same” for every coefficient. A multiple regression coefficient without that clause is easily read as a simple difference.
• Report units. “4.71” means nothing without “cubic feet per inch”.
• Put the interval next to the estimate, and give exact p-values rather than “n.s.”.
• For logistic, Poisson, negative binomial and ordinal models, report the exponentiated coefficient (odds ratio or rate ratio) with its interval, and say which category or level is the reference.
• Describe the diagnostics in the methods section: for linear regression, how you checked the form of the model, the spread of the residuals and influential points; for counts, the dispersion check; for an ordinal model, how you assessed proportional odds.
• Say how the predictors were chosen, before seeing the data or not. A model whose predictors were selected by their p-values needs that disclosed, because its intervals are too narrow.

Common mistakes

Reading a coefficient as a causal effect. A regression on observational data describes associations, adjusted for the predictors included. Unmeasured variables can still drive both a predictor and the outcome. Causal language needs a design that supports it, such as randomisation, not a regression. Confounding variables shows what adjustment can and cannot remove, and correlation vs causation works through a case where a third variable creates the whole association.

Forgetting what the coefficient is adjusted for. Adding or removing a predictor changes the meaning of every other coefficient. Two models with different predictors answer different questions, and their coefficients for the same variable are not comparable.

Choosing predictors by p-value. Stepwise selection and keeping only significant predictors produce coefficients biased away from zero and intervals that are too narrow. Choose predictors from the question and prior knowledge, before seeing the outcome.

Treating a high R-squared as proof of a good model. The worked example has an R-squared of 0.948 and still has a strong reason to suspect the wrong shape. Look at the residual plot.

Ignoring the design. Repeated measurements, clusters, matched pairs and censored follow-up all break the independence or completeness a standard regression assumes. Model them, do not ignore them.

Using Poisson for overdispersed counts. The point estimates may look reasonable while the intervals are far too narrow. Check the dispersion every time.

Treating “not significant” as “no association”. A coefficient whose interval includes zero and also includes associations large enough to matter is inconclusive, not null. Look at the interval, not only the p-value.

When regression is not the first step

Regression is the right tool when the question is about how an outcome relates to several predictors at once. Simpler questions have simpler answers:

• One continuous outcome, two groups and nothing to adjust for: a t-test answers the question directly. Student’s version is mathematically a linear regression with one binary predictor; Welch’s version, the safer default, relaxes the equal-variance assumption that the regression makes.
• Three or more groups, nothing to adjust for: a one-way ANOVA, which is again a linear regression with one categorical predictor.
• An association between two categorical variables: a chi-square or Fisher’s exact test.
• Agreement between two measurement methods: a Bland-Altman analysis, not a regression of one method on the other.

Fitting the regression version of these is never wrong, and it becomes the better choice as soon as there is a second predictor or a covariate to adjust for.

Rerunning the example

Below the article you will find the R script exactly as the engine generated it, with the seed and the bundled data path, and the output it printed. The script reads its data from a data/ folder; the data file is offered for download with it. On the same package versions it reproduces the regression summary (coefficients, standard errors, t statistics, p-values and R-squared) and the sensitivity figure. The confidence intervals and the diagnostics come from the engine itself, and these calls on the same model reproduce them: confint(fit) for the intervals, shapiro.test(residuals(fit)) for the normality check, cooks.distance(fit) compared with 4 / 31 for the influence flag, and boxplot.stats(trees$Volume)$out for the extreme value, where fit is lm(Volume ~ Girth + Height, data = trees).

The guide to reproducible statistical analysis explains what a reviewer needs in order to rerun a result like this one.

# Results

- Test: Linear regression
- n: 31
- statistic: 17.82
- p-value: < 0.001
- effect size: R-squared = 0.948
- 95% CI: [4.167, 5.249]
- seed: 20261022

## Assumption checks
- normality (Shapiro-Wilk): passed=TRUE (no evidence against normality)

## Power (sensitivity, not observed power)
- minimum detectable effect at n=31: Cohen's f2 = 0.3466
- to detect a conventional medium effect: 68 total to detect f2=0.15
- alpha = 0.05, target power = 0.8

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Linear regression  |  engine 0.30.0  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261022)
.data_path <- "data/regression-analysis-guide.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("Volume", "Girth", "Height")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- summary(lm(Volume ~ Girth + Height, data = data))
print(result)

# Sensitivity: minimum detectable effect at n=31 (power=0.8, alpha=0.05). NOT observed power.
power <- pwr::pwr.f2.test(u = 2, v = 28, sig.level = 0.05, power = 0.8)
print(power)