Regression analysis – method notes #

The decisions behind the regression analysis module – which estimator produces each quantity, how its interval is built, where a threshold comes from and what happens when a fit cannot run – and the numbers that validate them, computed on the example dataset (cafe.csv, its spec beside it). Sections follow the manual’s.

Variable selection #

Multinomial and two-part count outcomes get neither mediation nor moderation. Causal mediation has no model for a nominal outcome, and the product-of-coefficients fallback is not valid for a response with three or more unordered levels – path b would be one coefficient per non-baseline category, not one effect. The zero-inflated and hurdle families transmit an effect through two equations at once, so there is no single indirect path to decompose and no single simple slope to probe. In both cases the buckets hide, any picks are cleared, and the card points at the plain Poisson or negative binomial family. A multinomial outcome loses simple slopes for the same reason – one contrast per non-baseline category – while an ordinal one keeps both analyses: mediate() handles a polr outcome per category, and the proportional-odds model has one slope to probe.

Mediators must be numeric. Every mediator model on the Simple tab – the product-of-coefficients path and mediate()'s mediator model alike – is an lm(), which does not fail on a factor: it silently regresses on the integer level codes and returns NA standard errors. A categorical mediator is therefore rejected by name before the run rather than fitted on a scale that means nothing.

Causal mediation replaces the product of coefficients on every non-linear outcome. Path a is a linear coefficient and path b a log-odds or log-rate one, so their product is on no scale. mediation::mediate simulates, for every case, the outcome with the mediator set to its value under each level of the predictor and averages the differences, so ACME, ADE and the total effect come out on the response scale with quasi-Bayesian intervals and p-values from the same draws, their number the bootstrap replications setting. An ordinal outcome goes through mediate.order with a nonparametric bootstrap and percentile intervals, per outcome category, because no averaged effect or proportion mediated is defined for an ordered multi-category response.

The treatment contrast is mean ∓ 1 SD, not 0 → 1. mediate() contrasts control.value = 0 against treat.value = 1 by default. On a logistic or count model the response-scale effect is not constant – a one-unit step near the middle of the data moves the predicted probability far more than the same step in the flat tail – and 0 → 1 puts both counterfactuals outside the data for any uncentred variable (age in years, income, a test score), where the curve is flat and every effect collapses toward zero. Pinning the contrast to the range the data covers is what keeps the ACME readable; the price is that the number is per 2 SD of the predictor, which the note under the table states. A binary numeric predictor is contrasted between its two observed values and a categorical one between its first two levels, and the same contrast serves the conditional indirect effects and the sensitivity sweep.

The indirect effect’s interval is bias-corrected, without acceleration, and falls back to the percentile interval. An indirect effect is a product of two coefficients, and products of estimates have skewed sampling distributions, for which the plain percentile interval sits off-centre and rejects the null too readily. The correction shifts the quantile probabilities by z₀ = Φ⁻¹(share of replicates below the estimate); when every replicate falls on one side of the observed value z₀ is infinite and the plain percentile interval is used instead. Acceleration (BCa) is left out because its jackknife would cost one refit per case per path. One bc_ci serves the mediation, total indirect, conditional indirect and path-analysis bootstraps, so every bias-corrected interval on the module’s cards is the same construction.

Intervals on fewer than 1,000 replicates are disclosed, not refused. At B = 100 a single replicate below the point estimate drives z₀ to −2.33 and collapses the bias-corrected interval onto the smallest replicate. The conventional floor is 1,000, which is also mediate()'s default for its draws, so one threshold serves the bootstrap, the quasi-Bayesian and the percentile branch; WebR is too slow to impose it, so the card names the count and the procedure that ran and asks for Bootstrap replications to be raised before the bounds are reported. A resample that fails to fit is dropped and counted rather than hidden – an interval resting on 15 of 2,000 draws says so under the table.

Proportion mediated is withheld when the total effect is not significant. The ratio $a\times b/c$ has c in the denominator, and a c not distinguishable from zero at the configured alpha makes it unbounded – tiny denominators produce wildly fluctuating “proportions” – so the row is replaced by a note rather than left blank, which would read as an oversight. A ratio outside (0, 1) is still printed, with a note, because the bare number gives no signal of what happened: opposite signs on the indirect and total effects (inconsistent, or suppression, mediation) put it below 0, and an indirect effect larger than the total puts it above 1. mediate.order defines no proportion mediated, so nothing is suppressed on an ordinal card.

Several mediators share one joint outcome model and one resampling loop. The outcome model is Y ~ predictors + M₁ + … + Mₖ + covariates, so each path b controls for the co-mediators (parallel multiple mediation) and no mediator → mediator path exists – a serial chain is the Advanced tab’s. The total indirect effect per predictor is resampled as one quantity rather than assembled from the specific effects’ intervals, which would have no known coverage; for a linear outcome it equals c − c′. On a linear outcome one bootstrap loop resamples the data once per iteration and refits every mediator, moderator and conditional-indirect model on that same resample, which is what keeps a run with several predictors, mediators and moderators tractable in the browser.

Simple slopes are read off the full model. Each slope is the predictor’s conditional effect in the model the coefficient table reports – all predictors, all moderators and every predictor × moderator interaction – with the focal moderator held at the probe value and any other numeric moderator at its mean, so it is adjusted for the co-predictors rather than taken from a one-predictor side model.

Percentile probes are offered because mean − 1 SD can leave the data. Mean ∓ 1 SD is a sound probe on a symmetric moderator, but on a skewed one it can fall below the smallest value anyone in the sample has, and the slope reported there is an extrapolation of the fitted line rather than a description of a subgroup. The 16th, 50th and 84th percentiles – PROCESS’s default – are by construction values the data reaches. Whichever set is chosen, a probe outside the observed range is marked with a dagger and the range is named, on both tabs, so an extrapolated slope cannot be read as an observed one.

Johnson–Neyman boundaries outside the observed range are not reported. The boundaries are the moderator values at which the conditional effect’s t equals the critical value, solved from base R and MASS alone. A root beyond the data is a crossing of the extrapolated line, not a finding, so only the roots inside the observed range are printed, and when none is, the row says whether the effect is significant everywhere in range or nowhere. Which side of the boundaries the region lies on is read off the significance verdicts at the grid’s endpoints rather than the quadratic’s leading sign, so the table and the plot’s shading agree by construction; the share of cases in the region is counted over the actual data, not the plotted grid, since a region that sounds wide can hold few members.

The moderated-mediation row is Hayes’s index only where the conditional indirect effect is linear in the moderator. With one moderated path, or a two-level moderator, the indirect effect changes at one rate across the moderator and its slope is the index of moderated mediation – a single number whose interval tests whether the mediation is moderated at all. When both paths carry the moderator the effect is quadratic in it, and on a non-linear outcome a response-scale ACME does not change at a single rate; R reports whether the index is constant, and where it is not the row is relabelled as the change in the indirect effect per unit of the moderator – the rate at the moderator’s mean for the quadratic case, the average change between the outer probes for the ACME – with a note saying which quantity it is.

The sensitivity analysis is medsens, defined for a linear and a probit equation. ρ is the correlation between the errors of the mediator and outcome equations – the dependence sequential ignorability sets to zero – and medsens re-solves the ACME across a grid of ρ in steps of 0.1, its own default (halving the step doubles a run that is already the slowest thing an opt-in tick adds), so the breakdown point is reported to that precision and a curve that never reaches zero says so up to |ρ| = 0.9. A binomial outcome is refitted with a probit link, the only binary form the method has; a linear outcome’s equation stays as it is. The breakdown ρ and the two variance shares are invariant to the scale of either equation and carry over from the refit unchanged, but the refit’s own effect is for a one-unit change at a fixed covariate profile – not the card’s contrast – so its magnitude and curve are withheld there; on a linear outcome medsens’s effect is exactly the a × b above, so both are shown. With a probit link the effect under treatment and under control are different quantities, so those paths get two arms.

A categorical mediator’s probit refit is opt-in because it costs about a minute per path. On the Advanced tab a path through a categorical mediator was fitted with a logistic mediator equation, and medsens needs a probit one, so the equation is refitted before the sweep – about a minute per path in WebR – and only the scale-free half of the result, the breakdown ρ and the variance shares, is shown, since the refit’s effect is not the one reported above. Three shapes are named rather than skipped silently: a binary mediator under a binary outcome, for which medsens has no probit-on-both-sides method; the same mediator with the option off; and a chain through more than one mediator, which the method does not cover.

