Changelog

Source: NEWS.md

DMAR 1.0.0

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

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 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 -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 whose 100(1 - 2 alpha)% Fisher’s 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. conf_limits_ncf() 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_ncf_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 conf_limits_ncf(), 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.
• conf_limits_ncf() and conf_limits_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 (conf_limits_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 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 . 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, 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_nct(), moments_ncf(), 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 conf_limits_ncf() and conf_limits_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.