Added maihda_describe() for pre-model sample, strata, missingness, weighting, and outcome summaries, with print and plot methods.
Added pcv_importance() for Shapley and dominance-based PCV attribution, with optional bootstrap intervals.
Added a predict() method for maihda_model objects.
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 assumes that strata do not diverge over time. 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.
maihda_vpc_response() now also returns the exact count-scale VPC of Leckie et al. (2020) for Poisson and negative-binomial lme4 fits whose only random effect is the stratum intercept; summary(response_vpc = TRUE) attaches it.
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.
calculate_pcv() and maihda() gained estimation = "REML", which compares REML-type variances for binomial and Poisson lme4 fits; their maximum-likelihood PCV runs high with few or small strata.
A brms crossed-dimensions fit now reports its realised strata (Gelman 2005) in decomposition$realised and a second print() block: the VPC, shares, variances and median odds ratio of the effects of the strata and levels in the data, which depend little on the prior when each stratum has enough data. The few-levels note points to the block.
A brms fit with a context now reports its realised partition (Gelman 2005) in summary()$context$realised and a second print() block: the variances, VPC and context shares of the effects of the strata and context levels in the data, which depend little on the prior when each context and stratum has enough data (variances only for Poisson and negative-binomial fits). A two-model analysis prints the realised VPC and context share, a crossed-dimensions fit's realised-strata block its VPC, and the few-levels note points to them.
print() on a PCV built from maximum-likelihood variances (a glmer, wemix or ordinal fit) now says so; the numbers are unchanged.
print() on a two-model maihda() analysis and on a longitudinal PCV now names the variances the PCV compares, as calculate_pcv() does, where no note already names them; a longitudinal PCV records ml_variances. The numbers are unchanged.
print() on stepwise_pcv(), pcv_importance() and compare_maihda_groups() results built from maximum-likelihood variances (a glmer, wemix or ordinal fit, or an lmer fit made with REML = FALSE) now names them as calculate_pcv() does and notes that they run low with few or small strata, so the PCV runs high; the results record ml_variances. The numbers are unchanged.
print() on a maihda(decomposition = "crossed-dimensions") analysis fitted with engine = "brms", and on its summary(), now names the stratum dimensions with fewer than 10 levels: the variance of such a dimension depends on the prior on its standard deviation, and the VPC and the shares follow it. The level counts are in decomposition$dim_levels; the numbers are unchanged.
print() on a maihda() analysis fitted with engine = "brms", and on the summary() of a brms fit, now names a context with fewer than 10 levels: its variance depends on the prior on its standard deviation, and the context share and the VPC depend on it. The level counts are in the contextual partition's context_levels; the numbers are unchanged.
predict_maihda(allow_new_levels = TRUE) on a brms fit now refuses newdata missing a grouping column the prediction uses, such as the context column, as lme4 does; brms filled the column with NA and drew a new effect for every row.
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 before 2.1-0 leaves out of logLik() for a glmer fit with nAGQ > 1 (a fit made by lme4 2.1-0 or later already holds it), 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.
The prediction-deviation panels draw 95% bars for lme4 fits from the joint covariance of the fixed and random effects: a stratum's bar is no longer the average of its rows' SEs, and bars of fits with several grouping terms apply the sparse factor's permutation that lme4's own predict(se.fit = TRUE) leaves out (off several-fold there).
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.
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.
print() on a longitudinal PCV with a rank-deficient linear stratum block now names the time at which that block has no between-stratum variance, near which the PCV is not interpretable; the numbers are unchanged.
Gaussian growth models fitted by lme4 with a linear stratum slope are refitted from two more starting values, and a refit replaces lmer's fit when its REML criterion (deviance for a maximum-likelihood fit) is lower by more than 1e-4. lmer's default optimiser could stop on the wrong side of the stratum block's boundary; pass start to fit from your own values only.
Gaussian growth models fitted by lme4 with quadratic or higher growth, with or without a stratum slope, are refitted from lmer's fit of the same model on time centred at the middle of its range, and, when that refit and lmer's fit differ by more than 0.01, from a start with both covariance blocks diagonal; a refit replaces lmer's fit when its criterion is lower by more than 1e-4. lmer's default optimiser stopped above the REML optimum in 10 to 20 per cent of such fits, recorded as converged.
summary(df_method = "bootstrap") on a growth model fitted by lme4 refits each null-restricted model as fit_maihda() fits the model, with the same refits from further starting values, and keeps the lower of that fit and the former refit. Refitted without them, the null model of a quadratic growth fit stopped above its optimum in 14 of 110 refits on a balanced design, by up to 1.95 REML units.
summary(bootstrap = TRUE) and summary(df_method = "bootstrap") on a growth model fitted by lme4 refit each draw from the fitted model's estimates and then as fit_maihda() fits the model, from its further starting values and, for glmer, across a variance bound, at two to five times the time. Refitted once, 12 to 19 of 100 draws of a Poisson growth model stopped short of their maximum, the VPC(t) band's upper limit ran up to 25 per cent low at early times, and draws stopped beside a variance bound gave bootstrap statistics up to 53 from standard errors about 60 times too small.
Gaussian growth models fitted by lme4 whose time span is not between 1 and 10 are fitted on time divided by a power of ten (days over five years by 1000), with a message; results are reported on the original time scale, while the stored lme4 fit is on the internal time. Counted in fine units, lmer stopped far from the optimum, mostly recorded as converged: in days maihda_long_data's ML PCV was 0.92 / 0.35 (0.96 / 0.92 in waves) and maihda_ic() was off by 78 and 87 log-likelihood units, and quadratic growth failed from a span of about 12.
Growth models fitted by glmer (Poisson, binomial, negative binomial) are fitted on the internally scaled time too, so their fixed-effect standard errors no longer depend on the unit of time.
glmer growth models with a linear stratum slope are refitted from two more starts, and any glmer growth fit with a block's intercept variance at zero is refitted from the other side of that bound by bobyqa and, when that refit is kept, once more (across the bound again if it stops at one), at about three times the fitting time; glmer stopped more than 0.01 below the maximum in about a third of 120 Poisson growth fits, by up to 32 log-likelihood units.
A glmer growth fit stopped with a block's intercept variance at zero where the likelihood still rises is recorded as not converged, the blocks in diagnostics$boundary_traps.
glmer fits with a random slope outside growth models, such as (1 + x | stratum), whose fit stops with a block's intercept variance at zero are refitted once from the other side of that bound, by glmer's second stage with bobyqa; on simulated cross-sectional designs 39 of 200 fits had stopped there more than 0.01 below the maximum, by up to 23 log-likelihood units.
The convergence record flags such a stop in any glmer fit with nAGQ of 1 or more, not only in a growth fit, when the likelihood rises with the signs of the intercept's covariances reversed; the sentence names the block's terms.
glmer fits with a random slope outside growth models are also searched with each random-effect block's terms re-expressed as for lmer fits and the fixed effects centred and scaled, from lme4's starting values of the variance parameters and from the fit's estimates, and refitted from the lower end. With a slope on a calendar or birth year, glmer's own second stage stopped more than 0.01 short of the maximum in 37 of 133 simulated fits, by up to 37 in -2 log-likelihood, often with the slope's variance at zero, where lme4 calls the fit singular; 27 were recorded as converged. Fitting takes about one and a half times as long, up to about four times.
glmer growth models, and their bootstrap and null-model refits, are searched the same way and refitted from the lower end when it is at least 0.001 lower in -2 log-likelihood: about a fifth more fitting time, about twice as long where a refit is made. Quadratic growth fits had stopped more than 0.01 below the maximum in 47 of 240 simulated fits, by up to 21, 11 of them singular and recorded as converged.
A glmer fit outside growth models that either refit replaced is refitted once more from its estimates by the second stage of the call's own control (Nelder-Mead by default) and so records that optimiser, as glmer's own fits do: recorded as bobyqa, lme4::refit() refitted such a fit by bobyqa with minqa's default first step and ended more than 0.01 above it on 18 of 43 simulated fits (it failed on 15), by up to 188 in -2 log-likelihood.
glmer growth fits that a refit replaced, and their bootstrap and null-model refits, are refitted once more the same way and so record Nelder-Mead: recorded as bobyqa, lme4::refit() stopped at minqa's limit of evaluations on 19 of 87 responses simulated from 29 such fits, most within 0.001 of the maximum, where a bootstrap draw's refit keeps such a stop and counts it as not converged (2 of 20 draws of one fit, with bobyqa's warning), and it ended up to 1.8 above 2 of the 29 fits. On a simulated response, lme4::refit() of such a fit, which lme4::bootMer() now runs, stops short of the maximum more often (80 of the 87, against 13); the bootstraps of summary() reach it either way.
The last pass of a glmer.nb() fit, which searches the dispersion again, is kept only when its -2 log-likelihood at its own dispersion is not higher: compared at the new dispersion alone, a growth fit's pass that ended 1.36 higher was kept.
The searches of glmer fits with a random slope outside growth models run their deviance function to a relative tolerance of 1e-10, and every glmer search, growth models' included, takes a fit when that function at its estimates lies within 1e-3 (relative) of lme4's -2 log-likelihood: the 1e-5 check at lme4's tolerance had left some simulated binomial null models with strata far apart on the logit scale unsearched, up to 13.6 below the maximum and some recorded as converged. The singular-fit record's search does the same.
A glmer fit outside growth models stopped where the variance of a term followed by two or more others in its random-effect block is at zero, such as the intercept of (1 + x + z | stratum), is refitted from across that bound along the direction in which the likelihood rises fastest instead of with the signs of the term's covariances reversed; on simulated designs the reflection had left 10 of 57 such fits more than 0.01 below the maximum, by up to 1.78 log-likelihood units, 5 of them recorded as converged.
The convergence record measures each such stop along the direction in which the likelihood rises fastest as well as by the reflection, which alone missed stops of fits up to 1.03 log-likelihood units short; for a term followed by two or more others its sentence says the likelihood rises when that variance and its covariances move off zero.
The convergence record measures a stop at a variance bound at steps of 0.01 and 0.1 as well as 0.001, and by a Newton step in the block's variance parameters with its terms centred and scaled; the step of 0.001 alone called 30 of 65 simulated lme4 fits stopped there more than 0.01 below the maximum converged, most with a slope covariate in large units such as age in years or a covariate in thousands. Its sentence says what showed the likelihood rising and names a bound only for a column at one: a slope's entry on a covariate in thousands lies within lme4's 1e-4 of the bound at its natural size. boundary_traps gains the falls across the bound and reflected, and at_bound.
lmer fits with a random slope outside growth models are refitted from the other side of a bound their fit stops at and from a start carried back from their random-effect terms centred and scaled, keeping the lowest REML criterion, at about three times the fitting time; lmer stopped more than 0.01 above the optimum in 64 of 480 simulated fits with a covariate on a scale of about 1, 121 of 180 with age in years and 55 of 60 with a birth year.
The convergence record flags a stop at such a bound in an lmer fit outside growth models too.
lmer fits with a random slope outside growth models take that start with each random-effect block's terms re-expressed on their own, a term in a block without an intercept, as in (0 + x | ctx) or (1 + x || stratum), divided by its root mean square; the fit kept and lmer's own fit are also refitted from a variance at zero where the criterion still falls off zero. Of 570 simulated fits, 57 had stopped more than 0.01 above the REML optimum, by up to 12.8, all recorded as converged, and all but one had a calendar or birth-year slope in such a block; fitting takes about four times as long as lmer's own fit, up to about nine times where lmer's own fit stops with a variance at zero.
The convergence record flags an lmer fit outside growth models that lme4 calls singular, which lme4 does not check, when a search of the same model reaches a REML criterion 0.01 or more below it, in diagnostics$search_gain.
It searches an lmer fit with a random slope outside growth models that lme4 does not check because of its size, by default one of 10,000 rows or more or with 10 or more variance parameters, in the same way, unless fit_maihda() refitted it; in simulations such fits made from a start had been recorded as converged by up to 1,236 above the REML optimum.
It flags a glmer fit with a random slope and nAGQ of 1 or more outside growth models the same way, searched from the fit with its random-effect terms and fixed effects centred and scaled, when the deviance falls 0.01 or more; of 119 simulated such fits lme4 called singular, with a slope on a birth year, an age, an N(0, 1) covariate or an income in currency units, 37 were 0.01 to 23 short of the maximum, and 36 of them were recorded converged.
A glmer fit with a random slope and nAGQ of 1 or more outside growth models that lme4 does not check because of its size, by default one of 10,000 rows or more or with 20 or more parameters, fixed effects included, and does not call singular, is measured as one whose checks fired, and counts as not converged 0.1 standard errors or more from a stationary point, with diagnostics$lme4_unchecked; where that cannot be measured it is searched as a singular one is. In simulations such fits made from a start on 10,000 to 20,016 rows had been recorded as converged 0.017 to 14.3 log-likelihood units below the maximum.
The rank-deficiency verdict for a linear stratum block, behind PCV_slope's NA and the longitudinal PCV's note, no longer depends on the unit of time.
The longitudinal PCV of engine = "brms" fits -- PCV_intercept, PCV_slope and the time-specific PCV of maihda(decomposition = "longitudinal") -- now compares posterior medians of the per-draw between-stratum variances instead of variances read off the posterior-mean covariance blocks. With few strata the mean-based PCV ran low (0.73 for a true 0.83 at 12 strata, 0.80 from medians; 0.81 and 0.83 at 24), and from medians it still runs low on average when strata are few. The stored blocks Sigma_stratum_null / Sigma_stratum_adjusted remain posterior means; the new variance_summary and Sigma_summary fields say which is which, and print() notes it.
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.
The two-model PCV of engine = "brms" fits -- calculate_pcv(), maihda(), stepwise_pcv(), pcv_importance(), compare_maihda_groups() -- now compares the posterior medians of the between-stratum variances instead of their posterior means. With few strata the mean-based PCV ran low (0.58 for a true 0.84 at 8 strata, 0.74 from medians) and changed from one sampler seed to the next.
For engine = "brms" fits, maihda_mor() and the variance-components table of summary() -- with everything read from it: maihda_table()'s variance and SD rows, compare_maihda_groups()'s variance columns, stepwise_pcv()'s Context_Variance, the VPC plots -- now report the posterior median of each variance instead of its posterior mean, as the PCV does, with totals and proportions computed from those rows. Random-effect variances and MORs fall, most with few strata. A longitudinal fit's table is 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.
fit_maihda(engine = "ordinal") now compares a Laplace fit that ends at or near a between-stratum variance of zero with the boundary fit, made by ordinal::clmm()'s general routine (clmm.control(useMatrix = TRUE)), and returns the boundary fit unless its log-likelihood is lower by more than 1e-9 of its size: strata with the same outcome proportions, or a fixed part that spans the strata, no longer return clmm()'s starting variance of 1, and the boundary fit returned has a variance of zero, is flagged singular and has threshold standard errors.
fit_maihda(engine = "ordinal") now fits a model with fixed effects on its design with the columns centred and scaled and returns the clmm fit of the formula at that maximum: a covariate far from zero or on a large or tiny scale (a birth year, currency units, a squared age) could leave ordinal::clmm() at its starting values or short of the maximum. A fit given a start is made as the formula stands.
Standard errors of engine = "ordinal" fits now come from the Hessian of that rescaled fit, carried over exactly; clmm()'s own, taken by differences with a fixed step, were wrong for a covariate far from zero or on a large scale (for a birth year between 0.01 and 1.8 times the right value on simulated data).
attr(maihda_sparse_data, "truth") now holds what the 36 drawn stratum effects carry, and the effects themselves: VPC 0.24 on y and 0.29 on event, interaction share 41%, where it held the generator's targets (0.26, 0.31, 40%). ?maihda_sparse_data and the sparse-intersections vignette quote these, and the vignette adds lme4's share with 1,000 individuals in every stratum (28% and 35%).
calculate_pvc(), deprecated since 0.2.1, is removed; use calculate_pcv().
plot_comparison(), deprecated since 0.1.11, is removed; use plot() on the compare_maihda() result.
plot_group_comparison(), deprecated since 0.1.11, is removed; use plot() on the compare_maihda_groups() result.
decomposition = "cross-classified", deprecated since 0.1.11, is removed: maihda() and compare_maihda_groups() stop with an error that points to "crossed-dimensions" and to context =.
Reading the deprecated pvc element of a calculate_pcv() result with $ or [[ now warns; read pcv.
summary(), maihda_proportional_odds_test(), fit_maihda(), predict_maihda() and calculate_pcv(): simulation figures and notes on earlier behaviour are removed, the guidance is unchanged.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.?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.glmer()'s convergence and "nearly unidentifiable" warnings on the strata's own main effects; they come from a covariate far from zero or in large units, such as age in years, and from lme4's tolerance.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.?calculate_pcv, ?maihda and the PCV print output no longer describe estimation = "ML" as a "correction-free" comparison: maximum likelihood understates the adjusted model's between-stratum variance by more than the null model's, so it overstates the PCV when strata are few. Estimates are unchanged.?maihda and the cross-classified and planning vignettes no longer recommend engine = "brms" for a crossed-dimensions decomposition whose dimensions have few levels. Under the default prior the VPC and the additive share of such a fit are overstated, most with two- to four-level dimensions; the two-model decomposition, or a prior with a sensitivity check, is the alternative.?maihda, ?fit_maihda, the README and the cross-classified vignette no longer recommend engine = "brms" for a context with few levels: under the default prior its variance comes out high for context effects smaller than that prior's scale, most with two or three levels, and the context share and the VPC depend on it, while lme4 may estimate it at zero. The sparse-intersections vignette now says that its shares depend on the priors it uses.?maihda_sparse_data no longer calls the brms interval calibrated, and the planning and longitudinal vignettes no longer recommend engine = "brms" for sparse cells or few waves: with cells this small its estimate depends more on the prior than on the data, and its interval shows how little they say.?fit_maihda no longer present a rank-deficient stratum block as evidence that the strata do not diverge: with few strata it is common when they do. They say what stratum_slope = FALSE does to diverging strata, and the vignette that the AIC comparison has little power with few strata.lme4 fits put the interaction variance at the boundary: they put a dimension's variance at zero (two in the binary fit), and the interaction variance is positive in both.?pcv_importance no longer calls the exact methods deterministic: the bootstrap = TRUE intervals and engine = "brms" fits are random draws and need set.seed().?maihda_sparse_data no longer drops the text after each % or shows Markdown markup, and no longer says the vignette recovers the true interaction share.?plot_prediction_deviation_panels, ?stepwise_pcv and ?theme_maihda no longer show Markdown markup.lme4 cannot give the shares an interval (bootstrap = TRUE gives one) or that the brms prior keeps the standard deviations off zero, and it prints the fits' variances instead of their singular flags.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).lme4::bootMer() and confint(method = "boot") run on lme4 fits made with control = lmerControl() / glmerControl(), and on glmer fits that a refit replaced and Nelder-Mead finished with its own first steps: the call gives the control as lme4::lmerControl() / lme4::glmerControl(), not as a name bound in the formula's environment. They still stop on a glmer fit last refitted by bobyqa or with first steps of the package's own.
fit_maihda() now gives an lme4 fit the refits and convergence checks it gives the same model with its random-effect blocks written in another order. Where lme4 sorts the blocks by their number of levels and two have as many levels, as (x || stratum) before a factor with more does, the standardised refit start, the lmer convergence re-check, the searches of singular lmer and glmer fits and the trap record's Newton step were skipped silently, so that 12 of 12 simulated fits with a slope on a birth year stopped 7 to 44 log-likelihood units below the optimum, 9 of them recorded converged; the glmer re-check was taken in lme4's own units, where a slope in large units gives none, and the singular verdict by lme4's own test.
set.seed() now reproduces engine = "brms" fits whose chains run in parallel (cores > 1): an unseeded fit takes its sampler seed from R's random-number stream, as a serial fit does, instead of each chain's process drawing its own.
print() on a longitudinal PCV now shows a variance below 1e-4, such as a slope variance with time counted in days, to four significant digits instead of as 0.0000.
The rank-deficiency note of print() on a longitudinal PCV no longer offers stratum_slope = FALSE as the remedy: it says that such a block is common with few strata even when the strata diverge, and that the constraint then moves their slope differences to the person level.
fit_maihda(engine = "wemix") on a binomial model now says that the between-stratum variance was estimated at zero when WeMix::mix() stops there, instead of passing on its Matrix seems negative semi-definite error.
fit_maihda(engine = "wemix") on a binomial model now says what happened when WeMix::mix() stops because more than half of the stratum effects are exactly zero at a likelihood evaluation, instead of passing on its non-conformable matrix dimensions error.
A fit_maihda(engine = "ordinal") fit made with ordinal's ucminf optimiser is no longer always recorded as not converged: ucminf()'s return is read by its own codes, and where it says only that the search could go no further, the fit counts as converged when the gradient and Hessian clmm() stores with it put it within a tenth of a standard error of a stationary point.
A glmer fit whose optimiser reports success now counts as converged when lme4's convergence checks fire but the fit lies within a tenth of a standard error of a stationary point, measured with the fixed effects centred and scaled: lme4 takes its checks in the covariates' own units, where a calendar year, an income or an age in years on a few thousand rows can trip them at the maximum, so such fits were recorded as not converged and compare_maihda_groups() named their groups. print() says what lme4 warned of and how far the fit is.
A glmer fit with a random slope, as in (1 + x | stratum), is now judged with its random-effect terms centred and scaled as well as its fixed effects: lme4 takes its checks in the units of the slope's covariate, so fits with a slope on a birth year were recorded as not converged at the maximum (6 of 12 simulated fits).
That measure and the record of a stop at a variance bound run their deviance function to a relative tolerance of 1e-10 and take a glmer fit whose -2 log-likelihood lies within 1e-3 (relative) of the function: with lme4 before 2.1-0, their 1e-5 check at lme4's tolerance left 17 of 135 simulated 0 / 1 fits with strata far apart on the logit scale and a birth year unmeasured and recorded as not converged, all within 0.07 standard errors of the maximum, and 7 of 26 stops at a bound beside a block with a large variance unmeasured.
A glmer or lmer fit with a diag() random-effect block is judged with that block's terms divided by their root mean square, as the same model written with || is: glmer fits with a diag() slope on a birth year were recorded as not converged at the maximum (4 of 8 simulated Poisson fits), and an lmer fit with such a block kept lme4's verdict.
An lmer fit is judged the same way, with its random-effect terms centred and scaled: lme4 takes its checks in the units of a random slope's covariate, so fits with a slope on a birth year, an age in years or an income in thousands were recorded as not converged at the REML optimum (17 of 60, 47 of 180 and 24 of 60 simulated fits).
In an lmer fit, or a glmer fit with nAGQ of 1 or more, outside growth models, a random-effect block with a slope is judged singular when a diagonal entry of its factor with the terms centred and scaled is below 0.01: lme4's test in the formula's units read a slope on a birth year or an income in currency units as singular in blocks that are not at the boundary (14 of the 88 simulated birth-year lmer fits lme4 calls singular at the REML optimum). lme4's own verdict is kept in diagnostics$lme4_singular where the two differ, and print() says when lme4 reports a singular fit that the package does not.
compare_maihda_groups() now gives a failure reason shared by several groups once, in one warning for failed fits and one for failed PCV decompositions, names a one-level stratum dimension as the cause only for the groups where it is, and keeps the first line of an error message in status.
The fit-quality warning of compare_maihda_groups() no longer says that lme4 reported the problem on the wemix, ordinal and brms engines or that a singular fit always has a between-stratum variance of zero, explains only what it lists, and lists a brms group under a largest R-hat above 1.01 or divergent transitions.
The fit-quality warning of compare_maihda_groups() now names a group for a singular fit only when its between-stratum variance is estimated at zero; a zero interaction, dimension, context or other variance component no longer names it.
The fit-quality warning of compare_maihda_groups() now also names a group whose adjusted model did not converge (on brms: has a largest R-hat above 1.01 or divergent transitions) and points to its PCV.
The few-context-levels warning of compare_maihda_groups() and the context_few_levels note of maihda_describe() no longer advise engine = "brms"; they say what each engine does with so few levels: an lme4 variance is imprecise and at zero when the levels differ little, a brms one depends on the prior on its standard deviation.
The empty_strata note of maihda_describe() no longer says sparse designs are where engine = "brms" is most useful; it says that small strata are shrunk toward the mean and their effects are imprecise with either engine.
The singular-fit note printed for a decomposition = "crossed-dimensions" fit, and shown in the Shiny app, now reads a zero stratum term as the interaction and a zero dimension term as part of the additive variance, instead of calling the between-stratum variance unreliable or unaffected.
print() of a maihda_interactions() table now says that the stratum standard errors collapse only when the term the interaction is read from has a variance of zero, not for a zero in a dimension, context or other term, and calls it the interaction variance for a crossed-dimensions fit; the new attribute interaction_at_boundary records it.
The "Fit diagnostics" card of the Shiny app now shows one bullet per diagnostic, its sentences as running text and the convergence messages listed under it, instead of one bullet per printed line.
predict_maihda(allow_new_levels = TRUE) on a brms fit now predicts a row with a missing grouping value -- a stratum dimension, a supplied stratum, a context or a longitudinal id -- at a zero effect for that term, as the other engines do, instead of drawing a new effect for it on every call.
summary() on a binomial engine = "wemix" fit with a small between-stratum variance no longer passes on WeMix's diag(V) had non-positive or NA entries warning, and neither does anything that calls it, such as tidy(), maihda_table(), plot(), maihda() and compare_maihda_groups(); estimates are unchanged.
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.
predict_maihda(newdata = ) on an lme4 fit now gives each row its own stratum's random effect. lme4 2.0.1 matched a character copy of an integer grouping column -- the stratum predict_maihda() rebuilds from the dimension columns -- to the wrong strata whenever its text order differs from the numeric one ("10" before "2"): 0.85 off on a 20-stratum fit at the default settings. plot_prediction_deviation_panels(data = ) did the same for a character stratum.
predict_maihda(allow_new_levels = TRUE) on an lme4 fit now predicts a grouping column missing in every newdata row -- one row missing its stratum dimension, a supplied NA stratum, a context or longitudinal id that is a logical or numeric NA throughout -- at a zero effect for that term instead of failing with "Invalid grouping factor specification".
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 __.
integration-tests CI job for optional non-brms backends and expanded regression coverage for PCV, weighting, AUC, offsets, longitudinal models, and crossed-dimensions models.test-audit-2026-07-13.R for the latest regression fixes.refit() on incomplete cases, the Laplace deviance at the converged step, bootMer()'s control lookup, where glmer stops a quadratic growth fit) check either version.glmer.nb() fit) hold on every CI platform.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.maihda_analysis objects can show the null model, adjusted model, or both.stepwise_pcv().calculate_pvc() to calculate_pcv(). The old name and $pvc result element remain as deprecated aliases.trials() and weights() when formulas are rebuilt."upset" plot type and maihda_upset_size() for compact intersection displays.theme_maihda() as the default theme for package plots.select = "deviation" to retain the most extreme strata when plots are truncated.only_flagged to BLUP plots and ensured flagged strata are not hidden by display limits."risk_vs_effect" and "ternary" plot types and the optional ggtern dependency.allow_new_levels.NA rather than NaN for undefined crossed-dimensions shares.wemix and ordinal predictions..maihda_dim_ column names used for automatic binning.maihda_table() combines fitted REML variance rows with an ML-based PCV.cli, with automatic plain-text fallback.var_between_adjusted_ml to grouped comparisons.maihda() now computes and reports intersectional interaction diagnostics by default.tidy() and glance() methods for models, summaries, and analyses.id, time, and time_degree arguments, including VPC and PCV trajectories.ordinal and brms engines.sampling_weights, with wemix and brms support.context argument."cross-classified" decomposition to "crossed-dimensions"; the old name remains a deprecated alias.maihda() workflow for null/adjusted model fitting and PCV reporting.maihda_country_data, maihda_sparse_data, and maihda_long_data example datasets.maihda_table() for model summaries and ranked-strata exports.maihda_ic() for AIC/BIC and WAIC/LOOIC comparisons.stepwise_pcv() for binary outcomes.maihda_interactions() and plot highlighting for flagged strata.plot() generic; older plotting helpers are deprecated.wemix information-criterion tables.weights, subset, and offset arguments.predict_maihda(type = "strata") respect newdata.compare_maihda_groups() for VPC/ICC and variance comparisons across groups.plot_group_comparison() for forest and variance-composition plots.log(1 + 1 / mu).stepwise_pcv().make_strata().