Doubled roles, constant columns and count outcomes are checked before R sees them. A name in two of the five buckets would collapse silently in the formula, so it is rejected by name; a column with fewer than two distinct values among the complete cases produces either a raw R error or NA coefficients, depending on the family, so it is rejected too. glm() refuses a negative count response outright and only warns on a non-integer one – a warning R never surfaces – so the first is a validation error and the second a warning that the likelihood is a quasi-likelihood; pscl treats the same conditions, and a response with no zeros, as errors, which is why the two-part families state them as preconditions.

The parameter count follows the family. A single-equation GLM fits the design columns plus an intercept, but nnet::multinom fits an equation per non-baseline outcome level – (J − 1)(p + 1) coefficients for a J-level outcome – and MASS::polr trades the intercept for J − 1 thresholds, p + J − 1. Charging either the single-equation figure would let a wildly over-parameterized fit past the guard that blocks a run with no more complete cases than parameters, and under-fire the events-per-parameter warning by the same factor; the two-part families add their zero equation’s columns and intercept, and θ for a negative binomial count part. Each dependent variable in a batch is checked against its own complete-case count. Regularized methods are exempt – fitting more predictors than observations is what ridge, LASSO and the elastic net are for – but their penalty is chosen by cross-validation, which fails inside cv.glmnet and ordinalNetTune with a raw R error below 4 complete cases, or when an outcome category is too small to leave every training fold non-degenerate; both are checked first and reported as validation messages, the second at 3 cases per category.

Events per parameter below 10 warn without blocking. In a logistic or count model the binding constraint is the number of events, not of rows: 200 rows with 8 positive cases and 4 predictors pass any row-count rule and have 2 events per parameter, enough for large, unstable coefficients with intervals that mean little. Peduzzi et al. (1996) put the floor at 10, so the check warns below it and still computes the estimate, since blocking would remove capability. The events are the rarest outcome category for the classification families and the summed count for the count families, both over the rows that survive listwise deletion; the parameters are the family-aware slope count, without intercepts; and regularized estimation skips the check, being designed for exactly this regime.

Model setup #

A numeric outcome is accepted for multinomial. The proportional-odds assumption is tested on outcomes coded 1, 2, 3 as often as on labelled ones, and the multinomial model is its standard fallback, so refusing a numeric dependent variable would block the fallback exactly where it is needed. The module factors the outcome itself, and the reference outcome selector reads its levels as it reads text ones.

The two-part families withhold the ANOVA table and the influence statistics. A pscl fit’s coefficients span two equations while its model matrix spans one, so a sequential or partial test on it is not well defined, and neither is a leverage. The other withdrawals – regularized estimation, model comparison, path mode – are of features built on a single-equation fit, and mediation and moderation are withheld for the reason under Variable selection. The residual scatter keeps a third scale, Pearson residuals against the expected count, because a two-part fit has neither a deviance residual nor a single linear predictor to plot against.

An MM fit reports no AIC, BIC or likelihood-ratio test. robustbase::lmrob maximizes no likelihood, so the omnibus is a joint Wald χ² of the slopes on the fit’s own covariance matrix, and the R² is the robust one – the same ratio taken on the fitted case weights rather than on raw sums of squares, hence not comparable with an OLS R² on the same data – with the robust residual scale in place of the residual standard error. The standardized betas come from an lmrob refit on scaled data, so the β column keeps the meaning it has everywhere else.

A robust fit gets no Q–Q plot, no Shapiro test and no goodness-of-fit block. lmrob defines no standardized residual to draw the plot from, and normality is not what an MM fit defends; the goodness-of-fit tests are withdrawn because no candidate’s reference distribution survives the down-weighting. Influence is routed rather than gated: the case weights are the estimator’s own influence report – the mean weight, its quartiles, the smallest, the rows rejected at robustbase’s 0.1/n cutoff and the rows counting for less than half a case, each read against the package’s own alarm limits – so they are printed in place of Cook’s distance rather than a diagnostic computed after the fact.

The robust and quantile families are withdrawn from path mode and refused by model comparison. In path mode only the outcome equation would be MM-estimated or quantile-fitted, and every path through it would be part robust and part not; the comparison search ranks candidates on AIC and BIC, which neither family has. Under the robust family mediation and moderation are inherited rather than gated – the indirect effect’s bootstrap refits through the robust fit, while the mediator equation stays OLS under the module’s rule that a mediator equation follows its own target’s type. Under the quantile family both report a note instead of a table, on separate grounds: an indirect effect is a product of coefficients that decomposes an effect on a conditional mean, and the quantile of a sum is not the sum of the quantiles; and a simple slope’s region of significance is read off a coefficient covariance matrix, which quantreg defines none for.

The quantile grid is fitted whole, and the default is the quartiles. quantreg::rq solves a linear program, and the whole grid costs about what one point does, so a single τ saves nothing. A median-only fit is largely redundant with the robust family, the more efficient estimator of a central line; the claim that justifies quantile regression – an effect that varies across the distribution – is invisible at a single τ, which is why the default is the quartiles rather than the median alone. Standardized betas, collinearity, the diagnostic and comparison blocks cannot be “all τ at once”, so they are computed at the reported quantile, and each says which τ it came from.

R¹ is a one-quantile measure with no likelihood beside it. Koenker & Machado’s R¹ is the model’s check-function loss against the same loss under an intercept-only fit at the same τ: it partitions no variance and is comparable neither across τ nor with an OLS R². The F beside it is quantreg’s own analysis of deviance against that intercept-only fit, printed only when the nid covariance inverts, since anova.rq returns nothing where it cannot; a quantile fit maximizes no likelihood, so no AIC, BIC or likelihood-ratio test stands beside them.

The reference outcome is hidden under regularized estimation. glmnet’s multinomial fit is symmetric – one equation per outcome level with no baseline – so a baseline choice has nothing to act on. Under the classic fit the same choice is handed to model comparison, so that its per-outcome model-averaged tables read against the baseline picked rather than against nnet::multinom’s default first level.

Firth’s ANOVA table is Wald, and so are its intervals. The likelihood-ratio route would refit each reduced model through a fitter that does not accept the penalty, so the ANOVA table reports Wald chi-square tests, which have the same basis as the coefficients printed beside them; the confidence intervals are Wald intervals on the penalized fit, matching the basis the rest of the GLM output uses. Everything else – diagnostics, influence, goodness of fit, ROC, mediation, moderation, the β column, the model-comparison card – takes the classic path, since the brglm2 fit is an ordinary glm object with a penalized score.

Group lasso is offered for the linear, binomial and Poisson families only. grpreg fits gaussian, binomial and poisson; ordinal fits go through ordinalNet and negative binomial through mpath, neither of which has a grouped penalty, so switching to one of those families falls back to plain LASSO rather than erroring, and the option is restored when a family that offers it is picked again. With no categorical predictor every group holds a single column and the fit is the LASSO.

lambda.1se is the default. It is glmnet’s own convention, and it answers the question a penalty is usually reached for – which predictors matter: the two cross-validated errors differ by less than the noise in either, so the rule gives away no measurable accuracy and buys a simpler, more stable selected set, a term kept only when it clearly earns its place. lambda.min is the choice when prediction accuracy rather than parsimony is the goal. Its label is criterion-neutral because the criterion differs by family: cv.glmnet and cv.grpreg minimize the cross-validated error or deviance, ordinalNetTune and mpath’s tuning maximize the cross-validated log-likelihood.

A manual λ snaps to the fitted path. Every backend fits along its own λ grid, and coef(s = ) would interpolate between grid points while the cross-validation error and the deviance ratio are read at the nearest one, so the value typed is replaced by the nearest grid point and the λ (selected) row reports that. The grid is on each method’s own scale – a different order of magnitude for each – which is what the Lambda range row is for. Inside a cross-validation fold or a split the same value snaps to the fold’s own path: passing it as lambda = would ask glmnet for a one-point path at a value the fold’s data may put far outside its grid.

The selection rule is re-applied inside every fold and split; the stability refits keep the full-data λ. Cross-validating a penalized fit re-runs the λ selection inside each training fold – reusing the full-data λ would carry tuning that saw the held-out rows into every fold and give back the optimism the cross-validation exists to measure – and multi-sample splitting re-tunes it on the selection half of each split, since a λ read off the full sample would carry the inference half into the selection and void the split; a manual λ is the one exception, a fixed constant tuned at no level. Stability selection instead refits every half-sample at the λ the full-data fit selected, because a selection frequency is a statement about that one model, on a λ grid truncated at that point – every backend reaches it by warm start from λ_max, so the prefix gives the identical solution at a fraction of the cost. Every fold deal and half-sample draw runs under the random seed setting, so the selected λ and the frequencies are reproducible.

Additional statistics #

