                 Changes in version 1.0.0 (2026-09-21)                  

First public release of DMAR (pronounced "Dee-Mar," for "Design,
Measurement, and Analysis in R"), a greatly expanded reimagining of the
MBESS package.

Every Example Runs, and Nothing Writes or Seeds Behind the User's Back

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

An Adversarial Quality Control Pass

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

One Prefix for Every Confidence Interval

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

Fitting Change Models, Linear and Nonlinear

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

Nonlinear Growth Curve Simulators

  - Four random-coefficients nonlinear change simulators join
    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.

One Word for the Equivalence Family

  - The equivalence tests are named for the question they answer:
    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.

Sample Size Planning for Correlation Equivalence

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

The Examples Teach From the Package's Own Data

  - The help-page examples now use DMAR's own datasets rather than the
    base-R teaching sets: Dunnett-style comparisons run on
    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.
  - Spell checking is now wired in: DESCRIPTION declares Language: en-US
    and inst/WORDLIST carries the package's technical vocabulary, so
    devtools::spell_check() runs clean.

A Quieter, Faster Test Suite and ggplot2 4.0 Compatibility

  - 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.
  - The test suite runs its files in parallel (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.

One Return Schema Across the Sensitivity Family

  - Every 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.
  - The term names that departed from the schema moved to it, loudly
    (the package is unreleased; no old name survives).
    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.
  - Four members scaled their non-coverage rows by 100 while the rest of
    the family reported proportions: 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().
  - The family contract test caught a real defect while it was being
    written: the 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.
  - Members that reported only width summaries now also report their
    estimand: 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).
  - The power siblings joined the sweep: 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().

The ICC Sample Size Planner Now Plans the Average-of-k Forms

  - 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.
  - The claim in ?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.

The Noncentral F Clamp Warning Speaks the Caller's Language

  - When an observed F falls below the 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.
  - A genuine failure of the inner root finding is no longer a bare
    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.

design_consequences() Now Uses the Exact Noncentral t When df Is Finite

  - With finite 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.

The Robust Standardized Mean Difference Scales Correctly at Every Trim

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

The Robust Standardized Mean Difference Interval Now Uses the Yuen-Welch Degrees of Freedom

  - 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].

A Consistent Confidence Interval for Lin's CCC

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

The Fleiss Kappa Test Uses the Corrected Null Variance

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

One-Sided Dunnett Intervals Now Contain Their Point Estimates

  - 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() Declines to Estimate an Inestimable Epsilon

  - When an effect's numerator degrees of freedom exceed n - 1, the
    covariance matrix of the effect's orthonormal contrasts is singular
    and the Greenhouse-Geisser formula returns an artifact bounded above
    by (n - 1)/df regardless of the population epsilon (even under exact
    sphericity, where the population value is 1); the Huynh-Feldt value
    derived from it is equally meaningless. 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() Now Requires at Least Four Observations

  - 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.
  - The same unguarded 1/(N - 3) sat behind 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.

The AIPE Sensitivity Family Reports Coverage as Proportions

  - Three members of the 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.

The MBCO Likelihood Ratio Is Now Branch-Deterministic

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

tidy() and glance() Speak DMAR's Names, on Every Table

  - The broom verbs now return DMAR's native column names: 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.
  - The wide tables answer the verbs too. 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.

One Bootstrap Vocabulary, Three New Bootstrap Intervals

  - 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.
  - Failed bootstrap replications are dropped, never fatal: 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.
  - No help page wraps its examples in \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.

The CFA Family: One General Function, Two Convenience Wrappers

  - 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 Beyond the Exact-Invariance Ladder

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

Item Response Theory, as the Categorical Factor Model

  - 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

  - 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 inference via model-based constrained optimization

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

API and Documentation Consistency

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

Ten Help Pages Now State Exactly What the Code Computes

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 Second Documentation Pass: Quoted Values Now Match What Runs

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.

  - Quoted numbers now match what the examples print. 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.
  - The 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.
  - Smaller contract corrections. 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 and noninferiority for linear contrasts

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

Mixed-effects R-squared

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

Output, tidiers, and contrasts

  - Publication-ready output helpers: 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()).
  - The broom summary now spans the power-based sample size planners.
    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().
  - More estimators now route through the tidy 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.

Correctness

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

New measurement functions

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

ANCOVA multiple comparisons

  - Bryant--Paulson multiple comparisons for ANCOVA: a new family for
    simultaneous inference on covariate-adjusted means when the
    covariates are random. 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.

The remaining textbook critical-value tables

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

Pairwise comparisons without the usual assumptions

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

API and Documentation Consistency

  - The Type I error rate is 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.

Multiple-comparison critical values computed exactly

  - 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).
  - The 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.
  - A new vignette, "Reproducing the Textbook Critical-Value Tables,"
    walks the 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.

Bug fixes

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

Highlight features

  - 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.
  - Auxiliary variables in the FIML family. 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.
  - broom-style 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.
  - Performance: inner-loop fast paths for the convert_R2_* family cut
    iterative AIPE / sensitivity planning calls by roughly 3x;
    representative ss_aipe_R2() calls go from ~2.7s to ~0.9s.
  - Numerical-correctness tests: a dedicated test file compares marquee
    CI / variance / planner functions to MBESS and to closed-form
    derivations from the foundational papers, so silent numerical
    regressions are caught immediately.

Data sets

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

Migrating from MBESS

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.

What is new in DMAR (beyond MBESS)

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.

Internal changes

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

Authorship

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 kkelley@nd.edu; please put "DMAR" in the subject
line.
