Changelog

Source: NEWS.md

MAIHDA 0.3.0

New features

• plot(type = "upset") now honours order_by, which it previously accepted and silently ignored. order_by = "predicted_desc" gives the ranked caterpillar of type = "predicted" with the UpSet category matrix in place of the long text labels. The default is unchanged (intersection size, largest first), and the panel caption now reports the order actually drawn.

• order_by gained "size" (largest stratum first) on both type = "predicted" and type = "upset".

• summary() on a maihda analysis gained which = "adjusted", matching tidy(). It returns the null model by default, whose fixed effects are the intercept and covariates only because the strata dimensions are its random-effect grouping – easy to misread as the dimension main effects having gone missing. print() now names which of the two models it is showing.

• summary() now reports fixed-effect standard errors, Wald statistics, two-sided p-values and intervals for the lme4, WeMix and ordinal engines, so tidy(component = "fixed") returns the full broom shape – std.error, statistic, df, p.value, conf.low, conf.high – instead of estimates beside all-NA columns. The interval level follows summary(conf_level = ).

• A longitudinal summary() now also reports the intercept and slope VPCs of Bell et al. (2024, eq. 5) as $longitudinal$vpc_intercept and $vpc_slope. These exclude the occasion-level residual from the denominator, so they measure the between-stratum share of between-INDIVIDUAL variation in where a trajectory starts and how fast it changes, and are unaffected by measurement noise; the headline VPC keeps the residual and answers the discriminatory-accuracy question. vpc_slope is the only scale-free summary of how intersectional the rates of change are, and it pairs with the PCV_slope the decomposition already reported. Report the headline VPC unless you mean the trajectory question – ?summary.maihda_model and the longitudinal vignette set out the contrast.

• fit_maihda() and maihda() gained stratum_slope for longitudinal fits. stratum_slope = FALSE fits (time | id) + (1 | stratum): individuals keep their growth curves but strata differ in level only, so the between-stratum variance is constant over time and the decomposition reports PCV_intercept alone. It is the remedy when the stratum slope variance sits at the singularity boundary, which few strata, few occasions per stratum, or irregular measurement times routinely produce. The VPC still varies with time, through the person-level slope variance and the residual.

• maihda_interactions() gained scale = "response", which reports Evans et al.’s (2024, sec. 2.5.1) pi^B_j = pi_j - pi^A_j – a stratum’s total predicted outcome minus the outcome implied by its additive main effects alone – so a logistic fit’s interaction reads in percentage points of probability rather than log-odds, and a rope set with it is a smallest interaction of interest in those same units. The map from the BLUP is strictly increasing through zero, so flagged, direction and the p-values are identical on both scales and the interval is the exact image of the link-scale one rather than a simulated or delta-method approximation; se is dropped, the resulting interval being asymmetric about the estimate. An identity-link fit returns the same numbers either way, and a crossed-dimensions fit uses its dimension random effects as the additive baseline. print() on a non-identity-link fit now also shows that quantity beside the link-scale estimate, so the log-odds figure is never the only number on screen; the column is named prob_diff, count_diff or score_diff for what the model’s response scale actually is, which the link alone does not settle (a cumulative fit is logit-linked but scores categories rather than returning probabilities). The returned columns are unchanged, and scale = "response" remains how to obtain the quantity with its interval as data.

• fit_maihda() models now carry likelihood-adequacy diagnostics that print() and summary() surface alongside the singular/convergence caveats: Poisson/negative-binomial overdispersion (Pearson ratio) and zero inflation, stratum random-effect non-normality, and longitudinal residual autocorrelation. A well-specified model stays silent.

• fit_maihda() and maihda() gained count_approximation, selecting the log-link count level-1 variance behind the VPC from Nakagawa, Johnson & Schielzeth (2017): "lognormal" (default, unchanged), "delta" or "trigamma". The three agree above a marginal count of about 2 and diverge by a factor of six below it. Recorded on the fit, so summary(), the bootstrap intervals and the longitudinal VPC(t) all use the same one.

• summary() on a count model now reports which approximation produced the level-1 variance and the marginal count it was evaluated at, as count_vpc and in the printed output, and warns when that count (after any negative-binomial overdispersion) is at or below 2.

• print() on a maihda() analysis and on a summary() now reports the auto-binning recipe whenever make_strata() discretised a numeric stratum dimension: the cut-points used and the labels they carry. Those cut-points are quantiles of the analytic sample, so they define the strata and therefore the estimand, yet they were announced only by a message() at fit time – long gone by the time a saved model is printed in a later session. summary() objects now also carry them as strata_autobin_info.

• summary(df_method = "bootstrap") refers lme4 fixed-effect p-values and intervals to a null-restricted parametric bootstrap, the reference to use for a GLMM, whose Wald z is anticonservative for terms constant within a stratum. It costs n_boot refits per tested fixed-effect coefficient, and the intercept, having no null model to simulate from, is NA. tidy(component = "fixed") carries it through.

• maihda_proportional_odds_test() tests the proportional-odds assumption of a cumulative (clmm) fit by parametric bootstrap under the fitted model, redrawing the stratum random effects each replicate. Opt-in: every replicate refits two clm() models.

• maihda_describe() gained weights, the precision weights of fit_maihda(), written the same way (a bare column name) so one call describes and fits the same sample. Rows with a missing, zero, negative or non-finite weight leave the analytic sample as the engines drop them, and on a binomial model the weights also supply the denominator of a proportion response. It is the last argument, so positional calls written against 0.2.1 are unaffected.

• summary(bootstrap = TRUE) on a longitudinal lme4 fit now returns intervals for the two trajectory VPCs, as $longitudinal$vpc_intercept_ci and $vpc_slope_ci, reusing the refits the VPC(t) band already pays for. print() shows them and trajectory_vpc_method records the basis.

API changes

• fit_maihda(engine = "brms"), maihda() and compare_maihda_groups() now refuse a response addition term other than trials() and weights() – se(), rate(), cens(), trunc(), mi() and the like – which the VPC, predictions and summaries did not model: rate(expo) evaluated the count VPC at the rate per unit of exposure (0.18 where the same model gives 0.43), se(s) reported a VPC of 1, and cens() fits had censored values summarised as exact. Write an exposure as + offset(log(expo)): the same model for a Poisson, while for a negative binomial brms’s rate() also scales the shape by the exposure and the offset does not.

• A stratum the strata table holds but the fit never used – every one of its rows left the analytic sample, through a missing covariate for instance – is now refused by predict_maihda() like any other unseen stratum, and predicted at a zero random effect under allow_new_levels = TRUE. The wemix and ordinal engines returned that zero-effect prediction silently by default; lme4 and brms raised their own engine errors, which the package’s directed message replaces.

• predict_maihda() on the training rows now returns one value per analytic row for a fit made with na.action = na.exclude, matching nrow(object$data) and every other predict_maihda() return, rather than base R’s vector over the original input rows with an NA at each dropped position. The padded vector could not be bound to object$data, which is what produced the misalignment fixed below. Unchanged under the default na.action = na.omit.

• maihda_ic() reports NA criteria, with estimator "ML (glmer scale family: no likelihood)", for a glmer fit of a family with a scale parameter: a Gaussian with a non-identity link, Gamma, inverse Gaussian. lme4’s logLik() for those fits is not the marginal likelihood – it exceeded the achievable maximum by up to 15.4 on test fits – and nested AIC differences were off by 4 to 45.

• maihda_ic() adds back the saturated log-likelihood that lme4 leaves out of logLik() for a glmer fit with nAGQ > 1, so a quadrature fit’s criteria are the complete likelihood and compare with a Laplace fit’s. The two differed by 1351.6 AIC units on a 720-row Poisson fit; Bernoulli outcomes, whose term is zero, are unchanged.

• fit_maihda() now refuses family = "negbinomial" with nAGQ > 1 on the lme4 engine. lme4::glmer.nb() estimates theta on that incomplete likelihood and silently returned theta 0.087 instead of 1.62 and a stratum SD of 0; use nAGQ = 1, or fix theta with family = MASS::negative.binomial(theta).

• maihda_ic(), and the criteria compare_maihda() appends, now count only the estimated parameters of an lme4 negative-binomial fit with a fixed theta (family = MASS::negative.binomial(theta)): df is one lower, and AIC 2 and BIC log(n) lower, than AIC() and BIC() on the fitted model, which count theta for every negative-binomial family. glm() counts that family the same way; family = "negbinomial" fits estimate theta and are unchanged.

• plot_prediction_deviation_panels() on an engine = "brms" ordinal fit now draws the posterior-mean category probabilities from fitted(), the ones predict_maihda() uses, instead of tabulating responses simulated by predict(). The expected-score panel no longer changes between calls, and a category that no draw happened to produce no longer gets probability 0.

• summary() on a maihda() analysis now errors when passed df_method, bootstrap, n_boot, conf_level, response_vpc or seed instead of silently ignoring them. Those summaries are computed by maihda(); set them there, or summarise the fitted model directly with summary(x$model_adjusted, df_method = "bootstrap").