A Type I table exists only where something sequential can be read. Types II and III come from car::Anova; Type I is the sequential split car does not offer, so it comes from stats::anova, whose GLM table names its columns differently and opens on an intercept-only row – both mapped onto car’s shape so one shaping block reads either. polr and multinom have no single-model anova() method to sequence, anova.lmrob is a two-model comparison test rather than a table, and a robust covariance replaces the sum-of-squares partition a sequential split is read off – car offers Types II and III on one, nothing sequential – so the option is hidden in those four cases and a stale pick is coerced to Type II. Type III is free of the parametrization only where the factor contrasts sum to zero, so the model is refitted under contr.sum first, and a robust covariance is recomputed on that refit rather than carried across two coefficient bases.

The robust covariance is offered where sandwich can line a score up with the model matrix. vcovHC and NeweyWest both read estfun() against the model matrix, so they line up only where the fit has one score contribution per coefficient and per row – lm and the single-equation glm families: linear, binomial, Poisson and negative binomial. polr and multinom carry threshold and per-contrast rows the model matrix has no columns for, a two-part fit spans two equations, and a penalized fit has no covariance to correct in the first place; the robust and quantile families keep their own estimators’ covariance. A Firth fit is an ordinary glm object and would line up, but estfun.glm returns the unpenalized score, so a sandwich on it would describe an estimator nobody fitted.

HC3 is the heteroscedasticity-consistent form, and model-based stays the default. Of White’s family the module uses HC3, MacKinnon and White’s (1985) finite-sample form, which discounts high-leverage rows the most. It is not the default because when the residual variance really is constant it gives away a little power against the textbook errors, and the diagnostics on the same card are what says which case the data are in. The omnibus follows the coefficients: an OLS F or a likelihood-ratio χ² printed above robust coefficients would contradict the note under them, so the model test becomes lmtest::waldtest on the same matrix – an F on a linear fit, a χ² on a GLM – and the ANOVA table puts its tests on it too, which is why it then reports Wald tests with no sums of squares, none being defined on that basis, and why the pseudo-R², read off the likelihoods, is the one part of a GLM card that cannot follow.

Newey–West runs with automatic bandwidth, AR(1) prewhitening and the finite-sample adjustment. The call is sandwich::NeweyWest(m, adjust = TRUE): the bandwidth is Newey and West’s (1994) automatic lag selection rather than a number the user would have to justify, the AR(1) prewhitening of the residuals is the function’s default and is left on, and adjust = TRUE adds the finite-sample correction. The bandwidth chosen is reported as a number of lags (bwNeweyWest, floored) so the note under the table says what ran. It reads the rows in stored order as the time order – the same assumption the autocorrelation test on the same card already makes – which is why the manual reserves it for time-ordered rows.

The reference distribution follows the family, and the rows a correction missed are named. The corrected interval uses a t on the residual df where the table reports a t – the linear fit – and the normal where it reports a Wald z – the GLMs – so each interval still agrees with the p beside it, and the exponentiated interval is the same interval on the ratio scale, so it moves too. The corrected matrix is read by term name rather than position, because an aliased term keeps its coefficient row and has no robust variance, and reading positionally would shift every SE after it onto the wrong predictor; such rows keep the model-based SE, statistic, p and CI and the note names them. When the matrix cannot be formed at all the table stays on the model’s own SEs and the note names the correction that failed, so the card never claims a correction that did not run.

A GLM’s β is the user’s convention, with SDs from the model matrix. There is no single accepted standardization for a GLM, because the outcome is 0/1 or a count with no SD on the link scale to divide by. Predictor only is b × SD(x), the change in the link-scale outcome per SD of the predictor; latent variable divides further by SD(y*), with Var(y*) = Var(η) + π²/3 under the logit (Long 1997), which puts β on a linear model’s footing at the cost of assuming the threshold model. Only the logit carries that fixed residual variance, so the latent option is hidden for the count families and, should it reach a non-logit link, is reassigned to predictor-only before the note is rendered – the note always names the convention that produced the numbers. A multinomial fit gets predictor-only within each outcome equation, since a baseline-category model has one latent scale per contrast rather than one for the model, and the two-part families qualify on the count families’ grounds. The SDs are taken from the fitted model matrix rather than the raw data, so a dummy and an interaction column are standardized on the same footing as a plain numeric predictor.

A linear β is a refit on standardized data, with a factor contrast left 0/1. The column comes from refitting the model on z-scored numeric columns – the outcome and the numeric predictors – through the family’s own estimator, so a robust or quantile fit’s β stays a robust or quantile quantity rather than an OLS one wearing the same table’s label. A categorical contrast stays a 0/1 column, so a numeric predictor’s β is $b\times \text{SD}(x)/\text{SD}(y)$ and a category’s is b ⁄ SD(y) – the group difference in outcome-SD units – and the note under the table says the two are not on one footing. The analytic β interval is confint on that refit rather than a rescaling of the raw interval, because a standardized interaction column is the product of z-scores, not the z-score of the product, so β is a second quantity with an interval of its own.

β’s bootstrap interval is bias-corrected and refits on every resample. A β is the coefficient over standard deviations all estimated from the same sample; an interval built by confint on the standardized refit treats those SDs as fixed constants, which forces it symmetric about β whatever the shape of the sampling distribution, and no HC correction reaches it – a sandwich fixes a Wald interval’s width, never its shape, while a ratio’s sampling distribution skews. Each replicate resamples the cases, refits and re-standardizes on its own resample – standardizing once and resampling the standardized model would hold the same SDs fixed – with the two conventions the point estimates use: the z-scored refit for linear, robust and quantile, b × SD(x) off the replicate’s own model matrix for a GLM, over SD(y*) where the latent scale was asked for. The interval is the same bc_ci as the mediation bootstrap’s, at the bootstrap replications setting under the bootstrap seed; a row fewer than 10 finite replicates reached keeps the interval it already carried and is named, and when no row reaches that many the card says the basis stayed analytic. B and its interval are untouched, so one row never carries two inference bases. The tick is offered for linear, robust, quantile, binomial, Poisson and negative binomial fits: polr and multinom refit by iterative optimization and can fail to converge on a draw that loses an outcome level, the two-part families do the same and have a β per equation with no single model matrix to take the SDs from, and a penalized fit reports no β at all.

Reading results – classic regression #

R²’s interval inverts the model’s own F test. A confidence interval for R² has no closed form; the one on the card is the noncentral-F interval, obtained by solving for the noncentrality λ at which the observed F sits at each tail probability of the configured level (Steiger & Fouladi 1997) and mapping each limit back with $R^2=\lambda /(\lambda +n)$ – the same test the F row reports, read backwards – so the interval and the p-value beside it can never disagree; a bound the search cannot reach is clamped to 0 or 1. f² is printed as a row rather than left to the reader because Cohen’s 0.02, 0.15 and 0.35 are f² benchmarks, and an R² read against them without the conversion is read against the wrong scale; the verdict beside R² reads on those benchmarks converted – 0.02, 0.13 and 0.26 of variance explained.

McFadden’s R² reads on McFadden’s bands, Nagelkerke’s on Cohen’s, and adjusted R² on its shrinkage. McFadden’s own guidance (1977) puts 0.2–0.4 at “an excellent fit”, so the verdict reads Weak fit below 0.1, Acceptable fit to 0.2 and Excellent fit from there, and nothing higher exists to award; an OLS reader waiting for 0.7 would never see it. Nagelkerke rescales Cox & Snell onto a 0–1 maximum, which matches an OLS R²’s range, so its verdict reads on Cohen’s R² bands – 0.02, 0.13, 0.26 – not McFadden’s. Adjusted R² is judged by the gap between R² and itself, the shrinkage: under 0.01 minimal, under 0.05 acceptable, above it notable, since a gap that size says the model’s terms are fitting the sample’s noise.

Type II is the default. It is car’s default and the more powerful choice when there is no interaction, which is the common case; it tests a main effect ignoring the interactions it appears in, which is interpretable only when the interaction is assumed away, so a model built around an interaction wants Type III – SPSS’s default and the applied literature’s – where each main effect is read in the interaction’s presence at the centre of the other terms. Type I answers what each term adds in the order entered and is right only when that order encodes a hierarchy. The overall model test is the same under every scheme; only the per-term decomposition differs.

β leads on the forest plot, with intervals from the standardized refit. A forest plot shares one axis across every term, and B is measured in whatever units each predictor happens to use: a predictor measured in thousands sets the axis and every other interval collapses to a point. β puts the terms on one scale, so the plot opens on β wherever the refit produced intervals, and the switch is there for the effects in their real units; the table above the plot is always on the B scale. The β intervals come from the refit on standardized data rather than from rescaling the raw ones because, with interactions in the model, a standardized interaction column is the product of z-scores and not the z-score of the product, so the two are not the same interval.

