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().