A patch release for the bootstrap intervals of the reliability family, prompted by a report from an analysis of the Big Five Inventory constructs (2026-10-14), in which every BCa interval for categorical omega failed with "estimated adjustment 'a' is NA" and no interval said how many replications it rested on.
Every bootstrap interval in the reliability family now says how
many replications it conditions on. reliability_alpha(),
reliability_kr20(), reliability_omega(),
reliability_omega_categorical(), and reliability() through them
record, in a bootstrap attribute, how many replications were
attempted, how many returned a finite value, how many were set aside
and why (an error from the refit, or no finite value, as from a fit
that did not converge), and the distinct error messages with their
counts. The standard error is the standard deviation of the returned
replicates and the percentile and BCa limits are their quantiles,
which the help pages now state; the printed table shows the counts
in one line, and a single warning states them whenever a replication
was set aside. Through 1.0.0 a failed replication was dropped
silently. No existing row or attribute changes.
The BCa interval no longer fails when a jackknife refit fails.
With B below the number of cases, boot estimates the BCa
acceleration from jackknife influence values, and one leave-one-out
refit that returned no value made the whole vector NA and stopped
the call. The influence values are now computed in the package under
boot's own rule (regression on the resampling frequencies when at
least as many replications as cases returned a value, the jackknife
otherwise), the leave-one-out refits that fail are counted and set
aside, and the acceleration is estimated from the rest; the
jackknife counts and the acceleration join the bootstrap
attribute. When the acceleration or the bias correction cannot be
estimated at all, the percentile interval is returned in its place
with a warning, and the ci_method attribute reads
"percentile (BCa acceleration could not be estimated: k jackknife replicates failed)", so the substitution is on the record. On a run
where nothing fails the interval and the acceleration equal boot's
to machine precision, in both regimes, so no interval that 1.0.0
computed changes.
reliability_omega() states which intervals apply to which
denominator. The closed forms, the Wald intervals, and the profile
likelihood are derived under the model implied ratio and apply to
denominator = "model_implied" only; the bootstrap methods and
"none" apply to both. The help page now says so in one place, and
a request for a combination that does not apply, such as
ci_method = "mlr" with the default observed denominator, stops
with a message naming the denominator that does.
First public release of DMAR (pronounced "Dee-Mar," for "Design, Measurement, and Analysis in R"), a greatly expanded reimagining of the MBESS package, published on CRAN on 2026-09-21.
Every line of code in every help page example now runs. Earlier
drafts carried a slow call (a bootstrap interval, a Monte Carlo
planner, a refit) as a comment so that a reader could see the
syntax without the check paying for it; the CRAN review of 2026-09-04
asked that no example line be commented out, and none is. Where a
call was slow, it runs with a small replication count and a comment
naming the count a reported analysis deserves (the defaults are
unchanged); where a page demonstrated the same thing twice, it now
demonstrates it once. The package still ships no \donttest{}
and no \dontrun{}, and a parse-based detector in the test suite
fails on any comment line that would parse as code.
The sensitivity family no longer writes a file unless asked to, and
never to a default location. The save argument is gone from
<<SAVE_LIST>>; the single switch is filename, which defaults to
NULL. Supplying a path writes the per-replication results there
exactly as save = TRUE used to; a throwaway run points it at
tempfile(fileext = ".csv").
A supplied seed is now set through withr::local_seed() (withr,
which has no dependencies of its own, joins Imports), by way of one
internal helper that every seeded function calls. The behavior is
what it was, the user's generator state is restored when the
function returns, but the package no longer touches
.Random.seed or the global environment to do it, and no function
uses <<-: the handlers that count clamped or non-converged
replications keep their counts in a local environment.
The six Monte Carlo functions that silenced every warning for their
duration with options(warn = -1) now muffle only the warnings
their loops are known to raise (lavaan convergence and variance
warnings on borderline replicates, the noncentral clamp); any other
warning reaches the user.
The Bryant-Paulson distribution function integrates over the
covariate shrinkage factor after the change of variables
delta = 1 - u^2, under which the Beta weight of Bryant and Paulson
(1976, Equation 12) is smooth on the unit interval for every number
of covariates instead of carrying an integrable singularity at the
endpoint where the distribution puts most of its mass. The two
forms agree to about 1e-10 (checked on a grid of 288 quantile,
covariate, group, and degrees of freedom combinations), the Table 1
reproduction in the tests is unchanged, and a critical value now
costs about a twentieth of what it did, which is what lets every
qbryant_paulson(), cv_bryant_paulson(), and ci_c_ancova_bp()
example run live.
Two smaller changes from the same pass. mlmr(), mlmr_mv(), and
average_variance_extracted() resample through
lavaan::lavBootstrap(), the entry point that replaced
bootstrapLavaan() in lavaan 0.7-2 (the package already requires
that version); the results are identical. The five sensitivity
functions whose ... was never forwarded anywhere
(ss_aipe_sc_sensitivity(), ss_aipe_smd_sensitivity(),
ss_aipe_sm_sensitivity(), ss_aipe_R2_sensitivity(), and
ss_power_R2_sensitivity()) no longer accept it, so a misspelled
or retired argument such as save = TRUE now stops with an error
instead of being swallowed. ss_aipe_cv() with assurance
supplied reports its noncentral accuracy note once, with the
combined count, where it used to report it twice; and
ss_aipe_cv_sensitivity() and ss_aipe_sc_ancova_sensitivity()
report that note once per run, with the number of replications it
arose in, rather than once per replication.
A package-wide audit (numerical verification against complex-step
derivatives and Monte Carlo, documentation contracts run rather
than read, and mechanical convention gates) closed out with these
fixes. analysis_of_change() now fits an intercept-only
polynomial (order = 0) under both methods; it previously stopped
inside stats::poly(), and its no-fit message for the polynomial
now counts occasions instead of suggesting start, which the
polynomial rejects.
The equivalence sensitivity siblings speak their parents'
language: ss_aipe_equivalence_smd_sensitivity() takes
delta_lower / delta_upper and
ss_aipe_equivalence_r_sensitivity() takes rho_lower /
rho_upper, positive magnitudes with the upper bound required and
the lower defaulting to symmetric, exactly as in
equivalence_smd() and equivalence_r(). The signed
equivalence_lower / equivalence_upper spellings, and their
silent default of plus or minus 0.20, are gone.
Every AIPE planner that reports an expected interval width names
the row ci_width_expected and reports the expected full width:
ss_aipe_equivalence_smd() and ss_aipe_mixed_effects()
previously reported the half-width as ci_half_width_expected
while ss_aipe_equivalence_r() reported the full width, a trap
for anyone extracting the row programmatically.
The equivalence trio attaches the conf_level attribute the ci_*
family always carried, and the equivalence_c() verdict label is
spelled "Noninferior only", matching the solid noninferior the
rest of its page uses. equivalence_r() now validates raw data
(matching lengths, at least 4 complete pairs) the way its
summary-statistics path always did. convert_r_Z() and
convert_Z_r() insist on a single value, as documented; vector
input previously recycled the term column into duplicated rows.
A follow-up ruling extended the same single-value guard to the
whole scalar conversion family (convert_R2_f() and its three
siblings, convert_delta_lambda() and its inverse, and
convert_z_normal()), with domain checks where a map's algebra
ends (an R2 at or past 1, a nonpositive group size), so no
conversion can recycle vector input into duplicated rows.
Documentation corrections from the same audit: the four nonlinear
simulators state that the first-order delta method behind
reliability and reliability_by_occasion can drift from the
realized variance ratio when random variances are large, and
document the schedule attribute; the Richards example comment
places the delta = 3 inflection at about 63% of total change, per
its own formula; the mixed-effects examples on
?analysis_of_change now run live.
The four multiple-comparison interval functions moved into the
ci_* family: ci_dunnett, ci_scheffe, ci_games_howell, and
ci_tukey_kramer replace dunnett_ci, scheffe_ci,
games_howell_ci, and tukey_kramer_ci. The suffix forms were the
only four exports naming a confidence interval outside the family
prefix, and each function now pairs with its critical-value sibling
(cv_dunnett, cv_scheffe, cv_tukey_hsd). The package is
unreleased, so the old names are gone rather than aliased.
The noncentral-distribution toolkit is named by one rule
(2026-08-21): nc_ and then the distribution, so ci_nc_t(),
ci_nc_F(), and ci_nc_chisq() replace conf_limits_nct(),
conf_limits_ncf(), and conf_limits_nc_chisq() (the snake-cased
descendants of MBESS's conf.limits.nct() and kin), and
moments_nc_t() and moments_nc_F() replace moments_nct() and
moments_ncf() to match moments_nc_chisq(). The intervals join the
ci_* family because a confidence interval, on the noncentrality
parameter that the effect size intervals are built from, is what
they compute, and "confidence limits" was vocabulary that existed
only in MBESS. Splitting nc from the distribution token retires
the one inherited exception (nct, ncf against nc_chisq) and
restores the capital \emph{F} the package writes everywhere else
(F_value, convert_F_chisq). The noncentral marker stays because
the noncentrality parameter is the estimand. plot_ci() keeps its
name: it is a plot first. The old names are gone rather than
aliased; each page names its earlier self in the description.
The correlation intervals now carry the names of their estimands,
matching the rest of the correlation family (ss_aipe_r,
ss_power_r, var_r, expected_r, equivalence_r): ci_r() is the
confidence interval for the simple Pearson correlation (renamed from
the MBESS-heritage ci_cc(), with the estimate row renamed from
est_cor to r), and ci_R() is the interval for the population
multiple correlation coefficient (its former lowercase alias
spelling is gone). Cross-references across the package were audited
against this distinction and every one now points at the interval it
meant. Because the two natural file names differ only by case, which
case-insensitive filesystems refuse, both functions live in
R/ci_correlation.R and share the ci_correlation help page.
Kish's design effect function is design_effect(); the abbreviation
deft() is gone as a function name, while the returned table keeps
its design_effect and deft rows, deft being the standard term
for the standard error inflation factor, the square root of the
design effect.
The remaining exported alias pairs were resolved to single names.
Fisher's Z is written with a capital Z throughout the package,
because it is the variance-stabilizing transform of a correlation
and not a z-score, so the converts are convert_r_Z() and
convert_Z_r() only (the lowercase spellings are gone) and the
prose and math were swept to match. The limits of agreement
function is limits_of_agreement(), with loa() kept as its short
alias and bland_altman_loa() gone (the method is named for what
it computes; Bland and Altman are credited in the references).
reliability_omega_categorical() and covmat_from_cfa() are the
only names for those functions; the reliability_omega_c() and
covmat_from_cfm() aliases are gone. The dataset alias bindings
HS_Data and Prime_Time are gone as well: every dataset carries
exactly one name, the documented snake_case one, and unlike the
alias bindings the canonical names work with data().
The vocabulary of the longitudinal simulators names what the rows
are: the identifier column identifies units (persons, animals,
trees, classrooms), n counts units, and the former group column
is population, one level per data generating parameter vector.
Measurement schedules can now be unit-specific, in every simulator:
in place of a shared target_times grid, time_range = c(lower, upper) draws each unit's own measurement times uniformly between
the bounds, with occasions fixing the number of times per unit or
c(min, max) letting it vary, so designs such as age at testing
rather than grade at testing are simulated directly. The polynomial
simulator additionally requires every unit's count of occasions to
reach P + 1, so each simulated trajectory identifies the
polynomial it came from.
analysis_of_change() closes the loop the simulators open: it fits
any of the four nonlinear change models, or a polynomial of any
order, to longitudinal data. The default two-stage method fits one
curve per unit, from that unit's data alone, with data-driven
starting values, and summarizes the unit-level parameters by their
mean, their standard deviation and variance across units (the
individual differences), and the standard error of the mean;
method = "mixed" instead estimates the proper random-coefficients
model simultaneously, lme4::lmer() for the polynomial and
nlme::nlme() for the nonlinear curves started at the two-stage
estimates, so sd_units becomes a variance component purged of
estimation noise. The help page positions the function candidly
against nlme::lmList(), nlme::nlsList(), and the self-starting
SS* curves: what it adds is the landmark parameterizations, the
simulator match, and one tidy interface across linear and nonlinear
change. With a single trajectory (or
id = NULL) it reduces to one nonlinear least squares fit reported
with asymptotic standard errors, so N = 1 is one unit's change.
Non-converging units are dropped with a single counted warning and
the effective count travels as the n_used attribute; the full
matrix of unit-level estimates rides along for plotting or as
starting values for a simultaneous nlme::nlme() fit, and the help
page states plainly that the between-unit spread of estimates
includes estimation noise, which a variance-component model
separates.
simulate_longitudinal_polynomial(), sharing its interface (several
populations of units, between-unit parameter variances and
correlations, level-one error
by variance or by target reliability, error correlation structures,
assessment-time jitter) and its long-format return that feeds
plot_trajectories() and nonlinear mixed-model fitters:
simulate_longitudinal_negative_exponential() (asymptotic
regression), simulate_longitudinal_logistic(),
simulate_longitudinal_gompertz(), and
simulate_longitudinal_richards(), whose shape parameter delta
subsumes the logistic (delta = 1) and the Gompertz (delta -> 0) as
special cases, relations the tests pin exactly. The
parameterizations are those of Kelley (2005, dissertation; 2008,
Methodology), in which every parameter is a landmark of the change
process (floor, ceiling, moment of fastest change, curvature,
inflection height) and the intercept-shifting zeta frees the lower
asymptote from zero; the Richards application follows Guo, Cheng,
and Kelley (2016). A new vignette, "Nonlinear Growth Curves and the
Meaning of Their Parameters", illustrates all four curves, the
Richards unification, and the polynomial comparison: a ninth-order
polynomial needs ten uninterpretable coefficients to track a four
parameter Gompertz inside the data and still collapses the moment it
extrapolates.equivalence_smd(), equivalence_r(), and equivalence_c()
replace tost_smd(), tost_r(), and tost_c(), and the AIPE
planners follow (ss_aipe_equivalence_smd(),
ss_aipe_equivalence_r(), with their sensitivity siblings). The
whole family now speaks the same word as the planning and plotting
surfaces that always did (power_equivalence_md(),
ss_power_equivalence_c(), plot_equivalence()). The procedure is
unchanged, Schuirmann's two one-sided tests, and each help page
says so by that name; the reported TOST \emph{p}-value keeps its
standard p_tost term. The old names are gone rather than aliased.ss_aipe_equivalence_r() plans the sample size for an equivalence question
about a Pearson correlation: the smallest \emph{N} whose
100(1 - 2 alpha)% Fisher's \emph{Z} interval (the interval
equivalence_r() and ci_r() invert) has expected width at or below the
target, with an optional Monte Carlo assurance correction. At the
conservative default planning value of 0 the answer has a closed
form, and the tests anchor the search to it. Its Monte Carlo
sibling ss_aipe_equivalence_r_sensitivity() follows the family API and
reports, alongside the family's width and coverage summaries, the
realized proportion of equivalence verdicts inside the chosen
bounds. This fills the gap noted when the equivalence family was
reviewed: the SMD had an AIPE planner for its equivalence interval
and the correlation did not.depression_bdi with its wait list control, the studentized range
procedures on test_market's six panels, the unequal-variance and
rank methods on drinks_trial's skewed outcome, the factorial and
generalized eta squared examples on pygmalion's manipulated
treatment crossed with measured grade, and the correlation,
regression, reliability, and multivariate examples on the
holzinger_swineford battery, whose second-form tests carry real
missingness that now powers the FIML demonstrations in mlmr() and
mlmr_mv() in place of artificially punched holes. Examples that
fit lme4 models keep lme4::sleepstudy, since those pages require
lme4 regardless. The ecvi() example and tests likewise moved from
lavaan's copy of the 1939 data to the package's own
holzinger_swineford. Numeric claims in example comments were
recomputed against the new output throughout.Language: en-US and inst/WORDLIST carries the package's
technical vocabulary, so devtools::spell_check() runs clean.plot_equivalence() no longer passes the fatten argument of
ggplot2::geom_pointrange(), which ggplot2 4.0.0 deprecates. The
point size is set through the size aesthetic at the value that
draws the same figure on every supported ggplot2 (>= 3.4.0), so
building the plot under ggplot2 4.0 no longer raises a deprecation
warning.reliability_omega() and reliability_omega_categorical() now
validate their inputs before the courtesy message that a default
call reports no confidence interval, so a call that is about to
fail with an informative error no longer receives advice first.Config/testthat/parallel: true, with the slowest files scheduled
first through Config/testthat/start-first; the testthat floor
moves to 3.2.0). The heaviest Monte Carlo checks were recalibrated
to smaller replication counts whose assertions still hold with room
to spare, duplicated planner calls across neighboring tests were
consolidated into single calls, and tests no longer leak messages,
notes, or printed tables into the run's output. A full
devtools::test() on the maintainer's machine dropped from about
20 minutes to about 5, with no warnings and nothing skipped.ss_aipe_*_sensitivity() member now returns the family's
documented schema: mean_ci_width / median_ci_width / sd_ci_width
for the realized interval widths, pct_ci_less_w for the proportion
of intervals at or below the planning width, pct_ci_miss_low /
pct_ci_miss_high / total_type_I_error for the empirical
non-coverage (proportions on the 0 to 1 scale, so the total is the
sum of its tails), mean_X / median_X / sd_X for the member's
own estimand, and input echoes named for their unit (total_N or
n_per_group for the evaluated size, true_X, estimated_X,
width, conf_level, and, when one was supplied, assurance).
The core term vector lives once, as the internal registry constant
.SS_AIPE_SENS_CORE_TERMS in R/dmar_tidiers.R, and a family-wide
contract test asserts every member against it. Every member's
@return now lists exactly the rows the function returns.ss_aipe_R2_sensitivity()
drops the lowercase _r2 suffixes and pct_less_w (mean_r2 is
mean_R2, mean_ci_width_r2 is mean_ci_width, the realized-limit
rows are mean_lower_limit and siblings, the one-sided widths are
mean_ci_width_lower / mean_ci_width_upper); the four term-matching
reads inside ss_aipe_R2()'s assurance search were updated with it.
ss_aipe_sc_sensitivity(), ss_aipe_sc_ancova_sensitivity(), and
ss_aipe_sm_sensitivity() retire mean_full_width and
pct_Width_obs_narrower / pct_width_obs_narrower;
ss_aipe_smd_sensitivity() retires pct_less_desired;
ss_aipe_c_ancova_sensitivity() retires width_narrower and
mean_width_obs (its standard error comparison is now
mean_se_ratio); ss_aipe_rmsea_sensitivity() retires rmsea_pop,
desired_width, mean_width, and, most importantly, an output row
named assurance that actually reported the realized width-attainment
proportion, now pct_ci_less_w so the name no longer collides with
the planning input; ss_aipe_sem_path_sensitivity() retires
width_less_than_desired and the type_I_err* trio.ss_aipe_sm_sensitivity(),
ss_aipe_smd_sensitivity(), ss_aipe_sc_ancova_sensitivity() (both
divisor branches), and ss_aipe_c_ancova_sensitivity() (which also
scaled its width-attainment row). All are now proportions on the
0 to 1 scale, completing the sweep that earlier fixed
ss_aipe_sc_sensitivity(), ss_aipe_cv_sensitivity(), and
ss_aipe_reg_coef_sensitivity().divisor = "s_anova" branch of
ss_aipe_sc_ancova_sensitivity() read the confidence limits out of
ci_sc_ancova() by row position ([2, 2] and [4, 2]) when that
function returns three rows, so the branch treated the point estimate
as the lower limit and NA as the upper: its realized widths were NA
and its "Type I error" rows were nonsense. The limits are now read by
term name, as the s_ancova branch always did, and the branch's
summaries are meaningful for the first time.mean_psi and siblings in the three contrast members,
mean_sm, mean_smd, mean_rmsea, and mean_path, computed from
the same replications the widths come from. The echo rows the schema
calls for were added where missing (ss_aipe_R2_sensitivity(),
ss_aipe_cv_sensitivity(), ss_aipe_reg_coef_sensitivity() and its
rc / src wrappers, the contrast members, and
ss_seq_c_sensitivity(), which now echoes half_width, true_psi,
true_sigma, alpha_level, and m0).ss_power_R2_sensitivity() and
ss_power_reg_coef_sensitivity() echo their planning inputs (p,
true_R2 / true_b_j, estimated_R2 / estimated_b_j,
desired_power, NA when a size was specified directly, and
alpha_level), and their realized-R^2 rows follow the
meaningful-capital rule (mean_R2, and in the omnibus member
mean_F / F_crit). tidy() and glance() for the
dmar_ss_power_sensitivity class are unchanged; the echoes ride
along in glance().ss_aipe_icc() accepted and documented the six Shrout-Fleiss type
labels but planned every one of them on the single-rater scale, so
average-of-k plans were materially undersized. The planning value and
target width are now interpreted on the scale of the requested form:
an average-of-k value is mapped to the single-rater scale through the
inverse Spearman-Brown relation and each candidate confidence limit
is mapped back (the convention var_icc() uses), and type is now
validated, so an unrecognized label fails instead of silently
planning ICC(1,1). Average-of-k recommendations move: at
rho = .70, k = 3, and width = .20, planning for ICC(1,k) now
recommends n = 110 (realized mean width .20 across 2,000 Monte Carlo
replications of the F-based interval) where the undersized plan
recommended n = 69 (realized width .26); at width = .10 the
recommendation is n = 421, and at rho = .90, width = .10 it is
n = 52. Single-rater plans are unchanged. The planned form travels
as the icc_type attribute on the returned table.ss_aipe_icc_sensitivity() had the matching defect on the generator
side: it treated true_rho as the single-rater ICC no matter the
type. For the average-of-k forms it now simulates data whose
population ICC at the average-of-k level equals true_rho, so the
realized estimates, widths, and coverage refer to the form being
planned; type is validated the same way.?ss_aipe_icc that the assurance correction
over-recommends sample size by 25 to 40 subjects was backwards. At
the page's own condition (rho = .70, k = 3, width = .20,
assurance = .80) the recommended n = 79 delivers an empirical
assurance of about .77, and n = 81 is the smallest sample size that
reaches .80 (10,000 Monte Carlo replications). The Details section
now says so and points to ss_aipe_icc_sensitivity() for checking a
strict assurance target.alpha_lower critical value of
the central F-distribution, the lower noncentrality limit is 0 by
construction; that is a normal consequence of a small observed
effect, not a failure. ci_nc_F() still warns once in that
case, but the message now states the consequence for the interval
(the lower confidence limit is 0) instead of describing achieved
tail probabilities, and every function that builds its interval by
inverting the noncentral F through it (ci_snr(), ci_srsnr(),
ci_pvaf(), ci_R2() with fixed predictors, ci_eta_squared(),
ci_eta_squared_partial(), ci_omega_squared(),
ci_eta_squared_generalized() with the parametric method, and
ci_mahalanobis()) restates it for its own effect size, for example
"the lower confidence limit on the signal-to-noise ratio is 0",
at most once per call. Previously ci_snr() and ci_srsnr()
surfaced the inner wording, which pointed users to a prob_greater
column those functions do not return. The warning carries the
condition class dmar_nc_F_clamp, which the iterative callers
(ss_aipe_R2(), ss_aipe_omega_squared(), factorial_anova(),
simple_effects_AB()) now match by class when deduplicating.uniroot() message: the error now names the function that was
called, reports the F-statistic, degrees of freedom, and tail
probabilities involved, and says what to try.df, the significance lens of design_consequences()
(power, type_s_error, exaggeration_ratio) is now computed from
the noncentral t distribution of the test statistic, the exact
distribution when the standard error is estimated from the data and
the same sampling model the precision lens already used; the
exaggeration ratio integrates the truncated normal moments over the
chi distribution of the estimated standard error. The previous code
evaluated a location-shifted central t, the known-se approximation
behind Gelman and Carlin's retrodesign(), so the help page promised
the noncentral t while the code delivered the approximation.
Small-df results move: at n = 5 per group with a true effect of
d = 1, power is now 0.2863 (was 0.2469; a two-million replication
Monte Carlo of the design gives 0.2860 with simulation standard
error 0.0003, and base R's power.t.test(strict = TRUE) agrees
exactly), the Type S error is 0.00129 (was 0.00937), and the
exaggeration ratio is 1.658 (was 1.915). The differences fade as
df grows (negligible by about 60 per group), and the
df = Inf normal case is unchanged.smd_trimmed() applied the reciprocal of the Algina-Keselman-Penfield
scaling constant at every trimming proportion other than the 0.20
default: the internal constant computed the Winsorized variance of a
standard normal correctly and then returned 1 / sqrt(win_var) where
the definition calls for sqrt(win_var), and a hard-coded 0.642 at
trim = 0.20 masked the error at the default while making the
estimate discontinuous there (trim = 0.1999 returned 2.43 times the
value trim = 0.20 returned). The constant is now sqrt(win_var) at
every trim and the hard-coded branch is gone. At the default the
constant moves from the rounded 0.642 to its exact value 0.6419398,
so estimates at trim = 0.20 change by less than one part in ten
thousand; at any other trim the correction is substantial (previous
estimates were 1.47 times too large at trim = 0.10, 4.99 times too
large at trim = 0.30, and 17.99 times too large at trim = 0.40).
The help page identity is corrected to match: 0.642 is
SD(X_W) / SD(X), that is sqrt(Var(X_W) / Var(X)), for a standard
normal Winsorized at 0.20, not Var(X) / Var(X_W).smd_trimmed()'s help page promised a confidence interval on the
Yuen-Welch degrees of freedom, but the code inverted the noncentral
t at h_1 + h_2 - 2 degrees of freedom with noncentrality
d_R * sqrt(h_1 * h_2 / (h_1 + h_2)); the Yuen-Welch value was
computed, returned in the df_yuen row, and never used. The
interval now follows the construction in Keselman, Algina, Lix,
Wilcox, and Deering (2008): Yuen's t-statistic on the
trimmed-mean difference (their Equation 8) is referred to a
noncentral t distribution with the Yuen-Welch approximate degrees
of freedom (their Equation 9), and the noncentrality limits are
rescaled to the d_R metric. Feeding the summary statistics of the
paper's worked example (their Tables 1 and 3) through this
construction reproduces the printed robust effect size intervals
[0.31, 3.37] and [0.10, 1.11] (p. 119) to the precision the rounded
published inputs support, while h_1 + h_2 - 2 degrees of freedom
give [0.40, 3.31] for the first: the paper's intervals use the
Yuen-Welch value. The former noncentrality also understated the
observed statistic (by about 16% at the default trim under equal
Winsorized variances), so intervals move even where the two degrees
of freedom nearly agree: on the help page example the interval is
now [-0.88, 0.08] where the old construction gave [-0.97, 0.17].lin_ccc() now builds its confidence interval on Lin's (1989)
z-transformed standard error, with the correction noted in Lin
(2000), and the former default method = "king_chinchilli" is
removed. The removed variance was not a consistent estimator of the
sampling variance of the z-transformed CCC: when the two means and
variances are equal it inflates the correct asymptotic variance by
exactly the square of (1 + rho^2) / (1 - rho^2), already about
2.8-fold at a CCC of .5 and several hundred-fold at .95 (and by a
comparable factor otherwise), so its intervals were nearly vacuous.
The help page's own first example returned CCC = .928 with CI
[-0.977, 1.000]; the same example now returns [0.871, 0.960],
matching the independent DescTools::CCC() z-transform interval
to ten decimals, and the Lin interval's simulated coverage is .952
at nominal .95 (bivariate normal, rho = .5, n = 50, where the
removed default covered .982 with intervals 61% wider). The removed
formula was also not King and Chinchilli's (2001) estimator, so
keeping it under that name would have credited a wrong formula to
real authors; a correct King-Chinchilli variance may return in a
later release. Calls that request method = "king_chinchilli" now
fail loudly.fleiss_kappa() computed its z statistic and p-value from a
null variance that treats the category marginals as known constants,
not from the estimated-marginals null variance of Fleiss, Nee, and
Landis (1979), the correction of the standard errors in Fleiss
(1971), even though the help page credited that paper. The two
expressions coincide at uniform marginals, which is how the error
hid, and diverge as the marginals skew: at N = 1000, m = 5, and
marginals (.85, .10, .05), the old null standard deviation was 3.8
times the empirical one (.0303 versus .0080). Because the
known-marginals expression is never smaller than the corrected one
(they are equal only at uniform marginals), every affected z was
understated and every affected p-value overstated; the test was
conservative, never anticonservative. On the help page's own Fleiss
(1971) Table 1 example the z statistic moves from 15.64 to 17.65,
now matching irr::kappam.fleiss to ten decimals. The point
estimate and the confidence interval do not move: kappa is
unchanged, and the se, lower_limit, and upper_limit columns
come from the Gwet (2008) linearization, which was correct all
along.ci_dunnett(alternative = "less") formed its simultaneous upper
bounds with the signed critical value that cv_dunnett() reports for
that alternative, a negative lower-tail quantile, which placed every
bound at diff - |d| * se: below the point estimate it was supposed
to bound, and in contradiction with the adjusted p-values printed
beside it (on the PlantGrowth example both "less" rows reported an
upper limit below zero, a rejection, beside adjusted p-values of
0.162 and 0.989, no rejection). Simultaneous coverage of the true
differences measured 0.00115 against the nominal .95 in a balanced
null design with four groups and n = 10 per group. The critical
value's magnitude is now applied on the side the alternative
dictates, so the "less" bound is diff + |d| * se; coverage in the
same design measures 0.94685, and the PlantGrowth "less" upper
limits move from -0.928 and -0.063 to 0.186 and 1.051, agreeing with
multcomp::glht()'s one-sided limits to about four decimal places
(the remaining daylight is multcomp's simulated-quantile error) and
with the adjusted p-values on every row. The
alternative = "greater" interval already used the positive
critical value and is unchanged, as are the two-sided intervals and
all adjusted p-values.anova_within_two_way()
previously reported these artifacts as estimates. The
Greenhouse-Geisser and Huynh-Feldt rows for such an effect now carry
NA in epsilon, the adjusted degrees of freedom, and the p-value,
with one warning naming the condition; car::Anova likewise
declines to report the corrections when the effect's error matrix is
singular. The unadjusted and lower-bound rows are unchanged: the
lower bound 1/df is Geisser and Greenhouse's a priori bound, not an
estimate, and remains valid however few the subjects (Maxwell,
Delaney, & Kelley, 2027, Chapters 11 and 13).ci_r() accepted any sample size and went quietly wrong below
n = 4: at n = 3 it returned the vacuous interval [-1, 1] (the
Fisher's Z variance 1/(n - 3) is infinite there), and at n = 2 it
returned NaN limits. The function now stops with an informative
error for n < 4, and the n = 3 message explains why no usable
interval exists at that boundary.ci_method = "fisher" in
reliability_alpha(), reliability_omega(), and
reliability_kr20(); the Fisher interval there now stops for N < 4
with the same explanation. The package's other Fisher's Z consumers
(var_r(), equivalence_r(), meta_r(), ss_power_r(),
correlations_test(), correction_for_attenuation(),
expected_r()) already validated their sample sizes.ss_aipe_*_sensitivity family scaled some or
all of their coverage and width-attainment rows by 100 while the
rest of the family reported proportions. ss_aipe_sc_sensitivity()
reported type_I_error_upper and type_I_error_lower as
percentages; ss_aipe_cv_sensitivity() reported pct_ci_less_w,
pct_ci_miss_low, and pct_ci_miss_high as percentages beside a
total_type_I_error that was already a proportion, so its own rows
did not add up; and ss_aipe_reg_coef_sensitivity() (inherited by
ss_aipe_rc_sensitivity() and ss_aipe_src_sensitivity()) scaled
all four. These rows are now proportions on the 0 to 1 scale
everywhere, so total_type_I_error equals the sum of
pct_ci_miss_low and pct_ci_miss_high, and a realized Type I
error compares directly to 1 - conf_level with no per-function
rescaling. Term names are unchanged; only the scale moved.mediation_mbco()'s null hypothesis for an indirect effect is a
union of branches (the product is zero when any factor is), and the
constrained search could converge to a worse-fitting branch on some
platforms and OpenMx builds, inflating the likelihood ratio
statistic. When every constrained algebra is a pure product of free
parameters, the branches are now also fit directly as ordinary
unconstrained models with one factor fixed to zero, and the reported
statistic is defined by the best-fitting branch on every platform.
Non-product constraints keep the multi-start constrained search.p_value, se, ci_lower, ci_upper, p_adjusted,
conf_level, std_estimate, R2, adj_R2, df_residual, and so
on, replacing broom's dotted p.value / std.error / conf.low
across every tidier in the package. One naming system now covers
every DMAR surface; a pipeline that feeds a tool expecting broom's
dotted schema renames the columns at that boundary.content_validity_index(),
dmacs(), measurement_invariance(), and measurement_alignment()
gained bespoke tidy()/glance() pairs (items, ladder rungs, and
groups as terms; the scale and model level summaries as the one-row
glance), and a generic wide branch in the default methods serves
htmt(), average_variance_extracted(), and any future wide table
(term from the label columns, estimate from the first numeric
column, remaining columns passed through). The output vignette's
promise that the verbs answer everywhere is now literally true.B is the number of bootstrap replications everywhere. mlmr()
and mlmr_mv() rename boot_R, R2_mixed_effects() and
ci_eta_squared_generalized() rename R, and
krippendorff_alpha() renames n_boot; the effective-count row is
now B_used. Old argument names fail loudly.cohen_kappa() and fleiss_kappa() gain bootstrap intervals
(ci_method = "percentile" or "bca", B = 10000, seed),
redeeming the kappa page's own advice that the Wald interval can
have poor small-sample coverage (Blackman & Koval, 2000; Zapf,
Castell, Morawietz, & Karch, 2016). Subjects are resampled, table
input is expanded to the equivalent pairs, and the se, z_value,
and p_value columns keep their asymptotic definitions; only the
interval changes.average_variance_extracted() gains a percentile bootstrap
(ci_method = "percentile", B = 1000, refitting the model per
replication) and now always returns ci_lower / ci_upper columns
(NA under ci_method = "none"). No interval is possible from
loadings alone, and the function says so.krippendorff_alpha() no longer bootstraps by default
(boot = FALSE), matching the package-wide rule that no analysis
runs a bootstrap unless asked; its verbal reliability benchmarks
were removed in favor of reporting the coefficient with its
interval.mediate() no longer errors when a degenerate resample returns no
indirect effect, and htmt() no longer aborts when a resample
leaves a block's average within-construct correlation nonpositive.
Both drop the replication, warn once with the count, and stop only
when fewer than 100 replications survive.\donttest{} or \dontrun{}
any more. Under R CMD check --as-cran, one \donttest{} block
anywhere makes R run the whole example corpus twice, so the wrapper
cost time rather than saving it. The examples that were slow are
now fast. Replication counts in the cheap demonstrations were
lowered, with a comment naming what a reported analysis deserves.
No example runs a bootstrap confidence interval: those calls
are carried as commented code, so the syntax is still on the page
for a reader who wants it, and the surrounding prose explains the
interval and when to ask for it. The randomization tests keep their
resampling, since permuting the data is what those functions do and
their examples cost a tenth of a second. Anything else expensive is
commented out the same way, each passage introduced by a sentence
saying what it does and why it is not run. The data behind every
example are unchanged. The slowest help page now takes 0.38
seconds, and all 302 together take 8 seconds.cfa_1() is now the one factor special case of cfa_k(), a
convenience wrapper rather than an independent implementation, and
the new cfa_2() is the two factor sibling. cfa_1() only
requires the data (and, when the data hold more than the items, a
vector of item names); cfa_2() takes the items of each factor as
factor_1 and factor_2. Both forward everything to cfa_k(), so
there is now a single fitting implementation, one output schema
(with confidence interval columns), and one set of fit index
choices; under a robust estimator every CFA surface now reports the
robust index versions. The former cfa_1() extras moved or
dissolved: composite reliability with the observed denominator is
reliability_omega(denominator = "observed"), the omega output
mode is the omega_f1 row of the standard table, and the term
names follow the cfa_k() convention (lambda_f1_<item>,
psi_f1_<item>, omega_f1).data and S are now separate arguments across the CFA family.
Raw data are passed as data and a covariance matrix as S (with
N); an argument never means both. cfa_k() refuses a square
symmetric matrix passed as data, naming the fix. cfa_1() keeps
its conveniences: unnamed input is auto-named y1, y2, ..., and
items defaults to every column.measurement_invariance() now fits any measurement model, not only
the one-factor case. It accepts model as lavaan syntax or as a
named list mapping factors to their items (the items argument
remains as the one-factor convenience), along with ordered,
missing for full information maximum likelihood, group_partial
for the partial invariance case, and parameterization. With ordered
indicators the ladder itself changes: thresholds carry the location
information, so a thresholds rung is fitted between configural and
metric, following Wu and Estabrook (2016) and Millsap and Yun-Tein
(2004). Constraining loadings before thresholds, as the continuous
ladder does, tests the wrong hypothesis for ordered items. The chi
square difference test is the scaled one whenever the estimator is
robust or the data are ordered, and the fitted lavaan objects are
returned on a "fits" attribute so a caller can run score tests and
partial invariance refits without paying for the ladder twice.
dmacs() reports the dMACS effect size of measurement
noninvariance (Nye & Drasgow, 2011): the expected difference between
two groups' measurement equations for an item, integrated over the
focal group's latent distribution and standardized by the pooled item
standard deviation. A score test says a loading or intercept differs
detectably; dMACS says whether the difference is large enough to
matter. The defining integral has a closed form under a normal latent
variable, so no numerical integration is used, and the two agree to
machine precision in the tests. Accepts a fitted multiple group
lavaan model or the parameters a paper reports. The help page states
plainly that the index is meaningless from a configural fit, where
the groups share no metric.
measurement_alignment() implements the alignment method of
Asparouhov and Muthen (2014). Exact invariance essentially never holds
across many groups, which leaves the ladder stalled at configural and
group comparison blocked. Alignment estimates the group factor means
and variances that make the measurement parameters as nearly invariant
as possible, minimizing a simplicity function whose fourth-root
component loss tolerates a few large differences and punishes many
small ones. Both the fixed and free identifications are available, the
optimizer runs from several starts because the surface has local
minima, and the number of distinct optima found is reported so a
solution is never presented as unique when it is not.
irt_grm() fits Samejima's (1969) graded response model. It does
so the way the rest of the package works, by fitting the categorical
factor analysis model with lavaan and converting the solution:
discrimination a_i = lambda_i / sqrt(1 - lambda_i^2) and location
b_ik = tau_ik / lambda_i. The two models are the
same model in different parameterizations (Takane & de Leeuw, 1987),
so this adds item response theory without a second estimation engine
and without leaving the factor analytic tradition. Both the normal
ogive and the logistic metric are available.
irt_information() and plot_irt_information() give the item
and test information functions and the standard error of the latent
trait, SE(theta) = 1 / sqrt(I(theta)). This is precision as a
function of where the respondent sits on the trait, which is what a
single reliability coefficient cannot express: a scale can have
excellent omega and still measure poorly over the range a study cares
about. Verified against the closed form for the dichotomous case,
which the graded model must reproduce exactly.
content_validity_index() computes the item level content
validity index from a panel of expert relevance ratings, with the
modified kappa that corrects it for chance agreement (Polit, Beck, &
Owen, 2007), Lawshe's (1975) content validity ratio, and the scale
level S-CVI/Ave and S-CVI/UA. Each I-CVI carries an exact binomial
confidence interval, because an index computed from five or six
experts is a proportion with real uncertainty and reporting it as a
bare point estimate overstates what a small panel establishes. This
is evidence about the items, gathered before any data are collected.mediation_mbco() implements the model-based constrained
optimization (MBCO) procedure of Tofighi and Kelley (2020,
Psychological Methods): a likelihood ratio test of any smooth
function of path coefficients (an indirect effect, a total effect, a
contrast of two indirect effects) formed by refitting the mediation
model subject to the nonlinear constraint that the function equals
zero. The model is specified in lavaan syntax (observed or latent
variables; parallel or sequential mediators) and fit in OpenMx,
whose optimizers support the nonlinear equality constraints the
null model requires. The function enumerates the total, direct,
total indirect, and every specific indirect pathway from x to
y, reports each with a delta method standard error and a profile
likelihood, Monte Carlo, or Wald confidence interval, tests each
with the MBCO likelihood ratio statistic and its p-value, and
reports AIC and BIC differences and the change in each endogenous
R-squared under every null model. Because the null set of a product
constraint is a union of surfaces, each null model is refit from
several starting configurations and the best feasible solution is
kept, so the reported statistic reflects the globally best-fitting
null model rather than a local branch. Accepts raw data or the
summary statistics a paper reports (covariance matrix, means, and
sample size), which reproduce the raw-data analysis exactly; the
tests replicate the published empirical example from its Table 1
moments. A group argument turns the model into a multiple-group
SEM in which every effect is estimated per group and its
between-group difference is tested, which is moderated mediation
with a categorical moderator. A moderator argument probes
continuous moderation stated in the syntax with : interaction
terms (or precomputed product columns): each moderated pathway
effect becomes a symbolically derived polynomial in the moderator,
reported as conditional effects at probe values (the mean and one
standard deviation either side by default), the index of moderated
mediation (Hayes, 2015) tested by likelihood ratio rather than
bootstrap, and, when a pathway is moderated in several places, a
joint constancy test whose null model imposes several nonlinear
constraints at once. hypotheses accepts a named list whose
multi-expression elements are likewise tested jointly. Constrained
null models are started on every branch of the constraint's null
set (including exact conditional-coefficient branch starts for
probed effects) and each solution is polished by a warm restart, so
the reported statistic reflects the best-fitting feasible null
model. Structural guardrails refuse pathway enumeration for
nonrecursive (feedback) structures and warn on binary endogenous
variables, on interaction terms no declared moderator accounts for,
and on interactions missing their main effect (the principle of
marginality).plot_mediation_mbco() draws the conditional effects of a
moderated mediation_mbco() analysis as curves over the
moderator's observed range, one per moderated pathway effect, with
a pointwise Monte Carlo confidence band, the probed values marked,
a dashed zero line whose band crossings estimate the Johnson-Neyman
boundaries, and a rug of the observed moderator values so
extrapolation is visible. The band's polynomial is the same stored
quantity the table probes, evaluated over the whole range; the help
page is explicit that the band is pointwise and that boundary
locations read from it are estimates, with the table's constancy
test as the formal companion. Returns a plain ggplot object in the
Okabe-Ito palette.var_ete() computes the variance of the estimated treatment
effect at selected covariate values in a two-group ANCOVA with
heterogeneity of regression and a random covariate (Li, McLouth, &
Delaney, 2020), the reimplementation of MBESS::var.ete(); tested
against the MBESS reference on every branch and dogfooded in the
Pygmalion vignette.
tidy() and glance() answer on every DMAR result. Two default
methods on the dmar_tbl class are the floor under the package: any
result table answers tidy() with the broom-shaped long view (term,
estimate) and glance() with the one-row wide view, repeated terms
disambiguated rather than dropped. The family methods sit ahead of the
defaults and keep their richer views, and two families join them: the
ss_aipe_* planners carry a new dmar_ss_aipe class whose tidy()
reports the planned size beside the desired width it was planned
against (term, estimate, width), through the same size registry
the power planners use and a width registry beside it, and whose
glance() widens the echoed inputs. ss_power_pcm(), the one power
planner missing its family class, joins dmar_ss_power, so its
tidy()/glance() no longer error, and
ci_eta_squared_generalized() joins its siblings in dmar_ci_anova.
One sample size vocabulary across every planner. A planner's answer
row is now named for what it counts: necessary_N (total),
necessary_n_per_group, necessary_n_per_cell, or
necessary_n_clusters; a user-fixed size echoes as specified_*; and
total_N appears only as the implied-total companion beside a per-unit
answer, never as an answer's name. Before this sweep the same word
meant different things in different planners: sample_size was the
total in ss_aipe_R2() and its seven siblings but the per-group size
in ss_aipe_c() and its three, ss_aipe_smd() said
sample_size_per_group, ss_power_pcm() carried the MBESS-era ss_c
and ss_t pair (two rows for one number in a balanced design, now one
branch-named row), ss_power_factorial_ancova() used one bare
n_per_cell for both its planned and its user-supplied branch, and the
cluster planners' answers had no prefix at all. Sensitivity outputs
echo the evaluated size under its bare unit name (total_N,
n_per_group), and the six sensitivity functions whose specified_N
argument actually meant a per-group size now call it n_per_group
(ss_aipe_c/sc/sc_ancova/smd/equivalence_smd/pcm_sensitivity). Two term
names in the cluster family that contained literal spaces are
underscored. The tidy()/glance() size recognizer no longer accepts the
legacy names, so a stray old producer fails a test instead of slipping
through, and a 216-table characterization grid recorded before the
sweep reproduces identically after it, so only names moved, never
values. The package is unreleased; old names fail loudly.
The number of predictors is p everywhere. ci_rc(),
ci_src(), ss_aipe_rc(), ss_aipe_src(), ss_power_rc(), and
plot_R2() renamed J to p, and ci_R() renamed K to p,
matching ci_R2(),
ci_reg_coef(), and the rest of the regression family. J remains
only in the partial and semipartial correlation family, where it
counts the variables partialed out (a different quantity; J = 0 is
the simple correlation). The package is unreleased, so the old
names fail loudly rather than being aliased.
ci_rc() and ci_src() are documented as the thin wrappers they
are around ci_reg_coef(), the general engine, with titles that
distinguish the unstandardized, standardized, and general cases.
Every function help-page title is now AP title case with no trailing period, matching base R convention (previously the package mixed sentence case and title case).
A numerical audit compared every displayed formula and quoted magnitude on these pages against the quantity the function computes. In each case the computation was correct and the page was not, so no returned value changes; the documentation now matches the code, and new tests recompute each corrected expression independently and assert agreement with the function output.
expected_partial_r() displayed the expectation with every index
one lower than the Olkin-Pratt formula under the n to n - J
substitution the code makes; as printed, the formula can exceed 1, an
impossible value for the expectation of a correlation. The page now
prints E[r] = rho 2F1(1/2, 1/2; (n - J + 1)/2; rho^2)
Gamma((n - J)/2)^2 / (Gamma((n - J - 1)/2) Gamma((n - J + 1)/2)),
which matches the function to machine precision. The bias prose also
described E[r] - rho while the returned bias column is rho - E[r];
the sign convention now matches the column, and the quoted magnitudes
are the true +0.00624 (rho = 0.4, n = 30, J = 2) and +0.00888
(J = 10).expected_r() quoted bias magnitudes with the wrong sign and
size (-0.024 and -0.006); the true values are +0.021 (rho = 0.5,
n = 10) and +0.0065 (n = 30), as the page's own example already
reported.icc_lmer() printed the variance of the Bonett (2002)
L-transform as 1 / (2 (n - 2)), omitting the k / (k - 1) factor the
code applies; the page now prints k / (2 (k - 1) (n - 2)). A reader
hand-building the interval from the old page at k = 2 would have an
SE too small by a factor of sqrt(2) and roughly .84 coverage instead
of .95.anova_within_two_way() printed the lower-bound epsilon as
1 / (df - 1); the code computes 1 / df, the attainable infimum, with
df the effect's numerator degrees of freedom.loa() printed a symmetric central t interval for the CIs on the
limits of agreement; the code computes the Carkeet (2015) exact
noncentral t interval, which is asymmetric about the sample LoA. The
page now shows the noncentral form: quantiles of the noncentral t
with noncentrality parameter k sqrt(n) for the upper LoA and its
negative for the lower.variance_components_mls() printed a symmetric interval with
unsquared constants and a -MS_b MS_w / n cross term that appears
nowhere in the method; the code computes the genuine Burdick-Graybill
(1992, equations 2.4.1--2.4.5) form, and the page now prints V_L and
V_U with the squared constants and the G_12 / H_12 cross terms.ss_aipe_mixed_effects() printed Var(beta-hat) with a trailing
design-effect factor the code does not apply; for the
cluster-mean-centered level-1 slope the variance is
sigma2_y (1 - rho_I) / (N sigma2_x), and the page now explains why
the design effect does not enter.ci_smd() pointed paired designs to ci_smd_c(), which is the
interval for Glass's estimator (two independent groups, control group
SD) and takes no correlation between paired measurements, so
following the pointer reproduced essentially the independent groups
interval. The page now states plainly that a paired-design SMD
interval is not currently provided.power_fisher_exact() attributed the alternative distribution to
Wallenius; the code computes Fisher's noncentral hypergeometric, the
conditional distribution of one binomial count given the total of two
independent binomials, which is what conditioning on the margins of
the 2 x 2 table produces. The prose and references now cite Fisher
(1935) and Fog's (2008) sampling-methods paper.ci_scheffe() said the default returns a - 1 pairwise contrasts;
it returns all a (a - 1) / 2 of them, as the sibling pages already
stated.A continuation of the audit above, covering another two dozen help pages. As before, the computations were correct and the pages were not; no returned value changes, and where a corrected number is load-bearing a test now recomputes it independently.
ci_rmsea()'s
90 percent example described its upper limit as landing below 0.05;
the limit is 0.052, just above the Browne and Cudeck close fit
threshold, and the page now draws the conclusion that follows (close
fit is not established, even though the point estimate sits below the
threshold). obrien_test()'s Hunter example quoted p = .2595, which
matches neither the quoted F = 1.29 nor the exact statistic; the page
now quotes p = .260, as computed. combine_p()'s Edgington sum for
the Raudenbush example prints 7, not "near 6.9". ci_c()'s
hypertension example described 24 subjects across group sizes of 4,
6, 5, and 5; it says 20.mediation_mbco() separates the published memory-example
likelihood ratios from what its example computes. Run from the
rounded Table 1 summary statistics, the two null branches give
LRT = 71.31 (the best-fitting branch, the statistic reported) and
179.02; the 72.54 and 175.77 the page previously quoted are the
full-precision-moments values, and 175.77 is now attributed to
Tofighi and Kelley (2020) as the published value, from the
worse-fitting branch on which their optimizer stopped. The page notes
that the difference comes from running from the published (rounded)
summary statistics.ss_power_indirect_effect()'s example no longer equates the
approximation with the simulation benchmark. The joint significance
approximation returns necessary_N = 65 for a = b = .39 at power
.80; raw-data simulation puts the power at 65 nearer .77 and reaches
.80 near N = 70, in line with the somewhat larger requirement in
Fritz and MacKinnon's (2007) simulation-based table, and the example
comment now says exactly that. The raw-data validation test's
relative tolerance is widened from 0.03 to 0.06 with the reason
recorded in place: the approximation sits about 3 percent above the
simulation at those settings on any seed, so the old tolerance failed
on some seeds for the approximation gap alone.prime_time_achievement's corporation summaries are recomputed
from the derived corp_id key. The bare corp column merges the
two corporations that share code 2400, so the page overstated the
largest corporation: corporations have 15 to 756 students (median
117), not 15 to 808 (median 124), and the three-level null model
decomposition is 16.29 / 22.72 / 240.43 with a corporation ICC of
0.058 (previously 16.50 / 22.70 / 240.43 and 0.059). The data set's
tests now anchor the corrected values.bessel_errors' expected counts imply a normal sigma near 0.22,
not 0.2. A least squares fit of the half-normal bin probabilities
to Bessel's expected frequencies gives 0.216; at 0.2 the first bin
alone is off by about eight observations.which_width contract of ss_aipe_omega_squared(),
ss_aipe_icc(), and ss_aipe_partial_r() now says what the code
does. The pages described "Lower" / "Upper" as a one-sided
half-width; both settings interpret width as half the full width
and return the same sample size. The pages now say so, note that an
asymmetric interval's realized half-widths differ from each other
and from half the full width, and state that a genuinely one-sided
width target is not currently offered.ci_mahalanobis() documents the tail arguments its backend
accepts. The page promised that supplying alpha_lower and
alpha_upper recomputes conf_level; that call errors, because the
tails pass straight through to ci_nc_F(), which refuses a
non-NULL conf_level beside them. The page now directs users to set
conf_level = NULL and supply both alphas, and a test pins the
contract.ci_nc_F() and ci_nc_chisq() state the
monotonicity each search actually exploits. The pages claimed each
tail probability is strictly decreasing in the noncentrality
parameter; the lower-limit condition works on the upper tail, which
is strictly increasing. Each page now gives the direction per tail
and notes that both roots are located on the decreasing lower-tail
scale.ss_aipe_equivalence_smd() described an
expected half-width target of omega where the code targets omega / 2
(full width omega), and a 1-row return where the table has 5 rows;
both now match, and the return's stale term names are corrected.
The stale sample_size return term on the ss_aipe_c(),
ss_aipe_c_ancova(), and ss_aipe_sm() pages now reads
necessary_n_per_group, necessary_n_per_group, and necessary_N.
ci_eta_squared()'s N is documented as the total number of
observations (for aovlist fits, one more than the sum of all
effect and residual degrees of freedom), matching its Details.
mlmr()'s adj_R2 is documented as using the complete cases
(N_complete), which the code deliberately uses, rather than the
lavaan reported N. ci_c_ancova_bp()'s contrast_type entry
dropped inline width expressions that omitted the 1/sqrt(2) and
sqrt(1/n) factors and now points to the correct display in Details.
ci_sc()'s note no longer contains a sentence directing users to
pass the error variance where the argument (and the rest of the
page) wants the error standard deviation. The multiplicative
diff_size illustration on the ss_aipe_crd_* pages was a verbatim
copy of the additive one; it now works the multiplicative form's
arithmetic (cluster size 25 with diff_size = c(0.8, 1, 1.2) gives
sizes 20, 25, 30, recycled across clusters).equivalence_c() performs the two one-sided tests procedure and the
companion noninferiority test for a linear contrast of group means
against bounds stated in raw units of the response, with one pooled
error term, a summary-statistic and a direct (estimate, SE, df)
interface, a benchmark argument for comparisons against a known
constant, and a five-way verdict (equivalent, superior, inferior,
non-inferior only, inconclusive). Reproduces
emmeans::test(..., side = "equivalence" | "noninferiority") to
machine precision.power_equivalence_c() computes the exact TOST power for a
contrast by integration over the chi distribution of the estimated
error standard deviation, generalizing power_equivalence_md() to
arbitrary weights and unequal group sizes, and the noninferiority
power in closed noncentral t form.ss_power_equivalence_c() finds the smallest per-group sample
size whose exact TOST (or noninferiority) power reaches a target,
the declaration-probability counterpart of the width-targeting
ss_aipe_c().ss_seq_c() and ss_seq_c_sensitivity() implement the
purely sequential fixed-width confidence interval procedure for a
contrast of Chattopadhyay, Bandyopadhyay, Kelley, and Padalunkal
(2025), with cost-optimal allocation across groups; the sensitivity
sibling verifies first-order efficiency and near-nominal coverage by
Monte Carlo.plot_equivalence() draws contrast estimates and their
intervals against the equivalence region, colored by the equivalence_c()
verdict.R2_mixed_effects() reports the Nakagawa and Schielzeth (2013)
marginal and conditional R-squared for a fitted lme4 or nlme
mixed model. It agrees to machine precision with r2mlm and with the
direct Johnson quadratic form across the model classes the tests pin:
random-intercept, correlated random-slope, split (||) random-slope,
multi-slope, and nested two-level models. On split random-slope models
performance/insight omits the random-slope variance from its
marginal and conditional R-squared, so DMAR matches r2mlm and the
Johnson form there and, by design, not performance::r2_nakagawa; the
two agree on random-intercept and correlated-slope models.R2_mixed_effects_decomposition() implements the Rights and Sterba
(2019) integrative framework of mixed-effects model R-squared measures: the full family of
total, within-cluster, and between-cluster measures from a complete
five-source decomposition of the outcome variance, matching the
authors' r2mlm reference implementation to machine precision across
the same tested model classes.knit_print.dmar_tbl() renders
DMAR results as formatted tables in knitted documents, as_kable()
produces a knitr::kable view, and results_sentence() writes an
APA-style "estimate, CI" sentence from any interval-carrying result.broom support (tidy() / glance()) added for the elementary
tests (welch_t(), summary_t_test(), contrast_test()) and the
simultaneous-comparison intervals (ci_dunnett(), ci_tukey_kramer(),
ci_scheffe()).tidy() and glance() summarize every closed-form ss_power_*
planner that reports one size and one power, from the effect size
planners to the ANOVA, ANCOVA, contrast, cluster, and mediation
designs, reporting the design's planning unit (per group, per cell,
per subject, or per cluster, or the total for the one-way ANOVA) beside
its power. The lookup that resolves the size row learned the
design-specific names (n_per_cell, n_subjects, J_per_arm,
ss_per_group), so the meaningful names are kept rather than
homogenized. The Monte Carlo sensitivity siblings
ss_power_R2_sensitivity() and ss_power_reg_coef_sensitivity() carry
a dmar_ss_power_sensitivity class whose tidy() places the empirical
and analytic power side by side.ss_power_composite_ancova_2group() plans the per-group sample size for
composite power in a two-group ANCOVA: the probability that the group
effect, the covariate effect, and the group by covariate interaction
are all significant in the same study, the quantity a design must be
planned against when its conclusion needs more than one result to hold
at once. Composite power is not the product of the marginal powers,
because every test divides by the same error estimate and the tests are
positively dependent even when the effects are orthogonal. A one
dimensional integral over the chi square distribution of that estimate
evaluates the composite deterministically, with no simulation, and the
plot() method draws the population effects the plan rests on. This
opens the new composite power family.ss_power_composite_factorial_ancova() and
ss_power_composite_factorial_anova() carry composite power to any
balanced factorial design: name any set of main effects and interactions
in effects, each with its Cohen's f or partial eta squared, and the
planner returns the per-cell sample size at which all of them are jointly
significant. The effects are noncentral F tests sharing one error
estimate, so the same shared-error integral applies; a single effect
reproduces ss_power_factorial_anova() (or ss_power_factorial_ancova()
with a covariate) exactly, and the two-effect composite matches a direct
simulation. The ANOVA function is the no-covariate case named directly and
admits no covariate. With no covariate the composite is exact to quadrature
precision. The effects can be stated as sizes (Cohen's f or partial eta
squared) or as a full array of population cell means with a common within-cell
standard deviation, from which each effect's f is read off the analysis of
variance decomposition of the means. The plot() method draws the purported
population values: the cell-mean pattern (with error bars of one within-cell
SD) when means were given, or the effect sizes annotated with their marginal
power otherwise, with the composite power in the subtitle either way.ss_power_composite_factorial_ancova_het() carries composite power to a
factorial ANCOVA whose covariate slope differs across the cells. The average
slope (the covariate main effect) and the factor by covariate slope
heterogeneity are then testable effects that can join the composite alongside
the factorial mean effects. The full model fits the means, the covariate, and
every factor by covariate slope, so the residual has N minus twice the cells
degrees of freedom; a "mean", "covariate", or "slope" effect is named in
effects and sized by a Cohen's f or read from population values (cell
means and a covariate outcome correlation per cell, with a common within-cell
SD). The one-factor two-level case reproduces ss_power_composite_ancova_2group() to
machine precision, and a three-effect composite matches a direct simulation of
the heterogeneous-slope model. Kept separate from the common-slope
ss_power_composite_factorial_ancova() for ease of use. The plot() method
draws the population regression line in each cell, so heterogeneous slopes read
as lines of different angle.ss_power_composite_ancova() and ss_power_composite_anova() are the
general entry points to the composite power family, covering a one-way design
with any number of groups as well as any factorial arrangement. The ANCOVA
planner takes a slopes argument, "homogeneous" for one common covariate
slope or "heterogeneous" to let the slope differ across cells and make the
covariate and slope-heterogeneity effects testable; its heterogeneous one-way
case is the a-group generalization of the two-group
ss_power_composite_ancova_2group(), and its two-level case reproduces it
exactly. The ANOVA planner is the no-covariate design named directly, with a
covariate-free interface that points to the ANCOVA version when a covariate is
present. Both accept effect sizes or population values (cell means with a
common within-cell SD) and forward to the factorial planners, so the same
broom tidy()/glance() summaries and plot() figures apply.ss_power_composite_sem() and ss_aipe_composite_sem() carry
composite sample size planning to structural equation models, for both
design goals. The researcher states the population as a fully fixed lavaan
model (or its cov_sem() covariance matrix), labels the parameters of
interest in the free analysis model (structural paths, loadings,
covariances, or :=-defined quantities such as an indirect effect), and
the planner finds the smallest N at which every labeled parameter is
statistically significant in the same study with the desired probability
(ss_power_composite_sem()), or at which every confidence interval is
sufficiently narrow, in expectation or with a stated assurance for the
joint event (ss_aipe_composite_sem()). No closed form covers a set of
dependent SEM estimates, so both planners run an a priori Monte Carlo
simulation (Muthén & Muthén, 2002; Maxwell, Kelley, & Rausch, 2008):
data are drawn from the population covariance matrix, the analysis model
is fit G times per candidate N, and a search seeded by the analytic
Wald approximation brackets and bisects to the smallest integer meeting
the goals; each reported power or proportion travels with its simulation
standard error, and a seed argument makes a plan reproducible. The
power planner joins the dmar_ss_power broom family; marginal
power_<label> and width_within_desired_<label> rows show which
parameter binds the design. Both planners handle a population mean
structure, so latent growth curve targets such as the slope factor's
mean are planned the same way: cov_sem() now also returns mu_theta,
the model implied means of the observed variables, a mu argument
supplies means beside a hand-built Sigma, and the simulated data carry
those means whenever the analysis model has a mean structure. The
"Composite Sample Size Planning for SEM" vignette
(vignette("composite_sem_planning")) walks through the full workflow
for a mediation model with observed variables and for a linear latent
growth curve, planning both composite power and joint accuracy.cohen_h() returns Cohen's h, the effect size for the difference
between two proportions on the arcsine (variance-stabilizing) scale,
h = 2 asin(sqrt(p1)) minus 2 asin(sqrt(p2)). It is signed and is the
proportion analogue of the standardized mean difference (smd()), so a given
h carries the same detectability wherever the proportions sit, which a raw
difference does not.contrast_adjusted() tests an arbitrary contrast among
covariate-adjusted cell means in a factorial ANCOVA, matching
emmeans::contrast().dmar_tbl display layer,
and literature-synonym @concept tags (Cohen's d, Cronbach's alpha,
Cohen's U3, ...) make the marquee estimators searchable by their
eponyms.ss_power_smd() and ss_aipe_smd() echo their user-supplied
planning inputs as rows of the returned table, so the assumptions a
design was planned under travel with the result as tidy data rather
than a printed footer. The supposed effect is labeled supposed_smd
to make clear it is a value the researcher posits (a minimally
important effect or a value believed to be true in the population),
not a sample estimate; power planning also echoes desired_power,
alpha_level, and a numeric tails (2 or 1), and AIPE echoes
width (and assurance when supplied). The value column stays
numeric.reliability_omega_h() now places the observed total-variance
denominator on the same maximum-likelihood (N-divisor) metric as the
fitted loadings, matching MBESS::ci.reliability(type = "hierarchical")
and semTools::compRelSEM(obs.var = TRUE).DMAR is the modern, more general reimplementation of MBESS that reflects how the methods are now used across quantitative psychology, sociology, education, management, marketing, and information systems. MBESS remains stable on CRAN; DMAR is the recommended path forward for new users.
reliability_alpha() and reliability_omega() handle missing
data by full information maximum likelihood, with auxiliary
variables. The measurement family thereby catches up with mlmr(),
whose vignette makes the case against listwise deletion that these
functions previously ignored. Two new arguments, also passed through
by reliability(): missing = c("listwise", "fiml"), defaulting to
listwise deletion so no existing result changes, and aux, a
character vector naming auxiliary columns of data that enter as
saturated correlates (Graham, 2003): correlated freely with each
other and with every item's residual, never loading on the factor.
Supplying aux implies missing = "fiml". The analytic alpha
applies the classical formula to the FIML estimate of the item
covariance matrix; the model based estimators fit with lavaan's
missing = "ml"; robust omega's observed total variance comes from
the FIML covariance matrix. MBESS documents an aux argument on
ci.reliability() that no longer runs (upstream API drift in
semTools), so DMAR implements the saturated correlates model
directly in lavaan syntax rather than depending on it. The returned
table now always carries an N_complete row beside N, making the
cost of listwise deletion visible at a glance, and the treatment is
recorded in missing and aux attributes. Interval methods that
cannot be made correct under FIML (the complete-data closed forms,
ADF, the profile likelihood) are errors, never silent fallbacks;
the bootstrap resamples partially observed rows and refits by FIML.
Under MAR missingness driven by an auxiliary, the test suite shows
FIML-with-auxiliary reducing the bias of listwise deletion roughly
38-fold across 200 replications, and a hand-specified saturated
correlates model in lavaan reproduces the estimates exactly.
reliability_omega() gains a denominator argument selecting
how the total variance in the denominator of coefficient omega is
estimated: "observed" (the default: robust omega, the variance of
the composite estimated directly from the data, the coefficient
Kelley and Pornprasertmanit, 2016, call hierarchical omega and
MBESS::ci.reliability() calls type "hierarchical") or
"model_implied" (the textbook form). The two definitions coincide
in the population when the single-factor model is correctly
specified; only robust omega retains its interpretation as the
proportion of the variance of the composite actually computed when
the model is misspecified, which is why it is the default. The help
page states the properties of each choice and the distinction from
the bifactor omega-hierarchical of Zinbarg, Revelle, Yovel, and Li
(2005). The returned object records the choice in a denominator
attribute.
cfa_k() accepts ordered-categorical items (ordered = TRUE or
a vector of item names; raw data required, each factor all ordered
or all continuous). The model is fit by WLSMV to polychoric
correlations with thresholds in the theta parameterization, which
keeps every defined measurement quantity available, with robust
standard errors. Because a sum score of ordered items lives on the
metric of the observed categories rather than the latent response
metric, each ordered factor's omega is the Green and Yang (2009)
categorical sum score omega computed from the same fit; the
substitution is announced in a message, recorded per factor in an
omega_metric attribute, and the delta method interval columns are
NA for those rows (use reliability_omega_categorical() for a
bootstrap interval). A single-factor ordered cfa_k() reproduces
reliability_omega_categorical() to 1e-6 in the regression tests
despite the different parameterizations, which pins the
parameterization invariance of the computation.
reliability_alpha() gains an estimator argument carrying the
two routes to coefficient alpha, which were briefly two functions.
estimator = "analytic" (the default) is the classical closed-form
equation applied to the observed covariance matrix, the number a hand
calculation produces; estimator = "model_implied" is the reliability
implied by the tau-equivalent (equal loadings) single-factor model fit
by maximum likelihood, which brings a delta method standard error, a
robust (Satorra-Bentler) variant, the profile likelihood interval, and
a testable fit of the equal-loadings claim. Users of MBESS will know
them as ci.reliability(type = "alpha") and type = "alpha-cfa");
the point estimates match those to 1e-6 in the regression tests. The
value "model_implied" is deliberately the word reliability_omega()
already uses for the analogous choice, so one vocabulary covers both
coefficients, and the reliability() wrapper gains a matching
estimator argument beside its denominator.
This replaces the short-lived reliability_alpha_analytic(), whose
name was inherited MBESS jargon that described neither the estimand
nor the contrast. Merging the two also removes a defect their
separation created: reliability_alpha(ci_method = "likelihood") used
to report the classical point estimate beside an interval profiling
the model implied coefficient, two different quantities, so on a
misspecified model the interval could exclude the estimate printed
above it (on the psych bfi A scale, estimate 0.431 against an
interval of [0.591, 0.636]). Each interval method now belongs to the
estimator that can supply it, and asking for one the chosen estimator
cannot give is an error naming the estimator to use instead.
Profile likelihood confidence intervals join the closed forms:
ci_method = "likelihood" in reliability_omega() (model implied
denominator) and reliability_alpha() (profiling the
tau-equivalent, CFA-based alpha). The interval is the set of
population values not rejected by the likelihood ratio test,
computed by refitting the single-factor model under a nonlinear
constraint on the model implied reliability. It respects [0, 1], is
not forced to be symmetric, works from raw data or a covariance
matrix, and matches MBESS::ci.reliability(interval.type = "ll")
to the fourth decimal in the regression tests.
No bootstrap runs unless the user requests one, anywhere in the
reliability family. Robust omega and categorical omega, whose
confidence intervals are bootstrap based, report the point estimate
by default with a message naming the exact call that produces the
recommended interval (percentile or BCa for robust omega; BCa for
categorical omega). The model implied omega keeps its closed-form
robust ML interval as the default, and alpha and KR-20 keep their
closed forms. When a bootstrap is requested, B = 10000
replications is the default.
reliability_omega_h() was removed. Its coefficient is exactly
reliability_omega(denominator = "observed"), and the "h" (for
"hierarchical") described no hierarchy in the single-factor model
the function fits while inviting confusion with the bifactor
omega-hierarchical. The reliability() wrapper drops
type = "omega_h" and instead forwards a denominator argument to
reliability_omega(). The documentation refers to the coefficient
as robust omega, records the original motivation (model misfit as
minor common factors, with the coefficient isolating the general
factor's variance against the observed composite variance) and the
authors' retrospective preference for a name stating the behavior,
states the qualifications the word robust requires (robust to
misspecification of the total variance only; distinct from
outlier-robust estimation and from robust standard errors), and
notes the design principle shared with categorical omega: in both,
the total variance in the denominator is not taken from the fitted
factor model.
reliability_omega_categorical() is the categorical omega
function's full name, and its only name (an earlier
reliability_omega_c() alias was removed). The reliability()
wrapper's canonical type is "omega_categorical" with "omega_c"
accepted as a shorthand. The coefficient attribute is now
"omega_categorical".
The reliability vignette was rewritten around the framing of Kelley and Pornprasertmanit (2016): choosing the coefficient (a claim about the measurement model and the composite being scored) and choosing the interval (an empirical performance question their Monte Carlo studies answered, which is what the family's defaults encode). One running example walks alpha versus omega, the two omega denominators under a minor-factor contamination, categorical omega under same versus differing threshold patterns, and coefficient H as a different composite rather than a different assumption.
cohen_kappa() now accepts a published k x k frequency table
(table =) in place of raw rater vectors, and custom weight
matrices in either scaling: agreement weights (diagonal 1) or
Cohen's (1968) ratio-scaled disagreement weights
(weight_scaling = "disagreement", zero diagonal, invariant to
positive rescaling, converted internally via w = 1 - v / max(v)).
Asymmetric weight matrices are supported for validity designs where
the two directions of a confusion carry different costs. The help
page replicates Cohen's (1968) Table 1 analyses in full: unweighted
kappa .492, weighted kappa .348 under his disagreement weights,
.574 with the 6 and 1 weights interchanged, his asymmetric
computer-diagnosis validity example (.353), and his Formula 10 and
13 standard errors (.0901, .0916, z = 3.80), which the examples
reproduce for the historical record while the function reports the
Fleiss, Cohen, and Everitt (1969) standard error that superseded
them. When both raters are supplied as factors with the same level
set, their level order is now respected instead of alphabetical
sorting, which previously could silently misalign ordinal categories
with linear, quadratic, or custom weights. Every result also carries
a cells attribute holding the per-cell detail in the form of
Cohen's Table 1: observed proportion, chance-expected proportion,
and the weights (both scalings when disagreement weights were
supplied), one row per cell of the confusion matrix.
diagnosis_agreement ships Cohen's (1968) Table 1 as a data
set in its original layout (Judge B in rows, Judge A in columns):
one row per cell with the frequency, Cohen's ratio-scaled
disagreement weight, the observed proportion, and the
chance-expected proportion. The reconstruction is verified against
every quantity computed from the table in the paper.
A weighted kappa vignette works Cohen's (1968) illustration in
full on the diagnosis_agreement data: unweighted and weighted
kappa, both weight scalings, his Formula 10 and 13 standard errors
beside the Fleiss-Cohen-Everitt interval, linear and quadratic
weighting with the weighted-kappa-equals-r identity under equal
marginals, and the asymmetric validity example, including a neutral
working of how the printed weight display's orientation relates to
the published values.
cfa_k() fits a confirmatory factor analysis model with one or
more factors, each specified by naming its indicators. The
measurement structure is specified by describing what is constrained
(equal_loading, equal_intercept, equal_error, each a single
value or per-factor), and the function names the classical structure
the description implies (congeneric, essentially tau-equivalent,
tau-equivalent, essentially parallel, or parallel; Graham, 2006) in
the printed header and the "model" attribute. The table reports
every estimate with a confidence interval, and per factor coefficient
omega, the average variance extracted, and coefficient H as lavaan
defined parameters with delta method standard errors and intervals
(no semTools involvement). output = "measurement" gathers the
measurement properties, the latent correlations, and htmt() per
factor pair; output = "fit" hands back the lavaan object so two
descriptor fits feed lavaan::lavTestLRT() directly. The RMSEA
interval level is reported explicitly as the rmsea_ci_level row.
plot_cfa_k() displays the item-level loadings, error variances,
or intercepts of a cfa_k() fit with confidence intervals, one panel
per factor, with a dashed reference line that shows the equated value
(or, for free estimates, the informal "one common value" anchor), so
the equality questions behind the classical structures can be seen
before they are tested.
bifactor_indices() computes the bifactor dimensionality and
reliability indices (ECV, omega, omega hierarchical and hierarchical
subscale, PUC, and coefficient H) from a fitted bifactor lavaan model,
with a guard that flags improper (Heywood) solutions (Rodriguez, Reise, &
Haviland, 2016).
simple_structure() quantifies Thurstonian simple structure in a
loading matrix: Hofmann item complexity, the hyperplane proportion, and
pure/complex item counts.
ecvi() gives the Browne and Cudeck (1989) expected cross-validation
index for a covariance-structure model, with a confidence interval derived
from the noncentral chi square (ci_nc_chisq()); accepts a
lavaan fit or a published fit table.
common_method_single_factor() and common_method_marker()
implement the single-common-factor (Harman) screen and the Lindell and
Whitney (2001) marker-variable adjustment for common method variance.
The Harman screen is implemented factor analytically: a one-factor
model is fit by maximum likelihood (stats::factanal()) and the
reported proportion of variance is the common factor's, not a
principal component's.
cv_bryant_paulson() gives the simultaneous
critical value (the ANCOVA member of the cv_* family, reducing to
sqrt(2) * cv_tukey_hsd() when there are no covariates), and
ci_c_ancova_bp() places simultaneous (familywise) confidence
intervals on contrasts of adjusted means, the familywise counterpart
of the per-comparison ci_c_ancova(). Both are computed exactly from
the Bryant--Paulson generalized studentized range distribution
(qbryant_paulson() / pbryant_paulson() / dbryant_paulson()),
evaluated as a Beta mixture of ptukey() rather than read from a
table. Three vignettes cover the critical values, an end-to-end
ANCOVA workflow, and a simulation confirming exact familywise error
control. Implements Bryant and Paulson (1976) and Bryant and Bruvold
(1980).regions_of_significance(): the values of the covariate at which
two groups differ significantly when the within-group regression
slopes are not equal (heterogeneity of regression), where the group
difference is a function of the covariate rather than a single
number. The boundaries solve the quadratic that sets the squared
group difference against its sampling variance at the critical value,
the Johnson and Neyman (1936) procedure; with more than two groups
the calculation is carried out for every pair.
plot_regions_of_significance() draws the estimated difference
across the covariate with the confidence band the region is read
from, so the plot is the decision rule.cv_f(), cv_chisq(), and cv_bonferroni_f() complete the cv_*
family's coverage of the Maxwell, Delaney, and Kelley (2027) Appendix.
The family had cv_t() and cv_z() but not the F or chi square
counterparts, and the Bonferroni F table had no function at all.
cv_f() covers Appendix Table A.2, cv_bonferroni_f() Table A.3, and
cv_chisq() Table A.9; the tests assert each against the printed
values.cv_f() and cv_chisq() default to alternative = "greater" rather
than the "not_equal" that cv_t() and cv_z() use. Neither
distribution is symmetric, and both are used one-sided in the upper
tail for the tests they serve: a restricted model fits worse than a
full one, so evidence against a restriction is a large F, never a
small one. Both tails remain available, through
alternative = "not_equal" or through alpha_lower and alpha_upper,
for an interval on a variance or on a ratio of variances. Both accept a
noncentral parameter, as cv_t() does.cv_bonferroni_f() reports the per-comparison rate alpha / C back in
its area_greater column, which is the whole of what the adjustment
does. Its help page separates it from the rank-sum procedure of
dunn_test(): Dunn (1961) is the Bonferroni procedure, Dunn (1964) is
the nonparametric one.ci_games_howell(): simultaneous confidence intervals for all
pairwise comparisons when homogeneity of variance is not assumed. Every
other all-pairs procedure in the package pools the within-group
variances into MS_W, so none of them is robust when that assumption
fails. Games-Howell uses a separate error term and a
Welch-Satterthwaite degrees of freedom for each pair, then takes its
critical value from the studentized range, following Maxwell, Delaney,
and Kelley (2027, Chapter 5, Equations 5.13 and 5.14). It is the
heterogeneity-robust counterpart of ci_tukey_kramer() and handles
unequal n as a matter of course. With two groups it is exactly Welch's
t test, which the tests assert against t.test(var.equal = FALSE).
Implements Games and Howell (1976).dunn_test(): Dunn's rank-sum test of all pairwise differences,
the follow-up to a significant Kruskal-Wallis test. It ranks all
observations together, as the omnibus test does, and uses the variance
of the ranks implied by the Kruskal-Wallis null (with the tie
correction), so it stays coherent with the omnibus result in a way that
running a Mann-Whitney test on each pair does not. The method
argument passes the pairwise p-values to p.adjust(). Implements
Dunn (1964). Note this is the nonparametric Dunn procedure, not the
Bonferroni procedure of Dunn (1961) that Maxwell, Delaney, and Kelley
(2027) call Dunn's procedure; the help page disambiguates the two.randomization_test() and randomization_test_paired():
randomization (permutation) tests for two independent groups and for
paired observations. The two-group test refers the observed statistic
to its distribution over reassignments of the observed scores to the
groups, so the p-value needs no assumption about the population's
shape, and the result travels with the effect sizes the test only
screens for: the mean difference with a randomization-based interval,
the standardized mean difference with a noncentral t interval, the
common language effect size, and Cliff's delta. The paired test
treats the within-pair sign of each difference as the randomization
mechanism, enumerating all 2^n sign patterns exactly for small
samples and sampling them otherwise. Implements the logic of Fisher
(1935) as developed by Edgington and Onghena (2007).
plot_randomization_test() displays the randomization distribution
the p-value is read from, with the observed statistic marked.alpha_level everywhere. Twenty-two
functions took a bare alpha while twenty-eight already took
alpha_level, and the bare name was carrying three incompatible
meanings inside one package: the Type I error rate in the cv_*,
equivalence, TOST, sequential, and Fisher-exact families; coefficient
alpha in var_alpha(); and ggplot2 transparency in
plot_trajectories(). A reader who learned one meaning would misread
the others. The Type I error uses are now alpha_level throughout, so
fifty functions share one unambiguous name, and the argument that
echoes it in a result table is named to match. The other two keep
alpha, where it is unmistakable: var_alpha(alpha = ) is the
coefficient the function is named for and parallels
var_omega(omega = ), and transparency is the universal ggplot2
convention in a plotting call where no error rate appears. Since the
package is unreleased, the old name fails loudly rather than being
aliased. alpha_lower and alpha_upper, which name the two tails of
an asymmetric critical value, are unchanged.cv_smm() and cv_dunnett() are now deterministic. Both formerly
obtained their multivariate-t quantile from mvtnorm::qmvt(), a
Monte Carlo integrator whose sampling error reached a few hundredths
at alpha_level = .01, enough to move the second decimal of a tabled
critical value. They now evaluate the quantile by exact numerical
quadrature and root finding: the studentized maximum modulus
factorizes into a single one-dimensional integral over the shared
scale (the m statistics are independent given it), and the
balanced-design Dunnett statistic, with its constant correlation
1/2, factorizes through a one-factor representation into two nested
one-dimensional integrals. The returned values are reproducible to the
solver tolerance, reproduce the published Dunnett and studentized
maximum modulus tables more closely than the Monte Carlo path did, and
no longer depend on mvtnorm. Both functions now accept df = Inf
(the known-variance normal limit).seed argument of cv_smm() and cv_dunnett() is removed: the
computation is no longer random, so there is nothing to seed. Call
sites that passed seed should drop it.ci_dunnett() adjusted p-values are now exact. They were
computed from mvtnorm::pmvt(), a Monte Carlo integrator (and, when
\pkg{mvtnorm} was absent, from a conservative Sidak-Bonferroni
fallback). They now use the same deterministic one-factor integral as
the critical value, so the reported p_adjusted is reproducible and
no longer depends on \pkg{mvtnorm}. The shared numeric engine lives in
R/dunnett_internals.R.cv_* family reproducing the Appendix critical-value tables
and validates them by simulation.qbryant_paulson() / pbryant_paulson() / cv_bryant_paulson()
are now accurate at small error df. They obtained the studentized-range
part of the distribution from stats::ptukey(), whose algorithm loses
accuracy at small df (at nu = 3, k = 20 by ~3e-4 in probability, enough
to move the critical value by ~0.2, and more at nu = 2). For nu < 7 the
studentized-range distribution is now evaluated directly, without
ptukey, by integrating the probability integral of the range against
the chi squared error density. The range CDF is splined and cached per
group count; the knots are placed at a fixed spacing of about 0.003 in the
range argument, fine enough that the monotone (Fritsch and Carlson)
interpolation error stays below 1e-9 and the returned critical value
matches the direct integral to about seven figures, while the small-df
path still costs a fraction of a second per group count. The functions
reproduce Bryant and Paulson's (1976) Table 1 exactly, to the two decimal
places tabled, over the whole of its range: both tail areas, all three
covariate counts, every tabled group count, and every tabled error degrees
of freedom from nu = 2 to nu = 120, which is 1188 critical values in all.
Two entries, q_.01;2,8,3 = 23.165013 and q_.01;2,20,4 = 19.745008, sit
about 1e-5 above the point where the second decimal turns over, so they
round to 23.17 and 19.75 while the 1976 table rounds them down; both were
confirmed to fourteen significant figures by two independent high-order
quadrature engines that share no code with the package. A large-scale
simulation of the statistic confirms the computed values independently.
Values for nu >= 7 are unchanged.adjusted_means() extracts the full adjusted (least-squares)
means table. From an lm/aov factorial or ANCOVA fit it returns
every cell's covariate-adjusted mean with SE and CI, or marginal
means over chosen factors via by, with weights = "equal" (the
population marginal means of Searle, Speed, and Milliken, 1980) or
"proportional" (observed frequencies). The reference grid follows
the model's own terms, so transformed covariates and interactions
are handled exactly, and nonestimable cells in rank-deficient
designs error plainly by name. Validated against emmeans 2.0.3 at
1e-10 across fourteen fits (unbalanced factorials, multiple and
transformed covariates, character and ordered predictors, weighted
fits, an empty-cell design) for both weightings, with the
comparisons pinned in tests and re-run live in
tools/oracle_checks.R; also pinned to ancova()'s adjusted means
and to contrast_adjusted() at 1e-10.
loa() reports the Carkeet (2015) pair intervals by default.
The exact CIs for the limits of agreement considered as a pair, the
construction Carkeet recommends for most uses, replace the per-limit
intervals as the default; method = "individual" keeps the
one-sided tolerance-factor form. Both constructions are anchored in
tests to the paper's printed coefficients and worked example, and a
misattributed comment (the Bland and Altman 1999 approximate SE,
previously credited to Carkeet) is corrected.
Sixteen packages leave Suggests. DMAR's test oracles (MBESS, emmeans, semTools, metafor, mirt, sirt, multcomp, performance, r2mlm, irr, irrCAC, psych, BayesFactor, gsl) are no longer dependencies: every live comparison was replaced by its pinned value, with provenance comments naming package, version, and date, and the live comparisons themselves moved to tools/oracle_checks.R, which re-runs all of them against the installed oracles at release time. BiasedUrn is replaced by a self-contained Fisher noncentral hypergeometric density (agreement 1.2e-13), kableExtra by knitr-only HTML and LaTeX table construction, and AMCP by shipping the depression_bdi data directly. Suggests drops from thirty packages to fourteen, none of them needing compiled system libraries beyond what the remaining features genuinely use.
The Bayes factor functions take informed priors, stated either
way. The prior on the standardized effect can now be an informed
Cauchy, moving prior_location off zero (Gronau, Ly, & Wagenmakers,
2020), or a normal with prior_mean and prior_sd for beliefs
stated as moments; a Cauchy has no mean and no variance, so moment
beliefs could not previously be expressed at all. The two families
are exclusive, and the help pages explain what the Cauchy scale
fixes (the quartiles: half the prior mass within one scale of the
location) and the exact identity linking the families: a Cauchy is a
normal prior whose variance is itself uncertain, and the test suite
uses that identity to cross-validate the two code paths against each
other with no external oracle. Every result now carries the full
posterior of the effect in a "posterior" attribute (a data frame
of delta and density), so any posterior probability can be
computed, not only the reported ones; the prior is echoed in
prior_location and prior_scale rows with the family recorded as
an attribute. Default calls reproduce the previous JZS results
exactly.
The Bayes factor functions accept summary statistics. The JZS
Bayes factor depends on the data only through the t statistic and
the sample sizes, so bayes_one_sample_t(), bayes_paired_t(), and
bayes_independent_t() now take the summary statistics a paper
reports (mean, sd, and n; mean_diff, sd_diff, and n; or
mean_1, sd_1, n_1, mean_2, sd_2, n_2) as readily as raw
data, following the same exactly-one-path rule as the rest of the
package. The summary form is exact, not an approximation, and the
help pages show how a standardized effect size enters (a d of 0.5
is mean_1 = 0.5, mean_2 = 0 with unit standard deviations). Tests
pin the raw and summary forms to each other at machine precision.
The directional test family speaks snake_case. The canonical
alternative values of ci_dunnett(), power_fisher_exact(),
randomization_test(), randomization_test_paired(),
summary_t_test(), and welch_t() are now "two_sided", "less",
and "greater", with the base-R spelling "two.sided" accepted as
an alias so existing calls keep working; anything stored on a
returned object carries the underscore form. The cv_* critical
value family keeps its deliberately wider synonym vocabulary with
"not_equal" canonical, unchanged.
var_omega_squared() returned a variance that did not shrink with
the sample size. The worker multiplied the wrong combination of
terms, and the consequence was not a small bias: the returned variance
converged to a positive constant instead of to zero, so its error grew
with N. Against a 200,000-replication Monte Carlo of the sampling
distribution it was too large by a factor of 1.4 at the smallest cell
tested and by 115 at the largest, and N times the variance climbed from
0.8 at N = 30 to 51 at N = 5100 where the true value settles near 0.31.
Any standard error or Wald interval built on it was badly inflated, and
worse in larger samples. The variance is now the delta-method transfer
of Fleishman's (1980, Eq. 22) exact variance of the unbiased estimator
of the signal-to-noise ratio, carried to the omega squared scale by his
Eq. 8 with Jacobian (1 - omega^2)^2, which tracks the Monte Carlo to
within a few percent across every cell tested and to 1.00 by N = 300.
Two regression tests now guard it: agreement with a fixed-seed Monte
Carlo, and the consistency check that N times the variance stays
bounded, which the old formula would have failed. The help page also no
longer attributes an omega squared variance to Fleishman, who gives one
only for the signal-to-noise ratio and says explicitly that the
correlation ratio has an interval and a median but not a variance.
Generalized eta squared now classifies interactions the way its
sources do. Olejnik and Algina (2003, Eq. 5) and Bakeman (2005)
define an effect as a measured source of variance when any factor
in its term is measured, so a measured-by-manipulated interaction
belongs in the denominator automatically. eta_squared_generalized()
and ci_eta_squared_generalized() had instead included only the
effects explicitly listed in observed, a rule the help page
misattributed to Bakeman (2005), who recommends the opposite.
Listing a factor in observed now also places every interaction
containing it in the denominator, so observed = "c" reproduces the
papers' worked examples as written; an explicit interaction label is
still honored as given. Values change only for model-interface calls
on designs with measured-by-manipulated interactions, where the
corrected values are smaller. Tests now anchor the Keppel and Kirk
worked values printed in Olejnik and Algina (2003).
signal_to_noise_R2()'s nonlinear estimator now follows Muirhead's
(1985) Equation 10 exactly: the c/Y correction is added to the
untruncated linear estimate before truncation at zero, rather than
after the linear estimate was itself truncated. The two orders agree
whenever the sample R-squared is not tiny; for sample R-squared below
p/(N-1) the reported nonlinear value is now smaller, as the paper
defines it.
In ss_aipe_pcm()'s assurance search, the critical t now tracks
each candidate sample size's degrees of freedom instead of reusing
the value frozen from the expected-width search. The frozen t was
very slightly conservative; the correction moves nothing by as much
as one subject, and every documented example returns the same
sample size as before.
The ci_c_ancova() example carried the wrong sum of squares for its
covariate. The Maxwell, Delaney, and Kelley Chapter 9 depression
example needs the within-groups sum of squares of the pretest, which
the chapter data give as 752.5 (and which the sibling
ci_sc_ancova() example already used); the example instead passed
313.37, the within-groups sum of squares the covariate explains in
the posttest, a different row of the same ANCOVA decomposition. The
interval moves only in the third decimal, but the two pages now agree
with each other and with the chapter data, and a test recomputes
752.5 from the data (via the AMCP package) so the constant cannot
drift again.
ci_sc() now works when the noncentrality parameter is supplied
directly. The ncp branch set the noncentrality but never the
standardized contrast, so the documented call errored on the closing
table construction (three terms against two values). The branch now
inverts the same relationship the means and psi branches use
(the noncentrality is the standardized contrast divided by
sqrt(sum(c_weights^2 / n))), so the three parameterizations of one
effect return identical tables; that equivalence and the ncp path's
coverage are now tested.
ci_c() and ci_sc() now recycle a scalar n by the number of
contrast weights. The scalar was recycled by length(means), which
is zero on the psi-only path (and, for ci_sc(), the ncp-only
path), so those documented calls stopped on the n / c_weights
length check. No result changes on the means path, where the two
lengths agree; a means vector whose length differs from
c_weights is now rejected with a clear error instead of being
silently recycled.
ss_aipe_cv() now accepts the documented mu / sigma
parameterization, deriving C_of_V = sigma / mu; previously neither
formal was read and the call stopped on "argument is of length
zero". Supplying C_of_V together with mu or sigma is rejected
as conflicting. The mu / sigma path returns exactly the sample
size of the equivalent C_of_V call (for mu = 10, sigma = 1,
width = .1, conf_level = .99, both give N = 20). The example
prose on the help page said the population coefficient of variation
was .25 while the calls passed .1 (which plans N = 20, where .25
plans N = 100); the prose now matches the calls, and a mu /
sigma example was added.
cfa_1() with missing = "ml" now estimates the mean structure
with free item intercepts. The bare lavaan::lavaan() interface it
calls turns the mean structure on for FIML but leaves int.ov.free
at FALSE, so every item intercept was fixed to zero and the
loadings absorbed the item means; on any data not centered at zero
the FIML estimates were wrong (the function's own example, simulated
with mean zero, could not show it). The fit now matches
lavaan::cfa(..., missing = "ml") exactly. Listwise fits are
unaffected.
ss_power_split_plot_anova() now uses the correct between-subjects
noncentrality. The between-subjects test used the per-group sample
size n where its noncentrality should use the total sample size
N = n a, so the noncentrality was too small by the factor a (the
number of between-subjects groups), halved for two groups. This
understated between-subjects power and, planning in reverse,
overstated the necessary per-group sample size; the within-subjects
and interaction tests were unaffected. The corrected between-subjects
test now agrees with a split-plot aov() simulation and, for two
groups, reproduces the two-level treatment test of
ss_power_mixed_effects() exactly (the between-subjects
F(1, .) is the two-level treatment t squared), the identity that
ties the two planners together. Both a value anchor and the
cross-planner identity are now tested.
cv_smm() and cv_dunnett() now work at a large error df. Both
integrated the chi squared error variate over (0, Inf), and once df
is large that density is a narrow bump far from the origin, so the
adaptive rule sampled its way past the mass and returned zero: at
df = 200 the integral evaluated to exactly 0 for every candidate
quantile, which left the root finder with no bracket and the functions
stopped with "f() values at end points not of opposite sign." The
limits are now the extreme quantiles of that variate, which puts the
quadrature on the mass at any df. Both functions gained a large-df
test; the values they already returned are unchanged.
cv_f() no longer returns NaN when the numerator df is infinite.
Supplying ncp = 0 explicitly sends qf() and pf() down their
noncentral algorithm, which does not admit an infinite numerator df,
so the infinite-numerator column of Appendix Table A.2 came back
NaN. Both functions now pass ncp only when it is nonzero, so the
central algorithm serves the central case. cv_chisq() takes the same
guard.
ancova() now reports the covariate-adjusted omnibus F for the
treatment effect. It previously entered the treatment before the
covariate and read the sequential (Type I) sum of squares, which is
the unadjusted treatment F; the covariate is now entered first so
the treatment's sequential sum of squares is its adjusted sum of
squares.
ss_aipe_sm() corrects the noncentrality used in the assurance
branch (it was the reciprocal sm / sqrt(n) instead of
sm * sqrt(n)), which previously left the assurance target
unreachable and pinned the search to a bracket endpoint.
ci_c() now honors the documented df_error argument; supplying it
previously left the error degrees of freedom undefined and errored.
tidy() and glance() on ss_power_r() and ss_power_smd() now
return the user-supplied sample size on the realized-power path
(previously NA).
ss_aipe_reliability() now computes confidence intervals when the
default interval = TRUE is used (a string-only equality check had
skipped the interval for the logical default).
ss_aipe_reliability(type = "Factor Analytic") works again. The path
read cfa_1()'s legacy list layout ($factor_loadings,
$parameter_cov), which the reworked cfa_1() no longer returns, so it
errored; coefficient omega and its delta method interval now route through
the maintained reliability internals (.omega_fit_cfa(),
.ci_omega_delta()).
mlmr() and mlmr_mv() are the new lm-like front end to
full information maximum likelihood (FIML) regression. The
univariate mlmr() mirrors the lm() API (formula interface,
coef / vcov / confint / summary / anova / predict /
update S3 methods, profile / Wald / bootstrap CIs); the
multivariate sibling mlmr_mv() takes cbind(y1, y2) ~ ... and
models the joint distribution of correlated outcomes, with the
residual covariance among outcomes estimated as part of the fit.
See vignette("mlmr", package = "DMAR") for when the
FIML route actually buys you something over lm() + listwise
deletion.mlmr() and
mlmr_mv() gained an auxiliary argument that brings variables
related to the missingness or to the incomplete outcome into the
model as saturated correlates (Graham, 2003): each auxiliary is
correlated with the outcome residual, every predictor, and each
other auxiliary, but never enters as a predictor, so the focal
regression coefficients keep their meaning. This is the inclusive
analysis strategy (Collins, Schafer, & Kam, 2001); it leaves
complete-data estimates unchanged and, under MAR, recovers
information that listwise deletion discards. The differential
missingness scenario in vignette("mlmr", package = "DMAR") works
through when it helps and when randomization already protects the
estimand.tidy() and glance() integration via the
generics package: a uniform interface to the broom ecosystem
for mlmr, mlmr_mv, cfa_1, the reliability family, the long-
format CI family, the ANOVA effect size CI family, and the power
planner family. purrr::map_dfr(fits, generics::tidy) now works
across DMAR outputs.convert_R2_*
family cut iterative AIPE / sensitivity planning calls by
roughly 3x; representative ss_aipe_R2() calls go from ~2.7s
to ~0.9s.New data set test_market: a small balanced ANCOVA example (sales
by promotion type, adjusted for a baseline covariate) used to
illustrate the Bryant--Paulson simultaneous intervals in
vignette("bryant_paulson_ancova", package = "DMAR") and
ci_c_ancova_bp().
New data set drinks_trial: nine-month follow-up
drinks-per-week outcomes for the N = 88 homeless
alcohol-dependent participants in Smith, Meyers, and Delaney's
(1998) randomized trial of the Community Reinforcement Approach
at the Salvation Army Adult Rehabilitation Center in
Albuquerque, New Mexico. Two consecutive cohorts: Cohort 1
compared Standard, CRA, and CRA + Disulfiram with cell sizes
17, 15, 19; Cohort 2 dropped the disulfiram cell after Cohort 1
results and compared Standard against CRA with cell sizes 20,
17. The outcome ships in raw form (heavily right-skewed, range
0 to 624.6 drinks per week) and on the log10 scale used in the
published analyses to recover approximate normality. Reproduced
in Maxwell, Delaney, and Kelley (2027, Designing Experiments
and Analyzing Data: A Model Comparison Perspective, 4th ed.,
Routledge), Chapter 3, Section 3.10.4.
New data set bessel_errors: Friedrich Wilhelm Bessel's (1818)
9-bin grouped frequency distribution of the absolute errors of
300 stellar position observations made by British Astronomer
Royal James Bradley at the Greenwich Observatory between 1750
and 1762. Both the observed and expected (normal-model)
frequency columns sum to 300, matching Maxwell, Delaney, and
Kelley (2027, 4th ed.), Table 1.4. Documented as a worked
example for approximating moments from grouped frequency data
(frequency-weighted means and variances using bin midpoints)
and for plotting empirical-versus-theoretical frequency
comparisons. Ships in the original grouped form Bessel
reported; individual error values are not extant.
New data set prime_time_achievement (also accessible via the
short alias Prime_Time): the full Indiana Prime
Time third grade achievement evaluation file (Lapsley, Daytner,
Kelley, and Maxwell, 2002, ERIC ED466679), built from the
original Indiana Department of Education SPSS system file.
10,927 students nested in 586 classrooms in 163 schools in 61
school corporations (district x region combinations) in 9
educational service regions on 113 variables. Includes ISTEP+
NCE composites (the criterion in the published HLM analyses),
Gates-MacGinitie and AANCE test scores, NPA cognitive ability
scores, classroom enrollment, pupil to teacher ratio, Prime Time
aide indicator and status, school and corporation demographics
and finance, and three derived unique cluster identifiers
(corp_id, school_id, class_id) that respect the nesting
irregularities in the source file (one corporation ID spans two
regions). Original Indiana DOE variable names and the original
six-category race coding are preserved verbatim; the SPSS
variable labels are retained as a label attribute on every
column. The SPSS Select Cases artifact FILTER_$ and the 888
"not applicable" codes have been dropped or recoded to NA.
Includes documented examples that show level-1, level-2, and
level-3 lmer fits mapped onto the multilevel framework used
in Lapsley et al. (2002) and in Finch, Bolin, and Kelley
(2019, Multilevel Modeling Using R, 2nd ed., CRC Press,
chapters 3, 4, 6, 9, 10).
New data set holzinger_swineford (also accessible via the
short alias HS_Data): the complete Holzinger and
Swineford (1939) factor analysis data, 301 pupils on 26 ability
tests from the Pasteur (n = 156) and Grant-White (n = 145)
elementary schools in Chicago. Variable names follow the MBESS
convention (e.g., t1_visual_perception through t26_flags) so
scripts written against MBESS::HS port over with a single rename
of the object. The values are the corrected version of the data,
identical to MBESS::HS as of MBESS 4.9.3 and to
psychTools::holzinger.raw. The documentation discusses the bi-
factor study design, the five ability blocks (spatial, verbal,
mental speed, memory, reasoning), the Joreskog (1969) 9 test
subset, and the silent post-4.6.0 MBESS correction that is the
source of values still found in sem::HS.data and
OpenMx::HS.ability.data.
New data set pygmalion: the teacher-expectancy data from
Rosenthal and Jacobson's (1968) Pygmalion in the Classroom,
310 elementary school pupils in grades 1 to 6, of whom 64 were
randomly designated to their teachers as likely intellectual
"bloomers" and 246 served as controls, with pretest and follow-up
IQ. This is the classic benchmark for analysis of covariance with
heterogeneity of regression, and the running example for that
topic in Maxwell, Delaney, and Kelley (Designing Experiments and
Analyzing Data, Chapter 9). The within-group slopes (0.778 for
controls, 0.969 for bloomers), the pooled residual variance
(175.3251), and the covariate variance (348.91) reproduce the
worked example for the variance of the estimated treatment effect
at selected covariate values (cf. MBESS::var.ete). The same
numbers ship with the book's data companion AMCP as
chapter_9_exercise_15; here the experimental condition is a
labeled factor (Control, Bloomer) and the columns use DMAR's
descriptive names, with no measured value altered. Pairs with
ancova() and is documented in
vignette("pygmalion", package = "DMAR").
The largest visible change is a uniform snake_case interface for
both function names and argument names.
Argument names that were dot.case in MBESS are snake_case in
DMAR. The most common renames a user will hit when porting a
script:
| MBESS | DMAR |
|----------------------|-------------------|
| conf.level | conf_level |
| Random.Predictors | random_predictors |
| Specified.N | specified_N |
| alpha.lower | alpha_lower |
| alpha.upper | alpha_upper |
| degrees.of.freedom | Use df or the explicit df_1 / df_2 |
| Group.1, Mean.1 | group_1, mean_1 |
The rule of thumb: change every . between words to _, and
lowercase the non-statistical prefix words. Capitals are preserved
when the capital is statistically meaningful: R2 (squared
multiple correlation), N (sample size), S (a covariance
matrix), Lambda (a factor-loadings matrix), F_value (an
F-statistic).
aipe_smd() is renamed to ss_aipe_smd() to align it with the
rest of the ss_aipe_* family.
Every estimation, inference, and planning function returns a
data.frame with term and value columns. (Plotting functions
return a ggplot object, and a few utilities return their natural
type.)
Return objects in MBESS were sometimes named lists, sometimes
data.frames, and sometimes vectors; the unified tidy return makes the
package compose cleanly with the rest of the modern R ecosystem.
Confidence intervals and sample size annotations travel with every
effect size plot. plot_smd(), plot_ci(), and plot_R2()
default to show_ci = TRUE and show_n = TRUE.
There are no deprecation shims. Old MBESS argument names will throw "unused argument" errors at the call site, which is intentional; the fix is mechanical, and a silent-forwarding shim would hide it.
DMAR adds 40+ functions across families that were absent or under-developed in MBESS. New families include:
Mediation, both halves: mediate() analyzes the simple mediation
model (optional covariates) with percentile bootstrap, BCa, Monte
Carlo, and Sobel intervals for the indirect effect, and
ss_power_indirect_effect() plans the study by joint significance
with exact noncentral t component powers (Sobel power for
comparison), validated against raw-data simulation.
mediation_mbco() extends the inference side to arbitrary
mediation structures through the likelihood ratio model comparison
framework (see its own section above).
Measurement invariance: measurement_invariance() fits the
configural, metric, scalar, and strict multi-group ladder for any
measurement model, with a thresholds rung for ordered indicators,
and reports the full comparison table (fit, the likelihood ratio
test per step, and delta CFI / delta RMSEA), pinned to direct lavaan
fits in the tests. See the entry at the top of this file for the
arguments and the returned "fits" attribute.
Construct validity: htmt() (the Henseler heterotrait-monotrait
ratio with an optional bootstrap upper bound) and
average_variance_extracted() (Fornell-Larcker AVE from a lavaan fit
or standardized loadings). Composite reliability needs no new
function; it is coefficient omega (reliability_omega()).
Clinical and behavioral endpoints: responder_analysis() (per-group
responder proportions with Wilson intervals, the Newcombe risk
difference interval, the number needed to treat, an omnibus chi
square, and a threshold sweep) and ci_proportion(), the package's
Wilson score interval for a single proportion.
A meta-analysis family. meta_es() pools any effect sizes given
sampling variances (REML between-study variance by default, with
Paule-Mandel, DerSimonian-Laird, and fixed effect options, the
Hartung-Knapp adjustment on by default, Q-profile confidence
intervals for tau-squared mapped to I-squared, and a prediction
interval always reported); meta_smd() (exact Hedges correction by
default) and meta_r() (Fisher's Z pooling, with optional per-study
attenuation corrections) are the metric front ends. combine_p()
provides the four classical combined significance tests and
meta_contrast() the Rosenthal and Rubin contrast among effect
sizes; plot_forest() draws the studies, the pool, and the
prediction interval. Validated against metafor, and against
Raudenbush (1984) line by line in the teacher expectancy vignette;
the 19 effect sizes of that synthesis ship as teacher_expectancy.
Factorial ANCOVA power: ss_power_factorial_ancova() extends the
factorial ANOVA planner to baseline covariates (error variance
scaled by 1 - R2, one error df per covariate), with the complete
2 x 4 x 3 worked example, planning, simulation, Type III analysis,
interaction plots, and focused complex comparisons, in the
ancova_2x4x3_power vignette.
A new vignette, "Power and Precision for the One-Way ANOVA: A
Model Comparison Perspective," works the one-way design from the
comparison of a full model and a restricted model: the omnibus
power and sample size through ss_power_one_way_anova(), a planned
contrast through ss_power_contrast(), effect size confidence
intervals through ci_omega_squared(), ci_pvaf(), ci_snr(),
and ci_srsnr(), the model comparison made literal with mlmr(),
the Type S and Type M consequences of the design through
design_consequences(), and accuracy in parameter estimation for a
contrast through ss_aipe_c(), following Maxwell, Delaney, and
Kelley (2027).
A Bayesian t family with the probability statement front and
center: bayes_one_sample_t(), bayes_paired_t(), and
bayes_independent_t() report the JZS posterior of the
standardized effect (median, mean, credible interval, and
P(delta > 0 | data)) with the default Bayes factor as a secondary
row, computed by exact quadrature and validated against the
BayesFactor package.
correction_for_attenuation() (renamed from its working name
correct_attenuation()) documents and demonstrates the latent
variable route it approximates: when item-level data exist, prefer
the two-factor model's latent correlation to the plug-in formula.
Moments of the noncentral distributions: moments_nc_t(),
moments_nc_F(), and moments_nc_chisq() return the mean, variance,
standard deviation, skewness, and excess kurtosis of the noncentral
t, F, and chi square distributions, with NA for moments whose
degrees of freedom conditions fail. The noncentral t mean is the
quantity behind the upward bias of the standardized mean difference.
A random-coefficients polynomial growth simulator:
simulate_longitudinal_polynomial() generates longitudinal data of
any polynomial order (order 0 is a flat line) for one or several
groups, ties the level-one error to a target measurement reliability
(reported per occasion), allows assessment-time jitter around the
nominal schedule, and supports structured level-one error covariance
(AR(1), compound symmetry, Toeplitz, heteroscedastic, or a full
matrix). It is the Monte Carlo companion to ss_power_pcm().
ss_power_pcm() now plans power for any polynomial change
coefficient (intercept, linear, quadratic, cubic, ...) through a
trend argument, carrying the general Raudenbush and Liu (2001)
sampling variance; the linear default reproduces the National Youth
Survey benchmark.
orthogonal_polynomial() returns orthogonal polynomial trend
contrast weights stored levels-by-trends (ready for contrasts() and
lm()) and printed in the Maxwell, Delaney, and Kelley Table A.10
layout with a trailing sum-of-squared-weights column.
unbiased_R2() gives the Olkin and Pratt (1958) exactly unbiased
estimator of the population squared multiple correlation alongside
the Ezekiel (1930) adjusted estimator (the adj.r.squared of
summary.lm).
Design consequences: design_consequences() reports what a chosen
design delivers under both of the package's lenses, the significance
lens (power, the Type S wrong-sign error rate, and the Type M
exaggeration ratio, after Gelman and Carlin, 2014, computed exactly
rather than by their simulation) and the precision lens (the expected,
median, and standard deviation of the realized confidence interval
width, and the probability the realized interval beats a target width,
the closed-form versions of the ss_aipe_*_sensitivity() Monte Carlo
terms). Accepts a standard error directly or derives it from sd and
per-group n. Every ss_power_* and ss_aipe_* help page points to
it.
The Spearman correction for attenuation: correction_for_attenuation()
disattenuates a correlation for measurement error in either or both
variables, with a confidence interval when the sample size is
supplied, using reliabilities from the reliability_* family.
New parameterization conversions: convert_d_r() / convert_r_d()
(standardized mean difference and point-biserial correlation, with an
unequal-group factor) and convert_d_or() / convert_or_d() (the
Hasselblad and Hedges logistic link to the odds ratio).
convert_F_chisq() and convert_chisq_F() move a test statistic
between the F and chi square metrics, returning the converted
statistic itself (a value, not a p-value). Two conversions are
offered through one argument, df_denominator. The default,
df_denominator = Inf, treats the F's error variance as known and
returns the standard scaling, df_numerator * F_value (and back,
chi_square / df); an F is a chi square whose error variance is
estimated rather than known, and infinite denominator degrees of
freedom is the case where it is known. A finite df_denominator
instead returns the chi square value with the same upper-tail
probability (the same p-value) as the F, which is exact at any
denominator degrees of freedom; the two conversions agree as
df_denominator grows, so the scaling default is the large-sample
limit of the same family. The finite case is computed from ordinary
upper-tail p-values (no logarithms): the upper tail is used because
the lower-tail probability rounds to 1 in double precision by about
F = 500 at small denominator degrees of freedom, which would send the
result to infinity, whereas the upper-tail computation stays accurate
past F = 1e20. Each help page states the exact computation. The map
preserves the p-value but does not transport a noncentrality
parameter, so noncentral work belongs in ci_nc_F() and
ci_nc_chisq().
Display helpers extending the package's p-value convention to
objects DMAR does not produce: format_p(), print_anova(), and
print_summary().
DMAR no longer requires the GSL system library: the Gauss hypergeometric function behind the exact squared multiple correlation moments is computed in base R (and remains accurate as the squared multiple correlation approaches 1, where the previous route lost precision). Installation now has no system prerequisites.
Plot functions in the plot_* family: plot_smd(), plot_ci(),
plot_R2(), plot_trajectories(), plot_trajectories_fitted().
All built on ggplot2 and colored by default with a neutral,
colorblind-safe palette; each accepts a palette argument
("okabe_ito" or "tableau") and, where applicable, a colors
override for full manual control.
Plot colors come from base R. The plot_* family colors itself
with base R's Okabe-Ito colorblind-safe palette by default, with base
R's Tableau 10 available through each plot's palette argument. DMAR
defines no palette of its own and adds no color dependency; a user who
wants other colors adds an ordinary ggplot2 scale to the plot.
Ordinal and non-parametric effect sizes: cliff_delta(),
vargha_delaney_A(), probability_of_superiority_paired(),
cles(), proportion_of_superiority() (sometimes called Cohen's
U3).
Agreement and reliability extensions: bland_altman_loa(),
lin_ccc(), cohen_kappa(), fleiss_kappa(), gwet_ac(),
krippendorff_alpha(), icc(), icc_lmer(),
reliability_omega_categorical(), reliability_omega_h(),
reliability_H(), reliability_kr20(), plus the var_*
family of asymptotic variance utilities.
Within-subjects, mixed, and multivariate ANOVA: anova_within(),
anova_within_two_way(), mixed_anova(), manova_split_plot(),
mauchly_test(), epsilon_corrections(), pairwise_within(),
simple_effects_AB().
Sample size planning has been broadened, particularly the
ss_power_* family, which now covers between-subjects,
within-subjects, mixed, and multi-level designs across the
chapters of Maxwell, Delaney, and Kelley (2027).
Parameterization conversions in the convert_* family
(convert_R2_f, convert_f_R2, convert_lambda_R2,
convert_R2_lambda, convert_r_Z, convert_Z_r,
convert_delta_lambda, convert_lambda_delta,
convert_cor_cov). Every conversion is exact-invertible and the
inverse direction is shipped as a sibling function where it makes
sense.
A worked simulation study of the AIPE family reports expected and realized CI widths, realized coverage, and assurance-achievement rates across 10,000 Monte Carlo replications per cell. It establishes how the methods perform rather than how they are used, so it is maintained alongside the package rather than shipped in it.
ss_aipe_reliability() no longer embeds a vendored copy of MBESS
3.2.0's ci.reliability() for the Monte Carlo assurance search. The
interval at each candidate sample size now comes from DMAR's own
reliability internals (the single-factor delta method interval for
the factor analytic type; the van Zyl, Neudecker, and Nel closed
form for the normal theory type). A before-and-after grid covering
every model and interval type combination confirmed the computed
interval widths, and therefore every planned sample size, are
unchanged.
icc_lmer() now reads its grouping factor from lme4's grouping list,
so interaction groupings such as (1 | school:teacher) work; and
mlmr_mv gained the tidy() / glance() methods its documentation
promised (one row per coefficient per outcome, with a response
column).
The test suite runs with warnPartialMatchArgs,
warnPartialMatchDollar, and warnPartialMatchAttr enabled, so any
internal reliance on partial matching fails loudly, and the display
layer (the dmar_tbl print methods, format_p(), print_anova(),
print_summary()) is covered by testthat snapshots.
The seed argument defaults to NULL everywhere it is exposed,
meaning the function uses the caller's current RNG state and does not
seed; a user who wants reproducibility passes an explicit integer, and
a supplied seed restores the caller's RNG state on exit. The value
113 is used only in @examples and tests, never baked into a
function default.
The cluster-randomized design helpers under ss_aipe_crd* were
consolidated into a single set of shared internals
(R/ss_aipe_crd_internals.R); both the difference and
effect size families now call the same .find_*_crd_*()
back end.
ss_aipe_pcm() had a typo in the cubic change coefficient
constant (K_3 was coded as 1/1000800, about one tenth of the
correct 1/100800). Because the constant sits in the denominator of
the slope variance, this inflated the cubic variance, and hence the
resolved sample size, roughly tenfold for trend = "cubic". The
constant now matches the closed form
K_p = (p!)^2 / [(2p)! (2p+1)!] of Raudenbush and Liu (2001,
p. 392), and a regression test recovers it directly from an OLS fit.
The linear (1/12) and quadratic (1/720) constants were already
correct.
ss_aipe_pcm() now honors its documented contract that a user may
supply either error_variance or the converted variance
variance_true_minus_estimated_trend. Previously error_variance
was effectively required, and the consistency check between the two
used round(..., 3), which mishandled the small variances typical
of slope-change designs. The cross-check now uses all.equal(),
and supplying neither argument raises an informative error.
ss_aipe_pcm() now scales the within-subject contribution to the
slope variance by frequency^(2p), matching ss_power_pcm() and
Raudenbush and Liu (2001, p. 392). The change coefficient is a
per-unit-time rate, so its sampling variance is
error_variance * frequency^(2p) / (sum of the squared polynomial weights); the factor was previously omitted. Every Kelley and
Rausch (2011) benchmark uses frequency = 1, where the factor is 1,
so all tabled sample sizes are unchanged; the correction only
affects designs with frequency != 1, where it now agrees with the
direct OLS slope variance on the actual time grid.
ss_aipe_pcm_sensitivity() now simulates the same estimand its
planner targets: the between-group difference in mean slopes (the
group-by-time change parameter), using two independent groups, a
pooled standard error, and 2n - 2 degrees of freedom. The previous
simulator built a one-group mean-slope interval, whose width was
systematically 1/sqrt(2) of the planned target at every frequency.
The realized mean CI width now tracks the planned width. The
schema terms mean_slope / median_slope / sd_slope are renamed
to mean_slope_diff / median_slope_diff / sd_slope_diff to
reflect the difference estimand.
Ken Kelley (Department of Information Technology, Analytics, and Operations; Mendoza College of Business; University of Notre Dame) is the package author and maintainer. Bug reports and feature requests are welcomed by email to [email protected]; please put "DMAR" in the subject line.