A cut-point’s p-value gets no verdict. polr pins the latent scale’s zero by an identification convention, not by anything in the data, so “this threshold differs significantly from zero” is a statement about the convention; the interpretation column places the cut-point below or at-or-above the origin instead, and the estimate, SE, CI and p stay on screen for anyone who wants the numbers. Each row is labelled with the pair of categories it separates rather than polr’s internal coefficient name, so a boundary can be read without decoding.

Separation is read off the size of an estimate on its column’s own step, gated by a fitted probability at 1. A binomial, ordinal or multinomial fit is flagged when some fitted probability is within 10⁻⁵ of 1 and some coefficient exceeds 10 logits – a 22,000-fold change in the odds – per SD of a continuous predictor’s column, or across the two values of a two-valued one, such as a dummy; the terms that cross it are named. The step matters for the indicator: its SD is at most 0.5, and glm stops a diverged coefficient at about 17 to 20 logits, so per SD it would never reach the bar, and the commonest separation – one category holding a single outcome – would pass unnamed. The criterion is the estimate’s size and not its Wald test, because polr’s numerical Hessian can give a diverged term a standard error near zero rather than a huge one.

Validation. On the café members with an overall score, recommend on taste and an indicator for an overall of 10 – 28 members, every one a recommender – returns the indicator at 17.3 logits with a standard error of 1,201; it is flagged and named on the binomial, ordinal and multinomial fits, while recommend on overall alone, at 0.848, is not.

VIFs are read off a centred design, and off the model’s own covariance where there is one. A raw X*W interaction column is mechanically correlated with its parents – arithmetic, not a modelling problem – and an ordinary y ~ x1*x2 can report a VIF near 160 uncentred and near 1 once the parents are centred, so the numeric predictors are mean-centred and the design rebuilt from that frame before the auxiliary regressions run; the auxiliary regressions carry an intercept, so centring the design’s own columns would change nothing – only centring the parents takes the non-essential collinearity out of the interaction columns. The fitted model’s coefficients are untouched; only the diagnostic uses the centred design. Where the fit has a covariance matrix of its own, car::vif on a refit over the centred frame is preferred, since it reads the inflation off that matrix and so describes this model’s coefficients rather than an OLS design sharing its columns – the two differ once the model’s weights do; an aliased fit, whose covariance has nothing to slice, multinom, where car’s indexing returns arithmetic that is not a VIF, and the two-part families, whose covariance spans both equations while their model matrix spans one, fall back to the design-matrix VIF. A grouped categorical predictor is one GVIF row, car’s GVIF^(1/(2·Df)) squared so that it is comparable to a plain VIF; a column whose auxiliary R² passes 1 − 10⁻⁶ is aliased rather than collinear, and its row is marked not assessable instead of carrying a 13-digit artefact.

The Q–Q panel draws studentized residuals; Shapiro–Wilk reads the raw ones. Raw least-squares residuals have unequal variances – each one’s depends on its row’s leverage – and are mildly correlated, while the confidence band around the Q–Q line assumes an independent, identically distributed sample, which the internally studentized residuals (rstandard) approximate far better; R’s own plot.lm draws its Q–Q panel from them for the same reason. Shapiro–Wilk stays on the raw residuals because that is what “Shapiro–Wilk on the residuals” means everywhere it is reported. The two read slightly different quantities on purpose, and the caption says so. The studentized set travels with its own rank order, because dividing by the leverage term can reorder residuals, and a row at leverage 1 yields no finite value and leaves the sample.

Breusch–Godfrey is built on the model’s Pearson residuals, and Durbin–Watson is two-sided. lmtest::bgtest extracts X and y and refits by ordinary least squares, ignoring family and link, so on a GLM it would report a linear-probability diagnostic rather than this model’s; the card builds the same first-order LM statistic itself – the Pearson residuals regressed on the design plus their own lag, n·R² referred to χ² on 1 df – so it diagnoses the model that was fitted. dwtest defaults to the positive-only “greater” alternative, which never flags a DW above 2, so the module asks for the two-sided test and reads the direction off the statistic when it is significant. The two-part families join the GLM arm because they have Pearson residuals and a count design for the auxiliary regression but no Gaussian OLS residual for Durbin–Watson.

Shapiro–Wilk samples 5,000 residuals, and the plots thin to 5,000 in rank order. shapiro.test caps at 5,000 observations, so a larger model is judged on a random subsample of that size, drawn under the reproducibility seed and reported with the verdict. One <circle> per row stops rendering long before that on the residual panels, so each is thinned to 5,000 points taken in residual rank order – the same picture at lower resolution, the extremes still at the extremes – and the kept ranks and the true n travel with the points, because a Q–Q plotting position taken from the thinned length would place every point at the quantile of its rank among 5,000 rather than among n. Above the same N the card says that Shapiro–Wilk and Breusch–Pagan reject on deviations far too small to move the coefficients or their intervals, and that OLS inference is at its most robust in that range, so both verdicts are to be read from the plots. A response the model reproduces almost exactly leaves floating-point noise, which every test would read as a verdict on the fit, so residuals whose spread is under 10⁻⁸ of the fitted values’ scale are reported as degenerate and not tested.

The diagnostics read the assumption test significance level; the model’s own inference reads the main one. The residual diagnostics, the outlier test, the Brant test, the goodness-of-fit tests and the test of directed separation all hold “the data or the model are fine” as their null – normal, uncorrelated, homoscedastic residuals, no case more extreme than n explains, one odds ratio per predictor, nothing the linear terms missed, calibrated probabilities, Poisson dispersion, every implied independence – so a small p is the alarm and the stricter check is the higher threshold: had their verdicts followed the main level, a user lowering it to 0.01 for their hypotheses would have made every check more lenient, the opposite of the intent, and the assumption test significance level exists for that reversed-null case. It is the level the comparison module’s assumption panel and the normality tests read too, so a residual normality test here and a normality check there give one verdict from one p. The p and statistic cells of those tables highlight and star at the same level, so a row’s formatting agrees with its verdict however the two levels are set. The F test, the coefficients’ tests and the accuracy-over-no-information-rate test are inference for the model – significance is the result the user hopes for – and keep the main level.

Cook’s D carries two thresholds, and every threshold uses the model’s rank. The textbook D > 1 is intuitive and works for small samples, but in larger ones virtually nothing crosses it – every observation looks safe while several pull the model around – so the size-aware $4/(n-p)$, which is 4 over the residual df, is printed beside it: zero flags there is genuinely reassuring, many flags there but none past 1 means spread leverage with no dominant case, and a flag past 1 always warrants attention. The leverage cutoff is $2p/n$, DFFITS’s $2\sqrt{p/n}$ and COVRATIO’s band $1\pm 3p/n$, all with p the model’s rank rather than its coefficient count – an aliased coefficient is NA, not estimated – and n the rows the fit actually used, so a model with dropped terms or rows is judged against its real p and n. Each flag count lists at most ten cases, since a badly conditioned fit flags thousands and the list exists to point at a few worth opening; the count beside it carries the rest.

The outlier test is Bonferroni-corrected on the largest studentized residual, and one-step on a GLM. A fixed |t| > 3 count flags a fixed fraction of every clean normal sample, growing with n; the largest externally studentized residual (rstudent) referred to its own t with the p multiplied by the number of cases keeps the flag rate at alpha whatever the sample size, which is why the count was replaced by the test. On a GLM influence.measures takes its quantities from the fit’s IRLS weights – one-step approximations rather than exact leave-one-out refits – and rstudent is the approximate jackknife deviance residual, which has no exact t reference and so takes the normal one; a note under the table says so, since they remain the right diagnostics but not the exact quantities the linear versions are. The test is reported as not computed when no finite studentized residual exists or, on a linear fit, when fewer than one residual degree of freedom remains after leaving a case out.

The case-weight table reads robustbase’s own alarm limits. A row is rejected once lmrob gives it a weight under its outlier cutoff of 0.1/n, and the two warnings – a rejected share past the estimator’s limit, a mean weight under it – are read against robustbase’s warn.limit.reject and warn.limit.meanrw, both 0.5 by default, which is also the cutoff under which a row counts for less than half a case; the table applies the package’s own limits rather than inventing a threshold of its own.