• maihda_discriminatory_accuracy() now reads non-unit integral lme4 weights= on a binomial fit as trial counts, which is what ?glm documents them to be, so all four spellings of the same data give the same AUC and the same n_case / n_control. Previously aggregation was inferred from the response carrying a value strictly inside (0, 1), so individual records collapsed to frequency cells – every row all-success or all-failure – took the observation-level path: each cell contributed one case and one control at the same score, and the AUC was exactly 0.5. Non-integral weights cannot be counts and are unchanged. New binomial_weights argument forces either reading.

• The automatic ordinal proportional-odds caveat is gone. It compared two fixed-only clm() models and referred the nominal-effects LRT to a chi-squared distribution, but a conditional cumulative model with a normal random intercept is not in general proportional-odds once that intercept is marginalised away, so the null was false under the correct model and the flag rate grew with n (about a quarter of correct models at a 7% stratum VPC and n = 96,000). The statistic is still computed, as $diagnostics$adequacy$marginal_po_proxy carrying lrt/df/n_terms only – no p-value, no flag, previously $proportional_odds. Use maihda_proportional_odds_test() to test the assumption.

• maihda_discriminatory_accuracy() no longer rounds an aggregated-binomial proportion response into whole successes. When the response times its trial counts is not a whole number – a malformed binomial glmer warns about as “non-integer #successes” – the fractional case/control mass now enters the weighted AUC as it stands, with a warning; rounding invented observations, moving both the AUC and the reported n_case / n_control. Well-formed successes/trials and cbind(successes, failures) fits are unchanged.