Hosmer–Lemeshow’s bin count is the user’s, and le Cessie–van Houwelingen runs beside it. hoslem.test (ResourceSelection) runs at the count the control sets, clamped to 3–50 with 10 the default, because the verdict genuinely depends on it – the same model can pass at g = 10 and fail at g = 8 – and a count the module chose silently would be a verdict it chose silently. R reduces the count when the fitted values cannot fill the bins and says so only in a warning, so the effective count is read back as df + 2 and printed beside the requested one; below three distinct fitted values no count works, and the test is reported undefined naming the case – a model whose few predictors are all categorical. The le Cessie–van Houwelingen–Copas–Hosmer unweighted sum of squares is the binning-free complement, with no knob and no grouping arbitrariness, and it still answers where Hosmer–Lemeshow is undefined. Its variance d′(W − WX(X′WX)⁻¹X′W)d is evaluated the way rms:::residuals.lrm does it – the weighted residual sum of squares from regressing d = 1 − 2p on the design with weights p(1 − p), QR-stable and without an explicit inverse – ported rather than called, since rms::lrm crashes the WebR page and its gof is unreachable. A saturated design – one distinct fitted value per free parameter, or a null variance within 10⁻¹⁵ of its scale – leaves that projection nothing to weigh and is reported as saturated rather than as a z of noise over noise. The two are printed with an agreement note because they are sensitive to different departures: agreement is stronger evidence than either alone, and a disagreement is named with what each responds to instead of being left to read as a contradiction.

Accuracy is tested against the no-information rate, and precision and F1 are read off the confusion matrix. The accuracy of always predicting the largest observed level is the floor an argmax classifier has to clear, and “77% correct” is not a result until it does; the row tests the correct count against that rate with a one-sided exact binomial test – P(X ≥ correct) under Binomial(n, NIR), evaluated through the regularized incomplete beta so that the tail keeps its digits where a summed mass function would not – the comparison caret::confusionMatrix prints as Accuracy > NIR. R sends recall per level; precision is read along the confusion matrix’s rows in the app, because a level with high recall and low precision is one the model over-assigns rather than understands, and F1 is derived from the pair since neither producer ships it – pROC::coords drops an f1 entry from ret =, and the per-level block has only recall to start from; two zeros make F1 0 by convention, not undefined. Every figure in the block is resubstitution, and the note under it says so because the ROC card’s provenance note does not cover it.

The Brant test is built from the cumulative binary splits by hand. Brant’s (1990) statistic is not in MASS; the module refits the model as the K − 1 cumulative logits Y > j it implies, with glm.fit on the ordinal fit’s own design, assembles the joint covariance of the split slopes – Cov(b_j, b_l) = (X′W_jX)⁻¹(X′W_jlX)(X′W_lX)⁻¹ with W_jl = diag(π_l(1 − π_j)) – and tests the contrasts b_j − b_j+1 = 0 for each slope as a Wald χ², intercepts left free to differ: one row per predictor on K − 2 df and an omnibus over all of them. Two splits and one predictor are the minimum there is anything to compare, so a two-level outcome or an intercept-only model reports the test as not applicable, and a split that fails to converge fails the whole test, since a slope from an unconverged logit is not a slope to contrast.

The count-family χ² tests carry no verdict past a fifth of the fitted counts under 5. The deviance and Pearson χ² are referred to χ² on the residual df, an approximation that holds as the fitted counts grow – not as n does – so on sparse counts both p-values are arithmetic without a reference; Cochran’s rule, a fifth of the expected counts below 5, is the line, past which the card withholds the two verdicts and a note says neither p can be read as a statement about fit, with the share and the minimum fitted count printed in their own row. The overdispersion row is Cameron & Trivedi’s (1990) regression-based score test – the auxiliary regression of $({(y-\hat{\mu })}^2-y)/\hat{\mu }$ on μ̂ without an intercept, its t taken as the z – two-sided so that it speaks to underdispersion as well, since the dispersion ratio alone has no inferential reference; it needs only the Poisson fit and so appears for Poisson models only, where a negative binomial fit prints its θ and standard error instead, θ being what makes it a different family. The dispersion verdict’s bands, 0.8 and 1.2, are the conventional margins around a ratio whose expectation is 1.

The two-part comparison is ranked, not tested, and the hurdle’s third row is the one that carries information. The single-equation counterpart – the same count formula and distribution refitted with the zero equation switched off – is this model on the boundary of its parameter space rather than outside it, so a likelihood-ratio test does not have its stated null there and the Vuong test is for non-nested models, which these are not; ΔAIC and ΔBIC, the counterpart’s minus this model’s so that positive favours this model, rank the two without a claim they cannot support, and the note under the table says so where a reader would look for the p-value. The Model choice verdict reads ΔAIC at ±2, the conventional margin for a difference worth acting on; the Excess zeros verdict flags the counterpart at a tenth short of the observed zeros, below which the second equation is fitting zeros the count model already produces. Expected zeros are summed from each fit’s own probability of a zero – predict(type = "prob") for the two-part fit, dpois or dnbinom at the counterpart’s fitted means. A hurdle model’s zero equation is fitted on the zero/non-zero indicator itself, so its fitted probabilities sum to the observed zero count and the first two rows agree whatever the fit is like; the third row, how far the counterpart is from the same target, is the one to read, and a note on the card says so. The two-part θ’s standard error is the delta-method image of pscl’s log-θ estimate.

The threshold rules are pROC’s coords("best"), tied optima broken toward the corner, with PPV, NPV and accuracy re-derived at the effective prevalence. Youden’s J and closest-to-(0, 1) are best.method = "youden" and "closest.topleft"; the cost rule is Youden’s J under best.weights = c(cost ratio, prevalence), pROC’s own weighting, whose unweighted default is c(1, 0.5). pROC returns every tied-optimal row in incidental threshold order, so the pick is the row nearest the (0, 1) corner, the first on exact ties. coords()’ own ppv, npv and accuracy always use the raw sample prevalence and ignore the weights, so the three are re-derived from the chosen sensitivity and specificity at the prevalence in use – identical to pROC’s numbers when that is the sample rate. The cost rule reports both directions of the asymmetry, the ratio and its reciprocal, because the two reported directions are what covers the other side of 1 – which is why the input is read as at least 1 – and equal costs make them the same call, so one row is shown; a single cost ratio has no clean reading across K one-vs-rest curves, so the rule is withheld for multinomial and applied at every ordinal cutpoint under the same ratio.

Prevalence is one vector per outcome level, read three ways, and the Brier skill ignores it. The population shares are entered per level and each family reads what it needs off the one set: a binomial model takes the positive level’s share, a multinomial one-vs-rest curve its own class’s, an ordinal cutpoint P(Y > k) as the shares above it summed – the same partition its cumulative score is built from – so neither multiclass family has to leave the correction on the table for want of a scalar prevalence. The vector is named by factor level and renormalized over the levels the fit kept, so a level dropped before the fit takes its share out of the total with it; a vector R cannot match to the fitted levels leaves the sample rates standing, and the basis note is read off the rows’ own flags rather than the option, so it says which one ran. The Brier skill’s reference stays the sample base rate under an override, because the Brier being scored is the sample’s.

The operating-point intervals re-select the threshold inside every resample. The cutoff is chosen in-sample, so its metrics are optimistically biased point estimates, and holding the observed threshold fixed while resampling would price the metric’s noise and not the selection’s. Each of the bootstrap replications redraws the controls and the cases separately at their observed sizes, rebuilds the curve, re-runs every rule on it and records the seven metrics, and the interval is the percentile interval of those draws – one resample serving every rule, so the cost rule costs one loop rather than two; on the multiclass cards the per-class and per-cutpoint operating points are re-selected inside the same case-resampling draws that price the aggregates. These remain in-sample figures when the AUC is cross-validated, and the table’s note says so.

DeLong is offered where the curve is binary; the multiclass aggregates bootstrap cases, not pseudo-observations. DeLong’s interval is for one binary ROC, so it is offered for the binomial curve and each ordinal cutpoint and hidden for multinomial, whose aggregates need resampling. ci.auc(method = "bootstrap") on the stacked micro-average would draw n·K pseudo-observations and ignore the within-case correlation between a case’s K one-vs-rest labels, giving an anticonservative interval, so micro, macro and Hand-Till M are bootstrapped in one loop that resamples cases and recomputes all three per replicate – one shared draw, so the three intervals describe the same resamples – and a draw that loses a class is dropped whole, since multiclass.roc would then span K − 1 classes and the macro mean average over K − 1, a different estimand. Hand-Till M is pROC::multiclass.roc’s average pairwise AUC. On the ordinal card the per-cutpoint intervals follow the selector, but the mean cutpoint AUC, Somers’ D and tau-c are always bootstrapped by resampling rows and recomputing the K − 1 curves together, because the cutpoints share rows and a per-cutpoint bootstrap would discard their correlation; the note under the table says the mean row is bootstrapped whatever was picked. The plotted curve is thinned to one point per 1/2000 of the unit square, below a pixel at any plot size, before it crosses to the page.

The ordinal concordance rows rank on the expected score P·(1, …, K). Somers’ D and Kendall’s tau-c need one ranking score per case, and polr’s linear predictor would do – but ordinalNet parametrizes the cumulative model the other way round, so a score read off the engine would flip sign between the classic and regularized fits. The expected ordinal score, the probability matrix times 1…K, is rank-identical to polr’s linear predictor and carries no engine’s sign convention, so both fits report on the same scale in the same direction, and it is the score the cross-validated block pools. Somers’ D is D_xy, concordant minus discordant pairs over the pairs untied on the outcome – the form that satisfies $D=2\text{AUC}-1$ and so reads on the AUC rows’ scale – and tau-c rescales the same concordance excess by the shorter side of the table. Both are accumulated over blocks of tied scores in one pass rather than taken from cor(method = "kendall"), whose pairwise walk costs 1.3 s at n = 4,000 – once per bootstrap replicate.

The cross-validated AUC is the mean of the per-fold AUCs with an influence-curve interval; every other cross-validated row carries Monte-Carlo precision. Each repetition deals the positives and the negatives into k stratified folds, refits with each fold held out and scores it, so every case is predicted once out of fold; the estimand is the mean of the per-fold held-out AUCs rather than one AUC over the pooled scores, and its interval is the influence-curve interval of LeDell, Petersen & van der Laan (2015) – each case’s influence value centred on its own fold’s AUC, the second moment over the sample estimating the mean’s sampling variance – with the between-repetition variance of the estimate added, since the influence curve prices one fold deal and averaging several leaves only their spread to add. That is a confidence interval: it says how far the answer would move on another sample of this size. Before the split every cross-validated row carried the across-repetition t interval, which was easy to read as a confidence interval and is not one – it is the variation between repetitions of the fold deal on one fixed dataset, shrinking toward zero as repetitions are added and saying nothing about sampling the data: on a 300-case binomial fit, [0.819, 0.828] across repetitions against [0.776, 0.871] from the influence curve for the same estimate. The rows that keep it keep it for a stated reason – the micro-average and Hand-Till M are functionals the influence curve is not defined for, and Somers’ D, tau-c and the Brier are not AUCs – and the table keeps the two in separate columns so that the one is not read as the other. The folds are stratified so that every fold holds both classes, and a level with fewer cases than folds stops the block with a named reason; the fold deal is the only draw the template makes, seeded per repetition from the reproducibility seed, and a penalized fit pins its own inner folds. The cases-scored figure is the smallest over the completed repetitions, since a repetition that lost a fold scored fewer.

Reading results – regularized regression #

Covariates carry no penalty, and every backend standardizes internally. A covariate is a control the model is meant to hold whatever the penalty thinks of it, so its design columns get a penalty factor of 0 (penalty.factor in glmnet and mpath, penaltyFactors in ordinalNet, group “0” in grpreg, which reads it as unpenalized) and they are fitted at every λ; the card reports them as Always included, attaches no selection frequency to them, and the stability and post-selection schemes never count them as candidates. All four backends standardize the predictors before fitting, so the penalty falls on each equally whatever its units, and return the coefficients on the original scale – which is why the coefficient table reads like an unpenalized one while the group table’s norm has to be taken on SD-scaled coefficients. The intercept is never penalized by any backend, so it is exempt from every verdict column.

Cross-validation deals at most 10 folds, at least 3, stratified for the classification families, and its SE counts the folds that scored. The fold count is max(3, min(10, n)), since cv.glmnet errors at two folds and ordinalNetTune and cv.grpreg go no lower; for binomial, ordinal and multinomial fits the folds are dealt within each outcome category and capped at the rarest category’s count, because the packages deal at random and a category missing from a fold’s training set aborts the fit rather than degrading it. The negative binomial path has no cv.glmregNB that runs in WebR – it opens a socket cluster – so its folds are refitted by hand on the full path’s λ grid and scored by held-out log-likelihood; a fold that fails drops out, and the per-λ standard error is divided by the folds that scored rather than the folds requested, since a failed fold would otherwise shrink the SE and pull the 1-SE rule toward a less penalized λ. Every fold deal runs under the reproducibility seed.

The fit statistic is one row, the deviance ratio under the family’s name. It is the one quantity every backend reports: glmnet’s dev.ratio, which for a linear fit is R² and for a 0/1, nominal or ordinal outcome – whose saturated log-likelihood is 0 – is McFadden’s R², so the row is labelled for the family rather than printing the same number twice; ordinalNet reports $1-\ell /{\ell }_0$ outright as devPct and exposes no null deviance; mpath’s null and residual deviances are per λ because θ is re-estimated at each. grpreg’s first path point carries the forced-in covariate group unpenalized, so it is not the null model – the null deviance is refitted as an intercept-only glm to make its ratio comparable with glmnet’s.

Ridge’s shrinkage is read against a refit at λ = 0, which needs n > p + 1. glmnet’s ridge path stops far short of zero – lambda.min.ratio scales a large ridge λ_max – so its smallest λ is not the unpenalized fit, and the baseline is a separate glmnet(alpha = 0, lambda = 0) refit. The intercept counts against the residual degrees of freedom, so n = p + 1 already saturates that refit and every percentage would be measured against a perfect fit; below that, or when the refit fails to converge, the column stays empty and the note says why, which is precisely the regime ridge is chosen for. A negative ratio is a sign flip against the unpenalized estimate, real under ridge with collinear predictors, and is labelled as such rather than folded into a percentage; the 90% and 50% colour bands are conventional margins for “barely touched” and “halved”.

λ prints to four significant figures. A λ is a tuning value the user has to be able to read back and type into the manual field, and formatNumber at the precision setting floors anything under it to “< 0.001” – where a penalized path routinely selects; toPrecision(4) keeps the digits without pretending the setting applies.

The group table ranks by the SD-scaled L2 norm. The grouped penalty shrinks the L2 norm of a group’s standardized coefficients, so that norm is both the quantity the selection was decided on and the group’s effect magnitude; grpreg standardizes internally and returns coefficients on the original scale, so the norm is taken on the coefficients multiplied by their columns’ standard deviations – the raw norm would rank a predictor measured in milliseconds above one in seconds through its units alone. A group whose columns all carry a zero penalty factor is forced in and never a candidate, so it is left out of the Predictor groups selected count and carries no frequency.

The degrees-of-freedom row counts what its engine reports, and ridge gets none. glmnet’s df is the number of non-zero entries in β at that λ – for a multinomial fit the terms non-zero in at least one outcome equation, hence Non-zero terms (any outcome) – and for ridge it is the raw predictor count at every λ, an “effective number of parameters” that provably is not one, so the row is omitted there; grpreg’s is fractional and genuinely effective, counting each group’s effective parameters, so only its row is labelled that way. The CV curve carries the same count along its top under the same label.

Multi-sample splitting re-tunes λ per split, Bonferroni-corrects within it and aggregates by the quantile rule. Each of B splits (20–200, default 50) draws half the rows – stratified by outcome for the classification families – re-selects on that half with λ chosen by the user’s own rule on the half’s own cross-validation (the fold count at the engines’ floor of 3, since the split only has to pick a λ and B splits average away a coarse deal; the negative binomial grid cut to 20 points, 1.7 s a fit against 8 s on its own 100), and refits the selected terms unpenalized on the other half through the family’s fitter – lm, glm, MASS::glm.nb or polr, whose p is read off its t as a normal. A split whose selected model reaches the inference half’s size leaves no residual degrees of freedom and is redrawn, up to 20 times; a split that still cannot be fitted contributes p = 1 for every term, which keeps the aggregation over B rather than over a data-dependent subset and is why the note says every p-value is conservative by that much. Within a split the raw p-values are multiplied by the size of the model it selected; across splits each term’s B values are combined by Meinshausen, Meier & Bühlmann’s (2009) quantile rule – the smallest γ-quantile divided by γ over a grid from $\lceil 0.05B\rceil /B$ to $1-1/B$, scaled by $1-\log {\gamma }_{\min }$ for having searched over γ. The unpenalized estimate is a separate fit of the full-data selected terms on the whole sample, reported only when they number fewer than n − 1, and carries no test because its standard errors would be computed on the rows that chose the terms. The option is withheld for ridge, which selects nothing, and for multinomial, where glmnet fits every outcome equation symmetrically with no baseline while nnet::multinom, the only unpenalized nominal fitter, is baseline-category, so the refit’s coefficient would not be about the coefficient printed beside it.