• The maihda_table() intercept row now carries an interval (*_lower/*_upper, previously always NA): the summary’s Wald interval for the likelihood engines, the credible interval for brms. The variance and SD rows remain point estimates.

• The fixed_effects element of a summary() object gained statistic, df, p_value, lower and upper columns (and se for lme4), and summary(conf_level = ) now sets the brms credible-interval quantiles as well as the Wald ones. Code reading term/estimate/se is unaffected; code assuming the exact column set is not.

• Gaussian lme4 fixed-effect p-values and Wald intervals now use a t on containment (between-within) degrees of freedom, reported in a new df column, instead of a z. summary(df_method = "normal") gives the z. GLMM, WeMix, clmm and brms are unchanged.

• Gaussian PCV calculations now use each model’s fitted (REML) between-stratum variance by default. Use estimation = "ML" in calculate_pcv(), stepwise_pcv(), pcv_importance(), compare_maihda_groups(), or maihda() to restore the previous ML-refit behaviour.

• maihda() and compare_maihda_groups() now reject a stratum dimension written in transformed form in the fixed part – y ~ factor(edu) + gender + (1 | edu:gender), and likewise scale(), I(), poly() or a spline. Only a bare column name was ever recognised as a dimension’s main effect, so the transform was not removed when the null model was derived: the null already adjusted for that dimension, its between-stratum variance was deflated and the PCV was computed against the wrong baseline (in decomposition = "crossed-dimensions" the dimension entered as a fixed effect and a random intercept at once). Transform the column in data and write the bare name.

• maihda() and compare_maihda_groups() now also reject a fixed interaction between a covariate and a stratum dimension – y ~ age * gender + race + (1 | gender:race) – in every decomposition mode. Only the dimension’s bare main effect was removed when the null model was derived, so age:gender survived into the null, which then already adjusted for gender. With the main effect gone that term is a gender contrast scaled by age, so the PCV depended on age’s arbitrary origin (re-centring moved it 2.3 percentage points on maihda_health_data while the adjusted fits were identical to 1e-9); for a categorical covariate the null’s fixed part spanned the dimension’s main effect exactly, leaving the PCV 6.5 points off the value a dimension-free null gives. In decomposition = "crossed-dimensions" the dimension entered as a fixed effect and a random intercept at once, competing for the same contrast, so its additive variance was not identified as intended; in decomposition = "longitudinal" a user-written gender * wave put the dim:time term the adjusted growth model supplies into the null as well. Write the additive form.

• pcv_importance(method = "sequential") is soft-deprecated. Use stepwise_pcv() for an order-dependent path or method = "shapley" for order-invariant attribution.

• maihda_mor() on a crossed-dimensions fit now returns the mixture MOR over pairs of distinct strata instead of applying the independent-strata closed form to the summed variance. Reported values fall, most where the variance sits in the additive dimensions.

• fit_maihda(engine = "wemix") now rejects a max_iteration that is not a single whole number of at least 1.

• The Poisson and negative-binomial level-1 variance behind the VPC is now evaluated at a single mean count, log(1 + 1/mean(lambda)), as Stryhn et al. (2006) and Nakagawa, Johnson & Schielzeth (2017) define it, instead of averaging log(1 + 1/lambda_i) over rows. Null-model count VPCs are unchanged; adjusted, offset, weighted, and longitudinal count VPCs rise.

• A brms longitudinal fit’s trajectory VPCs ($longitudinal$vpc_intercept and $vpc_slope) are now the posterior median of the per-draw ratio, with a credible interval in vpc_intercept_ci / vpc_slope_ci, instead of a ratio of posterior-mean variance components reported without uncertainty. The between-stratum variance posterior is right-skewed whenever the strata are few, so the plug-in ran high – by 6.6% on the twelve strata of maihda_long_data (0.6145 against a posterior median of 0.5765), and by more as the strata get fewer. The interval it never reported spans [0.34, 0.82] on that same fit. lme4 point estimates are unchanged.

• plot_prediction_deviation_panels(type = "auto") now plots a wemix binomial fit’s stratum probabilities, as every other engine does, instead of log-odds labelled “Fitted Value”: stats::family() has no WeMixResults method, so the family read as NULL and the panel took its Gaussian default. type = "gaussian" reproduces the old panel.

• A two-level outcome written as an EXPRESSION is now detected as binary, so fit_maihda() with no family= picks binomial for I(ly > 0.9) ~ x where it previously picked gaussian and silently fitted a linear probability model; the same outcome stored as a column always picked binomial. Pass family = "gaussian" for the LPM.

Documentation

• maihda_ic() now states the predictive target of the Bayesian criteria: brms WAIC/LOOIC condition on the fitted random effects and so assess prediction of new observations within the represented strata, not generalisation to a new stratum (a leave-one-group-out cross-validation question), while the likelihood engines’ AIC/BIC are computed from the marginal likelihood.
• maihda_interactions() now documents that its default BH flags are a conservative screen: partial pooling deflates a truly-null stratum’s Wald tail (null z variance is about the shrinkage fraction), so the flags under-flag rather than exceed the nominal false-discovery rate. The Bayesian path is described as already partially pooled instead of “multiplicity-free”, and the flag description drops causal phrasing.
• The finding-interactions vignette no longer calls the default BH flags “false-discovery-rate controlled”. It now matches ?maihda_interactions: the flags are a conservative FDR screen that under-flags rather than over-flags, not an exact error-rate guarantee.
• maihda_interactions() now names what its numbers are on a non-identity link. The stratum BLUP is a log-odds departure for a logistic fit, so a flag – or its absence – is a claim about multiplicative interaction only: a model additive in log-odds is not additive in probabilities, so the main effects alone already generate probability-scale interaction. print() says so on such a fit, the result carries a link attribute, and ?maihda_interactions plus the binary-outcomes vignette set out both the new scale = "response" reading and the linear-probability-model route to an additive-scale model.
• The design-weighted documentation no longer describes the fixed-effect standard errors as design-consistent for complex surveys. sampling_weights takes one person-level weight column and cannot represent PSUs, sampling strata, higher-stage weights, finite-population corrections, or replicate weights, so it delivers population-weighted point estimates rather than general design-based inference.
• ?plot.maihda_model now describes the stratum intervals correctly. Its description called them “confidence intervals”; they are the random effect’s own conditional interval – Wald on the BLUP’s conditional SD for the likelihood engines, the posterior interval of u_j for brms – shifted onto the outcome scale by a fixed-effect prediction held at its point estimate. They answer whether a stratum departs from the dashed across-strata reference, not how precisely its outcome is predicted, and are wider than a full interval for that outcome would be, because holding the fixed effects fixed stops the fixed-effect uncertainty they carry from cancelling against the random effect it is correlated with. The plot titles and the maihda_table() footnote already said “conditional”.
• ?predict_maihda no longer calls the allow_new_levels = TRUE prediction a “population average”. Setting an unseen stratum’s random effect to zero gives the conditional prediction at u = 0, which equals the response-scale marginal mean only under an identity link: under a log link with stratum variance tau^2 the marginal mean is exp(tau^2/2) times larger – 38.6% on a fit with tau^2 = 0.65 – and a logit fit is attenuated towards 0.5. All four unseen- or missing-stratum error messages that offer the escape now say the same. Predictions themselves are unchanged.
• ?summary.maihda_model and ?maihda_proportional_odds_test no longer call their parametric-bootstrap p-values exact. The first said the df_method = "bootstrap" p-value is exact whenever (n_boot + 1) * alpha is a whole number, the second that correctly specified data are rejected at the nominal rate by construction. Both simulate from a fit estimated on the same data, so each is an approximation that no replicate count makes exact – a binomial fit with 4 strata rejected a true null about 14% of the time at a nominal 5% – and n_boot or n_sim sets only the Monte Carlo resolution. The interval’s agreement with the p-value is now described as algebraic rather than as a coverage guarantee, and more draws as converging on a fixed interval width rather than tightening it.
• ?maihda_proportional_odds_test no longer says that symmetric thresholds make the marginal slopes coincide and the fixed-only chi-squared statistic valid. Symmetry equates the marginal slopes at thresholds -c and +c where the location is zero but not elsewhere – -0.882 against -0.826 one unit away, at c = 1, a unit coefficient and stratum SD 1 – whether fitted thresholds look symmetric depends on how the covariates are coded, and observations that share a stratum stay dependent. With three categories cut at -1 and +1, a covariate symmetric about zero, 12 strata and a stratum VPC near 7%, the chi-squared reference still rejected 10% of datasets simulated from the correctly specified model at a nominal 5% with n = 96,000.

Performance

• maihda(decomposition = "crossed-dimensions") now fits the model once instead of twice. The preliminary pass that resolves the strata and family no longer refits the supplied formula only to discard it, which roughly halves the fitting cost (most noticeable for brms).

Bug fixes

• A brms trials() term that calls a function outside base R – your own, or one from an attached package such as coalesce() under library(dplyr) – no longer loses its trial counts: predict_maihda(scale = "response") returned expected success counts, the stratum summaries, plots, maihda_table() and maihda_describe() counted every row as one trial, and the AUC and the prediction-deviation panel refused the fit. A response calling such a function no longer leaves that panel with no row to score.

• calculate_pcv(), compare_maihda() and maihda_ic() now compare the trial counts of a brms y | trials(n) outcome, not just its successes, so fits of the same successes out of different trials are refused or flagged as different samples instead of compared as one.

• The same three functions now compare the weights of a brms y | weights(w) term, as they already compared lme4 weights =; brms fits with different weights passed as equally weighted.

• plot_prediction_deviation_panels() now draws every case of a case-level binomial panel. Its point shape showed the observed outcome, which a cbind() response, a proportion with trial weights, a brms trials() fit or data without the outcome does not have, and ggplot dropped each such point; those cases now take a fixed shape, and Bernoulli panels are unchanged.

• predict_maihda(newdata = ) and plot_prediction_deviation_panels(data = ) on an engine = "brms" fit no longer refuse rows carrying its two-level outcome as it was fitted – a factor, character labels or a logical. fit_maihda() hands brms that outcome recoded to 0/1, and brms checked supplied rows against the recoded column although no prediction reads it; the rows brms sees now leave it out.

• The binomial panel of plot_prediction_deviation_panels() now codes the outcome of supplied data as the model coded it, by label. It re-derived the 0/1 coding from the plotted rows, so the fitted rows handed back with their factor levels reversed had every residual scored against the wrong outcome and every “Wrong” / “Correct” swapped, and rows holding a single outcome value could not be coded, so they were scored as perfect fits.

• The binomial panel of plot_prediction_deviation_panels() now scores each row from its own successes and trials. Given data whose outcome is not one 0/1 value per row – a cbind() response or a proportion – it borrowed the fitted rows’ deviance residuals whenever data had as many rows, so reordering the rows moved each residual onto another row, and it gave every row of a subset a residual of 0. A 0/1 outcome with prior weights was scored as a single trial whatever its weight: five times too small at 25 trials a row, on the default panel too, and unlike the same model written with cbind(). A prior weight now counts as that many trials, as in R’s binomial family, and a row whose outcome or trial count is unknown is left out of the residual summaries and the labels rather than scored as a perfect fit.

• plot_prediction_deviation_panels(data = ) now weights each row of data by its own weight in the stratum summaries. The fitted rows’ weights were paired with the supplied rows by position whenever the counts matched, so reordering the fitted rows moved the weighted stratum means, new rows took the training rows’ weights, and a subset was silently unweighted. A row’s weight is read from its column in data, or from the fitted row with the same row name and prediction; if any row’s weight cannot be found, every row is weighted equally, with a warning. The weights of a MASS::polr() or bare ordinal::clmm() fit, which stats::weights() does not return, now weight its panel too; they were ignored.

• fit_maihda(engine = "wemix") no longer refuses a Bernoulli outcome written as an expression, such as I(ly > 0.9) ~ x, as though it were an aggregated binomial: maihda_analytic_response() returned NULL for every non-symbol response, so the is-binary test was a false negative. The fit now matches the precomputed-column spelling to the last bit, and maihda_describe() reports such an outcome as binomial rather than continuous. A character or 1/2-coded expression is still refused, now saying so and why: the 0/1 recoding reaches bare columns only.

• maihda_discriminatory_accuracy() on a brms fit whose outcome is an expression – an aggregated binomial such as round(raw) | trials(ntr) ~ ... – now reads the response the model was fitted on rather than the first variable named in the formula. brms keeps every raw input column beside the evaluated response in the frame it stores, so all.vars(formula)[1] found a real column holding different values and the wrong one was read silently. Measured on one 240-row fit, read with two spellings of that column: sitting 0.4 above the success counts it exceeds the trial count on each all-success row and stopped the call outright with an internal error about negative case/control mass, while sitting 0.4 below them it stayed in range and nothing complained – an AUC of 0.6587 against a correct 0.6562, with fractional totals of 909.6 cases and 1456.4 controls for data holding 1002 and 1364. The lme4 engine reads its response from the fit itself and was never affected; the wemix and ordinal engines refuse a non-symbol response at fit time, and a bare-symbol outcome returns the same values on every engine.

• plot(type = "obs_vs_shrunken") now puts the observed stratum means on the scale the model was fitted on when the response is a transformed expression. The observed outcome was read as the first variable named in the formula, taken raw out of the data, so a wemix fit of log(y) ~ x drew raw y on the x-axis against log-scale shrunken estimates on the y-axis – 3.539 where the fitted scale gives 1.100, and up to 3.0 out across the strata of that fit – which left the panel’s y = x diagonal, its only reference, meaningless. The identical model written ly ~ x with ly <- log(y) precomputed plotted the correct values, so two spellings of one fit disagreed while their coefficients and shrunken estimates were identical. The response expression is now evaluated as maihda_describe() already evaluated it, and reconstructed through model.frame() so it matches the lme4 route in class as well as value. Only wemix could reach this: lme4 and brms store a frame that carries the evaluated response, and the ordinal engine refuses a transformed response outright.

• The proportional-odds statistic of a cumulative (clmm) fit no longer depends on how a covariate was spelled in the formula. It was computed by refitting the original term labels against the model frame, which stores each variable already evaluated under its deparsed name, so y ~ log(x) sent R looking for a column x the frame does not have: log(x), I(x^2), scale(x), poly(x, 2) and factor(g) all lost the statistic, and maihda_proportional_odds_test() stopped with a message blaming a covariate-free model, while the identical fit written with the transformation precomputed kept it. Where a same-named object of the same length happened to be reachable from the package namespace, the statistic was silently computed from that object instead of from the fitted data. The refit now takes the fitted columns from the model matrix, reduced to full rank so a covariate entered twice or a constant one no longer costs the statistic either; the statistic is unchanged for the spellings that already worked, and a failed refit is reported separately from a covariate-free model.

• plot_prediction_deviation_panels(), and the prediction_deviation panel of plot(type = "all"), now work on a wemix fit whose response is an expression rather than a bare column. WeMix’s own predict() method names the response with as.name(form[[2]]), which is undefined for a call, so log(y) ~ ... stopped with “‘language’ object cannot be coerced to type ‘symbol’” and plot(type = "all") reported the panel as uncomputable and omitted it, while the same model spelled ly ~ ... over a pre-computed column drew normally: one fit, two spellings, a plot and an error. The panel now builds the linear predictor from the fit itself, as predict_maihda() already did for this engine, which also reuses the fitted transformation basis and factor coding; every case that drew before returns identical values, and a bare WeMix fit handed prediction data it cannot use for such a response now names the response and the way out instead of passing WeMix’s message on. lme4, brms and ordinal fits were never affected.

• predict_maihda(allow_new_levels = TRUE) now gives an unseen stratum combination the documented zero-random-effect prediction even when its label reads like an internal stratum id. Combinations are labelled from their dimension values while strata are numbered 1, 2, …, so a one-dimension value of "1" – or "12" under make_strata(sep = "") – named a fitted stratum, and lme4, WeMix, ordinal and brms all returned that stratum’s random effect instead: 1.508 against a correct 2.165 on one 20-stratum fit. Unseen combinations now carry generated ids checked against every id the model holds. The same collision defeated the check that a supplied stratum column agrees with the dimension columns beside it, including under the default allow_new_levels = FALSE, and such a row is now refused.

• summary(df_method = "bootstrap") now refers each fixed-effect coefficient to its own null rather than to its whole term’s. A term spanning several design columns – a factor with three or more levels, a polynomial, an interaction between factors – had every one of its columns deleted from the fit the reference draws were simulated from, while each row reported the p-value and interval of a single one of them. The deleted siblings’ effect moved into the variance components: on a binomial fit with a three-level dimension the reference for the c contrast came from a null whose stratum variance was 1.458, where the null that row states gives 0.299. The answer therefore depended on how the design was spelled – y ~ f and y ~ fb + fc are the same fitted model – and it no longer does: every spelling now simulates from the same restricted fit, and where the columns are also in the same order the two return identical numbers. The siblings are kept and re-estimated as the nuisance parameters they are, so the restriction is the one the coefficient names under the fitted contrasts. n_boot is now spent per tested coefficient rather than per term, so a k-level factor costs k - 1 blocks of refits; a model whose fixed terms each span one design column – every binary dimension, every continuous covariate – is unchanged.

• A model fitted with na.action = na.exclude no longer mixes predictions over the original input rows with summaries over the analytic ones. Base R pads predict(), fitted(), residuals() and weights(type = "prior") on such a fit back out to every input row, while object$data holds only the rows the engine kept, so on a 720-row frame with 717 analytic rows the stratum predictions stopped with “arguments imply differing number of rows”, maihda_discriminatory_accuracy() with “‘prob’ and ‘y’ must have the same length”, maihda_table() dropped its ranked-strata table through its error handler, plot(type = "all") left out four panels, and summary(bootstrap = TRUE) stopped with “All VPC bootstrap refits failed” because simulate() pads its draws the same way. Where the analytic count divided the input count R recycled instead, silently: stratum predictions were displaced by up to 0.016 on the probability scale and every stratum size was doubled. The zero-inflation adequacy check disappeared, and a growth count fit’s count_vpc reported a marginal count of 7.05 instead of 7.83. na.action = na.omit, the default, was never affected, and the fits themselves were always correct. Reachable by the argument, by an abbreviation of it, or from options(na.action = "na.exclude"); lme4 only, the ordinal engine having never carried the padding and WeMix refusing the argument.

• A binomial model, or an unweighted Poisson or negative-binomial model, fitted with na.action = na.exclude can be bootstrapped again. lme4 reads the fit’s prior weights itself when it simulates, and under na.exclude that vector is padded back out to the original input rows with an NA at each dropped position. The binomial simulator tests it with any(wts %% 1 != 0) and the Poisson and negative-binomial ones with any(wts != 1); the padded NA makes the binomial test NA for any whole-number weights and the other NA whenever all the real weights are 1, so the call stopped inside lme4 with “missing value where TRUE/FALSE needed” before any draw existed – and stripping the padding from the result, all the fix above could do, cannot repair a call that never returned. summary(bootstrap = TRUE), calculate_pcv(bootstrap = TRUE), the crossed-dimensions VPC and the longitudinal VPC(t) band were all dead on such a fit, and an aggregated cbind() binomial would have drawn its trial counts from the padded vector even without the NA. Gaussian fits and weighted Poisson and negative-binomial fits were never affected, and neither were summary(df_method = "bootstrap") and pcv_importance(), which simulate from complete-case refits. The draws now come back identical to the same fit under na.action = na.omit.

• compare_maihda_groups() now refuses an argument the engine does not take – weights, subset or offset on wemix, brms or ordinal, or sampling_weights together with weights – once, before any group is fitted. It returned a comparison table in which every group carried a “fit failed” status and a warning, repeating the refusal fit_maihda() makes up front. Five engine incompatibilities are refused the same way: wemix without sampling_weights, ordinal with sampling_weights, engine = "ordinal" on a non-ordinal outcome, an unsupported wemix family, and a family name no engine can resolve.

• A partial spelling of weights, subset or offset that the engine’s own fitting function cannot bind is now refused like the full name. offs = o on engine = "ordinal" fitted every row without the offset, warning only that an unknown control element was ignored; WeMix reported “unused argument” with the whole vector in the message, and brms only after compiling its Stan model. A spelling the engine does bind is untouched (w = 500 is brms’s warmup), as is the lme4 path.

• engine = "wemix" now rejects an offset() term in the formula, in fit_maihda(), maihda() and compare_maihda_groups(). WeMix::mix() leaves the offset out of the model it fits, so Gaussian and binomial fits returned the coefficients, variance components and VPC of the model without it (a VPC of 0.80 instead of 0.23 on one Gaussian fit), while predictions added the offset back. For a Gaussian outcome, fit the response minus the offset, which is the same model; for a binary outcome, use engine = "brms" with sampling_weights. Fits saved with an offset are not flagged; refit them.

• The binomial panel of plot(type = "prediction_deviation") and plot_prediction_deviation_panels() now draws an engine = "brms" y | trials(n) fit as per-trial probabilities. It drew brms’s expected success counts on its probability axis, with the stratum intervals clamped at 1, and it now ranks the strata by the deviance residual of their successes out of trials, which was 0 for every stratum. An lme4 cbind(successes, failures) fit no longer stops the panel with “obs_outcome must be size”, and plot(type = "all") no longer drops it.

• maihda_ic() now reports a delta between Poisson and negative-binomial fits of the same counts, and between binomial or cumulative fits that differ only in their link. It withheld the delta whenever the family or link differed, as the VPC and PCV must, although these likelihoods are on a common scale. A continuous family still needs the same family and link, and the printed comparability note now says so.

• predict_maihda() and predict() now refuse newdata for an lme4 fit that has both a formula offset() and an external offset =, as they already did for an external offset alone. The formula offset turned the check off, so predictions on new data silently left the external offset out.

• The ordinal surprise panel of plot(type = "prediction_deviation") and plot_prediction_deviation_panels() now finds each row’s observed category by its position among the model’s fitted categories, whatever their labels. It matched the label against probability columns named 1, 2, 3, so an engine = "ordinal" outcome coded low < mid < high was never scored and the panel was empty, while one coded 0 < 1 < 2 dropped its first category and scored the rest one category off, and one coded 3 < 2 < 1 had its end categories swapped. An engine = "brms" fit never found its response at all, which also gave every stratum of a brms binary-outcome panel a deviance residual of 0. A category named like one of the panel’s own columns (n, weight) is no longer overwritten by it, and a row whose category is not a fitted one is left out with a warning.

• plot(type = "context_vpc") now draws a crossed-dimensions fit that carries a context, showing the between-stratum variance as its additive dimension and interaction bars. It stopped with a message asking for the context = that had been given, and plot(type = "all") warned and left the panel out.

• plot(type = "effect_decomp") now leaves a contextual (context =) random effect, or any other grouping besides the stratum and dimension random effects, out of the decomposition, as plot(type = "predicted") already did. In crossed-dimensions mode each stratum’s context composition was drawn as part of the additive dimension component, and in both modes the global mean the deviations are measured from carried the row-weighted mean of the context effects.

• fit_maihda(), maihda() and compare_maihda_groups() now treat a partial spelling that the engine binds to subset, weights or offset (e.g. subs = keep, weig = w) like the full name. The engine bound the partial spelling while the package’s checks read the exact name, so engine = "ordinal" fitted a subset or weights it refuses, lme4 fits could detect the wrong family, strata or longitudinal time centring, and a partial contrasts name escaped the ordinal coding record. Supplying one of the three more than once, under any spellings, is now an error.

• engine = "wemix" now rejects WeMix::mix()’s center_grand and center_group. WeMix centred the covariates inside the fit while predictions, stratum tables, plots and binomial standard errors used the uncentred values; centre covariates in data instead. A saved fit that used them now refuses to predict, and a saved binomial maihda() analysis warns that its stored stratum standard errors and interaction tests are wrong.

• A partially named max_iteration passed to engine = "wemix" (e.g. max_iter = 0) is now validated like the full name instead of reaching WeMix unchecked.

• engine = "ordinal" and engine = "wemix" predictions now rebuild the fixed-effect design with the factor levels and contrast matrices the fit used. The rebuild took the contrasts in force at prediction time and every level declared on the stored data, so a fit coded by contrasts = (ordinal), by options(contrasts = ) at fit time or by a factor’s contrasts attribute failed with “missing column(s)”, or returned wrong values silently when a custom matrix’s columns were named like treatment coding’s; an ordered covariate with a declared but unobserved level was mis-coded even under the default options. This reached predict_maihda(), the stratum tables and plots, maihda_proportional_odds_test(), and for a binomial WeMix fit maihda_discriminatory_accuracy() (with the AUC that summary() and stepwise_pcv() report) and the summary() stratum standard errors. A newdata level declared on the data but absent from the fitted rows is now refused like any other unseen level, as lme4 refuses it.

• summary(df_method = "bootstrap") no longer warns that the n_boot = 199 its own help page recommends is too low. The low-replicate warning is written for a percentile interval, which the fixed-effect one is not: that interval is symmetric about the estimate at a critical value taken from a single |t*| order statistic, so what matters is how many draws lie beyond the cut-off, and that depends on the level rather than on n_boot alone. The check is now level-aware and separately worded, naming the rank: 99, 199 and 999 pass at the 10%, 5% and 1% levels, being the smallest counts that put ten draws beyond the cut-off, while 199 draws at conf_level = 0.99 (two beyond it), and any count too small to reach the level at all (an unbounded interval), still warn. The VPC, PCV, group-comparison and importance bootstraps keep the percentile warning unchanged.

• A summary(df_method = "bootstrap") interval now excludes zero exactly when its p-value is at most 1 - conf_level, at every level. The critical value’s rank was ceiling((1 - conf_level) * (n_boot + 1)) - 1, and 1 - 0.95 is stored just above 0.05 while 1 - 0.90 is stored just below 0.1, so wherever that product is a whole number the rule was p <= alpha at levels whose complement is stored above its decimal value (0.95, 0.99, 0.70, 0.85) but p < alpha at those stored at or below it (0.90, 0.80, 0.75, 0.50): with n_boot = 199, a p-value of exactly 0.05 excluded zero at 0.95 while one of exactly 0.10 did not at 0.90. Intervals at 0.95 and 0.99 are unchanged; the levels that took the strict rule move by one order statistic where the product is whole.

• An engine = "ordinal" formula whose right-hand side ends in a subtraction no longer refuses to fit. ordinal::clmm() reads a formula offset by position in the bar-free variables list, so the offset is moved ahead of the random-effect terms before the fit, but the relocation split the right-hand side on top-level + only – deliberately, so that nothing is promoted out of a subtraction. A right-hand side whose outermost operator is - therefore presented itself as one operand and was never rewritten, and the alignment backstop refused it. - 1 is what update() leaves behind for every no-intercept formula, so maihda(ord ~ 0 + x + a + b + offset(log(expo)) + (1 | a:b)) reached the engine as ... + (1 | stratum) + offset(log(expo)) - 1 and could not be fitted at all, though it is perfectly relocatable; it now fits, and to the same log-likelihood as the intercept-coded spelling. The refusal was never gratuitous – handed that formula as written, clmm() fits the stratum id as its offset, dropping the log-likelihood from -596.76 to -619.83 and inflating the between-stratum variance from 0.154 to 3.341 – so the relocation still declines to lift an offset out of a subtraction, where moving it would change the fixed design, and the alignment guard still refuses any formula it cannot align. That guard is a backstop rather than a spelling users meet: the strata resolution rewrites such a formula into its canonical + form before the engine sees it.

• An adjusted model written without an intercept no longer changes the derived null model. y ~ x + a + b + (1 | a:b) and y ~ 0 + x + a + b + (1 | a:b) are the same fit – with no intercept R codes the first factor as cell means, so the two designs span one column space and the two adjusted fits agreed to 3.5e-12 – but the null was derived with update(. ~ . - a - b), which in the second spelling took the grand mean out along with the dimensions. The null’s stratum random intercept then absorbed the outcome mean: on a 900-row Gaussian example the null between-stratum variance rose from 3.01 to 1941.78, the null VPC from 0.402 to 0.998 and the PCV from 0.600 to 0.999, silently. The same reduction feeds every mode, so the crossed-dimensions partition (dimension a’s additive variance 1.94 against 1940.64, additive share 0.650 against 0.999), the longitudinal PCV_slope (0.202 against 0.993), a binomial PCV (0.676 against 0.954) and every group’s PCV in compare_maihda_groups() moved with it. The grand mean is now restored whenever the reduction loses it and the source model carried it, reported once. A reduction that still spans the intercept another way – a surviving factor covariate coded as cell means, or a fixed part that collapses to 1 once the bars are stripped – keeps its parameterization, and a model that genuinely has no intercept in its span (numeric dimensions fitted through the origin) is left alone, because a grand mean added there would put a column in the null that the adjusted model does not have.

• A MAIHDA fit whose fixed part cannot represent the outcome’s mean now warns. With a 0 + or - 1 formula whose remaining terms are all numeric there is no column to carry the mean, so the stratum random intercept absorbs it and the between-stratum variance, VPC, MOR and PCV describe where the outcome sits rather than how it varies across strata: maihda(y ~ 0 + x + (1 | a:b)) reported a null VPC of 0.998 and a PCV of 0.999 in silence. The check reads the fitted design rather than the formula, so a - 1 spelling whose first factor is coded as cell means stays silent, and a cumulative model is exempt: its free thresholds are the intercepts, so it has no intercept column to miss and 0 + changes nothing there – clmm() returns the same log-likelihood and the same between-stratum variance either way.

• The longitudinal plot(type = "trajectories") now draws each stratum’s own fixed-part trajectory instead of one shared curve. The fixed part was evaluated once at the mean/modal covariate profile and added to every stratum’s random intercept and slope, but that profile replaced the stratum-defining dimensions too, so an adjusted growth model’s dimension main effects and dim:time interactions were held at the modal stratum’s values on every line: on the twelve strata of maihda_long_data the plotted baselines spanned 4.51-4.94 where the correct predictions span 4.06-5.97, a discrepancy of up to 2.32 over the time grid, and the only line drawn correctly was the modal stratum’s own. The dimensions now take each stratum’s values – the reconstructed tertile factor for an auto-binned numeric dimension, or that stratum’s mean when a hand-written formula enters one as a linear term – while every other covariate stays at the shared reference profile. A null growth model, whose fixed part carries no dimension terms, is bit-identical, as is a purely covariate-adjusted fit. A dimension reachable only through a transformed term (factor(gender) with no bare gender) cannot be set per stratum from the stored model frame and now warns rather than being frozen silently.

• An engine = "ordinal" fit whose formula carries an offset() term no longer fits the stratum column as its offset. ordinal::clmm() reads the offset by position in the bar-free variables list while indexing the frame it built from the barred one, so the two only line up while the offset precedes (1 | stratum) – and the resolved formula put it last, both through the (1 | var1:var2) shorthand and when written that way by hand. The mis-index landed on the integer stratum id: on a 900-row five-category fit the cut points moved from (-0.83, -0.01, 0.89, 1.79) to (4.85, 5.65, 6.54, 7.42), the between-stratum variance rose from 0.47 to 13.50 – the variance of the ids the random intercept had to cancel – and the log-likelihood fell 41.3 units, with no warning. maihda() errored outright instead, its adjusted formula putting a character dimension column at that position. The offset is now moved ahead of the random-effect terms for the clmm call – including a parenthesized one, (offset(off)) or (z + offset(off)), which terms() also treats as an offset and which silently dropped it instead – and a spelling that cannot be moved without changing the fixed design (an offset inside a -, * or : subtree after the random effect) is now refused with a message rather than fitted against the wrong column. An offset-free ordinal fit is bit-identical, and the lme4 and brms engines, which read the offset from their own model frame, are untouched.

• An engine = "brms" cumulative fit now honours the threshold (and link_disc) set on a brms::cumulative() family. The family was rebuilt from its link alone before handing it to brm(), which reset every other option to brms’s defaults, so family = brms::cumulative(threshold = "equidistant") fitted the flexible model and said nothing – a different model from the one asked for. The cut points themselves were always right: brms expands a constrained threshold into the full b_Intercept[1..K-1] vector, so the thresholds summary() reports are unaffected either way. A threshold brms cannot honour is now an error rather than a silent fallback. family = "ordinal" and maihda_cumulative() carry no such options and are unchanged.

• An engine = "ordinal" fit made with a structured threshold – threshold = "equidistant", "symmetric" or "symmetric2", passed through to ordinal::clmm() – no longer reads the free threshold parameters as if they were the cut points. clmm stores those two things in the same $alpha slot only under the default threshold = "flexible"; an equidistant five-category fit holds an anchor and a spacing there, so the response scale, the summary threshold table, the stratum response predictions, the proportional-odds bootstrap and the prediction-deviation panels all treated a five-category model as three-category. Rebuilt category probabilities were off by up to 0.51, more than a quarter of rows had an observed category with no column at all, the expected category score was capped near 3 instead of 5, and the bootstrap calibrated a five-category statistic against simulated three-category ones. The cut points now come from the fit’s $Theta, cut-point standard errors from the corresponding transform of the parameter covariance, and the number recovered is checked against the fitted category count. threshold = "flexible" is bit-identical.

• A longitudinal fit whose fixed part contains a transformed covariate – scale(z), log(z), poly(z, 2) – no longer fails with “object ‘z’ not found”. The count VPC(t) grid and the fixed trajectory are built from the stored model frame, which holds such a term only under its derived name and never as the raw z, so the term could not be re-evaluated there; it now falls back to the values the fit stored for it, per fitted row for the VPC(t) grid and at a representative value for the trajectory, as a formula offset() term already did. A grid carrying every raw variable is unchanged.

• Predictions from a fit whose fixed part contains a data-dependent transformation – scale(x), poly(x, 2), splines::ns(x, 3) – now evaluate it on the basis the model was fitted with. The lme4 fixed-effect, WeMix and ordinal prediction helpers rebuilt terms() from the bare formula, which carries no predvars, so the centre, scale or knots were recomputed from whatever rows were being predicted: predicting every row then subsetting disagreed with predicting the subset (by up to 1.1 on the ordinal link scale and 1.2 for a Gaussian WeMix fit here), and a one-row or otherwise constant grid divided by a zero standard deviation. A longitudinal count VPC(t) fixes every row to a single time, so a fit containing scale(<time>) reported VPC = NaN at every time. Untransformed fixed parts, and transformations that do not depend on the data such as log(x), I(x^2) and poly(x, 2, raw = TRUE), are bit-identical.

• summary(bootstrap = TRUE), summary(df_method = "bootstrap"), calculate_pcv(bootstrap = TRUE) and the longitudinal VPC band now run on an lme4 fit whose data carried missing values. lme4::refit() treats a simulated response with no na.action attribute as being on the pre-na.omit data and drops those rows a second time, so every draw failed and the call errored with “All … refits failed”; the bundled maihda_health_data has 208 such rows. A fit with no dropped rows is unchanged.

• summary(df_method = "bootstrap") now imposes the null it advertises for a term that is marginal to a higher-order term still in the model. The reference distribution was simulated from update(formula, . ~ . - <term>), but R’s marginality rules recode the surviving term to absorb the dropped one – for y ~ x * f the result is f + x:f, whose full dummy expansion spans exactly the original column space – so the coefficient under test was never restricted and the p-value was pinned near 0.5 whatever the effect size: it rejected at 0.0000 against a nominal 5% over 300 replicates, and had power 0.0000 against an effect the repaired test detects 97.5% of the time; the repaired reference rejects at 0.065 (200 replicates, Monte Carlo SE 0.017). The reduction is now verified against the fitted design and, where it restricts nothing, the term’s design columns are constrained directly. Affects the main effects of x * f and f * g, every main effect and two-way term under a three-way interaction, a nested f / g, and transformed terms such as poly(x, 2) under an interaction. An additive fixed part is unchanged.

• print() on an engine = "ordinal" fit no longer dumps the whole analytic data set. The clmm call embedded the data frame itself, and ordinal’s print/summary methods deparse it, so printing a 3,000-row fit ran to about 2,100 lines; deparse(getCall()) on the same fit echoed clmm’s source. The call now names the frame, as the lme4 and brms engines already did.

• A bootstrap PCV/VPC interval whose refit optimiser did not converge on more than half the contributing draws is now flagged interval_reliable = FALSE with an escalated warning, and print() discloses the retained non-converged count it already tracked. Such draws are still retained – the count is the optimiser’s own return code, not lme4’s false-positive gradient flag, and it fires on well under 1% of refits in practice – but the interval could previously clear the success-fraction gate unremarked at any non-convergence share, since those draws are finite and counted as successes. The flag and count reach calculate_pcv(), summary(), the crossed-dimensions, contextual and longitudinal decompositions, and pcv_importance().

• maihda() now warns that seed does not seed the brms sampler. It is the response-scale VPC simulation seed and is a formal argument, so it never reaches brms::brm() through ...: maihda(engine = "brms", seed = 1) gave a different posterior on every run with nothing to say so. Call set.seed() before maihda(), or use fit_maihda(seed = ), which has no such formal and does pass it through.

• maihda_describe() and the observed-vs-shrunken plot now read the denominator of a brms success | trials(n) outcome off the trials() term. Both took only the success counts, so the denominator defaulted to 1: a four-row 26-of-62 sample was reported as “26 events / 4 trials (650.0%)”, and the plot put mean success counts on the x-axis against predicted probabilities on the y-axis. The cbind(success, failure) form was always correct.

• A constant brms trial count (y | trials(20)) is now recycled to one value per row. The length-1 vector failed every caller’s length guard, so such a fit got unit prediction weights, predict_maihda(scale = "response") silently returned expected counts (9.7 to 23.2 out of 40 trials) where it documents a per-trial probability, and maihda_discriminatory_accuracy() errored on the aggregated response rather than computing the count-weighted AUC. A trial count supplied as a data column was always handled.

• maihda_describe() now reads the trial counts of an aggregated binomial supplied as weights = rather than as cbind(successes, failures), which it previously took to be 1 per row: a 12-row, 340-of-617 sample was reported as “6.420795 events / 12 trials (53.5%)” – a fractional event count over a row count – where the cbind() spelling of the same data gave “340 events / 617 trials (55.1%)”. It applies the same rule maihda_discriminatory_accuracy() does, so the two can no longer report different sample sizes for one model.

• Parametric bootstrap intervals on a weighted Gaussian lme4 fit now simulate the residual as sigma / sqrt(w_i). lme4::simulate() draws equal-variance noise for every row while refit() keeps the 1/w_i weights, which inflated the bootstrap residual variance by about mean(w) and pushed VPC, PCV, context, crossed-dimensions and longitudinal intervals off their own point estimates. Only all-unit weights were unaffected; a constant weight c != 1 was biased by c. Unweighted fits and every GLMM are unchanged, as is pcv_importance(), which takes no lme4 weights argument and so never fitted a weighted Gaussian model.

• The fixed-cell-interaction guard now sees a transformed dimension. y ~ factor(a) * b + (1 | a:b) passed the guard because factor(a) did not match the dimension name, so the fixed factor(a):b cell means survived into both derived formulas, saturated the strata and left the stratum variance unidentified (a degenerate Hessian, and a PCV of essentially zero).

• maihda_interactions() on a model carrying a transformed dimension now says so instead of reporting it as a null model, which was true in effect but named neither the cause nor the remedy.

• A singular fit is no longer reported as a convergence failure. lme4 files “boundary (singular) fit” among its own convergence messages, so a fit whose optimizer returned code 0 was recorded as converged = FALSE and printed twice – once as “Singular fit”, once under a “Convergence warnings reported by lme4” header – which also double-counted such groups in compare_maihda_groups()’s singular and non-converged warnings.

• The singular-fit report now names the random-effects block that is at the boundary and scopes its VPC/PCV caveat accordingly. It previously said the between-stratum variance “may be unreliable” for every singular fit, which is false when the boundary sits in a non-stratum block – routine in a longitudinal fit whose (time | id) block has no person-level slope variation.

• PCV_slope is now NA when the null model’s between-stratum block is rank-deficient (a perfect intercept-slope correlation), with a note saying why. The stratum variation has collapsed onto a single direction in (intercept, slope) space, so the slope variance is not a free parameter and the null and adjusted fits need not have collapsed onto the same direction – the ratio does not compare the same quantity before and after adjustment. It was not merely uncertain but explosive: across replicates of an irregular design with no true stratum slope variance, rank-deficient fits returned values from -196 to +1. The existing denominator guard could not catch it, because it scales a slope variance (outcome-units squared per time-unit squared) against the residual variance, so its relative threshold has no fixed meaning for that cell. PCV_intercept and PCV(t) are unaffected and still reported. print() also gained the adjusted-model boundary note the cross-sectional calculate_pcv() has carried all along.

• print() on a longitudinal maihda() now shows both growth fits’ diagnostics, labelled. Only the null model’s were shown, so a singular or non-converged adjusted fit – which pins the reported additive share – delivered its headline PCV in silence.

• print() on a cross-sectional two-model maihda() now shows fit diagnostics, which that branch printed for neither model: a non-converged fit, or one singular in a non-stratum block such as a context random intercept, produced its headline PCV in silence. A stratum-block singularity in the ADJUSTED fit is deliberately not banner-warned – it is the expected shape of an additive decomposition – and is reported instead by the boundary note below, so the diagnostics block stays a signal rather than firing on nearly every healthy analysis.

• The “PCV is pinned near 100%” boundary note now reaches print() on a maihda() analysis. calculate_pcv() has recorded adjusted_at_boundary all along and print.pcv_result() showed the caveat, but the analysis print formats its PCV inline rather than calling print(x$pcv), so a “100% additive” headline read off a singular adjusted fit carried no caveat anywhere the reader looks. Both methods now share one note.

• The per-group variance extractors behind the contextual and crossed-dimensions partitions now enforce the intercept-only contract the ordinary summary already applied. A contextual lme4 fit with (1 + x | stratum) previously returned a partition whose total silently dropped the slope variance and intercept-slope covariance, and a slope-only context (0 + x | site) produced an all-NA partition; the brms path accepted a slope-only sd_<group>__<slope> column, squaring slope draws as the group’s intercept variance.

• maihda(decomposition = "crossed-dimensions") now rejects a random slope on an allowed grouping factor, e.g. (1 + x | stratum), instead of silently rewriting it to the canonical intercept-only crossed formula. The builder previously checked only the grouping side of each random-effect bar, so the slope vanished without a warning through both maihda() and the per-group crossed workflow.

• predict_maihda(type = "strata") now rejects a non-longitudinal fit whose stratum random effects include random slopes, matching the summary guard. It previously returned the intercept column alone as a complete-looking table (estimate, SE, interval) – the stratum effect at zero of the slope variables, an extrapolation that can reorder or sign-flip strata. The brms extractor also no longer presents a slope-only (0 + x | stratum) block’s slope as the stratum effect. Longitudinal trajectory predictions and summaries are unchanged.

• fit_maihda() now rejects a longitudinal growth curve the observed times cannot identify: fewer than time_degree + 1 distinct measurement times, or no person measured at two distinct times. Such models previously fitted and returned arbitrary, optimizer-dependent slope variances flagged only as a singular fit. The check runs on the input and again on the analytic sample after row exclusions.

• maihda_ic() now reports WAIC/LOOIC as NA for a sampling-weighted brms fit, with estimator "Bayesian (weighted pseudo-posterior)", matching the wemix treatment: the weighted pointwise log-likelihoods of a pseudo-posterior are not log predictive densities and define no standard criterion.

• summary(), the VPC, and the PCV helpers now support a binomial model fitted with the complementary log-log (cloglog) link, using the extreme-value latent level-1 variance pi^2/6. Such a model previously fitted but then stopped with a “not implemented” error, even though the response-scale VPC already handled every binomial link.

• pcv_importance() and stepwise_pcv() now emit the same warning maihda() does when a numeric stratum dimension (a category code) enters the models as a single linear term rather than categorical main effects. The fitted models and the reported estimates are unchanged; only the previously missing warning was added.

• pcv_importance(bootstrap = TRUE) now warns and records n_boot_boundary when bootstrap draws are excluded because the null model’s between-stratum variance hit the zero boundary, so an attribution interval conditional on only a handful of surviving draws is no longer returned silently. calculate_pcv() already disclosed this.

• The crossed-dimensions, contextual, and longitudinal VPC bootstraps now count and report non-converged refits (n_boot_nonconverged and a warning), so their n_boot_ok no longer implies a convergence that was never checked. The longitudinal per-time bands also require a majority of finite draws to form, not only ten.

• fit_maihda(family = "binomial") no longer recodes a two-level proportion response (e.g. 0.25 and 0.5) to 0/1. A numeric response is Bernoulli only when its values are whole numbers, so a proportion of successes supplied with weights = trial counts stays an aggregated binomial fit instead of a silently different one.

• Bootstrap intervals (VPC, PCV, crossed and contextual decomposition, longitudinal VPC, pcv_importance()) now require a majority of the eligible refits to succeed, not only an absolute minimum of ten. An extreme failure rate makes the interval unavailable rather than returning one from a handful of draws; legitimately excluded draws (e.g. PCV boundary draws) are not counted against the success fraction.

• The crossed-dimensions MOR now returns 1, not NA, when at least half of the stratum pairs have zero contrast variance. The median odds ratio is exactly 1 there, but the mixture root search could not bracket it.

• Bootstrap draws whose refit optimizer did not converge are counted and reported (n_boot_nonconverged and a warning) instead of silently counted as successful refits, so n_boot_ok no longer implies convergence that was never checked.

• Integer weights= on a Bernoulli fit are no longer read as aggregated binomial trial counts by maihda_discriminatory_accuracy(). Aggregation is now inferred structurally from a cbind() matrix response or a proportion response, so precision weights no longer change the AUC estimand or inflate the reported case/control totals.

• wemix fits are no longer reported as converged unconditionally. WeMix returns the last iterate without warning when its optimisation is abandoned, so convergence is now judged from its own gradient criterion and reported as NA when no evidence is readable.

• Longitudinal count VPC trajectories now evaluate raw-time fixed-effect terms such as x:wave at the reporting time. Under internal time centering only the derived centered column was moved, so those terms stayed at each row’s own observed time.

• maihda_ic() now warns and omits the delta when the models differ in prior or sampling weights, which change the likelihood being maximised.

• stepwise_pcv() now reports Step_PCV as NA when the preceding step’s between-stratum variance is at the singularity boundary, instead of dividing by optimizer noise. The affected steps are listed in an "undefined_step_pcv" attribute and noted by print(); Total_PCV is unchanged.

• The response-scale VPC for a model with non-stratum random effects now estimates the total probability variance on the same basis as its numerator, removing an n_sim / (n_sim - 1) inflation that was largest at small n_sim.

• calculate_pcv(bootstrap = TRUE) now permutes each simulated response into each model’s own row order before refitting, so a bootstrap across two fits of the same observations in a different row order no longer corrupts the interval.

• Individual predictions now reject a supplied stratum that contradicts the intersectional dimension columns in the same newdata, instead of pairing one intersection’s fixed effects with another’s random effect. The check also covers dimension combinations the model never saw, and a row whose dimensions cannot be resolved no longer exempts the rest of the newdata.

• Predictions now check the numeric auto-bin ranges whether or not newdata supplies a stratum column. A row whose auto-binned numeric dimension falls outside the training range previously went unremarked when a stratum was supplied, silently pairing that stratum’s random effect with a combination the training bins cannot contain; it now warns. Existing results are unchanged, and the warning will become an error in a future release.

• lme4 convergence reporting now also consults the optimizer return code and message, so a fit whose optimizer stopped early (e.g. bobyqa hitting maxfun) is no longer reported as converged when lme4’s own gradient check happens to pass.

• brms individual predictions with allow_new_levels = TRUE now zero every unseen grouping level (context and longitudinal, not only stratum) to match lme4, instead of sampling unseen non-stratum effects from the random-effects distribution.

• Zeroing an unseen grouping level no longer overrides a caller-supplied re_formula or re.form. The requested scope is narrowed to drop the unseen terms rather than replaced, so re_formula = NA stays a fixed-effects-only prediction and a partial re_formula is not swapped for a different term.

• Longitudinal PCV now records estimation_used ("fitted", "ML", or "mixed"), so a boundary skip or failed ML refit that leaves a mixed REML/ML comparison is reported rather than mistaken for a clean fitted one.

• calculate_pcv(), compare_maihda(), and maihda_ic() now treat two fits to the same observations supplied in a different row order as the same analytic sample, instead of rejecting the reordered fit or warning about a differing sample. The comparison aligns on the row identifiers and matches the response, stratum partition, and weights after alignment.

• Longitudinal PCV calculations now return NA when the null growth variance is effectively zero. The result records this in null_at_boundary.

• fit_maihda() now rejects duplicate stratum intercepts, including equivalent spellings such as (1 | stratum) + (0 + 1 | stratum) and compound terms that repeat the intercept such as (1 | stratum) + (1 + x | stratum), while still allowing an uncorrelated slope such as (1 | stratum) + (0 + x | stratum).

• Time-dependent formula offsets are now re-evaluated at each time point in longitudinal predictions and count VPC trajectories.

• Failed or skipped ML refits are reported as estimation_used = "mixed"; brms comparisons are reported as "posterior".

• PCV boundary detection now covers lme4, wemix, and ordinal fits and is shared by calculate_pcv(), stepwise_pcv(), and pcv_importance().

• Singular adjusted fits are marked in PCV results and print output without treating a near-100% PCV as automatically erroneous.

• ML-refit failures now warn, retain their actual estimation basis, and no longer produce information-criterion deltas across mixed REML/ML fits.

• External offsets are retained in fitted-row predictions and included in family detection, response recoding, automatic binning, and longitudinal time centering.

• compare_maihda_groups() now computes shared numeric-strata cut points from the pooled analytic sample.

• compare_maihda_groups() and maihda() now exclude rows with missing offsets when selecting the family and calculating analytic sample sizes.

• Failed per-group PCV decompositions now set pcv_status = "failed" and report the affected groups.

• Subsetting a maihda_ic object now preserves the metadata required by its print method.

• maihda_ic() now warns and omits deltas when likelihood and Bayesian information criteria are mixed.

• Ordinal models now recheck the number of observed response categories after analytic-sample filtering.

• Longitudinal validation and time centering now use the transformation-aware analytic model frame.

• PCV result objects now retain and print the requested and actual estimation basis.

• fit_maihda(engine = "brms") now rejects unsupported lme4-style weights, subset, and offset arguments.

• Design-weighted brms workflows now complete derived null and adjusted fits without conflicting with the internal weight column.

• Design-weighted brms fits now keep the full pre-fit data as original_data, so maihda_describe() reports the original total, missingness, and excluded-weight counts rather than the analytic sample alone.

• Non-positive or non-finite Gaussian precision weights are dropped before fitting.

• compare_maihda_groups() and maihda() now apply the same non-positive/non-finite precision-weight exclusion before family/engine detection, shared-strata binning, and analytic sample-size checks.

• Individual predictions now reject missing strata unless allow_new_levels = TRUE.

• Count-family longitudinal VPCs now evaluate residual variance at each time point; the summary includes var_resid_t.

• Crossed-dimensions predictions now exclude contextual random effects from the stratum baseline.

• AUC calculations no longer treat lme4 precision weights as population frequencies.

• The AUC scope is now labelled "intersectional" when fixed effects contribute to the prediction.

• Crossed-dimensions response-scale VPCs now include additive dimension variances in the between-stratum component.

• Contextual binary models now report discriminatory accuracy and optional response-scale VPCs.

• brms convergence is reported as unknown unless an R-hat is available; a divergent-transition count alone is not evidence that the chains converged.

• Optional summaries and plots now warn when a requested component fails instead of silently omitting it.

• WAIC and PSIS-LOO reliability warnings are now passed through.

• Bootstrap requests below 200 replications now warn about unstable percentile intervals.

• Training predictions with an external offset now use the fitted model’s offset-aware prediction path.

• Crossed-dimensions models now reject additional random effects in the formula and direct users to context or the two-model decomposition.

• Sampling-weighted stepwise and importance analyses now exclude invalid weights before family detection.

• Automatic numeric-strata cut points are now based on the full analytic sample.

• PCV comparisons now require matching non-stratum random-effects structures.

• The brms count longitudinal VPC interval now propagates residual uncertainty by posterior draw.

• Contextual binary stepwise_pcv() results now include the AUC and MOR trajectory.

• Scalar PCV helpers now reject crossed-dimensions fits and direct users to the crossed-dimensions decomposition.

• Longitudinal ID/stratum checks now run on the analytic sample.

• The brms variance parser now supports grouping-variable names containing __.

Testing and CI

• Added an integration-tests CI job for optional non-brms backends and expanded regression coverage for PCV, weighting, AUC, offsets, longitudinal models, and crossed-dimensions models.
• Added test-audit-2026-07-13.R for the latest regression fixes.

MAIHDA 0.2.1

CRAN release: 2026-07-09

New features

• Added maihda_describe() for pre-model sample, strata, missingness, weighting, and outcome summaries, with print and plot methods.
• maihda() and compare_maihda_groups() now support stratified contextual analyses using group and context together.
• plot(type = "predicted") now orders strata by predicted value by default and accepts order_by.
• VPC plots for maihda_analysis objects can show the null model, adjusted model, or both.
• Added contextual models to stepwise_pcv().
• Added pcv_importance() for Shapley and dominance-based PCV attribution, with optional bootstrap intervals.

API changes

• Added a predict() method for maihda_model objects.
• Renamed calculate_pvc() to calculate_pcv(). The old name and $pvc result element remain as deprecated aliases.

Bug fixes

• Preserved brms response addition terms such as trials() and weights() when formulas are rebuilt.
• Normalized brms sampling weights after analytic-sample filtering.
• Rejected longitudinal IDs that are reused across strata.
• Corrected AUC handling for aggregated binomial data and fractional precision weights.
• Aligned AUC and MOR to the same intersectional model scope.
• Included non-stratum random effects in response-scale VPC calculations.
• Reported bootstrap boundary draws and rejected effectively zero PCV denominators.
• Fixed zero-row strata creation and aggregated-binomial formulas that use only the strata shorthand.
• Updated brms linear-predictor extraction for current brms releases.
• Fixed single-row ordinal predictions, longitudinal time centering, count residual variance, longitudinal PCV fitting, and higher-order slope PCV indexing.
• Removed whitespace padding from mixed-type stratum labels.

Improvements

• Warned when numeric stratum dimensions enter adjusted models as linear effects.
• Improved guidance for unsupported bootstrap engines and documented latent-scale PCV rescaling.
• Corrected documentation for bootstrap Monte Carlo error, plotting, group comparisons, interaction estimates, and fitted-sample descriptions.

MAIHDA 0.2.0

CRAN release: 2026-07-02

New features

• Added the "upset" plot type and maihda_upset_size() for compact intersection displays.
• Set theme_maihda() as the default theme for package plots.
• Added select = "deviation" to retain the most extreme strata when plots are truncated.
• Added only_flagged to BLUP plots and ensured flagged strata are not hidden by display limits.
• Changed the default interaction adjustment to Benjamini-Hochberg and added ROPE-based decisions.
• Removed the redundant "risk_vs_effect" and "ternary" plot types and the optional ggtern dependency.

Bug fixes

• Warned when model comparisons mix likelihood and Bayesian information-criterion scales.
• Rejected fixed interactions among stratum dimensions in standard MAIHDA decompositions.
• Restored discriminatory-accuracy summaries for brms aggregated-binomial fits.
• Limited ML-refit skipping to boundary stratum variances rather than any singular random effect.
• Weighted aggregated-binomial stratum predictions by trial count.
• Corrected brms ordinal stratum predictions to return expected category scores.
• Rejected scalar interaction diagnostics for longitudinal fits.
• Corrected Bayesian probability-of-direction values for negative effects.
• Standardized unseen-stratum prediction behaviour across engines and added allow_new_levels.
• Normalized brms aggregated-binomial predictions to per-trial probabilities.
• Returned NA rather than NaN for undefined crossed-dimensions shares.
• Included formula offsets in wemix and ordinal predictions.
• Selected families and engines from the analytic sample in high-level workflows.
• Reserved internal .maihda_dim_ column names used for automatic binning.
• Added a note when maihda_table() combines fitted REML variance rows with an ML-based PCV.

Improvements

• Added colour-aware console output through cli, with automatic plain-text fallback.
• Added var_between_adjusted_ml to grouped comparisons.

MAIHDA 0.1.11

CRAN release: 2026-06-18

New features

• maihda() now computes and reports intersectional interaction diagnostics by default.
• Added tidy() and glance() methods for models, summaries, and analyses.
• Added longitudinal growth-curve MAIHDA through the id, time, and time_degree arguments, including VPC and PCV trajectories.
• Added cumulative ordinal models through the ordinal and brms engines.
• Added negative-binomial models for overdispersed count outcomes.
• Added design-weighted models through sampling_weights, with wemix and brms support.
• Added contextual cross-classified MAIHDA through the context argument.
• Renamed the "cross-classified" decomposition to "crossed-dimensions"; the old name remains a deprecated alias.
• Added the high-level maihda() workflow for null/adjusted model fitting and PCV reporting.
• Added the maihda_country_data, maihda_sparse_data, and maihda_long_data example datasets.
• Added maihda_table() for model summaries and ranked-strata exports.
• Added automatic AUC and MOR summaries for binomial models and optional response-scale VPCs.
• Added maihda_ic() for AIC/BIC and WAIC/LOOIC comparisons.
• Added AUC and MOR trajectories to stepwise_pcv() for binary outcomes.
• Added maihda_interactions() and plot highlighting for flagged strata.

Improvements

• Restricted MOR to cumulative-logit and binary-logit models.
• Clarified the conditional interpretation of response-scale VPCs for adjusted models.
• Unified comparison and diagnostic plotting under the base plot() generic; older plotting helpers are deprecated.
• Added lme4 fit diagnostics to model output and grouped-comparison warnings.
• Added group plots for absolute between-stratum variance.
• Clarified PCV, VPC, weighting, supported-family, and plot terminology across the documentation and Shiny app.
• Removed checked-in rendered vignette HTML and corrected documentation dependencies and links.
• Removed package-install side effects from the health-data regeneration script.

Performance

• Replaced row-by-row stratum matching with vectorized, collision-safe key matching.

Bug fixes

• Applied analytic-sample identity checks to design-weighted fits.
• Evaluated longitudinal baseline PCV at the observed reference time.
• Rejected cross-sectional stratum summaries for longitudinal models.
• Reported analytic sample sizes for wemix information-criterion tables.
• Corrected Gaussian PCV comparisons by using ML refits for models with different fixed effects.
• Rejected VPC/ICC calculations for Gaussian models with non-identity links.
• Based binary-family detection and numeric-strata binning on the analytic sample.
• Fixed nested forwarding of data-masked weights, subset, and offset arguments.
• Preserved original response labels while evaluating subset expressions.
• Correctly sliced external weights, subsets, and offsets in grouped fits.
• Accounted for precision weights in Gaussian residual variance and stratum-level plot aggregation.
• Corrected effect-decomposition random effects and ordinal surprise scores.
• Aligned Shiny family detection and bootstrap VPC/ICC output with the core API.
• Added grouped-comparison checks for analytic sample size, populated strata, and shared stratum support.
• Corrected response-scale plotting for count models.
• Made predict_maihda(type = "strata") respect newdata.
• Required at least 10 bootstrap replications.
• Rejected PCV comparisons with different prior weights.
• Reported the event mapping when two-level outcomes are recoded to 0/1.
• Kept Shiny analyses available when the null between-stratum variance is zero.
• Accepted bare family functions in grouped comparisons.

Diagnostics

• Added maximum R-hat and divergent-transition counts for brms fits.
• Added successful-draw counts and Monte Carlo standard errors to bootstrap output.

MAIHDA 0.1.10

New features

• Added compare_maihda_groups() for VPC/ICC and variance comparisons across groups.
• Added plot_group_comparison() for forest and variance-composition plots.

Bug fixes

• Excluded failed parametric-bootstrap refits from VPC and PCV intervals.
• Corrected the Poisson residual-variance approximation to log(1 + 1 / mu).
• Omitted unsupported confidence intervals from prediction-deviation plots instead of drawing zero-width intervals.

MAIHDA 0.1.9

Bug fixes

• Clarified the interpretation of negative PCV values in the Shiny dashboard.
• Fixed the coverage workflow’s failure-artifact upload configuration.

MAIHDA 0.1.8

CRAN release: 2026-05-16

New features

• Added prediction-deviation, predicted-value, risk-versus-effect, and effect-decomposition plots.
• Added automatic tertile binning for numeric grouping variables with more than 10 unique values.
• Expanded the Shiny dashboard with the new plots and automatic-binning controls.
• Added automatic detection of binomial outcomes.

Bug fixes

• Corrected latent residual variance for binomial and ordinal VPC/ICC calculations.

MAIHDA 0.1.7

CRAN release: 2026-04-05

New features

• Added stepwise_pcv().
• Added the interactive Shiny dashboard.
• Improved bootstrap confidence-interval performance.
• Updated package tests, imports, dataset documentation, and the introductory vignette.

MAIHDA 0.1.0

CRAN release: 2026-04-03

Initial release

• Added intersectional-strata creation and multilevel model fitting with lme4 and brms.
• Added variance summaries, individual and stratum predictions, model comparison, and core plots.
• Added documentation, vignettes, and tests.
• Improved missing-value handling in make_strata().