Complementary pairs count a term over the λ region and price the bound with stabs. Half samples (Meinshausen & Bühlmann 2010) refit each of B draws (10–500, default 50) at the selected λ on the truncated grid the Model setup note describes, and count a term non-zero there – the same rule the Status column reports. Complementary pairs (Shah & Samworth 2013) draw a half and its complement – $\lfloor n/2\rfloor$ within each level under stratification, so the complement is the larger half where a count is odd – refit both, and count a term selected in a half when it is non-zero anywhere on the path down to the selected λ, because the bound is stated for the union over that region and the expected size of the union is what prices it; a pair whose halves do not both fit is dropped whole, the bound being defined on complete pairs. The candidate unit follows the penalty – a group with the forced-in groups excluded, one outcome equation’s term under multinomial, a design column otherwise – and the bound is stabs::stabsel_parameters(p, cutoff = π, q = q̂, B = pairs, sampling.type = "SS")$PFER under assumption = "none" for the assumption-free row and "r-concave" for the sharper one, with π clamped to 0.55–0.99. It is withdrawn, with the reason named, when fewer than two candidates exist or an average half selects them all, when π is at or below q̂ ⁄ p – a term then clears the threshold on chance alone, and the smallest π that would work is printed – or when stabs could not be loaded; the r-concave row alone can be undefined at a threshold and set size the package refuses.

The coefficient trace is not drawn for multinomial, and ridge’s says that nothing enters or leaves. The cross-validation curve needs only the per-λ scores, which every backend produces; the trace needs one coefficient matrix along the path, and a multinomial fit has one per outcome level, so no single trace exists to show. An elastic net set to α = 0 is ridge, and neither sets a coefficient to exactly zero, so their trace caption says no term enters or leaves instead of inviting the entry-order reading.

Model comparison #

The search window is validated on its smallest candidate, not its largest. A minimum above the maximum or above the number of selected predictors leaves nothing to enumerate, and the card would then report a comparison over zero models – with the effective number of models, $1/\sum w^2$, read off an empty vector as a bare infinity – so the window is refused before R sees it. The parameter check is the single-model one, but taken on the smallest candidate the window admits: intercept, covariates, moderator main effects and the cheapest minimum-count predictors, counted per family as the parameter count is. Every candidate carries that fixed block, so when it alone exhausts the complete cases nothing in the window can be fitted and the ranking would come back empty behind a wall of per-model failures; the full model may legitimately be saturated while smaller candidates are informative, and a saturated candidate already gets its own singularity report. A constant predictor is refused on the same principle: it would fail every candidate containing it, silently halving the search and filling the failed-models list with hundreds of entries in place of one message.

AICc ranks, k is the full parameter count, and an over-parameterized candidate gets ∞ rather than a clamped denominator. The ranking, the deltas and the Akaike weights key off one criterion, and it is AICc because it converges on AIC as the sample grows and pays the part of a parameter’s price AIC leaves out when it does not; k is attr(logLik(), "df") – the dispersion of a linear or negative binomial fit and the thresholds of an ordinal one included – rather than the length of the coefficient vector, since the correction is defined on every parameter estimated. Once n − k − 1 ≤ 0 the correction has no finite value, and clamping the denominator would print a finite, too-good AICc for exactly the candidate the sample cannot support, so such a candidate is emitted as ∞, sorts last and takes a weight of zero. When the fixed terms alone exhaust the degrees of freedom every candidate is ∞, ∞ − ∞ would turn every weight into NaN, and the run falls back to the uncorrected AIC for the ranking, the deltas and the weights and says so on the card. A candidate whose AIC is itself non-finite – a perfect fit with n ≤ k – cannot be ranked at all and is listed among the failed models instead.

Selection uncertainty is printed against its maximum and withheld for a lone candidate. $-\sum w\log w$ has no scale of its own – its maximum is log of the number of candidates, reached when every one holds the same weight – so the card prints the two side by side; the effective number of models, $1/\sum w^2$, is the same uncertainty in model units. A single candidate drives both to zero, where the pair reads as a defect rather than as the certainty it is, so they appear only past one candidate; and both are computed over the retained set, as the weights are.

Large searches retain to Δ ≤ 20 past 2,000 candidates, on every criterion the card ranks by. Holding every fitted model is what runs a large dredge out of the WASM heap before it finishes, so past a floor of 2,000 candidates – or the display maximum, when that is larger – a candidate is kept only while it is within 20 information units of the running best. Burnham & Anderson put a model more than 10 units off the best at essentially no support; the cap is twice that, and a model 20 units behind carries less than e⁻¹⁰ of the best model’s weight, so the retained set changes no reported quantity to any printed precision. The test is taken on AICc, AIC and BIC alike, so a model hopeless on one criterion but competitive on another is still held; a running minimum only falls, so a candidate within the cap of the final best is within it of every running best and pruning against the running one drops nothing the full run would keep; a new best re-tests everything held, and the floor outranks the cap – when too few clear it, the best fill the rest. The per-predictor model and sign counts are accumulated as the loop goes, so they describe the whole search after pruning, while the weights, the averaging, the confidence set and the two uncertainty figures read the retained set; the card names both counts when they diverge, which is exactly when the cap fired. Below the floor nothing is discarded, so a search small enough to report in full reads as it did before the cap existed.

Model-averaged standard errors are Buckland’s, and a term’s absence enters as its weight times the estimate. The full average keeps the zero of every candidate lacking the term without renormalizing, so it is the weighted sum over the models that contain it; the conditional average divides that sum by the term’s importance. The standard error follows Buckland, Burnham & Augustin (1997): $\sum w\sqrt{{\text{SE}}^2+{(\text{estimate}-\text{average})}^2}$ over the models containing the term, the between-model spread entering per model under the root, which is why the average has to be computed first and the error in a second pass. For the full average the models without the term contribute an estimate of 0 with variance 0, and since the weights sum to one across the whole set their mass is 1 − importance, so they add (1 − importance) × |average| to the error; the conditional error is the same sum divided by the importance. Intervals are normal on those errors at the configured level. Averaging is done on the link scale, as it must be, and the exponentiated columns apply exp() to the averaged estimate and its bounds for display alone; the intercept has no conditional average, since every candidate contains it.

Importance is read against the prior inclusion share, and its bands are descriptive. A summed Akaike weight has no fixed no-effect point: a predictor the search window puts into every candidate scores exactly 1.0 whatever the data say, and in a plain all-subsets search every predictor starts at 0.5 before any model is fitted. The share of candidates containing the predictor – its model count over the models fitted – is that baseline, so the interpretation reads the excess over it scaled by the room above it, $(\text{importance}-\text{share})/(1-\text{share})$, and labels a predictor in every candidate, or in none, as not distinguishable rather than as very high or very low. The cut-offs – 0.9, 0.7, 0.5 and 0.3 of the room – are descriptive bands with no published criterion behind them, and the card says so. The sign counts read the predictor’s own columns only – its main effect and any dummy expansion, a categorical predictor counting as positive when any level is – because an interaction column carries a different quantity, which the interaction table counts separately; both counts and the model count are accumulated over every candidate fitted, while the importance itself is summed over the retained set, so the two describe the whole search and its competitive part respectively.

The AUC comparison is DeLong’s paired test against the AICc-first model, adjusted across M − 1, over the models displayed. Each reported binomial candidate’s curve is rebuilt from its retained fitted probabilities, given the interval the ROC section’s picker holds, and tested against the top-ranked model’s curve with pROC::roc.test(paired = TRUE) over the same cases – DeLong’s test under the DeLong interval, the bootstrap test under the bootstrap one, the bootstrap arm seeded per model so that an interval does not depend on how many models were ranked ahead of it. The reference is the model AICc ranked first, not the one with the highest AUC, because the question the column answers is whether any other candidate discriminates differently from the one the criterion chose – and since that choice was made on these same data the comparison is post-selection and its p-values are exploratory, which the note under the table says. The raw p-values are corrected with the global adjustment method across the M − 1 comparisons the card shows, the display maximum bounding the curves rebuilt, the intervals computed and the family adjusted alike.

Mediation inside the search is Baron & Kenny with Sobel, averaged, and tested once. A bootstrap for every candidate would multiply the search’s cost by the replicate count, so each candidate’s paths come from sub-models: path a from a linear regression of the mediator on the candidate’s predictors plus the covariates and moderator main effects, paths b and c′ from the candidate’s own equation – interactions included, so the paths belong to the model that was ranked – with the mediator added. The indirect effect is a × b with Sobel’s error $\sqrt{a^2\text{SE}{(b)}^2+b^2\text{SE}{(a)}^2}$; across candidates it is averaged exactly as a coefficient is, the Sobel error standing in for the coefficient’s, and the table’s test is the full average against its unconditional error read as a normal deviate. The smallest p-value over the candidate set would have been a selected extreme, optimistic by construction and with nothing to adjust it against, which is why one test is reported per path rather than one per model. The normal approximation is anti-conservative in small samples and for a skewed product, so the card’s note sends a confirmatory analysis to the targeted mediation with its bootstrap interval.

The failed-models list is capped at 100 entries under a title that carries the true count. A systematically rank-deficient design fails every candidate, which over a large dredge is hundreds of thousands of list items; the accordion title already states the count, so the list stops at 100 and an “… and N more” line closes it. A candidate is listed with the R error that stopped its fit, with the aliased terms R dropped when its coefficients came back NA – it is excluded rather than ranked on a partial coefficient vector – or as a non-finite AIC when there was nothing to rank it by.

Path analysis (advanced mode) #

The formula expansion is capped at 256 terms. (a + b + … + t)^20 is a legal formula that expands to about a million terms, and the editor re-parses on every keystroke with a deduplication that is quadratic in the term count, so an unbounded expansion locks the tab before the keystroke can be undone. The expanding helpers give up past 256 terms and the editor raises a diagnostic instead; the limit is far above any model anyone fits deliberately.

Indirect paths are capped at 100, and the enumeration at 20,000 visits. The number of simple paths to an outcome grows combinatorially with how densely the model is connected, and each one costs a bootstrap loop on the R side and a row in the decomposition table, which stops being readable well before the computation stops being feasible; past 100 the depth-first walk stops and the run is refused with the count. A branch that dead-ends short of the outcome collects no path, so the path bound alone would let the walk explore exponentially many of them – a node reachable only through an interaction term absorbs every branch that reaches it – and the visits are bounded too, far above what a model at the path limit costs, so only a model already past rescuing is cut short, with a message that names the density rather than a count the truncated walk could not have found.

Every equation is fitted on the system’s one complete-case set, on its own. The cases are those complete on every variable of the system, covariates included, and every equation is fitted on exactly that set: coefficients multiplied along a path are comparable only when the equations they come from saw the same rows, and an N reported per equation would hide that a chain crosses several. Cycles are rejected by the editor, so the drawn model is always recursive, and for a recursive model equation-by-equation least squares and maximum likelihood are consistent – joint estimation buys efficiency and a global fit test, not validity – which is why the card sends a reader who wants simultaneous estimation to the SEM module rather than approximating it. A two-level mediator is fitted by binary logistic regression rather than refused, since a two-level factor has one contrast and its equation is a well-defined GLM; a mediator with more levels would be regressed on its integer codes, so it is refused by name.

Post-treatment variables are dropped from a path’s outcome model, and a missing treatment or moderator is added to it. The direct and indirect effects of a mediation analysis rest on comparing like with like – the same units with and without the mediator’s push – and if the outcome equation also holds fixed a second variable the treatment influences, part of the treatment’s effect is blocked before it is measured and the split between direct and indirect stops meaning what its name says. mediate() and the sensitivity sweep both read the drawn outcome equation as their outcome model, so every other mediator standing in it that is a descendant of the path’s treatment is dropped from that path’s model, together with any interaction it enters – for that path only, since another path’s treatment has other descendants. The treatment is added when the drawn model has no direct edge from it to the outcome, because causal mediation is defined on an outcome model that contains it; a moderator of an edge into the mediator is added as a main effect when it is absent from the outcome equation, because mediate() holds a moderator fixed through the outcome model’s frame and one outside it cannot be pinned at all. Each departure leaves a model the card’s outcome coefficient table does not show, which is why the path’s own model is tabled under it, and why the several paths on one card are read against several models. The product bootstrap reads its b coefficient off the drawn equation as it stands, so a product path’s sensitivity curve is computed on that model and hands the sweep no drop list – a curve through a different model would not qualify that b.

A downstream covariate is left out of an equation. Covariates join every equation’s right-hand side, so one the drawn model places downstream of an equation’s own target would be conditioned on as a mediator or a collider of that equation – a control that changes the coefficient it was meant to protect. Reachability over the drawn edges decides it per equation, the covariate stays in every equation it is not downstream of, and the card names the covariates and equations involved so that the covariate line at the top is not read as what every equation adjusted for.

A multi-level categorical node along a path is reported as not computed. A categorical node enters the design as dummy columns, so its bare name is never a coefficient name: a two-level node has one dummy, which is a single well-defined contrast and multiplies along the path, while three or more levels have no single coefficient to multiply and the product would be an arbitrary pick among them. The path is listed with the reason rather than dropped, because an absent path reads as no effect, and it contributes nothing to its origin’s total; a total is reported only for an origin with at least one usable indirect path, since elsewhere it would repeat the direct effect. A categorical origin’s direct effect and its paths are keyed differently – by coefficient name and by node name – and are joined by that mapping so that one predictor does not split across two rows each carrying half of its total.

Direct, indirect and total effects come out of one bootstrap, with the observed sample as its first row. Every resample refits the whole system and recomputes every product, direct coefficient and conditional effect from it, so the total effect – the origin’s direct coefficient plus every indirect path from it – is assembled replicate by replicate from the two storage matrices rather than from separately bootstrapped parts, its interval rests on the same draws its parts do, and a replicate that lost one part loses the total. Row 1 of every storage matrix holds the observed sample, which is what the bias correction centres on; the interval is the bias-corrected one of the Simple tab, computed only on more than ten completed draws and otherwise withheld with a reason, and the sparsest row of a table is what qualifies its completion note, since the note speaks for every interval in it. A direct effect exists for every non-interaction coefficient of the outcome equation, mediators included – an endogenous node has a direct effect on the outcome as much as an exogenous one, and its total is only right when both halves are counted. The standardized effect rescales by SD(origin) ⁄ SD(outcome) with a categorical contrast’s factor 1, the intermediate mediators’ SDs cancelling along a path; the moderated-edge footnote exists because an unconditional product through a moderated edge is the effect at moderator = 0, a point usually outside an uncentered moderator’s range, and a direct coefficient carrying an interaction inherits the same reading.

A two-way interaction on the outcome is read both ways, and higher-order terms are not probed. The builder assigns no roles, so a drawn Y ~ X*W says nothing about which variable moderates which; the simple slopes are read off the drawn outcome equation symmetrically – X at levels of W and W at levels of X – and both readings are reported, since the reader, not the model, knows which is the moderator. Conditional indirect effects are probed from two-way coefficients, so a three-way or larger interaction on an edge is listed and left out of the values shown rather than probed at an arbitrary setting of its other factors. A conditional effect over a categorical moderator is obtained by re-levelling and refitting the affected equations inside each resample, and a resample that misses a rare level is skipped for that probe rather than aborting the loop.

A suppressed intercept’s R² is about the origin. summary.lm switches to an about-the-origin definition of R² when the intercept is removed, so an equation written with - 1 or + 0 reports a number that is usually inflated and not comparable to the R² of an equation that keeps its intercept; the card attaches the note to that equation’s fit statistics rather than hiding the row.

Fisher’s C is computed over the basis set of the drawn DAG, exogenous pairs left out, and an untestable claim withholds it. The basis set (Shipley 2000) has one claim per pair of variables the model joins by no edge, conditioned on the parents of both; the pair’s response is its later member in a layered topological order, so the claim is tested with an equation the model has, and the test adds the other variable to a model of the response on the conditioning set – an exact F where that equation is Gaussian, a likelihood-ratio χ² for a logistic mediator or the outcome’s own family, so a multi-level predictor is tested whole rather than one contrast at a time. $C=-2\sum \ln p$ is χ² on twice the number of claims under the joint null; an exactly-zero p would underflow C to infinity, so the smallest positive double stands in for it, which keeps C finite and still far past every critical value. A pair of exogenous variables is left out and counted: the builder draws directed edges only, so a user who believes two variables the model does not explain covary has no way to say so, and a claim nothing can retract is not one the model makes – C is then the test of the claims listed and of nothing else. A claim whose equation would not fit, or whose added variable contributed no estimable term – a predictor aliased against the conditioning set adds no parameter, and pchisq on zero degrees of freedom returns p = 0, a rejection read off a test that never ran – withholds C rather than shrinking its degrees of freedom, since the statistic is defined over the whole basis set; and a saturated model yields no claim, so the section reports that rather than an empty table and cites nobody, since no test ran. C’s verdict and every claim’s cells read the assumption test significance level, as the other model checks do – the model is the null, so a lowered main level would loosen the test.

Missing data #

One complete-case set serves the whole card. The runner drops every row incomplete on any variable that entered – the outcome, predictors, covariates, moderators and the mediators too, though the main formula never uses a mediator – inside R and whatever the missing-data setting says, before the first fit. The alternative, each block on its own complete cases, would let the coefficient table, the mediation tables, the diagnostics and the ROC block each describe a different sample under one N; the price is the rows the main model alone could have used, and the card names them beside N as dropped as incomplete.