Package {DMAR}

Type:	Package
Title:	Design, Measurement, and Analysis in R (DMAR)
Version:	1.0.0
Date:	2026-09-07
Imports:	grDevices, MASS, generics, parallel, stats, utils, withr
Suggests:	boot, car, ggplot2 (≥ 3.4.0), ggrain, knitr, lavaan (≥ 0.7-2), lme4, lmerTest, mvtnorm, nlme, OpenMx, patchwork, reformulas, rmarkdown, testthat (≥ 3.2.0)
Description:	Methods for design, measurement, and analysis, with the aim of being user friendly yet methodologically sound. 'DMAR' (pronounced "Dee-Mar") implements many advanced and nonstandard methods and makes them available for straightforward use, with interfaces, defaults, and documentation that are consistent across the package and grounded in the methodological literature, in support of sound and reproducible results. The package emphasizes effect size estimation with confidence intervals; sample size planning through accuracy in parameter estimation (AIPE) and power analysis (including composite power for designs whose conclusions require several results to hold at once), with minimum risk, sequential, and equivalence frameworks; reliability, agreement, and measurement more broadly, from coefficient omega with confidence intervals to measurement invariance; factor analysis and structural equation modeling, in which constructs, latent variables measured by multiple indicators, are modeled directly, with confirmatory factor analysis, convergent and discriminant validity, and sample size planning for structural equation models; mediation analysis, from the simple mediation model with bootstrap intervals to likelihood ratio tests of arbitrary indirect effects by model-based constrained optimization (MBCO), with multiple groups and the probing of moderated mediation; equivalence and noninferiority testing; meta-analysis; repeated measures, multivariate, ANOVA, and ANCOVA designs; and inference grounded in model comparison throughout. Measurement is approached from a psychometric perspective, and although many of the methods grew up in human-centered research, they apply broadly across the empirical sciences. Much of what is implemented traces to the author's methodological work, interests, and collaborations. 'DMAR' is a more modern, more general, and greatly expanded reimagining of the 'MBESS' package (Kelley, 2007a, <doi:10.18637/jss.v020.i08>; 2007b, <doi:10.3758/BF03192993>), which has been on CRAN for more than two decades and remains available there in stable form. Most functions accept either raw data or the summary statistics typically reported in published articles, so an analysis can be reproduced from a paper without the original data, which is useful both for extending a published analysis and for meta-analytic work. The estimation, inference, and planning functions return one consistently formatted data frame per function that composes with the broader R ecosystem, and confidence intervals are reported alongside effect sizes throughout, as best practice recommends. Researchers who have data and a question but who are not R experts will find the package approachable, while methodologists gain access to advanced and nonstandard methods, including tables of critical values not available elsewhere.
Depends:	R (≥ 4.0.0)
URL:	https://kenkelley.org, https://yelleknek.github.io/DMAR/, https://github.com/yelleKneK/DMAR
BugReports:	https://github.com/yelleKneK/DMAR/issues
License:	GPL (≥ 3)
Encoding:	UTF-8
Language:	en-US
LazyData:	true
LazyDataCompression:	xz
VignetteBuilder:	knitr
Config/testthat/edition:	3
Config/testthat/parallel:	true
Config/testthat/start-first:	ss_aipe_crd_es, reliability_omega_categorical, cv_bryant_paulson, ci_dunnett, reliability_omega, ss_power_composite_sem, ss_aipe_composite_sem, mediation_mbco, ss_aipe_new_sensitivities
Config/roxygen2/version:	8.1.0
NeedsCompilation:	no
Packaged:	2026-09-07 23:48:51 UTC; kkelley
Author:	Ken Kelley [aut, cre]
Maintainer:	Ken Kelley <kkelley@nd.edu>
Repository:	CRAN
Date/Publication:	2026-09-21 20:50:14 UTC

Design, Measurement, and Analysis in R

Description

A modern R package for design, measurement, and analysis, with special strength in effect sizes, confidence intervals, size planning, reliability and agreement, mediation analysis, equivalence testing, meta-analysis, experimental and quasi-experimental designs, repeated measures, and model comparison-based inference. DMAR (pronounced “Dee-Mar”) is heavily methodological in nature, drawing on the psychometric and statistical traditions, and is aligned with the methodological and applied research program and interests of the author. Many aspects of the package traces to the author's methodological work and collaborations, including sample size planning via accuracy in parameter estimation (AIPE; Kelley & Maxwell, 2003; Kelley & Rausch, 2006; Maxwell, Kelley, & Rausch, 2008), the definition and communication of effect sizes (Kelley & Preacher, 2012; Preacher & Kelley, 2011), and the model comparison perspective of Maxwell, Delaney, and Kelley (2027). It aims to be methodologically sound and particularly well-suited to research in which the independent or dependent variables involve the person, across psychology, sociology, education, behavioral economics, management, marketing, and information systems.

Details

The package makes accessible to researchers a variety of methods that are easy to use, including from sample estimates or results reported in published articles: effect size estimation, confidence intervals for effect sizes, sample size planning, multivariate methods, factor analysis, and certain latent variable models. Particular strengths include sample size planning under several complementary frameworks: accuracy in parameter estimation (AIPE), power analysis, minimum risk, and equivalence. Most exported functions return a tidy data.frame with a term column and a numeric value column (some carry additional typed columns per term), and ggplot2 is used for the plotting functions. A few functions return the shape their task calls for instead: model fits such as mlmr return a richer list-like object with coef / vcov / confint methods, descriptives returns a list of summary tables, and a small number of scalar utilities such as skewness and kurtosis return a bare numeric. The interface is consistent, modern, and opinionated, and is designed for clarity and reproducibility.

DMAR builds heavily on the MBESS package (Kelley, 2007a, 2007b), which has been on CRAN for more than two decades. MBESS was originally framed for the behavioral, educational, and social sciences, but its use has grown well beyond that scope; DMAR is quite general, though especially well-aligned with human-centered research.

Function families. A user-facing tour:

Effect sizes

smd, smd_c, smd_trimmed, eta_squared, eta_squared_partial, eta_squared_generalized, omega_squared, omega_squared_partial, cohen_f, cles, cliff_delta, vargha_delaney_A, proportion_of_superiority, probability_of_superiority_paired, lin_ccc.

Confidence intervals on effect sizes

ci_smd, ci_smd_c, ci_R2, ci_R, ci_rc, ci_src, ci_eta_squared (and partial / generalized variants), ci_omega_squared, ci_pvaf, ci_snr, ci_srsnr, ci_mahalanobis, ci_eigenvalue, ci_cv, ci_sm, ci_reg_coef, ci_r, ci_rmsea.

Maximum likelihood regression

mlmr (univariate full information maximum likelihood (FIML), lm-like), mlmr_mv (multivariate FIML).

ANOVA and ANCOVA

ancova, anova_within_two_way, mixed_anova, manova_split_plot, simple_effects_AB, contrast_test, pairwise_within, mauchly_test, obrien_test.

Reliability and agreement

reliability, reliability_alpha, reliability_omega (with a model implied or observed total-variance denominator, and a reliability_omega_categorical for ordered items), reliability_kr20, reliability_H, cohen_kappa, fleiss_kappa, krippendorff_alpha, gwet_ac, limits_of_agreement.

Mediation

mediate (the simple mediation model with bootstrap, Monte Carlo, and Sobel intervals), mediation_mbco (likelihood ratio tests of arbitrary mediation effects by model-based constrained optimization, with multiple groups and moderated mediation probing), and plot_mediation_mbco (conditional effect curves with confidence bands).

Confirmatory factor and SEM tools

cfa_1, cov_sem, covmat_from_cfa, compare_cov_structures.

Sample size planning (AIPE)

ss_aipe_smd, ss_aipe_R2, ss_aipe_reg_coef, ss_aipe_partial_r, ss_aipe_omega_squared, ss_aipe_icc, ss_aipe_cv, ss_aipe_pcm, ss_aipe_rmsea, the cluster-randomized planners ss_aipe_crd_*, plus their Monte Carlo sensitivity companions ss_aipe_*_sensitivity.

Sample size planning (power)

ss_power_smd, ss_power_R2, ss_power_r, ss_power_reg_coef, ss_power_sem, ss_power_c, ss_power_c_ancova, ss_power_contrast, ss_power_factorial_anova, ss_power_split_plot_anova, ss_power_mixed_effects, ss_power_one_way_anova, ss_power_pcm, ss_power_rm_anova, ss_power_sc.

Critical values and tests

cv_t, cv_z, cv_smm, cv_scheffe, cv_tukey_hsd, cv_dunnett, ci_dunnett, ci_tukey_kramer, ci_scheffe, welch_t, summary_t_test, correlations_test, power_fisher_exact, randomization_test_paired, equivalence_smd, equivalence_r, power_equivalence_md.

Parameterization conversions

convert_R2_f / convert_f_R2, convert_R2_lambda / convert_lambda_R2, convert_delta_lambda / convert_lambda_delta, convert_r_Z / convert_Z_r, convert_cor_cov.

Visualization

plot_smd, plot_ci, plot_R2, plot_trajectories, plot_trajectories_fitted.

Multilevel and clustering

icc, icc_lmer, variance_components_mls, design_effect (Kish design effect), ss_aipe_crd_*.

Data sets

bessel_errors (Bessel's 1818 grouped distribution of Bradley's astronomical observation errors), diagnosis_agreement (Cohen's 1968 weighted kappa illustration), drinks_trial (Smith, Meyers, and Delaney's 1998 Community Reinforcement Approach drinking trial), holzinger_swineford (the 1939 factor analysis study), prime_time_achievement (the Indiana Prime Time third grade achievement evaluation), pygmalion (Rosenthal and Jacobson's 1968 teacher-expectancy data), teacher_expectancy (Raudenbush's 1984 meta-analysis of 18 teacher-expectancy experiments), and test_market (Bryant and Bruvold's 1980 controlled test-market experiment for ANCOVA with a random covariate).

Feedback. Bug reports, feature requests, and suggestions for new methods are welcomed by email to Ken Kelley kkelley@nd.edu (please put “DMAR” in the subject line). See https://kenkelley.org for Ken Kelley's web site, https://kenkelley.org/publications/ for related publications, and https://github.com/yelleKneK/DMAR for the project's GitHub page.

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007a). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K. (2007b). Methods for the behavioral, educational, and social sciences: An R package. Behavior Research Methods, 39(4), 979–984. doi:10.3758/BF03192993

Kelley, K., & Maxwell, S. E. (2003). Sample size for multiple regression: Obtaining regression coefficients that are accurate, not simply significant. Psychological Methods, 8(3), 305–321. doi:10.1037/1082-989X.8.3.305

Kelley, K., & Preacher, K. J. (2012). On effect size. Psychological Methods, 17(2), 137–152. doi:10.1037/a0028086

Kelley, K., & Rausch, J. R. (2006). Sample size planning for the standardized mean difference: Accuracy in parameter estimation via narrow confidence intervals. Psychological Methods, 11(4), 363–385. doi:10.1037/1082-989X.11.4.363

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge.

Maxwell, S. E., Kelley, K., & Rausch, J. R. (2008). Sample size planning for statistical power and accuracy in parameter estimation. Annual Review of Psychology, 59, 537–563. doi:10.1146/annurev.psych.59.103006.093735

Preacher, K. J., & Kelley, K. (2011). Effect size measures for mediation models: Quantitative strategies for communicating indirect effects. Psychological Methods, 16(2), 93–115. doi:10.1037/a0022658

See Also

Useful links:

• https://kenkelley.org

• https://yelleknek.github.io/DMAR/

• https://github.com/yelleKneK/DMAR

• Report bugs at https://github.com/yelleKneK/DMAR/issues

Marginal and Conditional R^2 for a Mixed-Effects Model

Description

Computes the Nakagawa and Schielzeth (2013) marginal and conditional coefficients of determination from a fitted mixed-effects model, returned in tidy long form. The marginal R^2 is the proportion of total variance explained by the fixed effects alone; the conditional R^2 is the proportion explained by the fixed and random effects together. The two quantities are the mixed-effects companions to the intraclass correlation returned by icc_lmer: where the ICC isolates the share of variance attributable to a single grouping factor, the marginal and conditional R^2 summarize how much of the outcome variance the fixed and random parts of the model account for.

Usage

R2_mixed_effects(
  model,
  conf_level = 0.95,
  ci_method = c("none", "boot"),
  B = 1000,
  seed = NULL,
  ...
)

Arguments

Details

Variance decomposition. Writing \sigma^2_f for the variance of the fixed-effect linear predictor, \sigma^2_r for the variance attributable to the random effects, and \sigma^2_\varepsilon for the residual variance,

R^2_{\mathrm{marginal}} \;=\; \frac{\sigma^2_f}{\sigma^2_f + \sigma^2_r + \sigma^2_\varepsilon}, \qquad R^2_{\mathrm{conditional}} \;=\; \frac{\sigma^2_f + \sigma^2_r}{\sigma^2_f + \sigma^2_r + \sigma^2_\varepsilon}.

The fixed-effect variance is \sigma^2_f = \mathrm{var}(\mathbf{X}\boldsymbol{\beta}), the variance of the fitted fixed-effect linear predictor across the observations. The residual variance is \sigma^2_\varepsilon = \mathrm{sigma}(\mathrm{model})^2.

Random-effect variance. For a random-intercept model the random-effect variance is the sum of the variance components read off VarCorr. For a model with random slopes the variance contributed by a random-effects term depends on the values of the associated covariates, so the sum of the diagonal variance components is not correct on its own. This function uses the Johnson (2014) extension: for each random-effects term with design matrix \mathbf{Z} and estimated covariance matrix \boldsymbol{\Sigma}, its contribution is the mean over the observations of the quadratic form \mathbf{z}_i^\top \boldsymbol{\Sigma}\, \mathbf{z}_i, that is, \tfrac{1}{n}\,\mathrm{tr}(\mathbf{Z}\boldsymbol{\Sigma}\mathbf{Z}^\top), and \sigma^2_r is the sum of these contributions across all random-effects terms. For a random-intercept term this reduces to the intercept variance component, so the two paths agree.

Scope. The decomposition here is the one appropriate for a Gaussian (identity-link) linear mixed model, which is what lmer and lme fit. Generalized linear mixed models introduce a distribution-specific variance term and are not handled by this function.

The bootstrap interval. The default ci_method = "none" reports the two point estimates alone, so the bootstrap is what to ask for when the marginal and conditional R^2 are to be reported with an interval and the refits it costs are affordable. With ci_method = "boot" the interval comes from a parametric bootstrap (bootMer): each of the B replicates (1000 by default) simulates a new response vector from the fitted model, refits the model, and recomputes the two R^2 values. The unit of resampling is therefore a whole simulated data set drawn from the estimated model, not a resampled set of cases. Only the percentile interval is offered: the limits are the empirical quantiles of the B bootstrap values (Efron & Tibshirani, 1993); there is no BCa or bootstrap standard error variant. Replicates whose refit fails are dropped, and the interval is computed from the replications that return a value. The default B = 1000 is adequate for the central quantiles a percentile interval uses; raising it tightens the Monte Carlo error of the reported limits. Bootstrap results vary from run to run; supply seed for reproducibility.

Value

A data.frame with rows "R2_marginal" and "R2_conditional" in the value column. When ci_method = "boot", lower- and upper-limit rows for each quantity are appended and the confidence level is carried on the object.

Author(s)

Ken Kelley kkelley@nd.edu

References

Efron, B., & Tibshirani, R. J. (1993). An introduction to the bootstrap. New York, NY: Chapman & Hall/CRC.

Johnson, P. C. D. (2014). Extension of Nakagawa & Schielzeth's R^2_{GLMM} to random slopes models. Methods in Ecology and Evolution, 5(9), 944–946. doi:10.1111/2041-210X.12225

Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R^2 from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133–142. doi:10.1111/j.2041-210x.2012.00261.x

See Also

icc_lmer, ss_power_mixed_effects, lmer, VarCorr

Other agreement and measurement: content_validity_index(), gwet_ac(), icc_lmer(), krippendorff_alpha(), limits_of_agreement(), lin_ccc(), variance_components_mls()

Other mixed models: R2_mixed_effects_decomposition(), icc_lmer(), manova_split_plot(), mixed_anova(), ss_aipe_mixed_effects(), ss_aipe_mixed_effects_sensitivity(), ss_power_mixed_effects(), ss_power_split_plot_anova()

Examples

fit <- lme4::lmer(Reaction ~ Days + (Days | Subject),
                  data = lme4::sleepstudy)

# Marginal R2 is the proportion of variance the fixed effects account for;
# conditional R2 adds what the random effects account for, so the gap
# between the two is what the subject-level terms buy.
R2_mixed_effects(fit)

# A parametric bootstrap percentile interval for both quantities. Each
# replication refits the model, so B = 20 keeps the example quick; a
# reported interval deserves the default B = 1000, and raising B
# further tightens the Monte Carlo error of the limits. The seed makes
# the limits reproducible and leaves the caller's generator state as
# it was.
R2_mixed_effects(fit, ci_method = "boot", B = 20, seed = 113)

R-Squared Measures for Mixed-Effects Models

Description

Computes the Rights and Sterba (2019) integrative framework of R-squared measures for a fitted two-level linear mixed-effects (multilevel) model. The model implied outcome variance is fully decomposed into five sources: variance due to level-1 predictors via fixed slopes (f_1), level-2 predictors via fixed slopes (f_2), predictors via random slope (co)variation (v), cluster-specific outcome means via random intercept variation (m), and level-1 residuals (\sigma^2). Proportions of the total, within-cluster, and between-cluster outcome variance attributable to combinations of these sources give the family of R-squared measures.

Usage

R2_mixed_effects_decomposition(model)

Arguments

Details

The companion R2_mixed_effects returns the Nakagawa and Schielzeth marginal and conditional R-squared; this function contains those two as special cases (total_f and total_fvm) within the fuller source decomposition.

The measure superscripts index the variance sources in the numerator and the subscripts index the outcome variance in the denominator: total_* use the total outcome variance, within_* the within-cluster variance (f_1 + v + \sigma^2), and between_* the between-cluster variance (f_2 + m). total_fvm is the omnibus measure (all explained sources over the total variance) and, for a random-intercept model, coincides with the Nakagawa and Schielzeth conditional R-squared computed by R2_mixed_effects; total_f coincides with their marginal R-squared. See Rights and Sterba (2019, Table 1) for the definitions.

The measures were derived under the assumption that the fitted model uses cluster-mean-centering of the level-1 predictors (with the cluster means entered as level-2 predictors). When that centering is not detected, only the total-variance measures are returned, matching the reference implementation.

Fitting the model requires lme4 (for a merMod fit) or nlme (for an lme fit) to be installed.

Value

A data.frame (dmar_tbl) with columns term and value. When the level-1 predictors are cluster-mean-centered, the full set of 12 measures is returned, named total_f1, total_f2, total_v, total_m, total_f, total_fv, total_fvm, within_f1, within_v, within_fv, between_f2, and between_m; otherwise the five total-variance measures total_f, total_v, total_m, total_fv, and total_fvm are returned (the within/between split requires cluster-mean-centering). The returned object carries the source-by-target variance decomposition in attr(x, "decomposition").

Author(s)

Ken Kelley kkelley@nd.edu

References

Rights, J. D., & Sterba, S. K. (2019). Quantifying explained variance in multilevel models: An integrative framework for defining R-squared measures. Psychological Methods, 24(3), 309–338. doi:10.1037/met0000184

Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R^2 from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133–142. doi:10.1111/j.2041-210x.2012.00261.x

See Also

R2_mixed_effects, icc_lmer

Other mixed models: R2_mixed_effects(), icc_lmer(), manova_split_plot(), mixed_anova(), ss_aipe_mixed_effects(), ss_aipe_mixed_effects_sensitivity(), ss_power_mixed_effects(), ss_power_split_plot_anova()

Examples

fit <- lme4::lmer(Reaction ~ Days + (Days | Subject),
                  data = lme4::sleepstudy)
R2_mixed_effects_decomposition(fit)

Adjusted Cell and Marginal Means From a Fitted Linear Model

Description

Given a fitted lm or aov object with one or more factors among its predictors, adjusted_means() returns the means the model actually compares, sometimes called least-squares means or estimated marginal means. By default the table has one row per cell of the crossed factor design, each cell's mean being the model's predicted response at that combination of factor levels with every covariate held at its sample mean (the adjusted cell means of an ANCOVA; for a model without covariates, the model-based cell means). Naming one or more factors in by instead returns the marginal means of those factors, formed by averaging the cell predictions over the remaining factors with either equal or frequency-proportional weights. Every mean is accompanied by its standard error and a t confidence interval on the model's residual degrees of freedom.

Usage

adjusted_means(
  model,
  by = NULL,
  weights = c("equal", "proportional"),
  conf_level = 0.95
)

Arguments

Details

The reference grid and adjusted cell means. The reference grid is the crossing of the model's factor levels, enumerated in the order the factors appear in the model formula with the first factor varying fastest (the order expand.grid produces). This is the same cell order contrast_adjusted expects, so contrast weights can be read off this table row by row. Every covariate enters the grid at its sample mean, and a transformed covariate is evaluated by applying the transformation to the mean of the raw variable: with log(x) in the formula the grid carries mean(x) and the model matrix applies log(), and a poly(x, 2) basis is evaluated at \bar{x}, matching predict on new data at the covariate mean. Writing L for the matrix whose rows are the design-matrix rows of the grid cells, the cell means are L \hat{\beta}, each standard error is the square root of the corresponding diagonal element of L \, \mathrm{vcov}(\hat{\beta}) \, L', and each interval is the t interval on the model's residual degrees of freedom.

Marginal means and the two weightings. With by, the cell predictions are averaged over the factors not named there, and the averaging happens in the coefficient map itself: the marginal mean's L row is the weighted average of its cells' rows, so the estimate and the standard error both follow from one linear function of the coefficients. weights = "equal" weights every combination of the averaged-over factors equally; this is the population marginal mean of Searle, Speed, and Milliken (1980), the mean for a population in which every cell is equally represented regardless of the sample's cell sizes. weights = "proportional" weights each averaged-over combination by its observed frequency (in a weighted fit, by its total prior weight), so the marginal mean targets a population whose margins are shaped like the sample's. With balanced data the two weightings coincide; with unbalanced data they generally differ, and the choice between them is a substantive question about the population of interest, not a technical one (Maxwell, Delaney, and Kelley, 2027, Chapter 7).

Nonestimable means. When the fitted design is rank deficient (for example an empty factorial cell), the model has no predicted value for the affected cell, and a marginal mean that averages over such a cell does not exist either. adjusted_means() refuses with an error naming the affected rows rather than reporting a value contaminated by lm's arbitrary zero for the aliased coefficient.

Scope. The function covers single-stratum lm and aov fits with a single response. Multi-stratum aovlist fits (within-subjects designs fit with an Error() term) are refused, because a within-subjects marginal mean takes its standard error from the matching error stratum, which this function does not compute. Factors must enter the model as variables in the data, not as conversions inside the formula: y ~ factor(g) + x is refused, so convert g in the data first.

Value

A data.frame (class dmar_tbl) with one row per cell of the reference grid or, with by, one row per combination of the named factors. The leading columns give the factor levels; the numeric columns are estimate (the adjusted mean), se (its standard error), and ci_lower / ci_upper (the t confidence limits). The residual degrees of freedom of the intervals are attached as the df_residual attribute and, when by is supplied, the weighting as the weights attribute. The stored values keep full precision; only the display rounds (see dmar_tbl).

Author(s)

Ken Kelley kkelley@nd.edu

References

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 7 on nonorthogonal factorial designs and Chapter 9 on designs with covariates.)

Searle, S. R., Speed, F. M., & Milliken, G. A. (1980). Population marginal means in the linear model: An alternative to least squares means. The American Statistician, 34(4), 216–221.

See Also

contrast_adjusted for a confidence interval on a single contrast of the adjusted cell means; ancova for the one-way ANCOVA table; ci_dunnett for simultaneous many-to-one comparisons.

Other hypothesis tests: ancova(), anova_within(), ci_dunnett(), ci_scheffe(), ci_tukey_kramer(), compare_cov_structures(), contrast_test(), correlations_test(), equivalence_r(), equivalence_smd(), factorial_anova(), manova_split_plot(), mauchly_test(), mixed_anova(), obrien_test(), pairwise_within(), randomization_test(), randomization_test_paired(), regions_of_significance(), simple_effects_AB(), summary_t_test(), welch_t()

Examples

# 1. Cell means of a 2 x 3 factorial (no covariate): end-of-study IQ in
#    the pygmalion expectancy experiment, grades 1 through 3. The grade
#    factor is created in the data, not inside the formula.
pyg <- subset(pygmalion, grade <= 3)
pyg$grade <- factor(pyg$grade)
fit <- lm(iq_8 ~ treatment * grade, data = pyg)
adjusted_means(fit)

# 2. Marginal means of grade, averaging the cell means over treatment.
adjusted_means(fit, by = "grade")

# 3. An ANCOVA: each adjusted mean holds the covariate, here the pretest
#    depression score, at its sample mean.
fit_ancova <- lm(bdi_post ~ condition + bdi_pre, data = depression_bdi)
adjusted_means(fit_ancova)

# 4. With unbalanced cells (few bloomers in every grade) the two
#    weightings answer different questions.
adjusted_means(fit, by = "treatment")
adjusted_means(fit, by = "treatment", weights = "proportional")

Analysis of Change: Fit Change Models to One or Many Trajectories

Description

Fits a change model to longitudinal data: any of the package's four nonlinear change models (negative exponential, logistic, Gompertz, Richards; the parameterizations of Kelley, 2005, 2008, exactly as generated by the simulate_longitudinal_*() simulators) or a polynomial of any order (the linear change model of simulate_longitudinal_polynomial). Two estimation methods are offered. The default, method = "two_stage", fits each unit's curve separately, using only that unit's data, and then summarizes the unit-level parameters: their mean, their standard deviation and variance across units (the individual differences), and the standard error of the mean. method = "mixed" fits the proper random-coefficients mixed-effects model simultaneously, with every parameter carrying a random effect, so the reported between-unit spread is a variance component purged of estimation noise. With a single trajectory (N = 1, or id = NULL) the two-stage method reduces to one least squares fit of that unit's change, reported with its standard errors.

Usage

analysis_of_change(
  data,
  id,
  time,
  outcome,
  model = c("negative_exponential", "logistic", "gompertz", "richards", "polynomial"),
  method = c("two_stage", "mixed"),
  order = 1L,
  start = NULL,
  maxiter = 500L
)

Arguments

Details

Choosing between the methods. The two-stage (curve-by-curve) route is the transparent classic: every unit's curve is inspectable, no unit's data influence another's fit, and it is the only method available for a single trajectory. Its known cost is that sd_units reflects the spread of estimates, which adds each fit's estimation noise to the true individual differences; with short or noisy trajectories it therefore overstates the population standard deviation. The mixed method estimates that between-unit variation as a variance component, separating it from level-one error, and borrows strength across units, at the price of a harder estimation problem (and, for the nonlinear curves, occasional convergence trouble; the fit is started at the two-stage estimates, and a failure suggests simplifying the model or falling back to two-stage). Under method = "mixed" the random effects are uncorrelated across parameters, matching the simulators' default.

Relation to existing tools. The two-stage method is the idea behind nlme::lmList() and nlme::nlsList(), and the mixed method wraps lme4::lmer() and nlme::nlme(); base R also ships self-starting curves (SSasymp, SSfpl, SSgompertz) in other parameterizations. What this function adds is the package's landmark parameterizations (the intercept-shifting \zeta floor, the Richards \delta; SSgompertz's a \exp(-b_2 b_3^x) answers no substantive question directly), the exact match to the simulate_longitudinal_*() simulators so design studies close the loop, one interface across linear and nonlinear change, and the package's tidy summary with failed fits dropped under a single counted warning.

Units whose two-stage fit does not converge are dropped with a single warning reporting how many, and the effective count is the "n_used" attribute. If no unit's fit converges the function stops. The polynomial fit is closed form and does not fail on any trajectory with at least P + 1 occasions.

Starting values. Unless start is supplied, each unit's nonlinear starting values are derived from that unit's data: the floor and the span from early and late observations, the inflection time from where the trajectory crosses the middle of its range, the curvature from the time the trajectory needs to travel the central half of its range, and, for the Richards model, a logistic start (delta = 1). Cleanly measured trajectories rarely need more; hard cases (very short series, strong decreasing curves, near-flat change) may need an explicit start.

Value

A data.frame with one row per model parameter and columns

term

The parameter name.

estimate

The mean of the unit-level estimates (two-stage; with one trajectory, that unit's estimate) or the fixed effect (mixed).

se

The standard error of estimate: the between-unit standard deviation over the square root of the number of fitted units (two-stage), the asymptotic standard error of the single fit (one trajectory), or the fixed-effect standard error (mixed).

sd_units

The between-unit standard deviation of the parameter: the spread of the unit-level estimates (two-stage; NA for a single trajectory) or the random-effect standard deviation (mixed).

var_units

The corresponding variance.

Attributes: "model", "method", "n_units" (units supplied), "n_used" (units whose fit converged; equal to "n_units" under method = "mixed"), "sigma" (the level-one residual standard deviation), and "per_unit_estimates" (a units-by-parameters matrix: the separate unit-level estimates under two-stage, or the unit-level predictions coef() under mixed, which are shrunken toward the fixed effects).

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2005). Estimating nonlinear change models in heterogeneous populations when class membership is unknown: Defining and developing the latent classification differential change model (Doctoral dissertation). University of Notre Dame.

Kelley, K. (2008). Nonlinear change models in populations with unobserved heterogeneity. Methodology, 4(3), 97–112.

Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-PLUS. Springer.

Richards, F. J. (1959). A flexible growth function for empirical use. Journal of Experimental Botany, 10(2), 290–301.

See Also

The generators simulate_longitudinal_negative_exponential, simulate_longitudinal_logistic, simulate_longitudinal_gompertz, simulate_longitudinal_richards, simulate_longitudinal_polynomial; plot_trajectories for plotting the data being fit; nlme::lmList(), nlme::nlsList(), lme4::lmer(), and nlme::nlme() for the engines and their generalizations.

Examples

# Simulate a Gompertz population with individual differences, then
# recover both the mean curve and the spread of its parameters.
set.seed(113)
d <- simulate_longitudinal_gompertz(
  n = 40, target_times = 0:10,
  fixed_parameters = c(alpha = 75, beta = 3, gamma = 0.55, zeta = 10),
  random_variances = c(alpha = 25, beta = 0.4, gamma = 0.005, zeta = 4),
  error_variance = 4
)
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "gompertz")

# The same data as a straight-line (order 1) polynomial: the linear
# model has nothing to say about floors, ceilings, or timing.
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "polynomial", order = 1)

# A single unit's trajectory: N = 1 is one fit, with asymptotic
# standard errors in place of between-unit spread.
one <- d[d$id == levels(d$id)[1], ]
analysis_of_change(one, id = NULL, time = "time", outcome = "y",
                   model = "gompertz")

# The proper mixed-effects model, fit simultaneously across units:
# nlme::nlme() for a nonlinear curve, lme4::lmer() for the polynomial
# (the latter requires lme4 to be installed).
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "gompertz", method = "mixed")
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "polynomial", order = 1,
                   method = "mixed")

Analysis of Covariance (ANCOVA)

Description

Fits a one-way analysis of covariance so the covariate-adjusted group comparison is available from a single call, without assembling the adjusted means, the omnibus test, and the effect sizes by hand. It returns the adjusted (covariate-corrected) group means with standard errors, the covariate-adjusted omnibus F for the group effect (Type III sums of squares), partial \eta^2 and partial \omega^2 with noncentral F confidence intervals, and a homogeneity-of-regression check, returned in one data.frame.

Usage

ancova(data, outcome, treatment, covariates, conf_level = 0.95)

Arguments

Details

Covariate-adjusted test (Type III sums of squares). The omnibus F tests the group effect after adjusting for the covariate(s), that is, the Type III sum of squares for the grouping factor. For a one-way ANCOVA (one grouping factor, with the covariate slopes held constant) the Type II and Type III sums of squares for the group effect coincide, and both equal the sequential sum of squares obtained with the covariate(s) entered first and the grouping factor last, which is how it is computed here; the value matches car::Anova(fit, type = 3). The choice of sum-of-squares type changes the result only in designs with more than one factor or with interactions among factors (Maxwell, Delaney, and Kelley, 2027, Chapter 7); for those, the two-way and mixed analyses report their sum-of-squares type and allow Type I, II, or III.

Adjusted means. The adjusted mean for treatment level j is the model-predicted response at X = \bar X (the covariate grand mean):

\hat \mu_j^{\mathrm{adj}} \;=\; \hat\mu_j - \sum_k \hat\beta_k (\bar X_{kj} - \bar X_k),

where \hat\beta_k is the within-cell slope on covariate k and \bar X_{kj}, \bar X_k are the per-cell and grand means of covariate k.

Homogeneity of regression. The model fit here holds the within-group covariate slopes \beta_k constant across groups. This is a property of the particular model being fit, not an assumption of analysis of covariance in general: it is a testable claim. Adding all group-by-covariate interactions gives an expanded model, and a model comparison F-test of the additive model against the expanded one (stats::anova) assesses whether the slopes differ across groups (Maxwell, Delaney, and Kelley, 2027, Chapter 9). A large F indicates the slopes are not constant, in which case the single adjusted comparison is not the whole story and the interaction model should be entertained directly. The check is reported in the F_homogeneity_of_regression rows.

Effect size CIs. Partial \eta^2 and partial \omega^2 use the noncentral F framework (ci_eta_squared_partial, ci_omega_squared).

Value

A data.frame (class dmar_tbl) with rows for the omnibus test (F_value, df_1, df_2, p_value), the sum-of-squares type used (sum_of_squares_type; 3 for Type III), the point estimates and confidence intervals of partial \eta^2 and partial \omega^2, the adjusted group means and their standard errors (one row per level of treatment), and the homogeneity-of-regression F-test. The result carries the dmar_tbl class, so it prints to 3 significant figures with whole numbers (such as the degrees of freedom) shown without a decimal part and p-values to 4 decimal places (a p-value below 0.0001 prints as “< 0.0001”); the stored values keep full precision. Control the display with print(x, digits = ) or globally with options(dmar.digits = ) (see dmar_tbl).

Author(s)

Ken Kelley kkelley@nd.edu

References

Huitema, B. E. (2011). The analysis of covariance and alternatives (2nd ed.). Wiley.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 9 on analysis of covariance and Chapter 7 on higher-order designs.)

See Also

ci_c_ancova, ci_sc_ancova, ss_aipe_sc_ancova, omega_squared_partial

Other hypothesis tests: adjusted_means(), anova_within(), ci_dunnett(), ci_scheffe(), ci_tukey_kramer(), compare_cov_structures(), contrast_test(), correlations_test(), equivalence_r(), equivalence_smd(), factorial_anova(), manova_split_plot(), mauchly_test(), mixed_anova(), obrien_test(), pairwise_within(), randomization_test(), randomization_test_paired(), regions_of_significance(), simple_effects_AB(), summary_t_test(), welch_t()

Examples

# 1. Compare the two groups in the Pygmalion data on eighth-grade IQ,
#    adjusting for pre-test IQ.
ancova(outcome = "iq_8", treatment = "treatment",
       covariates = "iq_pre", data = pygmalion)

# 2. The same comparison adjusting for two covariates (pre-test IQ and
#    grade).
ancova(outcome = "iq_8", treatment = "treatment",
       covariates = c("iq_pre", "grade"), data = pygmalion)

Likelihood Ratio Test for Nested Mlmr Fits

Description

Compares two or more nested mlmr fits with the likelihood ratio test, delegating the chi square computation to lavTestLRT. The models must be nested (every parameter in the more restricted model is also in the more general one) and must be fit to the same data with the same missing data handling.

Usage

## S3 method for class 'mlmr'
anova(object, ...)

Arguments

Details

"The same data" means the same observations, with the variables the fits share holding the same values. It does not mean the same number of complete cases. Under missing = "fiml" a legitimate nested comparison routinely has different complete-case counts: adding a predictor that is itself incompletely observed lowers the number of rows complete on every modeled variable, yet the observations, and the information each of them contributes to the likelihood, are unchanged. Comparing y ~ x1 with y ~ x1 + x2 when x2 has missing values is exactly the comparison full information maximum likelihood exists to support, and it is accepted. Fits run on genuinely different data are refused.

Each smaller model is refit as a constrained version of the largest one on the largest fit's data, with the slopes of its absent predictors fixed at zero. That keeps the comparison on a single joint observed variable set, which is what makes the chi square difference interpretable when the predictors are modeled as random (fixed_x = FALSE, the mlmr default).

Value

A data.frame of class anova reporting the degrees of freedom, AIC, BIC, log-likelihood, chi square test statistic, and p-value for each consecutive pairwise comparison.

Author(s)

Ken Kelley kkelley@nd.edu

Examples

# Both fits ask for the Wald interval and skip the effect size block,
# since the likelihood ratio test needs neither and each costs refits.
fit1 <- mlmr(t6_paragraph_comprehension ~ t5_general_information,
             data = holzinger_swineford, ci_method = "wald",
             effect_sizes = FALSE)
fit2 <- mlmr(t6_paragraph_comprehension ~ t5_general_information +
               t9_word_meaning,
             data = holzinger_swineford, ci_method = "wald",
             effect_sizes = FALSE)
anova(fit1, fit2)

Compare Nested Multivariate FIML Regression Fits

Description

Compares two or more nested mlmr_mv fits with the likelihood ratio test, delegating the chi square computation to lavTestLRT. The models must be nested (every predictor in the more restricted model is also in the more general one) and must be fit to the same outcomes and the same data with the same missing data handling. This is the multivariate counterpart of anova.mlmr: because every outcome is regressed on the shared predictor set, dropping a predictor drops its slope on every outcome, so the nesting constraint sets that slope to zero across all outcomes at once.

Usage

## S3 method for class 'mlmr_mv'
anova(object, ...)

Arguments

Details

"The same data" means the same observations, with the variables the fits share holding the same values. It does not mean the same number of complete cases. Under missing = "fiml" a legitimate nested comparison routinely has different complete-case counts: adding a predictor that is itself incompletely observed lowers the number of rows complete on every modeled variable, yet the observations, and the information each of them contributes to the likelihood, are unchanged. That comparison is accepted. Fits run on genuinely different data are refused.

Value

A data.frame of class anova reporting the degrees of freedom, AIC, BIC, log-likelihood, chi square test statistic, and p-value for each consecutive pairwise comparison.

Author(s)

Ken Kelley kkelley@nd.edu

Examples

# Both fits ask for the Wald interval and skip the effect size block,
# since the likelihood ratio test needs neither and each costs refits.
fit1 <- mlmr_mv(cbind(t6_paragraph_comprehension, t9_word_meaning) ~
                  t5_general_information,
                data = holzinger_swineford,
                ci_method = "wald", effect_sizes = FALSE)
fit2 <- mlmr_mv(cbind(t6_paragraph_comprehension, t9_word_meaning) ~
                  t5_general_information + t7_sentence,
                data = holzinger_swineford,
                ci_method = "wald", effect_sizes = FALSE)
anova(fit1, fit2)

One Way Within-Subjects ANOVA With Sphericity Diagnostics and Corrections

Description

Performs the univariate one-way within-subjects F test together with Mauchly's test of sphericity and the three standard \varepsilon-corrected p-values (Greenhouse-Geisser, Huynh-Feldt, and lower-bound). Returns everything in a single tidy data.frame so the user can decide which adjustment to report.

Usage

anova_within(x, id = NULL, time = NULL, outcome = NULL)

Arguments

Details

The unadjusted within-subjects F statistic is the same regardless of sphericity; corrections shrink the numerator and denominator degrees of freedom by a factor of \hat\varepsilon \in [1/(k - 1),\, 1], and the p-value is recomputed against the adjusted reference F distribution. When Mauchly's test rejects, prefer the Huynh-Feldt-corrected p-value (less conservative than Greenhouse-Geisser).

For multi-factor within-subjects designs or mixed designs, fit the model with stats::aov(... + Error(id/within)) or with lme4::lmer() directly.

Value

A data.frame with one row per reported F test: adjustment ("none", "Greenhouse-Geisser", "Huynh-Feldt", "lower_bound"), F_value, df_1, df_2, p_value, and epsilon (the correction factor used; NA for the unadjusted row). attr(<output>, "mauchly") contains the row from mauchly_test, and the partial \eta^2 is attached as attr(<output>, "partial_eta_squared").

Author(s)

Ken Kelley kkelley@nd.edu

References

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 11.)

See Also

mauchly_test, epsilon_corrections, aov

Other within-subjects analysis: anova_within_two_way(), epsilon_corrections(), mauchly_test(), pairwise_within(), plot_trajectories_fitted()

Other hypothesis tests: adjusted_means(), ancova(), ci_dunnett(), ci_scheffe(), ci_tukey_kramer(), compare_cov_structures(), contrast_test(), correlations_test(), equivalence_r(), equivalence_smd(), factorial_anova(), manova_split_plot(), mauchly_test(), mixed_anova(), obrien_test(), pairwise_within(), randomization_test(), randomization_test_paired(), regions_of_significance(), simple_effects_AB(), summary_t_test(), welch_t()

Examples

# Simulated within-subjects data with no real effect.
set.seed(113)
Y <- matrix(rnorm(20 * 4), nrow = 20)
anova_within(Y)

# Built-in within-subjects example: nlme::Orthodont (distance ~ age).
res <- anova_within(nlme::Orthodont,
                    id = "Subject", time = "age", outcome = "distance")
res
attr(res, "mauchly")
attr(res, "partial_eta_squared")

Two-Factor Within-Subjects ANOVA With Sphericity Adjustments

Description

Computes the full two-factor within-subjects ANOVA table, main effects A, B, and the A \times B interaction with each effect tested against its own residual stratum and with Greenhouse-Geisser, Huynh-Feldt, and lower-bound sphericity adjustments applied per effect. Returns a tidy long-form data.frame that composes with the rest of DMAR.

Usage

anova_within_two_way(data, outcome, factor_A, factor_B, subject)

Arguments

Details

Design. The two within-subjects factors A (with a levels) and B (with b levels) are fully crossed; every subject contributes a \cdot b observations. Each of the three fixed effects is tested against its own subject-by-effect residual stratum:

• F_A = MS_A / MS_{A:S}

• F_B = MS_B / MS_{B:S}

• F_{A:B} = MS_{A:B} / MS_{A:B:S}

Sphericity. Each effect's univariate F-ratio assumes sphericity of its corresponding subject-by-effect residual covariance matrix. Three adjustments are reported per effect: Greenhouse-Geisser (Greenhouse & Geisser, 1959), Huynh-Feldt (Huynh & Feldt, 1976), and the lower bound \epsilon = 1 / df, where df is the effect's numerator degrees of freedom; this is the smallest value \epsilon can attain, reached under maximal departure from sphericity.

Subjects needed to estimate epsilon. The Greenhouse-Geisser epsilon for an effect with q numerator degrees of freedom is estimated from the sample covariance matrix of q orthonormal contrasts among the effect's cell means, a different matrix for each effect (Maxwell, Delaney, & Kelley, 2027, Chapters 11 and 12). That matrix has rank at most n - 1, so when n - 1 < q it is necessarily singular; the same rank deficiency makes the multivariate approach to a within-subjects design mathematically impossible when n < a (Maxwell, Delaney, & Kelley, 2027, Chapter 13). The Greenhouse-Geisser formula still returns a number in that case, but the number is an artifact of the rank deficiency rather than an estimate: it cannot exceed (n - 1)/q no matter what the population epsilon is, even under exact sphericity, where the population value is 1. Rather than report a value the design cannot support, the function reports NA for the Greenhouse-Geisser and Huynh-Feldt rows of any effect with n - 1 < q and issues a single warning naming the condition (car::Anova likewise declines to report the corrections for an effect whose error matrix is singular). The unadjusted row and the lower-bound row remain: the lower bound 1/q is Geisser and Greenhouse's a priori bound on epsilon, valid no matter how badly sphericity is violated, and it requires no estimate of the covariance matrix (Maxwell, Delaney, & Kelley, 2027, Chapter 11).

Per-effect partial \eta^2. Computed as SS_\mathrm{effect} / (SS_\mathrm{effect} + SS_\mathrm{effect,\, error}) using the appropriate subject-by-effect residual sum of squares.

Balanced data assumed. The implementation assumes a fully balanced design (every subject observed once in every cell). When the design is unbalanced, the function errors and recommends a mixed-effects fit via lmer.

Sums of squares are unambiguous here. Because the design is balanced, the within-subjects factors are orthogonal and the Type I, Type II, and Type III sums of squares for each effect coincide. A Type toggle is therefore not meaningful, and the reported decomposition is unambiguous: the sum of squares attributed to each effect does not depend on the order in which terms enter the model (Maxwell, Delaney, & Kelley, 2027, Chapter 12).

Value

A data.frame with rows for each of the three effects (A, B, A:B) crossed with each sphericity adjustment (none, Greenhouse-Geisser, Huynh-Feldt, lower_bound). Columns: effect, adjustment, F_value, df_1, df_2, p_value, epsilon, partial_eta_squared. When an effect has too few subjects for its epsilon to be estimable (n - 1 smaller than the effect's numerator degrees of freedom; see Details), the Greenhouse-Geisser and Huynh-Feldt rows for that effect carry NA in epsilon, df_1, df_2, and p_value, and a single warning names the condition; the unadjusted and lower-bound rows are unaffected.

Author(s)

Ken Kelley kkelley@nd.edu

References

Greenhouse, S. W., & Geisser, S. (1959). On methods in the analysis of profile data. Psychometrika, 24(2), 95–112.

Huynh, H., & Feldt, L. S. (1976). Estimation of the Box correction for degrees of freedom from sample data in randomized block and split-plot designs. Journal of Educational Statistics, 1(1), 69–82.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapters 11–13.)

See Also

anova_within, mauchly_test, epsilon_corrections

Other within-subjects analysis: anova_within(), epsilon_corrections(), mauchly_test(), pairwise_within(), plot_trajectories_fitted()

Examples

# 1. Balanced 3x2 within-subjects design simulated for illustration.
set.seed(113)
n_sub <- 12
grid  <- expand.grid(subject = factor(1:n_sub),
                     A = factor(c("a1", "a2", "a3")),
                     B = factor(c("b1", "b2")))
grid$y <- with(grid,
  3 * (A == "a2") + 1 * (A == "a3") +
  2 * (B == "b2") + 1 * ((A == "a3") & (B == "b2")) +
  rnorm(nrow(grid), 0, 1) +
  rep(rnorm(n_sub, 0, 1.5), times = 6))
anova_within_two_way(grid, outcome = "y", factor_A = "A",
                     factor_B = "B", subject = "subject")

Average Variance Extracted (AVE)

Description

The average variance extracted (AVE) is the mean proportion of indicator variance a factor accounts for, a standard convergent validity summary for a reflective measurement block in confirmatory factor analysis and structural equation modeling. With standardized loadings \ell_j,

\mathrm{AVE} = \frac{1}{J} \sum_j \ell_j^2.

The quantity itself is elementary. In a model with cross-loadings, an item contributes its loading to the AVE of every factor it loads on; the same shared variance then counts toward each factor's summary, so compare AVE values across factors of such a model with that overlap in mind. Fornell and Larcker (1981) are credited for establishing it as a validity criterion: a construct shows convergent validity when its AVE reaches the conventional 0.50 (the construct explains at least half its indicators' variance), and the Fornell-Larcker discriminant criterion compares each construct's AVE with its squared correlations with the other constructs. AVE is closely related to composite reliability (omega); the modern complement on the discriminant side is htmt.

Usage

average_variance_extracted(
  fit = NULL,
  loadings = NULL,
  conf_level = 0.95,
  ci_method = c("none", "percentile"),
  B = 1000L,
  seed = NULL
)

Arguments

Details

The percentile bootstrap interval resamples the cases behind fit and refits the model once per replication, so its cost is B model fits. That refitting is why the examples below stop at the point estimates: even the smallest permitted B = 100 runs for several seconds on the two-factor model there. To obtain the interval, pass ci_method = "percentile" together with a seed, as in average_variance_extracted(fit, ci_method = "percentile", seed = 113), and keep the default B = 1000 or more for a reported analysis.

Value

A data.frame (class dmar_tbl) with one row per factor: factor (label), ave, and ci_lower / ci_upper (the percentile bootstrap limits; NA when ci_method = "none").

Author(s)

Ken Kelley kkelley@nd.edu

References

Efron, B., & Tibshirani, R. J. (1993). An introduction to the bootstrap. New York, NY: Chapman & Hall/CRC.

Fornell, C., & Larcker, D. F. (1981). Evaluating structural equation models with unobservable variables and measurement error. Journal of Marketing Research, 18(1), 39–50.

See Also

htmt for discriminant validity; reliability_omega for the composite reliability of the same block (coefficient omega is what the composite-reliability literature computes); cfa_1 to obtain the fit.

Other multivariate and latent variable methods: bifactor_indices(), cfa_1(), cfa_2(), cfa_k(), ci_eigenvalue(), common_method_marker(), common_method_single_factor(), dmacs(), ecvi(), htmt(), irt_grm(), irt_information(), measurement_alignment(), measurement_invariance(), procrustes_phi(), simple_structure()

Examples

# Directly from the standardized loadings a paper reports:
average_variance_extracted(loadings = c(.8, .7, .6))

# From a fitted model, one AVE per factor (requires lavaan).
data(holzinger_swineford)
fit <- lavaan::cfa(
  "verbal    =~ t6_paragraph_comprehension + t7_sentence +
                t9_word_meaning
   deduction =~ t20_deduction + t22_problem_reasoning +
                t23_series_completion",
  data = holzinger_swineford)
ave_tbl <- average_variance_extracted(fit)
ave_tbl

# The Fornell and Larcker (1981) discriminant criterion compares each
# factor's AVE with the squared correlation between the factors: a
# factor should account for more of its own indicators' variance than
# it shares with the other factor. Here the comparison favors verbal
# and goes against deduction, whose AVE falls below the shared
# variance. Fitting with cfa_k(..., output = "measurement") puts the
# AVE values and the latent correlations in one table.
lavaan::lavInspect(fit, "cor.lv")["verbal", "deduction"]^2

# The broom verbs: one row per factor.
generics::tidy(ave_tbl)
generics::glance(ave_tbl)

Bayesian Independent-Samples t Analysis

Description

The Bayesian counterpart of the two-sample (pooled-variance) t test. It reports the posterior of the standardized mean difference \delta = (\mu_1 - \mu_2)/\sigma under the default Jeffreys-Zellner-Siow (JZS) prior, a Cauchy prior on \delta with Jeffreys priors on the nuisance parameters (the variance), summarized by its median, mean, a credible interval, and the probability statement P(\delta > 0 \mid \mathrm{data}). The JZS default Bayes factor (Rouder, Speckman, Sun, Morey, & Iverson, 2009) is also reported. See bayes_one_sample_t for the model, the package's interpretive stance, and the computational details (exact quadrature, no Monte Carlo error).

Usage

bayes_independent_t(
  x = NULL,
  y = NULL,
  mean_1 = NULL,
  sd_1 = NULL,
  n_1 = NULL,
  mean_2 = NULL,
  sd_2 = NULL,
  n_2 = NULL,
  prior_location = 0,
  prior_scale = sqrt(2)/2,
  prior_mean = NULL,
  prior_sd = NULL,
  conf_level = 0.95
)

Arguments

Details

The likelihood of the pooled-variance t statistic given \delta is noncentral t with \mathit{df} = n_1 + n_2 - 2 and noncentrality \delta \sqrt{n_1 n_2 / (n_1 + n_2)}; equal variances are assumed, as in the standard JZS development. The raw-scale rows transform the \delta summaries through the pooled standard deviation.

Value

A data.frame (class dmar_tbl) with the same rows as bayes_one_sample_t, plus n_1 and n_2.

Specifying the prior. The default prior on the standardized effect \delta is the JZS Cauchy centered at zero. Its prior_scale r is not a standard deviation: a Cauchy has no mean and no variance (those integrals diverge), so beliefs stated as prior moments cannot be expressed through it. What the scale does fix is the quartiles: half the prior mass lies within \pm r of the location, so the default r = \sqrt{2}/2 says a 50 percent prior bet that |\delta| < 0.71. A directional prior keeps the Cauchy and moves prior_location (Gronau, Ly, & Wagenmakers, 2020). A researcher who thinks in prior moments instead sets prior_mean and prior_sd, which use a normal prior with exactly those moments; the two families are exclusive.

The families are also linked by an exact identity: a Cauchy with location \mu and scale r is a normal prior N(\mu, r^2/z^2) whose z is standard normal, that is, a normal prior whose variance you are not sure of. Choosing the Cauchy is therefore choosing a normal prior with built-in doubt about its own width, which is why its tails are heavier and its Bayes factors more conservative. A normal matched to the Cauchy's interquartile range has prior_sd = 1.4826 * prior_scale. The full posterior of \delta is returned in the "posterior" attribute as a data frame of delta and density, so any posterior probability, not only the reported ones, can be computed from it.

A standardized effect size enters through the summary form directly: an observed Cohen's d with group sizes n_1 and n_2 is mean_1 = d, mean_2 = 0, sd_1 = 1, sd_2 = 1, since d is the mean difference in pooled standard deviation units.

Author(s)

Ken Kelley kkelley@nd.edu

References

Gronau, Q. F., Ly, A., & Wagenmakers, E.-J. (2020). Informed Bayesian t-tests. The American Statistician, 74(2), 137–143. doi:10.1080/00031305.2018.1562983

Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin & Review, 16(2), 225–237. doi:10.3758/PBR.16.2.225

Jeffreys, H. (1961). Theory of probability (3rd ed.). Oxford University Press.

Zellner, A., & Siow, A. (1980). Posterior odds ratios for selected regression hypotheses. In J. M. Bernardo, M. H. DeGroot, D. V. Lindley, & A. F. M. Smith (Eds.), Bayesian statistics: Proceedings of the First International Meeting (pp. 585–603). University of Valencia Press.

See Also

bayes_one_sample_t and bayes_paired_t for the other designs; ci_smd and welch_t for the frequentist analyses of the same comparison.

Other Bayesian t analyses: bayes_one_sample_t(), bayes_paired_t()

Examples

set.seed(113)
g1 <- rnorm(35, 105, 15)
g2 <- rnorm(35, 100, 15)
bayes_independent_t(g1, g2)

# The probability that the effect is positive is read directly off the
# p_delta_positive row: a probability statement, not a p-value.

Bayesian One-Sample t Analysis

Description

The Bayesian counterpart of the one-sample t test. It reports the posterior distribution of the standardized effect \delta = (\mu - \mu_0)/\sigma under the default Jeffreys-Zellner-Siow (JZS) prior, a Cauchy prior on \delta with Jeffreys priors on the nuisance parameters (the variance), summarized by its median, mean, a credible interval, and the direct probability statement P(\delta > 0 \mid \mathrm{data}). The JZS default Bayes factor (Rouder, Speckman, Sun, Morey, & Iverson, 2009) is also reported.

Usage

bayes_one_sample_t(
  x = NULL,
  mu_0 = 0,
  mean = NULL,
  sd = NULL,
  n = NULL,
  prior_location = 0,
  prior_scale = sqrt(2)/2,
  prior_mean = NULL,
  prior_sd = NULL,
  conf_level = 0.95
)

Arguments

Details

The posterior is a probability statement about the parameter given the model, the prior, and the data, and that is how these functions are meant to be read. A Bayes factor is a different kind of claim, a comparison of how well two models predicted the data, and it leans harder on the prior; it is reported because it may be helpful for some questions. Neither replaces the estimation-first habits of the rest of the package; ci_sm and ci_smd remain the frequentist complements.

With a Jeffreys prior on (\mu_0, \sigma^2) and \delta \sim \mathrm{Cauchy}(0, r), all inference flows through the observed t statistic, whose likelihood given \delta is noncentral t with noncentrality \delta \sqrt{n}. The posterior of \delta is computed by quadrature (no Monte Carlo error) and the Bayes factor by the one-dimensional integral of that likelihood against the Cauchy prior, the exact JZS form. The raw-scale rows transform the \delta summaries through the sample standard deviation (a plug-in, as is conventional for reporting).

Value

A data.frame (class dmar_tbl) with the posterior summaries of \delta (delta_posterior_median, delta_posterior_mean, delta_lower, delta_upper, p_delta_positive), the same summaries mapped to the raw mean difference scale (raw_*), the Bayes factors (bf_10, bf_01), the observed t and df, the prior_scale, and n.

Specifying the prior. The default prior on the standardized effect \delta is the JZS Cauchy centered at zero. Its prior_scale r is not a standard deviation: a Cauchy has no mean and no variance (those integrals diverge), so beliefs stated as prior moments cannot be expressed through it. What the scale does fix is the quartiles: half the prior mass lies within \pm r of the location, so the default r = \sqrt{2}/2 says a 50 percent prior bet that |\delta| < 0.71. A directional prior keeps the Cauchy and moves prior_location (Gronau, Ly, & Wagenmakers, 2020). A researcher who thinks in prior moments instead sets prior_mean and prior_sd, which use a normal prior with exactly those moments; the two families are exclusive.

The families are also linked by an exact identity: a Cauchy with location \mu and scale r is a normal prior N(\mu, r^2/z^2) whose z is standard normal, that is, a normal prior whose variance you are not sure of. Choosing the Cauchy is therefore choosing a normal prior with built-in doubt about its own width, which is why its tails are heavier and its Bayes factors more conservative. A normal matched to the Cauchy's interquartile range has prior_sd = 1.4826 * prior_scale. The full posterior of \delta is returned in the "posterior" attribute as a data frame of delta and density, so any posterior probability, not only the reported ones, can be computed from it.

A standardized effect size enters through the summary form directly: an observed d relative to mu_0 = 0 is mean = d, sd = 1.

Author(s)

Ken Kelley kkelley@nd.edu

References

Gronau, Q. F., Ly, A., & Wagenmakers, E.-J. (2020). Informed Bayesian t-tests. The American Statistician, 74(2), 137–143. doi:10.1080/00031305.2018.1562983

Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin & Review, 16(2), 225–237. doi:10.3758/PBR.16.2.225

Jeffreys, H. (1961). Theory of probability (3rd ed.). Oxford University Press.

Zellner, A., & Siow, A. (1980). Posterior odds ratios for selected regression hypotheses. In J. M. Bernardo, M. H. DeGroot, D. V. Lindley, & A. F. M. Smith (Eds.), Bayesian statistics: Proceedings of the First International Meeting (pp. 585–603). University of Valencia Press.

See Also

bayes_paired_t and bayes_independent_t for the two-sample designs; ci_sm for the frequentist standardized mean.

Other Bayesian t analyses: bayes_independent_t(), bayes_paired_t()

Examples

set.seed(113)
x <- rnorm(40, mean = 0.4, sd = 1)
bayes_one_sample_t(x)

# Against a nonzero comparison value, with a wider prior.
bayes_one_sample_t(x, mu_0 = 0.1, prior_scale = 1)

Bayesian Paired-Samples t Analysis

Description

The Bayesian counterpart of the paired t test, the analysis of bayes_one_sample_t applied to the within-pair differences. It reports the posterior of the standardized difference \delta = \mu_D / \sigma_D under the default Jeffreys-Zellner-Siow (JZS) prior, a Cauchy prior on \delta with Jeffreys priors on the nuisance parameters (the variance), summarized by its median, mean, a credible interval, and the probability statement P(\delta > 0 \mid \mathrm{data}). The JZS default Bayes factor (Rouder, Speckman, Sun, Morey, & Iverson, 2009) is also reported. See bayes_one_sample_t for the model, the package's interpretive stance, and the computational details (exact quadrature, no Monte Carlo error).

Usage

bayes_paired_t(
  x = NULL,
  y = NULL,
  mean_diff = NULL,
  sd_diff = NULL,
  n = NULL,
  prior_location = 0,
  prior_scale = sqrt(2)/2,
  prior_mean = NULL,
  prior_sd = NULL,
  conf_level = 0.95
)

Arguments

Value

A data.frame (class dmar_tbl) with the same rows as bayes_one_sample_t, where \delta is the standardized within-pair difference and the raw rows are on the difference scale; n is the number of pairs.

Specifying the prior. The default prior on the standardized effect \delta is the JZS Cauchy centered at zero. Its prior_scale r is not a standard deviation: a Cauchy has no mean and no variance (those integrals diverge), so beliefs stated as prior moments cannot be expressed through it. What the scale does fix is the quartiles: half the prior mass lies within \pm r of the location, so the default r = \sqrt{2}/2 says a 50 percent prior bet that |\delta| < 0.71. A directional prior keeps the Cauchy and moves prior_location (Gronau, Ly, & Wagenmakers, 2020). A researcher who thinks in prior moments instead sets prior_mean and prior_sd, which use a normal prior with exactly those moments; the two families are exclusive.

The families are also linked by an exact identity: a Cauchy with location \mu and scale r is a normal prior N(\mu, r^2/z^2) whose z is standard normal, that is, a normal prior whose variance you are not sure of. Choosing the Cauchy is therefore choosing a normal prior with built-in doubt about its own width, which is why its tails are heavier and its Bayes factors more conservative. A normal matched to the Cauchy's interquartile range has prior_sd = 1.4826 * prior_scale. The full posterior of \delta is returned in the "posterior" attribute as a data frame of delta and density, so any posterior probability, not only the reported ones, can be computed from it.

Author(s)

Ken Kelley kkelley@nd.edu

References

Gronau, Q. F., Ly, A., & Wagenmakers, E.-J. (2020). Informed Bayesian t-tests. The American Statistician, 74(2), 137–143. doi:10.1080/00031305.2018.1562983

Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin & Review, 16(2), 225–237. doi:10.3758/PBR.16.2.225

Jeffreys, H. (1961). Theory of probability (3rd ed.). Oxford University Press.

Zellner, A., & Siow, A. (1980). Posterior odds ratios for selected regression hypotheses. In J. M. Bernardo, M. H. DeGroot, D. V. Lindley, & A. F. M. Smith (Eds.), Bayesian statistics: Proceedings of the First International Meeting (pp. 585–603). University of Valencia Press.

See Also

bayes_one_sample_t for the model and stance; bayes_independent_t for unpaired groups; probability_of_superiority_paired and randomization_test_paired for other paired analyses.

Other Bayesian t analyses: bayes_independent_t(), bayes_one_sample_t()

Examples

set.seed(113)
before <- rnorm(30, 100, 12)
after  <- before + rnorm(30, 3, 6)
bayes_paired_t(after, before)

Bessel's (1818) Grouped Frequency Distribution of Bradley's Astronomical Observation Errors

Description

The nine-bin grouped frequency distribution that Friedrich Wilhelm Bessel published in 1818 for the absolute errors of 300 stellar position observations made by British Astronomer Royal James Bradley at the Greenwich Observatory between 1750 and 1762. Bessel compared the empirical distribution of these errors to the normal distribution, providing one of the early empirical demonstrations that observational errors are approximately normally distributed, a position Gauss had developed on theoretical grounds a decade earlier. The data are reproduced from Maxwell, Delaney, and Kelley (2027, Designing Experiments and Analyzing Data: A Model Comparison Perspective, 4th ed., Routledge), Table 1.4.

Usage

bessel_errors

Format

A data frame with 9 observations on 6 variables, one row per bin of the grouped frequency distribution. The error magnitudes are in seconds of arc.

bin

Integer bin index, 1 through 9.

lower

Lower edge of the bin (inclusive), in seconds of arc.

upper

Upper edge of the bin (exclusive), in seconds of arc.

midpoint

Bin midpoint, (lower + upper) / 2, in seconds of arc. The conventional plug-in value when approximating moments from a grouped frequency distribution.

observed

Empirical frequency: the number of Bradley's 300 absolute errors that fell in the bin.

expected

Expected frequency under a normal distribution with mean 0 and standard deviation approximately 0.22 seconds of arc (a least squares fit to the expected counts gives 0.216). These are Bessel's own normal-model expectations as reproduced in Maxwell, Delaney, and Kelley (2027). Both the observed and expected columns sum to 300.

Details

This data set ships in the original grouped form Bessel reported. The 300 individual error values are not available; what Bessel published, and what is reproduced here, is the 9-bin frequency distribution. Computations that require the underlying continuous values must either be approximated from the bin midpoints (the usual weighted-moments approach, illustrated in the examples) or estimated parametrically by assuming a distributional form within each bin.

Historical context. Friedrich Wilhelm Bessel (1784–1846) was a German astronomer and mathematician best known to statisticians for the Bessel correction (n - 1 in the unbiased variance estimator) and the Bessel functions. The 1818 monograph that contains this frequency distribution is part of a much larger effort to produce a reference catalog of stellar positions, the Fundamenta astronomiae, derived from the observations of James Bradley (1693–1762), the third Astronomer Royal of Britain and a pioneering observational astronomer. Bradley's position measurements were the most accurate of his era; the observational errors are small (most under 0.5 seconds of arc) and approximately normally distributed.

Why this data set matters for measurement and analysis. Bessel's 1818 comparison is one of the earliest empirical demonstrations that observational error is approximately normal, complementing the theoretical case Gauss had made on independent grounds. It is also a clean worked example for the approximation of moments from grouped frequency data when only binned counts (rather than individual observations) are available, a common situation in published reports.

Approximating moments from grouped data. When the underlying continuous values are unavailable, the standard approach is to plug the bin midpoints in for the unknown individual values and form a weighted mean and weighted variance using the bin frequencies as weights. The frequency-weighted mean and variance are

\bar{x}_w = \frac{\sum_k f_k m_k}{\sum_k f_k}, \qquad s^2_w = \frac{\sum_k f_k (m_k - \bar{x}_w)^2}{(\sum_k f_k) - 1}

where f_k is the bin frequency and m_k is the bin midpoint. The examples below compute both, on the observed frequencies and on Bessel's normal-model expected frequencies. The two are close, consistent with Bessel's conclusion that the empirical and theoretical distributions agree.

Author(s)

Ken Kelley

Source

Bessel, F. W. (1818). Fundamenta astronomiae pro anno MDCCLV deducta ex observationibus viri incomparabilis James Bradley in specula astronomica Grenovicensi per annos 1750–1762 institutis [Foundations of astronomy for the year 1755, deduced from the observations of the incomparable man James Bradley at the Greenwich astronomical observatory during 1750–1762]. Friedrich Nicolovius.

Reproduced in Maxwell, Delaney, and Kelley (2027), Table 1.4.

References

Bessel, F. W. (1818). Fundamenta astronomiae pro anno MDCCLV deducta ex observationibus viri incomparabilis James Bradley in specula astronomica Grenovicensi per annos 1750–1762 institutis. Friedrich Nicolovius.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 1, Section 1.4.5.2, on the historical and empirical basis for the normal distribution.)

Kelley, K. (2026). DMAR: Methods for the design, measurement, and analysis of human-centered outcomes in R [R package]. https://github.com/yelleKneK/DMAR

Stigler, S. M. (1986). The history of statistics: The measurement of uncertainty before 1900. Belknap Press of Harvard University Press. (See Chapter 5 on the development of the normal distribution and the role of astronomical errors.)

Examples

data(bessel_errors)
bessel_errors

# Total frequencies (each column should sum to 300).
colSums(bessel_errors[, c("observed", "expected")])

# Weighted mean of the absolute error, approximating the
# individual observations by the bin midpoints. Because every
# bin lower edge is at or above zero, this is the mean
# absolute error rather than the mean error itself.
wmean_obs <- with(bessel_errors,
                  sum(observed * midpoint) / sum(observed))
wmean_exp <- with(bessel_errors,
                  sum(expected * midpoint) / sum(expected))
c(observed = wmean_obs, expected_under_normal = wmean_exp)

# Weighted variance and standard deviation of the absolute
# error, using bin midpoints as plug-in values for the
# individual observations.
wvar_obs <- with(bessel_errors,
                 sum(observed * (midpoint - wmean_obs)^2) /
                 (sum(observed) - 1))
c(weighted_variance = wvar_obs,
  weighted_sd       = sqrt(wvar_obs))

# Side-by-side bar plot of observed and expected counts.
# Bessel's normal-model expectation tracks the empirical
# distribution closely except in the long right tail, where
# Bradley's three largest errors (counts 3, 1, 1) exceed what
# the normal model predicts.
op <- par(mar = c(5, 4, 4, 2))
barplot(rbind(bessel_errors$observed, bessel_errors$expected),
        beside    = TRUE,
        names.arg = sprintf("%.1f-%.1f",
                            bessel_errors$lower,
                            bessel_errors$upper),
        legend.text = c("Observed (Bradley)",
                        "Expected under normal model"),
        args.legend = list(x = "topright", bty = "n"),
        xlab = "Absolute error (seconds of arc)",
        ylab = "Frequency",
        main = "Bessel (1818) on Bradley's 300 stellar positions")
par(op)

Bifactor Model Dimensionality and Reliability Indices

Description

Computes the indices used to judge whether a multidimensional scale is nonetheless unidimensional enough to score as a single total (Rodriguez, Reise, and Haviland, 2016): the explained common variance (ECV), coefficient omega and omega hierarchical (omega_H) for the general factor, omega hierarchical subscale (omega_HS) for each group factor, the percentage of uncontaminated correlations (PUC), and coefficient H, the construct reliability (maximal reliability) of Hancock and Mueller (2001). The input is a fitted bifactor model in which one general factor loads on every item and each item loads on exactly one orthogonal group factor.

Usage

bifactor_indices(fit, general = NULL)

Arguments

Details

Let \lambda^g_i be item i's standardized loading on the general factor, \lambda^s_i its loading on its group factor, and \theta_i its standardized residual variance. With orthogonal factors the overall ECV is \sum_i (\lambda^g_i)^2 / \sum_i [(\lambda^g_i)^2 + (\lambda^s_i)^2]; omega and omega_H share the total-score variance (\sum_i \lambda^g_i)^2 + \sum_g (\sum_{i \in g} \lambda^s_i)^2 + \sum_i \theta_i as denominator, with the general part (\sum_i \lambda^g_i)^2 in the numerator of omega_H. Each group factor's omega_HS uses the analogous numerator and that subscale's own total variance. PUC is one minus the share of item pairs that fall within the same group factor. Coefficient H is computed from \sum \lambda^2 / (1 - \lambda^2) on the relevant loadings (see reliability_H). High ECV and PUC with a high omega_H support scoring a single total; substantial subscale omega_HS argues for reporting subscales as well.

A bifactor model that is over-parameterized for the data frequently yields an improper (Heywood) solution – a standardized loading outside [-1, 1] or a negative residual variance – whose indices are not trustworthy. When that happens the indices are still returned but a warning is issued and the "improper" attribute is set to TRUE.

Value

A data.frame (class dmar_tbl) with one row per factor (the general factor first, then each group factor) and columns factor, ECV, omega, omega_H, omega_HS, PUC, and H. Quantities that do not apply to a row are NA (for example omega_H and PUC on a group factor). The "improper" attribute flags a Heywood solution.

Author(s)

Ken Kelley kkelley@nd.edu

References

Hancock, G. R., & Mueller, R. O. (2001). Rethinking construct reliability within latent variable systems. In R. Cudeck, S. du Toit, & D. Sörbom (Eds.), Structural equation modeling: Present and future (pp. 195–216). Scientific Software International.

Reise, S. P. (2012). The rediscovery of bifactor measurement models. Multivariate Behavioral Research, 47(5), 667–696. doi:10.1080/00273171.2012.715555

Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models: Calculating and interpreting statistical indices. Psychological Methods, 21(2), 137–150.

See Also

reliability_omega, reliability_H.

Other multivariate and latent variable methods: average_variance_extracted(), cfa_1(), cfa_2(), cfa_k(), ci_eigenvalue(), common_method_marker(), common_method_single_factor(), dmacs(), ecvi(), htmt(), irt_grm(), irt_information(), measurement_alignment(), measurement_invariance(), procrustes_phi(), simple_structure()

Examples

# Nine items: one general factor and three orthogonal group factors.
set.seed(113)
n <- 600
g <- rnorm(n); grp <- list(rnorm(n), rnorm(n), rnorm(n))
X <- vapply(1:9, function(i)
  0.5 * g + 0.5 * grp[[ceiling(i / 3)]] + sqrt(0.5) * rnorm(n), numeric(n))
colnames(X) <- paste0("x", 1:9)
model <- "g  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
          f1 =~ x1 + x2 + x3
          f2 =~ x4 + x5 + x6
          f3 =~ x7 + x8 + x9"
fit <- lavaan::cfa(model, data = as.data.frame(X),
                   orthogonal = TRUE, std.lv = TRUE)
bifactor_indices(fit)

The Bryant-Paulson Generalized Studentized Range Distribution

Description

Distribution function (pbryant_paulson), quantile/critical-value function (qbryant_paulson), and density (dbryant_paulson) for the Bryant-Paulson generalized studentized range, the sampling distribution of the studentized range of covariate-adjusted means in the analysis of covariance (ANCOVA) when the covariate(s) are random. These are the analysis-of-covariance analogues of ptukey / qtukey and supply the critical values needed for Tukey-Kramer-type simultaneous confidence intervals on (and tests of) contrasts of adjusted means.

Usage

pbryant_paulson(q, num_covariates, num_groups, df, lower_tail = TRUE, ...)

qbryant_paulson(prob, num_covariates, num_groups, df, lower_tail = TRUE, ...)

dbryant_paulson(q, num_covariates, num_groups, df, ...)

Arguments

Details

The statistic. In a balanced ANCOVA with k groups and p random covariates, let \hat\theta_i be the adjusted group means and \hat\sigma_{y \mid x} the square root of the ANCOVA error mean square (on \nu degrees of freedom). The Bryant-Paulson statistic is the studentized range of the adjusted means,

Q \;=\; \frac{\max_i \hat\theta_i - \min_i \hat\theta_i}{\hat\sigma_{y\mid x}\sqrt{K_1 - K_2}},

where K_1 - K_2 is the design constant that scales the variance of a single adjusted mean (for a one-way design with n per group, K_1 - K_2 = 1/n). Crucially, the studentizer uses only this “between-only” standard error: the extra sampling variability induced by having to estimate the covariate adjustment from random covariates is carried by the distribution of Q itself, not by a per-comparison standard-error correction. This is what distinguishes the procedure from naively applying Tukey's method to adjusted means.

The distribution. Bryant and Paulson (1976) give the exact CDF of Q_p in their Equation (17), a single integral over a variable that combines the \chi^2_\nu error estimate with a random covariate-shrinkage factor \delta that, by their Equations (11)–(12), has a \mathrm{Beta}((\nu+1)/2,\, p/2) distribution. Carrying out the error integral with the studentized-range routine ptukey reduces Equation (17) to the equivalent one-dimensional form

P(Q_p \le q) \;=\; \int_0^1 \mathrm{ptukey}\!\left(q\sqrt{\delta};\, k,\, \nu\right)\, f_{\mathrm{Beta}}\!\left(\delta;\, \tfrac{\nu+1}{2},\, \tfrac{p}{2}\right) d\delta,

which this package evaluates (the reduction is exact; see the source-code comments in ‘R/bryant_paulson.R’ for the one-line derivation from Bryant and Paulson's p. 634 conditioning argument). When p = 0 the factor \delta degenerates at 1 and Q_p is exactly the ordinary studentized range (Bryant and Paulson, 1976, Sec. 1), so the code short-circuits to ptukey. The integral is evaluated with integrate after the change of variables \delta = 1 - u^2, which removes the endpoint singularity of the Beta weight at \delta = 1 when p = 1 and makes the quadrature converge in a few subdivisions; qbryant_paulson inverts it with uniroot. Bryant and Bruvold (1980) later showed the same distribution and critical values remain valid when the covariates are not identically distributed across groups (their grouped-covariate model, Eq. 1.3), and added the Duncan multiple-range extension.

Accuracy at small df. The error integral is carried out with ptukey for \nu \ge 7, where it is accurate to about 10^{-9}. For \nu < 7 ptukey's algorithm loses accuracy (at \nu = 3, k = 20 its probability error reaches \approx 3\times10^{-4}, enough to move the critical value by about 0.2, and it is larger at \nu = 2), so the studentized-range distribution is instead evaluated directly, without ptukey, by integrating the probability integral of the range against the \chi^2_\nu error density. The small-\nu path costs a fraction of a second.

Validation. The implementation reproduces Bryant and Paulson's (1976) Table 1 exactly, to the two decimal places tabled, over the whole of its range: both tail areas (\alpha = .05 and \alpha = .01), all three covariate counts (p = 1, 2, 3), every tabled number of groups (k = 2, \ldots, 8, 10, 12, 16, 20), and every tabled error degrees of freedom (\nu = 2, \ldots, 8, 10, 12, 14, 16, 18, 20, 24, 30, 40, 60, 120), which is 1188 critical values in all. The corners of the table are included: q_{.01;\,1,2,2} = 19.09, q_{.01;\,3,20,2} = 73.01, q_{.01;\,3,20,3} = 33.13, and q_{.05;\,1,6,14} = 4.83 (the value used in the Bryant and Bruvold, 1980 worked example). Two of the 1188 entries, q_{.01;\,2,8,3} and q_{.01;\,2,20,4}, have exact values of 23.165013 and 19.745008, each roughly 10^{-5} above the 23.165 and 19.745 half-way points, so they round to 23.17 and 19.75; the 1976 table rounds them down, to 23.16 and 19.74. Both values were confirmed to fourteen significant figures by two independent high-order quadrature engines that share no code with the package implementation. A large-scale simulation of the Bryant and Paulson statistic confirms the computed values independently, and the Bryant and Bruvold (1980) Table 2 Duncan ranges are reproduced as well. See the package tests.

Value

Numeric vectors. pbryant_paulson returns cumulative (or upper-tail) probabilities, dbryant_paulson returns density values, and qbryant_paulson returns critical values (quantiles) of the generalized studentized range. Results are recycled to the length of the longest of q/prob and the parameter arguments.

Author(s)

Ken Kelley kkelley@nd.edu

References

Bryant, J. L., & Paulson, A. S. (1976). An extension of Tukey's method of multiple comparisons to experimental designs with random concomitant variables. Biometrika, 63, 631–638. doi:10.1093/biomet/63.3.631

Bryant, J. L., & Bruvold, N. T. (1980). Multiple comparison procedures in the analysis of covariance. Journal of the American Statistical Association, 75(372), 874–880. doi:10.2307/2287175

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 9.)

See Also

ci_c_ancova_bp for the simultaneous confidence intervals these critical values produce; ptukey and qtukey for the ordinary (fixed-covariate or no-covariate) studentized range.

Examples

# Critical value from the worked example of Bryant and Bruvold (1980):
# k = 6 panels, p = 1 covariate, nu = 14 error df, alpha = .05. The quantile
# is found by inverting the distribution function with uniroot; it is 4.83,
# the entry in Table 1 of Bryant and Paulson (1976).
qbryant_paulson(0.95, num_covariates = 1, num_groups = 6, df = 14)

# The ordinary Tukey value (ignoring that the covariate is random and
# estimated) is smaller, 4.64, so it yields intervals that are too narrow:
qtukey(0.95, nmeans = 6, df = 14)

# How much too narrow: the Bryant-Paulson area beyond the Tukey value is the
# familywise error rate Tukey's method actually delivers in this design.
# With no covariate the two distributions coincide, so the area is .0499,
# the nominal .05 up to the rounding of 4.64 itself. Each additional random
# covariate carries more estimation uncertainty, stretches the distribution
# to the right, and pushes the rate up, to .063, .076, and .091.
pbryant_paulson(4.64, num_covariates = 0:3, num_groups = 6, df = 14,
                lower_tail = FALSE)

# The p = 0 entry above is exactly the ordinary studentized range:
ptukey(4.64, nmeans = 6, df = 14, lower.tail = FALSE)

One Factor Confirmatory Factor Analysis Model

Description

Fits a single factor (congeneric by default) confirmatory factor analysis model to raw item data or a sample covariance matrix. This is the one factor special case of cfa_k: the function is a convenience wrapper that only requires the data (and, when the data hold more than the items, a vector of item names), builds the one factor specification, and forwards everything else to cfa_k(). The factor is named f1, so the rows of the returned table are lambda_f1_1, lambda_f1_2, ..., psi_f1_1, ..., and omega_f1, exactly as a one factor cfa_k() call would report them (the syntax column names the item behind each number).

Usage

cfa_1(
  data = NULL,
  items = NULL,
  S = NULL,
  N = NULL,
  equal_loading = FALSE,
  equal_error = FALSE,
  estimator = "ML",
  missing = "listwise",
  se = "standard",
  conf_level = 0.95,
  output = c("verbose", "measurement", "summary", "standardized", "fit"),
  ...
)

Arguments

Details

The model is identified by fixing the factor variance to 1 and estimating every loading. equal_loading and equal_error impose the classical measurement structures on the single factor, and the header of the printed table names the structure implied by the constraints. For the composite reliability coefficient with the observed total variance in the denominator, use reliability_omega with denominator = "observed"; for the model implied omega, the omega_f1 row of this function's output and reliability_omega agree.

Ordered categorical items are not supported here; use cfa_k, whose ordered argument fits WLSMV with the theta parameterization and reports the Green and Yang (2009) categorical sum score omega, or reliability_omega_categorical.

Value

The value of the corresponding cfa_k call: a data.frame (classes dmar_cfa_k, dmar_tbl) with one row per parameter (estimate, se, z_value, p_value, ci_lower, ci_upper) followed by the fit rows, or the alternative shapes selected by output (see ?cfa_k).

Author(s)

Ken Kelley kkelley@nd.edu

References

Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using structural equation modeling: An alternative to coefficient alpha. Psychometrika, 74(1), 155–167. doi:10.1007/s11336-008-9099-3

Kelley, K., & Pornprasertmanit, S. (2016). Confidence intervals for population reliability coefficients: Evaluation of methods, recommendations, and software for composite measures. Psychological Methods, 21(1), 69–92. doi:10.1037/a0040086

McDonald, R. P. (1999). Test theory: A unified treatment. Lawrence Erlbaum Associates.

See Also

cfa_k for the general function this wraps (factor analysis with any number of factors, intercept constraints, ordered categorical items, and the measurement output); cfa_2 for the two factor wrapper; reliability_omega for coefficient omega with confidence intervals.

Other multivariate and latent variable methods: average_variance_extracted(), bifactor_indices(), cfa_2(), cfa_k(), ci_eigenvalue(), common_method_marker(), common_method_single_factor(), dmacs(), ecvi(), htmt(), irt_grm(), irt_information(), measurement_alignment(), measurement_invariance(), procrustes_phi(), simple_structure()

Examples

set.seed(113)
f <- rnorm(200)
loadings <- c(0.5, 0.6, 0.65, 0.7, 0.8)
X <- sapply(loadings, function(l) l * f + rnorm(200, sd = sqrt(1 - l^2)))
colnames(X) <- paste0("y", 1:5)

# All columns are items, so the data are all that is needed.
cfa_1(X)

# Equal loadings (essentially tau-equivalent), named in the header.
cfa_1(X, equal_loading = TRUE)

# From a covariance matrix and sample size, as a paper reports them.
cfa_1(S = cov(X), N = 200)

# A subset of columns via \'items\'.
cfa_1(X, items = c("y1", "y2", "y3", "y4"))

Two Factor Confirmatory Factor Analysis Model

Description

Fits a two factor confirmatory factor analysis model to raw item data or a sample covariance matrix. This is the two factor special case of cfa_k: the function is a convenience wrapper that only requires the items of each factor, builds the two factor specification, and forwards everything else to cfa_k(). The factors are named f1 and f2, so the rows of the returned table are lambda_f1_1, lambda_f2_1, phi_f1_f2 (the factor correlation), omega_f1, omega_f2, and so on, exactly as a two factor cfa_k() call would report them (the syntax column names the item behind each number). To name the factors substantively, or for three or more factors, call cfa_k directly.

Usage

cfa_2(
  data = NULL,
  factor_1,
  factor_2,
  S = NULL,
  N = NULL,
  equal_loading = FALSE,
  equal_error = FALSE,
  correlated_factors = TRUE,
  estimator = "ML",
  missing = "listwise",
  se = "standard",
  conf_level = 0.95,
  output = c("verbose", "measurement", "summary", "standardized", "fit"),
  ...
)

Arguments

Details

Each factor is identified by fixing its variance to 1 and estimating every loading; each item loads on exactly one factor (simple structure). The factor correlation is estimated by default and correlated_factors = FALSE fixes it to zero. Per-factor constraint vectors use the factor names, for example equal_loading = c(f1 = TRUE, f2 = FALSE).

Value

The value of the corresponding cfa_k call: a data.frame (classes dmar_cfa_k, dmar_tbl) with one row per parameter (estimate, se, z_value, p_value, ci_lower, ci_upper) followed by the fit rows, or the alternative shapes selected by output (see ?cfa_k).

Author(s)

Ken Kelley kkelley@nd.edu

See Also

cfa_k for the general function this wraps; cfa_1 for the one factor wrapper; htmt and average_variance_extracted for the discriminant and convergent validity summaries the output = "measurement" table reports alongside omega.

Other multivariate and latent variable methods: average_variance_extracted(), bifactor_indices(), cfa_1(), cfa_k(), ci_eigenvalue(), common_method_marker(), common_method_single_factor(), dmacs(), ecvi(), htmt(), irt_grm(), irt_information(), measurement_alignment(), measurement_invariance(), procrustes_phi(), simple_structure()

Examples

data(holzinger_swineford)

# Two factors, each named by its items.
cfa_2(holzinger_swineford,
      factor_1 = c("t6_paragraph_comprehension", "t7_sentence",
                   "t9_word_meaning"),
      factor_2 = c("t20_deduction", "t22_problem_reasoning",
                   "t23_series_completion"))

# The measurement properties: omega, ave, and H per factor, the
# factor correlation, and the htmt ratio.
cfa_2(holzinger_swineford,
      factor_1 = c("t6_paragraph_comprehension", "t7_sentence",
                   "t9_word_meaning"),
      factor_2 = c("t20_deduction", "t22_problem_reasoning",
                   "t23_series_completion"),
      output = "measurement")

Multiple-Factor Confirmatory Factor Analysis Model

Description

Fits a confirmatory factor analysis model with one or more factors, where each factor is specified by naming its indicator variables and the measurement structure is specified by describing what is constrained (equal_loading, equal_intercept, equal_error) rather than by more technical terms. The function then reports which classical measurement structure the description implies (congeneric, essentially tau-equivalent, tau-equivalent, essentially parallel, or parallel), the parameter estimates with confidence intervals, fit information, and, per factor, coefficient omega, the average variance extracted (AVE), and coefficient H, each with a delta method standard error and confidence interval computed by lavaan from defined parameters (no additional packages are involved).

Usage

cfa_k(
  data = NULL,
  factors,
  S = NULL,
  N = NULL,
  M = NULL,
  equal_loading = FALSE,
  equal_intercept = FALSE,
  equal_error = FALSE,
  correlated_factors = TRUE,
  meanstructure = NULL,
  estimator = "ML",
  missing = "listwise",
  ordered = NULL,
  se = "standard",
  conf_level = 0.95,
  output = c("verbose", "measurement", "summary", "standardized", "fit"),
  ...
)

Arguments

Details

With ordered items, the model is fit by WLSMV to polychoric correlations with thresholds. A sum score of ordered items lives on the metric of the observed categories, not on the latent response metric of the polychoric loadings, so for a factor whose items are ordered the reported omega is the Green and Yang (2009) categorical sum score omega computed from the same fit; the substitution is announced in a message and recorded in the omega_metric attribute, and the delta method interval columns are NA for those rows because the delta method interval describes the latent response metric. The categorical sum score omega is the coefficient Kelley and Pornprasertmanit (2016) call categorical omega. AVE and coefficient H concern the latent response correlations themselves and are reported unchanged on that metric. For a bootstrap confidence interval on a categorical omega, use reliability_omega_categorical on the factor's items.

cfa_1 and cfa_2 are convenience wrappers around this function for the one and two factor cases: cfa_1() takes a vector of items and fits one factor over them, and cfa_2() takes the items of each of two factors. Both forward every argument here, so their results are this function's results, with the factors named f1 (and f2).

Describing the model instead of naming it. The classical measurement structures are nested patterns of within-factor equality constraints (Lord & Novick, 1968; Graham, 2006):

congeneric

loadings, intercepts, and error variances all free.

essentially tau-equivalent

equal loadings; intercepts and error variances free.

tau-equivalent

equal loadings and equal intercepts; error variances free.

essentially parallel

equal loadings and equal error variances; intercepts free.

parallel

equal loadings, equal intercepts, and equal error variances.

The caller states the constraints; the function reports the implied name, per factor, in the printed header and in the "model" attribute of the returned table. The distinction between tau-equivalent and essentially tau-equivalent (and between parallel and essentially parallel) lives entirely in the mean structure: the covariance structure of the two members of each pair is identical, so without intercepts in the model only the "essentially" form can be claimed. That is why equal_intercept requires raw data or M: covariances alone cannot speak to it. A constraint pattern outside the classical list (for example equal error variances with free loadings) is fit as requested and labeled descriptively, since it has no conventional name.

Identification fixes each factor variance to 1 (and each factor mean to 0 when the mean structure is modeled), so all loadings are estimated and within-factor equality constraints are meaningful. Two cfa_k() fits that differ only in descriptor settings are nested, so output = "fit" feeds lavaan::lavTestLRT() directly (with estimator = "MLR", lavaan applies the scaled difference test).

Measurement properties. For factor f with unstandardized loadings \lambda_j and error variances \psi_j (factor variance 1): coefficient omega = (\sum_j \lambda_j)^2 / ((\sum_j \lambda_j)^2 + \sum_j \psi_j) (McDonald, 1999), the reliability of the unit-weighted composite; the average variance extracted = J^{-1} \sum_j \lambda_j^2 / (\lambda_j^2 + \psi_j) (Fornell & Larcker, 1981), the mean proportion of item variance the factor accounts for; and coefficient H = (1 + (\sum_j \lambda_j^2/\psi_j)^{-1})^{-1} (Hancock & Mueller, 2001), the reliability of the optimally weighted composite, which no single item can drag below its value for any subset. All three are computed as lavaan defined parameters, so each carries a delta method standard error and a conf_level confidence interval in the same table as the model parameters.

Discriminant validity. Three complementary readings come from output = "measurement": (a) the latent correlation phi for a factor pair, with a confidence interval whose upper limit near 1 means the data cannot distinguish the two factors; (b) the Fornell and Larcker (1981) comparison, which asks whether each factor's ave exceeds the squared phi of its pairs (the factor should share more variance with its own items than with the other factor); and (c) for raw data, the model-free htmt ratio (Henseler, Ringle, & Sarstedt, 2015). The rows report the numbers and their uncertainty; the judgment is the researcher's.

Confidence interval conventions. Parameter rows (including omega, ave, H, and phi) use conf_level. The RMSEA interval follows its own convention: the rmsea_ci_level row records the level actually used (0.90, the conventional level for RMSEA, as in lavaan), and rmsea_ci_lower / rmsea_ci_upper are that interval.

Common row names under term: lambda_<factor>_<j> (loadings; lambda_<factor> when equated), psi_<factor>_<j> (error variances; psi_<factor> when equated), nu_<factor>_<j> (intercepts, when the mean structure is modeled; nu_<factor> when equated), phi_<factor> (factor variance, fixed to 1), phi_<factor1>_<factor2> (latent correlation), the per-factor defined parameters (loading_sum_<factor>, error_sum_<factor>, omega_<factor>, ave_<factor>, H_<factor>), and the fit rows chi_square, df, p_chi_square, cfi, tli, nnfi, rmsea, rmsea_ci_lower, rmsea_ci_upper, rmsea_ci_level, srmr, AIC, BIC, H0, H1.

Value

For output = "verbose" (default) and output = "measurement", a data.frame (classes dmar_cfa_k, dmar_tbl) with columns syntax, term, estimate, se, z_value, p_value, ci_lower, ci_upper. The "model" attribute is a named character vector giving, per factor, the implied classical structure; the printed header displays it. For output = "summary", the lavaan summary object; for "standardized", the standardized solution; for "fit", the lavaan fit object.

Author(s)

Ken Kelley kkelley@nd.edu

References

Browne, M. W. (1974). Generalized least squares estimators in the analysis of covariance structures. South African Statistical Journal, 8, 1–24.

Browne, M. W. (1984). Asymptotically distribution-free methods for the analysis of covariance structures. British Journal of Mathematical and Statistical Psychology, 37, 62–83.

Muthén, B. (1984). A general structural equation model with dichotomous, ordered categorical, and continuous latent variable indicators. Psychometrika, 49(1), 115–132.

Muthén, B., du Toit, S. H. C., & Spisic, D. (1997). Robust inference using weighted least squares and quadratic estimating equations in latent variable modeling with categorical and continuous outcomes. Unpublished technical report.

Satorra, A., & Bentler, P. M. (1994). Corrections to test statistics and standard errors in covariance structure analysis. In A. von Eye & C. C. Clogg (Eds.), Latent variables analysis: Applications for developmental research (pp. 399–419). Thousand Oaks, CA: Sage.

Fornell, C., & Larcker, D. F. (1981). Evaluating structural equation models with unobservable variables and measurement error. Journal of Marketing Research, 18(1), 39–50.

Graham, J. M. (2006). Congeneric and (essentially) tau-equivalent estimates of score reliability: What they are and how to use them. Educational and Psychological Measurement, 66(6), 930–944. doi:10.1177/0013164406288165

Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using structural equation modeling: An alternative to coefficient alpha. Psychometrika, 74(1), 155–167. doi:10.1007/s11336-008-9099-3

Hancock, G. R., & Mueller, R. O. (2001). Rethinking construct reliability within latent variable systems. In R. Cudeck, S. du Toit, & D. Sörbom (Eds.), Structural equation modeling: Present and future (pp. 195–216). Scientific Software International.

Henseler, J., Ringle, C. M., & Sarstedt, M. (2015). A new criterion for assessing discriminant validity in variance-based structural equation modeling. Journal of the Academy of Marketing Science, 43(1), 115–135. doi:10.1007/s11747-014-0403-8

Kelley, K., & Pornprasertmanit, S. (2016). Confidence intervals for population reliability coefficients: Evaluation of methods, recommendations, and software for composite measures. Psychological Methods, 21, 69–92. doi:10.1037/a0040086

Lord, F. M., & Novick, M. R. (1968). Statistical theories of mental test scores. Addison-Wesley.

McDonald, R. P. (1999). Test theory: A unified treatment. Erlbaum.

See Also

cfa_1 and cfa_2 for the one and two factor convenience wrappers; plot_cfa_k to display the estimates and the equality question visually; reliability_omega, reliability_H, average_variance_extracted, and htmt for the measurement properties as standalone functions; measurement_invariance for the across-group analog of these within-factor constraints; lavaan, lavTestLRT.

Other multivariate and latent variable methods: average_variance_extracted(), bifactor_indices(), cfa_1(), cfa_2(), ci_eigenvalue(), common_method_marker(), common_method_single_factor(), dmacs(), ecvi(), htmt(), irt_grm(), irt_information(), measurement_alignment(), measurement_invariance(), procrustes_phi(), simple_structure()

Examples

data(holzinger_swineford)
hs_factors <- list(
  verbal    = c("t6_paragraph_comprehension", "t7_sentence",
                "t9_word_meaning"),
  deduction = c("t20_deduction", "t22_problem_reasoning",
                "t23_series_completion"))

# Congeneric measurement model for both factors (the default:
# nothing is constrained, and the header names the structure).
cfa_k(holzinger_swineford, hs_factors)

# Equal loadings within every factor. Because only the covariance
# structure identifies this constraint, the implied structure is
# essentially tau-equivalent, and the header says so.
cfa_k(holzinger_swineford, hs_factors, equal_loading = TRUE)

# Measurement properties: omega, ave, and H per factor (each with a
# delta method standard error and confidence interval), the latent
# correlations, and the htmt ratios.
cfa_k(holzinger_swineford, hs_factors, output = "measurement")

# Descriptors can differ by factor: here the verbal loadings are
# equated and the deduction loadings are left free, and the header
# names each factor's structure separately.
cfa_k(holzinger_swineford, hs_factors,
      equal_loading = c(verbal = TRUE, deduction = FALSE))

# Equal loadings and intercepts (tau-equivalent), then also equal
# error variances (parallel). The mean structure is added because
# equal_intercept asks about it, so the nu terms join the table.
cfa_k(holzinger_swineford, hs_factors, equal_loading = TRUE,
      equal_intercept = TRUE)
cfa_k(holzinger_swineford, hs_factors, equal_loading = TRUE,
      equal_intercept = TRUE, equal_error = TRUE)

# Ordered categorical items: two correlated factors of three
# four-category items each, generated here. The model is fit by WLSMV
# to polychoric correlations, and each ordered factor's omega is
# reported on the categorical sum score metric (Green & Yang, 2009),
# which is why the interval columns of those two rows are NA.
set.seed(113)
n <- 200
eta_a <- rnorm(n)
eta_b <- 0.5 * eta_a + sqrt(1 - 0.5^2) * rnorm(n)
lambda <- c(0.6, 0.7, 0.8)
lat <- cbind(outer(eta_a, lambda), outer(eta_b, lambda)) +
  matrix(rnorm(n * 6), n, 6) %*% diag(sqrt(1 - c(lambda, lambda)^2))
likert <- as.data.frame(apply(lat, 2, function(x)
  as.integer(cut(x, breaks = c(-Inf, -1, 0, 1, Inf)))))
names(likert) <- paste0("item_", 1:6)
cfa_k(likert,
      list(scale_a = paste0("item_", 1:3),
           scale_b = paste0("item_", 4:6)),
      ordered = TRUE, output = "measurement")

# Does the equal-loadings description hold? Two fits that differ only
# in a descriptor are nested, so output = "fit" hands them straight to
# lavaan's likelihood ratio test.
fit_free  <- cfa_k(holzinger_swineford, hs_factors, output = "fit")
fit_equal <- cfa_k(holzinger_swineford, hs_factors,
                   equal_loading = TRUE, output = "fit")
lavaan::lavTestLRT(fit_free, fit_equal)

Confidence Interval for the Population Squared Multiple Correlation Coefficient

Description

Constructs a confidence interval for the population squared multiple correlation coefficient \rho^2 by inverting the sampling distribution of the sample R^2. The confidence interval is for the population value \rho^2; the required input is the corresponding sample value, the observed sample squared multiple correlation coefficient R^2 (or, equivalently, the observed F-statistic and degrees of freedom). The right choice of sampling distribution, and so the right interval, depends on whether the predictors are treated as random draws from a joint multivariate normal distribution (the default and the typical case in the behavioral, educational, and social sciences) or as fixed by design (planned dosing levels, factorial covariates, etc.). The function selects the sampling distribution via random_predictors and inverts the corresponding noncentral distribution; the construction is the regression analogue of a noncentral distribution based CI on a standardized effect size (Steiger & Fouladi, 1992; Kelley, 2007).

Usage

ci_R2(
  R2 = NULL,
  df_1 = NULL,
  df_2 = NULL,
  conf_level = 0.95,
  random_predictors = TRUE,
  F_value = NULL,
  N = NULL,
  p = NULL,
  alpha_lower = NULL,
  alpha_upper = NULL,
  tol = 1e-09
)

Arguments

Details

Fixed vs.\ random predictors. The two regression models give different sampling distributions for the sample R^2, and so different confidence intervals. Under fixed predictors the design matrix is treated as constant in hypothetical replications of the study, and the omnibus F-statistic F = (R^2 / p) / ((1 - R^2) / (N - p - 1)) follows a noncentral F with p and N - p - 1 degrees of freedom and noncentrality \lambda = N \rho^2 / (1 - \rho^2) (Cohen, 1988); the CI is obtained by inverting that distribution at the supplied confidence level (see ci_nc_F). Under random predictors the design matrix is itself a draw from a joint multivariate normal distribution and the unconditional sampling distribution of the sample R^2 is given by Lee (1971); ci_R2 uses the Lee (1971) bisection (the same construction Algina and Olejnik 2000 implemented in SAS) to invert that distribution. Gatsonis and Sampson (1989) document the comparison and show that treating random predictors as fixed tends to over-state precision (and so under-state the CI width); the discrepancy is modest at moderate to large N but non-trivial at small N with moderate-to-large effects. In the behavioral, educational, and social sciences predictor variables are almost always random, so the default is random_predictors = TRUE; pass random_predictors = FALSE for designs in which the predictor variables are fixed by design.

Value

A 3-row data.frame with columns term, value, prob_less, and prob_greater. The rows are ordered "lower_limit" (lower confidence limit on the population \rho^2), "R2" (the sample squared multiple correlation coefficient supplied by the user, the point estimate), and "upper_limit" (upper confidence limit on the population \rho^2), so the point estimate sits between its confidence limits. The prob_less and prob_greater columns report the achieved lower-tail and upper-tail error probabilities at each limit (they are NA for the "R2" estimate row). For random-predictor mode (random_predictors = TRUE) the limits are computed via the Lee (1971) bisection over the multiple-correlation sampling distribution; for fixed-predictor mode they are computed by inversion of the noncentral F distribution.

Author(s)

Ken Kelley kkelley@nd.edu

References

Algina, J. & Olejnik, S. (2000). Determining sample size for accurate estimation of the squared multiple correlation coefficient. Multivariate Behavioral Research, 35, 119–137. doi:10.1207/s15327906mbr3501_5

Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Hillsdale, NJ: Lawrence Erlbaum.

Gatsonis, C., & Sampson, A. R. (1989). Multiple correlation: Exact power and sample size calculations. Psychological Bulletin, 106(3), 516–524.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K. (2008). Sample size planning for the squared multiple correlation coefficient: Accuracy in parameter estimation via narrow confidence intervals. Multivariate Behavioral Research, 43, 524–555. doi:10.1080/00273170802490632

Lee, Y. S. (1971). Some results on the sampling distribution of the multiple correlation coefficient. Journal of the Royal Statistical Society, Series B, 33(1), 117–130.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on R^2 as a model comparison effect size.)

Smithson, M. (2003). Confidence intervals. Thousand Oaks, CA: Sage Publications.

Steiger, J. H., & Fouladi, R. T. (1992). R2: A computer program for interval estimation, power calculations, sample size estimation, and hypothesis testing in multiple regression. Behavior Research Methods, Instruments, & Computers, 24(4), 581–582. doi:10.3758/BF03203611

See Also

ss_aipe_R2, ci_nc_F

Other confidence intervals for effect sizes: ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# For random predictor variables.
ci_R2(R2 = .25, N = 100, p = 5, conf_level = .95, random_predictors = TRUE)

ci_R2(F_value = 6.266667, N = 100, p = 5, conf_level = .95, random_predictors = TRUE)

# For fixed predictor variables.
ci_R2(R2 = .25, N = 100, p = 5, conf_level = .95, random_predictors = FALSE)

ci_R2(F_value = 6.266667, N = 100, p = 5, conf_level = .95, random_predictors = FALSE)

# One sided confidence intervals when predictors are random.
ci_R2(R2 = .25, N = 100, p = 5, alpha_lower = .05, alpha_upper = 0,
      conf_level = NULL, random_predictors = TRUE)

ci_R2(R2 = .25, N = 100, p = 5, alpha_lower = 0, alpha_upper = .05,
      conf_level = NULL, random_predictors = TRUE)

# One sided confidence intervals when predictors are fixed.
ci_R2(R2 = .25, N = 100, p = 5, alpha_lower = .05, alpha_upper = 0,
      conf_level = NULL, random_predictors = FALSE)

ci_R2(R2 = .25, N = 100, p = 5, alpha_lower = 0, alpha_upper = .05,
      conf_level = NULL, random_predictors = FALSE)

Confidence Interval for a Contrast in a Fixed Effects ANOVA

Description

Computes the confidence interval for an unstandardized contrast of means in a fixed effects analysis of variance, so a focused comparison among groups (a pairwise difference or any weighted combination of the means) is reported with its precision and in the units of the response. Homogeneity of variance is assumed, as in the ANOVA on which s_anova is based.

Usage

ci_c(
  means = NULL,
  s_anova = NULL,
  c_weights = NULL,
  n = NULL,
  N = NULL,
  psi = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  df_error = NULL,
  ...
)

Arguments

Value

A 3-row data.frame with columns term and value. The term values are "lower_limit" (the lower confidence limit on the population contrast), "contrast" (the estimated unstandardized contrast), and "upper_limit" (the upper limit).

Note

Be sure to use the standard deviation and not the error variance for s_anova, not the square of this value (the error variance) which would come from the source table (i.e., use the root mean square error, not the mean square error).

Be sure to use fractional c-weights when doing complex contrasts (not integers) to specify c_weights. For example, in an ANCOVA of four groups, if the user wants to compare the mean of group 1 and 2 with the mean of group 3 and 4, c_weights should be specified as c(0.5, 0.5, -0.5, -0.5) rather than c(1, 1, -1, -1). Make sure the sum of the contrast weights is zero.

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Lai, K., & Kelley, K. (2012). Accuracy in parameter estimation for ANCOVA and ANOVA contrasts: Sample size planning via narrow confidence intervals. British Journal of Mathematical and Statistical Psychology, 65, 350–370. doi:10.1111/j.2044-8317.2011.02029.x

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 4 on individual comparisons of means.)

Steiger, J. H. (2004). Beyond the F Test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

ci_sc

Other confidence intervals for effect sizes: ci_R2(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Here is a four group example. Suppose that the means of groups 1--4 are 2, 4, 9,
# and 13, respectively. Further, let the error variance be .64 and thus the standard
# deviation would be .80 (note we use the standard deviation in the function, not the
# variance). The contrast of interest here is the average of groups 1 and 4 versus the
# average of groups 2 and 3.
ci_c(means = c(2, 4, 9, 13), s_anova = .80, c_weights = c(.5, -.5, -.5, .5),
n = c(3, 3, 3, 3), N = 12, conf_level = .95)

# Here is an example with two groups.
ci_c(means = c(1.6, 0), s_anova = .80, c_weights = c(1, -1),
n = c(10, 10), N = 20, conf_level = .95)

# An example given by Maxwell, Delaney, & Kelley (2027) :
# 20 subjects of mild hypertensives are assigned to one of four treatments: drug
# therapy, biofeedback, dietary modification, and a treatment combining all the
# three previous treatments. Subjects' blood pressure is measured two weeks
# after the termination of treatment. Now we want to form a 95% level
# confidence interval for the difference in blood pressure between subjects
# who received drug treatment and those who received biofeedback treatment

## Drug group's mean = 94; group size=4
## Biofeedback group's mean = 91; group size=6
## Diet group's mean = 92; group size=5
## Combination group's mean = 83; group size=5
## Mean Square Within (i.e., 'error_variance') = 67.375

ci_c(means = c(94, 91, 92, 83), s_anova = sqrt(67.375), c_weights = c(1, -1, 0, 0),
n = c(4, 6, 5, 5), N = 20, conf_level = .95)

Confidence Interval for an (Unstandardized) Contrast in ANCOVA With One Covariate

Description

Calculates the confidence interval for an unstandardized contrast in the one-covariate ANCOVA. Two procedures are available through the procedure argument. The default ("t") returns a single per-comparison interval based on the t distribution: it gives the correct (1 - conf_level) coverage for one contrast chosen in advance, and its standard error includes the (\sum c_i \bar X_i)^2 / SS_{\mathrm{within}(x)} term that accounts for the covariate separation between the groups in that one contrast. The "bryant_paulson" procedure instead returns Bryant–Paulson simultaneous (familywise) intervals over a whole family of contrasts of adjusted means; when selected, ci_c_ancova simply forwards its arguments to ci_c_ancova_bp and returns that result. See Details for which to use when.

Usage

ci_c_ancova(
  psi = NULL,
  adj_means = NULL,
  s_ancova = NULL,
  c_weights,
  n,
  cov_means,
  SSwithin_x,
  conf_level = 0.95,
  procedure = c("t", "bryant_paulson"),
  ...
)

Arguments

Details

Per-comparison versus simultaneous. The two procedures answer different questions and are not interchangeable. Use the default procedure = "t" when a single contrast was planned in advance: the interval has exact per-comparison coverage and its width reflects the covariate adjustment for that specific contrast through the (\sum c_i \bar X_i)^2 / SS_{\mathrm{within}(x)} term. Use procedure = "bryant_paulson" when several contrasts (for example all pairwise comparisons of adjusted means) are examined together and the coverage statement must hold simultaneously across the family: the Bryant–Paulson generalized studentized range supplies a larger critical value that controls the familywise error rate and, because the covariates are random, correctly absorbs the extra sampling uncertainty from estimating the covariate adjustment (which holds on average over the covariate distribution, so the per-contrast separation term is not added again). A per-comparison interval used as if it were simultaneous understates the family error rate; a simultaneous interval used for one planned contrast is wider than necessary. See ci_c_ancova_bp for the full description of the simultaneous procedure and its arguments.

Value

A 3-row data.frame with columns term and value (numeric). The term values are "lower_limit" (the lower confidence limit on the unstandardized ANCOVA contrast), "psi" (the unstandardized contrast point estimate), and "upper_limit" (the upper limit).

When procedure = "bryant_paulson", the return value is whatever ci_c_ancova_bp returns (a table with one row per contrast and columns contrast, estimate, lower_limit, and upper_limit).

Note

Be sure to use the standard deviation and not the error variance for s_ancova, not the square of this value which would come from the source table (i.e., do not use the variance of the error but rather use the square root).

If n receives a single number, that number is considered as the sample size per group. If n receives a vector, the vector is considered as the sample size of each group.

Be sure to use fractions not the integers to specify c_weights. For example, in an ANCOVA of four groups, if the user wants to compare the mean of group 1 and 2 with the mean of group 3 and 4, c_weights should be specified as c(0.5, 0.5, -0.5, -0.5) rather than c(1, 1, -1, -1). Make sure the sum of the contrast weights are zero.

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Lai, K., & Kelley, K. (2012). Accuracy in parameter estimation for ANCOVA and ANOVA contrasts: Sample size planning via narrow confidence intervals. British Journal of Mathematical and Statistical Psychology, 65, 350–370. doi:10.1111/j.2044-8317.2011.02029.x

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 9.)

See Also

ci_c,ci_sc_ancova, ci_c_ancova_bp for the Bryant–Paulson simultaneous procedure

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Maxwell, Delaney, & Kelley (2027) offer an example that 30 depressive
# individuals are randomly assigned to three groups, 10 in each, and ANCOVA
# is performed on the posttest scores using the participants' pretest
# scores as the covariate. The means of pretest scores of group 1 to 3 are
# 17, 17.7, and 17.4, respectively, and the adjusted means of groups 1 to 3
# are 7.5, 12, and 14, respectively. The error variance in ANCOVA is 29,
# and the sum of squares within groups from ANOVA on the covariate is 752.5.

# To obtain the confidence interval for adjusted mean of group 1 versus group 2:
ci_c_ancova(adj_means = c(7.5, 12, 14), s_ancova = sqrt(29),
            c_weights = c(1, -1, 0), n = 10,
            cov_means = c(17, 17.7, 17.4), SSwithin_x = 752.5)

# That interval is the right one for a single contrast planned in advance.
# For the family of all three pairwise comparisons of the adjusted means,
# with coverage that holds simultaneously across the family, select
# procedure = "bryant_paulson"; the call is forwarded to ci_c_ancova_bp().
# The simultaneous limits are wider, which is the price of the family
# statement: group 1 against group 2 runs from -10.61 to 1.61, where the
# single planned contrast above ran from -9.46 to 0.46.
ci_c_ancova(adj_means = c(7.5, 12, 14), s_ancova = sqrt(29), n = 10,
            procedure = "bryant_paulson")

Bryant–Paulson Simultaneous Confidence Intervals for Contrasts of Adjusted Means in ANCOVA

Description

Constructs Tukey–Kramer-type simultaneous confidence intervals on one or more contrasts of covariate-adjusted means in the analysis of covariance (ANCOVA) when the covariate(s) are random, using the Bryant–Paulson generalized studentized range (see bryant_paulson). Unlike the per-comparison interval of ci_c_ancova, these intervals control the familywise error rate over the whole family of comparisons and, through the Bryant–Paulson critical value, correctly account for the extra sampling uncertainty that comes from estimating the covariate adjustment from random covariates. Naively applying Tukey's method to adjusted means ignores that uncertainty and produces intervals that are too narrow (below-nominal coverage); a simulation study of that undercoverage is maintained alongside the package.

Usage

ci_c_ancova_bp(
  adj_means,
  s_ancova,
  c_weights = NULL,
  n,
  num_covariates = 1,
  df = NULL,
  conf_level = 0.95,
  contrast_type = c("pairwise", "allowance"),
  ...
)

Arguments

Details

The interval. For a contrast \psi = \sum_i c_i \theta_i of adjusted means, the simultaneous interval is

\hat\psi \;\pm\; q_{\alpha;\,p,k,\nu}\; \hat\sigma_{y\mid x}\; w(c),

where q_{\alpha;\,p,k,\nu} is the upper-\alpha Bryant–Paulson critical value (qbryant_paulson) and the width factor is w(c) = \tfrac{1}{\sqrt2}\sqrt{\sum_i c_i^2/n_i} for contrast_type = "pairwise" or w(c) = \tfrac12 \sum_i |c_i| \sqrt{1/n} for contrast_type = "allowance" (balanced n). For a pairwise difference with common n both reduce to q_{\alpha;\,p,k,\nu}\,\hat\sigma_{y\mid x}\sqrt{1/n}, reproducing the critical difference of Bryant and Bruvold (1980).

No per-comparison covariate term. By design the standard error here does not include the (\bar X_i - \bar X_j)^2 / SS_{\mathrm{within}(x)} term that appears in a single-comparison ANCOVA interval (ci_c_ancova). In the Bryant–Paulson framework the random-covariate uncertainty is carried by the (larger) critical value, which holds on average over the covariate distribution; adding the per-pair term as well would double-count it.

Unequal sample sizes. For unbalanced designs the "pairwise" standard error uses \sqrt{c_i^2/n_i} directly (the Tukey–Kramer generalization); coverage is then approximate but typically very close to nominal and slightly conservative.

Value

A data.frame (class dmar_tbl) with one row per contrast and columns contrast (a label), estimate (the contrast of adjusted means \hat\psi), lower_limit, and upper_limit. The Bryant–Paulson critical value used is stored in the "critical_value" attribute and the confidence level in the "conf_level" attribute (printed beneath the table).

Author(s)

Ken Kelley kkelley@nd.edu

References

Bryant, J. L., & Paulson, A. S. (1976). An extension of Tukey's method of multiple comparisons to experimental designs with random concomitant variables. Biometrika, 63, 631–638.

Bryant, J. L., & Bruvold, N. T. (1980). Multiple comparison procedures in the analysis of covariance. Journal of the American Statistical Association, 75(372), 874–880. doi:10.2307/2287175

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 9.)

See Also

bryant_paulson for the underlying distribution; ci_c_ancova for the per-comparison interval; ancova for an ANCOVA fit.

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Bryant and Bruvold (1980) worked example: 6 panels, 1 covariate, nu = 14,
# ANCOVA error MS = 0.01326. Here the design is a randomized block with
# s = 4 blocks, so the per-group "n" for the adjusted-mean SE is 4 and the
# error df (14) must be supplied directly. The multiplier for these
# intervals is a Bryant-Paulson quantile, which has no closed form; it is
# obtained by inverting the distribution function with a root search.
adj <- c(3.595, 3.619, 4.102, 4.515, 4.618, 4.876)
bp <- ci_c_ancova_bp(adj_means = adj, s_ancova = sqrt(0.01326),
                     n = 4, num_covariates = 1, df = 14)
bp
# The multiplier is 4.83 and every pairwise critical difference is 0.278,
# matching the paper; the 15 intervals hold jointly at the 95 percent level.
# The multiplier is kept on the result, so the critical difference can be
# rebuilt by hand as q * s_ancova * sqrt(1/n):
attr(bp, "critical_value")
attr(bp, "critical_value") * sqrt(0.01326) * sqrt(1 / 4)

# A complex contrast (say panels 1 and 2 against panels 3 through 6) is
# requested by passing its weights to c_weights, together with
# contrast_type = "allowance", the all-contrasts form of Eq. (2.4) of
# Bryant and Bruvold. That contrast of adjusted means is -0.921, with
# simultaneous limits of -1.199 and -0.643.
ci_c_ancova_bp(adj_means = adj, s_ancova = sqrt(0.01326),
               c_weights = c(0.5, 0.5, -0.25, -0.25, -0.25, -0.25),
               n = 4, num_covariates = 1, df = 14,
               contrast_type = "allowance")

Confidence Intervals for the Population Correlation and Multiple Correlation

Description

Two confidence intervals on the correlation scale share this page, named by the convention that lowercase r is the Pearson product-moment correlation between two variables and capital R is the multiple correlation between an outcome and a set of predictors.

ci_r() forms a confidence interval for the population correlation coefficient \rho. The confidence interval is for the population value \rho; the required input is the corresponding sample value, the observed sample correlation coefficient r. This approach assumes that the two variables on which the correlation is based are bivariate normally distributed (e.g., Hays, 1994, Chapter 14).

ci_R() constructs a confidence interval for the population multiple correlation coefficient \rho = \sqrt{\rho^2} from the sample multiple correlation coefficient (or, equivalently, from the observed F-statistic and degrees of freedom). The interval is obtained by inverting the sampling distribution of the sample R^2 and propagating the limits through the monotone (square root) transform.

The two estimands meet at a single predictor: the multiple correlation from a regression on one predictor is the absolute value of the Pearson correlation between the outcome and that predictor.

Usage

ci_r(r, n, conf_level = 0.95, alpha_lower = NULL, alpha_upper = NULL)

ci_R(
  R = NULL,
  df_1 = NULL,
  df_2 = NULL,
  conf_level = 0.95,
  random_predictors = TRUE,
  F_value = NULL,
  N = NULL,
  p = NULL,
  alpha_lower = NULL,
  alpha_upper = NULL,
  ...
)

Arguments

Details

The Pearson correlation interval (ci_r). This approach will not generally lead to a symmetric confidence interval. The function first transforms r into Z', forms a confidence interval for the population value (i.e., \zeta), and then transforms the confidence limits for \zeta into the scale of the correlation coefficient. The interval requires a sample size of at least 4. The variance of Z' is 1/(n - 3), which is infinite at n = 3; there the interval would be vacuous, covering [-1, 1] regardless of r, and for smaller n the variance is undefined. The function therefore stops with an error when n < 4.

Fixed vs. random predictors (ci_R). The two regression models give different sampling distributions for the sample R^2, and so different confidence intervals on \rho. Under fixed predictors the design matrix is treated as constant in hypothetical replications of the study, and the omnibus F-statistic follows a noncentral F with p and N - p - 1 degrees of freedom and noncentrality \lambda = N \rho^2 / (1 - \rho^2) (Cohen, 1988); the CI on \rho^2 is obtained by inverting that distribution and then taking the square root (see ci_nc_F). Under random predictors the design matrix is itself a draw from a joint multivariate normal distribution and the unconditional sampling distribution of the sample R^2 is given by Lee (1971); the same Lee bisection that ci_R2 uses for the random-predictor CI on \rho^2 is applied here and the limits are mapped to \rho. Gatsonis and Sampson (1989) document the comparison; in the behavioral, educational, and social sciences predictor variables are almost always random, so the default is random_predictors = TRUE. Pass random_predictors = FALSE for designs in which the predictor variables are fixed by design.

Value

ci_r() returns a 3-row data.frame with columns term and value. The term values are "lower_limit" (the lower confidence limit on the population correlation \rho), "r" (the observed sample correlation coefficient), and "upper_limit" (the upper limit on \rho).

ci_R() returns a 3-row data.frame with columns term, value, prob_less, and prob_greater. The rows are ordered "lower_limit", "R" (the sample multiple correlation coefficient supplied by the user, the point estimate), and "upper_limit", so the point estimate sits between its confidence limits. The lower and upper limits are the confidence limits on the population multiple correlation coefficient \rho (square roots of the corresponding limits on \rho^2). The prob_less and prob_greater columns report the achieved lower-tail and upper-tail error probabilities at each limit (they are NA for the "R" estimate row).

Note

The ci_r() confidence interval assumes that the two variables the correlation is based on are bivariate normal. See Hays (1994, Chapter 14) for details.

Author(s)

Ken Kelley kkelley@nd.edu

References

Algina, J. & Olejnik, S. (2000). Determining sample size for accurate estimation of the squared multiple correlation coefficient. Multivariate Behavioral Research, 35, 119–137. doi:10.1207/s15327906mbr3501_5

Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Hillsdale, NJ: Lawrence Erlbaum.

Gatsonis, C., & Sampson, A. R. (1989). Multiple correlation: Exact power and sample size calculations. Psychological Bulletin, 106(3), 516–524.

Hays, W. L. (1994). Statistics (5th ed.). Fort Worth, TX: Harcourt Brace College Publishers.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K. (2008). Sample size planning for the squared multiple correlation coefficient: Accuracy in parameter estimation via narrow confidence intervals. Multivariate Behavioral Research, 43, 524–555. doi:10.1080/00273170802490632

Kelley, K., & Maxwell, S. E. (2003). Sample size for multiple regression: Obtaining regression coefficients that are accurate, not simply significant. Psychological Methods, 8(3), 305–321. doi:10.1037/1082-989X.8.3.305

Lee, Y. S. (1971). Some results on the sampling distribution of the multiple correlation coefficient. Journal of the Royal Statistical Society, Series B, 33(1), 117–130.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on R^2 as a model comparison effect size.)

Smithson, M. (2003). Confidence intervals. Thousand Oaks, CA: Sage Publications.

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

Steiger, J. H., & Fouladi, R. T. (1992). R2: A computer program for interval estimation, power calculations, sample size estimation, and hypothesis testing in multiple regression. Behavior Research Methods, Instruments, & Computers, 24(4), 581–582. doi:10.3758/BF03203611

See Also

ci_R2, ss_aipe_r, ss_power_r, var_r, ss_aipe_R2, convert_r_Z, convert_Z_r, ci_nc_t

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Pearson correlation, from Hays. Suppose n = 100 and r = .35.
ci_r(r = .35, n = 100, conf_level = .95)

# Here is another way to enter the above example.
ci_r(r = .35, n = 100, conf_level = NULL,
     alpha_lower = .025, alpha_upper = .025)

# Here are examples of one-sided confidence intervals.
ci_r(r = .35, n = 100, conf_level = NULL, alpha_lower = 0, alpha_upper = .05)
ci_r(r = .35, n = 100, conf_level = NULL, alpha_lower = .05, alpha_upper = 0)

# Multiple correlation from a five-predictor regression.
ci_R(R = .7071, df_1 = 5, df_2 = 50, conf_level = .95,
     random_predictors = TRUE)

Confidence Interval for the Coefficient of Variation

Description

Computes the noncentral t-based confidence interval for the population coefficient of variation, the standard deviation relative to the mean, so variability can be reported on a scale that is free of the measurement units.

Usage

ci_cv(
  cv = NULL,
  mean = NULL,
  sd = NULL,
  n = NULL,
  data = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  ...
)

Arguments

Details

Uses the noncentral t-distribution to calculate the confidence interval for the population coefficient of variation.

Value

A 4-row data.frame with columns term, value, prob_less, and prob_greater. The rows are ordered so the two point estimates sit between the confidence limits: "lower_limit" (lower confidence limit on the coefficient of variation), "c_of_v" (the sample coefficient of variation), "c_of_v_unbiased" (the unbiased estimator), and "upper_limit" (upper confidence limit). The prob_less and prob_greater columns report the achieved tail probabilities of the noncentral t search at the limit values; they are NA for the point-estimate rows.

Author(s)

Ken Kelley kkelley@nd.edu

References

Johnson, N. L., & Welch, B. L. (1940). Applications of the non-central t-distribution. Biometrika, 31(3–4), 362–389. doi:10.1093/biomet/31.3-4.362

Kelley, K. (2007). Sample size planning for the coefficient of variation from the accuracy in parameter estimation approach. Behavior Research Methods, 39(4), 755–766. doi:10.3758/BF03192966

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3.)

McKay, A. T. (1932). Distribution of the coefficient of variation and the extended t distribution. Journal of the Royal Statistical Society, 95(4), 695–698.

See Also

cv

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

set.seed(113)
N <- 15
X <- rnorm(N, 5, 1)
mean.X <- mean(X)
sd.X <- var(X)^.5

ci_cv(mean = mean.X, sd = sd.X, n = N, alpha_lower = .025,
alpha_upper = .025, conf_level = NULL)
ci_cv(data = X, conf_level = .95)
ci_cv(cv = sd.X / mean.X, n = N, conf_level = .95)

Dunnett's Simultaneous Confidence Intervals Against a Control

Description

Computes Dunnett's (1955, 1964) simultaneous confidence intervals for the a - 1 comparisons of a - 1 treatment means against a single control mean, with family-wise coverage at the specified conf_level. Returns the result in tidy long form.

Usage

ci_dunnett(
  x,
  group = NULL,
  control = NULL,
  alternative = c("two_sided", "less", "greater"),
  conf_level = 0.95
)

Arguments

Details

Critical value. The two-sided Dunnett critical value d_{\alpha, a - 1, \nu}^{(2)} is obtained from the multivariate t distribution with a - 1 dimensions, common correlation 0.5 (the Dunnett correlation under balanced n; the function does not adjust for unequal n), and \nu error degrees of freedom. The function uses the existing cv_dunnett() critical value.

Adjusted p-values. Computed exactly from the same equicorrelated multivariate t distribution. The one common correlation 1/2 admits a one-factor representation, so the probability that all comparisons fall inside (or below) the observed statistic collapses to two nested one-dimensional integrals, evaluated by quadrature. The adjusted p-value is one minus that probability. The computation is deterministic (no Monte Carlo) and needs no additional package.

Value

A data.frame with one row per non-control level. Columns: contrast, mean_difference, se, t_statistic, lower_limit, upper_limit, p_adjusted.

Author(s)

Ken Kelley kkelley@nd.edu

References

Dunnett, C. W. (1955). A multiple comparison procedure for comparing several treatments with a control. Journal of the American Statistical Association, 50(272), 1096–1121.

Dunnett, C. W. (1964). New tables for multiple comparisons with a control. Biometrics, 20(3), 482–491.

Hsu, J. C. (1996). Multiple comparisons: Theory and methods. Chapman & Hall.

See Also

cv_dunnett, ci_tukey_kramer, ci_scheffe

Other hypothesis tests: adjusted_means(), ancova(), anova_within(), ci_scheffe(), ci_tukey_kramer(), compare_cov_structures(), contrast_test(), correlations_test(), equivalence_r(), equivalence_smd(), factorial_anova(), manova_split_plot(), mauchly_test(), mixed_anova(), obrien_test(), pairwise_within(), randomization_test(), randomization_test_paired(), regions_of_significance(), simple_effects_AB(), summary_t_test(), welch_t()

Examples

# 1. Compare the SSRI and placebo arms of the depression_bdi study
#    against the wait list control:
fit <- lm(bdi_post ~ condition, data = depression_bdi)
ci_dunnett(fit, control = "wait_list")

# 2. One-sided: a treatment that works pulls the posttest BDI down,
#    so the directional alternative is "less":
ci_dunnett(fit, control = "wait_list", alternative = "less")

Confidence Interval on the Largest Eigenvalue of a Sample Covariance Matrix

Description

Computes an asymptotic confidence interval on the largest population eigenvalue \lambda_1 of a population covariance matrix, given a sample covariance matrix from n observations on p variables under multivariate normality. Useful in principal-components analysis and dimension-reduction settings to gauge whether the largest eigenvalue is well-separated from the second.

Usage

ci_eigenvalue(cov_matrix, n = NULL, conf_level = 0.95, k = 1)

Arguments

Details

Asymptotic distribution. Under multivariate normality with eigenvalues \lambda_1 > \lambda_2 \ge \cdots \ge \lambda_p, the sample eigenvalues \hat\lambda_j are asymptotically independent and approximately normal with mean \lambda_j and variance 2 \lambda_j^2 / (n - 1) when the eigenvalues are simple (well-separated) (Anderson, 2003, Theorem 13.3.1; Muirhead, 1982, Section 9.7). The asymptotic CI is therefore

\hat\lambda_j \cdot \exp\!\left(\pm z_{1 - \alpha/2} \sqrt{\frac{2}{n - 1}}\right),

on the multiplicative scale (equivalently, a Wald CI on \log \lambda_j with variance 2/(n - 1)). The log scale is the natural variance-stabilizing transformation for an eigenvalue.

Caveats. The asymptotic CI assumes well-separated population eigenvalues. When the largest two eigenvalues are close, the sample eigenvalue exhibits a "repulsion" phenomenon and the CI is biased (typically too narrow). Diagnostic: if \hat\lambda_1 / \hat\lambda_2 is close to 1, the asymptotic CI should not be relied upon; a bootstrap is preferable.

Value

A 3-row data.frame with rows ordered "lower_limit", "eigenvalue" (the sample eigenvalue point estimate), and "upper_limit", so the point estimate sits between its confidence limits.

Author(s)

Ken Kelley kkelley@nd.edu

References

Anderson, T. W. (2003). An introduction to multivariate statistical analysis (4th ed.). Wiley. (See Chapter 13.)

Muirhead, R. J. (1982). Aspects of multivariate statistical theory. Wiley. (See Section 9.7.)

See Also

prcomp, eigen

Other multivariate and latent variable methods: average_variance_extracted(), bifactor_indices(), cfa_1(), cfa_2(), cfa_k(), common_method_marker(), common_method_single_factor(), dmacs(), ecvi(), htmt(), irt_grm(), irt_information(), measurement_alignment(), measurement_invariance(), procrustes_phi(), simple_structure()

Examples

# 1. From a data frame:
set.seed(113)
X <- data.frame(matrix(rnorm(200), nrow = 50))
ci_eigenvalue(X, k = 1)

# 2. From an explicit covariance matrix:
S <- cov(X)
ci_eigenvalue(S, n = nrow(X), k = 1)

Confidence Interval for Eta Squared (Effect Size for ANOVA)

Description

Computes the point estimate and an exact, noncentrality-based confidence interval for the population eta squared (\eta^2), the proportion of variance in the dependent variable accounted for by a fixed effect. Accepts either the raw ANOVA summary (F, effect df, error df, total N) or a fitted model object. Supports both between-subjects designs (aov / lm) and within-subjects / mixed designs (aovlist fits with an Error() term in the formula). For factorial and within-subjects designs the function returns one row per effect with the CI for partial \eta^2 computed against that effect's own error stratum.

Usage

ci_eta_squared(
  object = NULL,
  F_value = NULL,
  df_effect = NULL,
  df_error = NULL,
  N = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL
)

Arguments

Details

Point estimate. \hat{\eta}^2 = df_{\text{effect}} \cdot F / (df_{\text{effect}} \cdot F + df_{\text{error}}), which equals \mathit{SS}_{\text{effect}}/(\mathit{SS}_{\text{effect}} + \mathit{SS}_{\text{error}}). In a one-way ANOVA this is also \mathit{SS}_{\text{effect}}/\mathit{SS}_{\text{total}}. In a factorial design the same expression gives the per-effect partial \eta^2.

Confidence interval. The CI is constructed by Steiger's (2004) confidence interval transformation principle: a CI for the noncentrality parameter \lambda of the F distribution is obtained (via ci_nc_F) and then mapped through

\eta^2_{\text{bound}} = \frac{\lambda_{\text{bound}}}{\lambda_{\text{bound}} + N}.

This is the same transformation used by ci_pvaf and ci_omega_squared; the three functions share CI machinery and differ only in their sample point estimators. When the lower CI on \lambda is not identified (i.e., the observed F is below the one-sided critical value), the lower limit on \eta^2 is set to 0.

Designs supported.

• Between-subjects ANOVA: fitted aov/lm, or raw F/df/N.

• Within-subjects or mixed ANOVA: fitted aovlist. Each effect's CI is built from its own stratum's F and residual df; N is the total number of observations across all strata. The reported stratum column identifies which error term each row used.

Sums of squares in factorial designs. When a fitted model is supplied, F-values are read from anova() (single-stratum) or summary() (multi-stratum), both of which use Type I (sequential) sums of squares in base R. For balanced designs Types I, II, and III agree; for unbalanced designs they differ. If Type II or III F-values are required, compute them with e.g.\ car::Anova(object, type = 3) and pass the relevant F, degrees of freedom, and N into the raw-argument interface.

Value

A data.frame with one row per effect. Single-stratum fits and the raw interface return columns effect, eta_squared, lower_limit, upper_limit, F_value, df_effect, df_error, N. aovlist (within-subjects / mixed) fits additionally include a stratum column. With the raw-argument interface effect is "overall".

Author(s)

Ken Kelley kkelley@nd.edu

References

Fleishman, A. I. (1980). Confidence intervals for correlation ratios. Educational and Psychological Measurement, 40(3), 659–670.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Preacher, K. J. (2012). On effect size. Psychological Methods, 17, 137–152. doi:10.1037/a0028086

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on \eta^2, Chapter 7 on factorial designs, and Chapter 11 on generalized \eta^2 for within-subjects designs.)

Smithson, M. (2001). Correct confidence intervals for various regression effect sizes and parameters: The importance of noncentral distributions in computing intervals. Educational and Psychological Measurement, 61, 605–632. doi:10.1177/00131640121971392

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

eta_squared, ci_eta_squared_partial, ci_omega_squared, ci_pvaf, ci_nc_F

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# 1. Raw-argument interface. Bargman's (1970) example.
ci_eta_squared(F_value = 11.221, df_effect = 4, df_error = 50, N = 55)

# Same example with a 90% confidence interval.
ci_eta_squared(
  F_value = 11.221, df_effect = 4, df_error = 50, N = 55,
  conf_level = 0.90
)

# 2. One way ANOVA from a fitted model: mean IQ gain differs across
#        the six grades of the pygmalion data (N = 310).
fit_one <- aov(iq_gain ~ factor(grade), data = pygmalion)
ci_eta_squared(fit_one)

# 3. Two-factor ANOVA: partial eta squared per effect for the
#        manipulated expectancy treatment and the measured grade
#        classification (pygmalion data, N = 310). The treatment by
#        grade interaction is weak here (F = 1.19), so the additive
#        model is used.
fit_additive <- aov(iq_8 ~ treatment + factor(grade), data = pygmalion)
ci_eta_squared(fit_additive)

# 4. Within-subjects ANOVA. CI computed against the within-subjects
#        error stratum (the row reports which one via 'stratum').
set.seed(113)
n <- 20
rm_data <- data.frame(
  subject = factor(rep(seq_len(n), each = 3)),
  time    = factor(rep(c("Pre", "Mid", "Post"), n),
                   levels = c("Pre", "Mid", "Post")),
  y       = rnorm(n, sd = 1.5)[rep(seq_len(n), each = 3)] +
            0.7 * rep(1:3, n) + rnorm(n * 3, sd = 1.2)
)
fit_rm <- aov(y ~ time + Error(subject/time), data = rm_data)
ci_eta_squared(fit_rm)

Confidence Interval for Generalized Eta Squared (Approximate)

Description

Returns the point estimate of generalized eta squared (\eta^2_G; Olejnik & Algina, 2003) along with an optional confidence interval computed by one of two approximate methods. The default is to return only the point estimate (method = "none"), because both available CI methods are approximations whose coverage properties have not been broadly validated for this estimand and warrant independent evaluation before being used in substantive inference.

Usage

ci_eta_squared_generalized(
  object = NULL,
  observed = NULL,
  SS_effect = NULL,
  SS_observed = NULL,
  SS_error = NULL,
  F_effect = NULL,
  df_effect = NULL,
  F_observed = NULL,
  df_observed = NULL,
  df_error = NULL,
  N = NULL,
  method = c("none", "parametric", "bootstrap"),
  B = 10000L,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  seed = NULL
)

Arguments

Details

Why CI = "none" is the default. Confidence interval construction for \eta^2_G is not as settled as for partial \eta^2 or \omega^2, because the denominator mixes sums of squares from heterogeneous sources (the focal effect, one or more measured factors, and the error term). No noncentral F transformation maps the population noncentrality parameter directly to \eta^2_G. Both methods below are approximations and are exposed for exploration rather than as defaults.

method = "parametric". The function first obtains a confidence interval for the population noncentrality parameter \lambda of the focal effect's F-test via ci_nc_F. The NCP bounds are mapped through the partial-\eta^2 transformation \eta^2_{p,\text{bound}} = \lambda_{\text{bound}}/(\lambda_{\text{bound}} + N) (matching the convention used by ci_pvaf and ci_omega_squared), and then re-expressed as \eta^2_G bounds via

\eta^2_{G,\text{bound}} = \frac{r_{\text{bound}}}{r_{\text{bound}} + r_{\text{obs}} + 1},

where r_{\text{bound}} = \eta^2_{p,\text{bound}}/(1-\eta^2_{p,\text{bound}}) and r_{\text{obs}} = \sum \mathit{SS}_{\text{measured}}/\mathit{SS}_{\text{error}}. This treats the observed-factor sums of squares as fixed at their sample values, so the interval inherits whatever sampling variability those contribute. It has not been validated for coverage and should be treated as preliminary.

method = "bootstrap". A residual bootstrap from the fitted model: the resampling unit is a residual, drawn nonparametrically from the model's own residuals rather than from a fitted distribution. For each of the B replications, the response is regenerated as \hat{y}_i + \varepsilon^*_i where \varepsilon^* is sampled with replacement from the model's residuals; the model is refit; \eta^2_G is recomputed; and the percentile interval, the empirical quantiles of the B bootstrap values (Efron & Tibshirani, 1993), is reported. The percentile interval is the only bootstrap interval offered; there is no bias-corrected and accelerated (BCa) variant. Replicates whose refit fails are dropped, and the interval is computed from the replications that return a value. Bootstrap results vary from run to run; supply seed for reproducibility. Requires a single-stratum aov/lm fit. Not yet supported for aovlist (multi-stratum / within-subjects) fits, since the bootstrap needs to respect the subject-level correlation structure, which naive residual resampling does not. Use method = "parametric" for aovlist fits. Coverage has not been broadly validated for this estimand.

Within-subjects designs (aovlist). For multi-stratum fits the parametric CI uses each focal effect's stratum-specific F test and degrees of freedom. The denominator ratio is adjusted to reflect the full set of error strata: the implied \mathit{SS}_{\text{error}} for the focal effect is its own stratum's residual SS, while the r_{\text{obs}} term includes both measured-factor SS and the residual SS of all other strata. This generalizes the partial-\eta^2 CI machinery to the Bakeman (2005) denominator.

Value

A data.frame with one row per focal effect and the columns effect, eta_squared_generalized, lower_limit, upper_limit, and method. For aovlist fits a stratum column is also present. When method = "none", the limit columns contain NA.

Author(s)

Ken Kelley kkelley@nd.edu

References

Algina, J., Keselman, H. J., & Penfield, R. D. (2005). An alternative to Cohen's standardized mean difference effect size: A robust parameter and confidence interval in the two independent groups case. Psychological Methods, 10(3), 317–328. doi:10.1037/1082-989X.10.3.317

Bakeman, R. (2005). Recommended effect size statistics for repeated measures designs. Behavior Research Methods, 37(3), 379–384. doi:10.3758/BF03192707

Efron, B., & Tibshirani, R. J. (1993). An introduction to the bootstrap. New York, NY: Chapman & Hall/CRC.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Preacher, K. J. (2012). On effect size. Psychological Methods, 17, 137–152. doi:10.1037/a0028086

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on \eta^2, Chapter 7 on factorial designs, and Chapter 11 on generalized \eta^2 for within-subjects designs.)

Olejnik, S., & Algina, J. (2003). Generalized eta and omega squared statistics: Measures of effect size for some common research designs. Psychological Methods, 8(4), 434–447. doi:10.1037/1082-989X.8.4.434

Smithson, M. (2001). Correct confidence intervals for various regression effect sizes and parameters: The importance of noncentral distributions in computing intervals. Educational and Psychological Measurement, 61, 605–632. doi:10.1177/00131640121971392

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

eta_squared_generalized, ci_eta_squared

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# The pygmalion expectancy experiment: treatment is manipulated, while
# grade is a measured classification, so grade belongs in the denominator.
pyg <- pygmalion
pyg$grade <- factor(pyg$grade)
fit <- aov(iq_8 ~ treatment * grade, data = pyg)

# The default returns the point estimate and leaves the limits NA, because
# both interval methods are approximations.
ci_eta_squared_generalized(fit, observed = "grade")

# The parametric approximation maps a noncentral F interval for the focal
# effect through the observed sums of squares. Warnings mark the method as
# preliminary and report that the interaction's lower limit is clamped at
# 0; treat the limits accordingly.
ci_eta_squared_generalized(fit, observed = "grade", method = "parametric")

# The third option is a residual bootstrap, which refits the model once per
# replication. B = 1000 is the smallest count the function accepts and is
# what keeps this page quick; a reported interval deserves B = 10000 or
# more. The seed makes the interval reproducible and leaves the random
# number stream of the surrounding session as it was.
ci_eta_squared_generalized(fit, observed = "grade",
                           method = "bootstrap", B = 1000, seed = 113)

# Within-subjects ANOVA. The parametric CI uses each effect's own stratum.
set.seed(113)
n <- 20
rm_data <- data.frame(
  subject = factor(rep(seq_len(n), each = 3)),
  time    = factor(rep(c("Pre", "Mid", "Post"), n),
                   levels = c("Pre", "Mid", "Post")),
  y       = rnorm(n, sd = 1.5)[rep(seq_len(n), each = 3)] +
            0.7 * rep(1:3, n) + rnorm(n * 3, sd = 1.2)
)
fit_rm <- aov(y ~ time + Error(subject/time), data = rm_data)
ci_eta_squared_generalized(fit_rm, method = "parametric")

Confidence Interval for Partial Eta Squared (Effect Size for ANOVA)

Description

Computes the point estimate and an exact, noncentrality-based confidence interval for the population partial eta squared (\eta^2_p). Accepts either the raw ANOVA summary (F, effect df, error df, total N) or a fitted aov/lm/aovlist object, in which case the function returns one row per effect (with stratum identification for within-subjects fits).

Usage

ci_eta_squared_partial(
  object = NULL,
  F_value = NULL,
  df_effect = NULL,
  df_error = NULL,
  N = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL
)

Arguments

Details

This is the explicitly-named counterpart of ci_eta_squared. The two share point-estimate and CI machinery: in a one-way ANOVA partial \eta^2 coincides with \eta^2; in a factorial or within-subjects ANOVA both functions return the per-effect partial value computed against that effect's own error stratum. Use ci_eta_squared_partial when you want the function name to make the partial interpretation explicit.

Value

A data.frame with one row per effect. Single-stratum fits and the raw interface return columns effect, eta_squared_partial, lower_limit, upper_limit, F_value, df_effect, df_error, N. aovlist fits additionally include a stratum column. With the raw-argument interface effect is "overall".

Author(s)

Ken Kelley kkelley@nd.edu

References

Cohen, J. (1973). Eta-squared and partial eta-squared in fixed factor ANOVA designs. Educational and Psychological Measurement, 33(1), 107–112.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Preacher, K. J. (2012). On effect size. Psychological Methods, 17, 137–152. doi:10.1037/a0028086

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on \eta^2, Chapter 7 on factorial designs, and Chapter 11 on generalized \eta^2 for within-subjects designs.)

Smithson, M. (2001). Correct confidence intervals for various regression effect sizes and parameters: The importance of noncentral distributions in computing intervals. Educational and Psychological Measurement, 61, 605–632. doi:10.1177/00131640121971392

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

eta_squared_partial, ci_eta_squared

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Raw-argument interface.
ci_eta_squared_partial(F_value = 11.221, df_effect = 4,
                       df_error = 50, N = 55)

# Two-factor ANOVA: per-effect partial eta squared with CI for the
# manipulated expectancy treatment and the measured grade
# classification (pygmalion data, N = 310). The treatment by grade
# interaction is weak here (F = 1.19), so the additive model is used.
fit <- aov(iq_8 ~ treatment + factor(grade), data = pygmalion)
ci_eta_squared_partial(fit)

# Within-subjects ANOVA.
set.seed(113)
n <- 20
rm_data <- data.frame(
  subject = factor(rep(seq_len(n), each = 3)),
  time    = factor(rep(c("Pre", "Mid", "Post"), n),
                   levels = c("Pre", "Mid", "Post")),
  y       = rnorm(n, sd = 1.5)[rep(seq_len(n), each = 3)] +
            0.7 * rep(1:3, n) + rnorm(n * 3, sd = 1.2)
)
fit_rm <- aov(y ~ time + Error(subject/time), data = rm_data)
ci_eta_squared_partial(fit_rm)

Provides Games–Howell Simultaneous Confidence Intervals for All Pairwise Comparisons Without Assuming Homogeneity of Variance

Description

Provides Games–Howell Simultaneous Confidence Intervals for All Pairwise Comparisons Without Assuming Homogeneity of Variance

Usage

ci_games_howell(x, group = NULL, conf_level = 0.95)

Arguments

Details

Tukey's HSD and the Kramer modification for unequal n (ci_tukey_kramer) both pool the within-group variances into \mathit{MS}_W, so both assume homogeneity of variance. Neither is robust when that assumption fails. The Games–Howell procedure drops the assumption: it uses a separate error term for each pair and a Welch–Satterthwaite degrees of freedom for each pair, then takes its critical value from the studentized range.

For groups g and h, the standard error of the difference uses only those two groups' variances, and the degrees of freedom are

\mathit{df} = \frac{(s_g^2/n_g + s_h^2/n_h)^2}{s_g^4/[n_g^2(n_g-1)] + s_h^4/[n_h^2(n_h-1)]},

the same Welch–Satterthwaite expression that base R's t.test uses by default for two groups. A pair is declared different when the observed t exceeds q/\sqrt{2}, with q the studentized range critical value (cv_tukey_hsd) at the pair's degrees of freedom, so the interval for the difference of means is (\bar Y_g - \bar Y_h) \pm q_{\alpha;a,\mathit{df}}\sqrt{(s_g^2/n_g + s_h^2/n_h)/2}. Maxwell, Delaney, and Kelley (2027, Chapter 5) develop this as one of the two modifications of Tukey's HSD for heterogeneous variances (their Equations 5.13 and 5.14).

When to use it. Reach for Games–Howell when the group variances are not interchangeable and the design is between subjects. It is the heterogeneity-robust counterpart of ci_tukey_kramer and, like it, controls the family-wise error rate across all a(a-1)/2 pairs. It handles unequal n as a matter of course, so it does not need a separate unequal-n variant.

When something else is better. Dunnett (1980) found that Games–Howell becomes slightly liberal (the family-wise error rate runs somewhat above the nominal level) when the samples are small. Maxwell, Delaney, and Kelley (2027, Chapter 5) therefore recommend Games–Howell for larger samples and Dunnett's T3, which takes its critical value from the studentized maximum modulus (cv_smm) rather than the studentized range, when the groups have fewer than roughly 50 observations each. When the variances are in fact homogeneous, use ci_tukey_kramer instead: it pools the variances, so it has more error degrees of freedom and more power. When only treatments are compared to a single control, use ci_dunnett.

With a = 2 groups the procedure is exactly Welch's t test: the interval and the p-value equal those from t.test(..., var.equal = FALSE), because q_{\alpha;2,\mathit{df}} = \sqrt{2}\,t_{1-\alpha/2,\mathit{df}}.

Value

A data.frame with one row per pairwise comparison and columns contrast, mean_difference, se, df, q_statistic, lower_limit, upper_limit, and p_adjusted. The df column is the Welch–Satterthwaite degrees of freedom for that pair, which is why it varies from row to row. The table prints through the dmar_tbl display layer and works with tidy and glance (see dmar_tidiers).

Author(s)

Ken Kelley kkelley@nd.edu

References

Games, P. A., & Howell, J. F. (1976). Pairwise multiple comparison procedures with unequal n's and/or variances: A Monte Carlo study. Journal of Educational Statistics, 1(2), 113–125. doi:10.2307/1164979

Dunnett, C. W. (1980). Pairwise multiple comparisons in the unequal variance case. Journal of the American Statistical Association, 75(372), 796–800. doi:10.2307/2287161

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 5 on the multiple-comparisons problem, where the modifications of Tukey's HSD for unequal n and unequal variances are developed.)

See Also

ci_tukey_kramer for the homogeneity-assuming counterpart, ci_dunnett for many-to-one comparisons, ci_scheffe for arbitrary contrasts, cv_tukey_hsd for the critical value it uses, and dmar_tidiers for the tidy methods.

Examples

# The raw drinks_per_week outcome of the drinks_trial data: the
# Standard arm's variance is more than three times that of either CRA
# arm, so the pooled error term Tukey's HSD relies on is questionable.
ci_games_howell(drinks_trial$drinks_per_week, drinks_trial$treatment)

# A fitted one-way model may be passed instead of the two vectors.
ci_games_howell(aov(drinks_per_week ~ treatment, data = drinks_trial))

# With two groups the procedure is Welch's t test, so the limits agree;
# the trial's two enrollment cohorts give a two-group comparison.
ci_games_howell(drinks_trial$drinks_per_week, drinks_trial$cohort)$lower_limit
-t.test(drinks_per_week ~ cohort, data = drinks_trial)$conf.int[2]

Confidence Interval for the Squared Mahalanobis Distance

Description

Computes the squared Mahalanobis distance D^2 together with an exact confidence interval for the population squared distance \Delta^2, obtained by inverting Hotelling's T^2 statistic through its (noncentral) F-distribution as in Reiser (2001). Both the one-sample setting (mean vector against a hypothesized population mean) and the two-sample setting (between-groups distance from discriminant analysis or multivariate group comparison) are supported, and either raw data or a pre-computed D^2 with sample sizes can be supplied.

Usage

ci_mahalanobis(
  D2 = NULL,
  group_1 = NULL,
  group_2 = NULL,
  mu_0 = NULL,
  n_1 = NULL,
  n_2 = NULL,
  p = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  ...
)

Arguments

Details

Definition. For a p-vector \mathbf{x} drawn from a multivariate normal with mean \boldsymbol{\mu} and covariance \boldsymbol{\Sigma}, Mahalanobis's (1936) squared distance from a reference vector \boldsymbol{\mu}_0 is

\Delta^2 = (\boldsymbol{\mu} - \boldsymbol{\mu}_0)^\top \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu} - \boldsymbol{\mu}_0).

In the two-sample case the population distance between groups is \Delta^2 = (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^\top \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2), assuming a common covariance. The corresponding sample estimates plug the sample means and the sample (or pooled) covariance into the same quadratic form.

Link to Hotelling's T^2. Hotelling's (1931) T^2 statistic is T^2 = n D^2 (one sample) or T^2 = \{n_1 n_2 / (n_1 + n_2)\} D^2 (two samples). Under multivariate normality

\frac{n_1 + n_2 - p - 1}{(n_1 + n_2 - 2)\,p}\,T^2 \sim F'\!\left(p,\, n_1 + n_2 - p - 1,\; \lambda = \frac{n_1 n_2}{n_1 + n_2}\,\Delta^2\right)

in the two-sample case, and analogously \{(n - p)/[(n-1)p]\}\,T^2 \sim F'(p, n-p, n\Delta^2) in the one-sample case (see Anderson, 2003, Section 5.2).

Confidence interval. The CI on \Delta^2 is obtained by inverting these distributional results (Reiser, 2001): a CI on the noncentrality parameter \lambda is constructed via ci_nc_F and then mapped back to \Delta^2 by \Delta^2 = \lambda\,(n_1 + n_2)/(n_1 n_2) (two sample) or \Delta^2 = \lambda / n (one sample). When the observed F is below the lower-tail critical value of the central F-distribution at the requested confidence level, the lower CI on \lambda (and hence on \Delta^2) is clamped to zero, in keeping with the ci_nc_F convention.

Bias. The plug-in estimator D^2 is upward biased for \Delta^2; the CI from this function is exact for \Delta^2 under multivariate normality and reflects the bias structure correctly, but the point estimate reported is the standard plug-in D^2.

Value

A one-row data.frame with columns sample_type ("one-sample" or "two-sample"), D2 (point estimate of the squared distance), lower_limit and upper_limit (the confidence limits on the population squared distance \Delta^2), F_value, df_1, df_2, n_1, n_2 (NA in one-sample mode), and p.

Author(s)

Ken Kelley kkelley@nd.edu

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis (4th ed.). Wiley.

Hotelling, H. (1931). The generalization of Student's ratio. The Annals of Mathematical Statistics, 2(3), 360–378.

Mahalanobis, P. C. (1936). On the generalized distance in statistics. Proceedings of the National Institute of Sciences of India, 2(1), 49–55.

Reiser, B. (2001). Confidence intervals for the Mahalanobis distance. Communications in Statistics–Simulation and Computation, 30(1), 37–45. doi:10.1081/SAC-100001856

Steiger, J. H., & Fouladi, R. T. (1997). Noncentrality interval estimation and the evaluation of statistical methods. In L. L. Harlow, S. A. Mulaik, & J. H. Steiger (Eds.), What if there were no significance tests? (pp. 221–257). Mahwah, NJ: Lawrence Erlbaum.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

See Also

ci_nc_F, ci_smd, ci_R2

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Two-sample distance between the two schools of the Holzinger and
#     Swineford (1939) study on a four-test cognitive battery.
battery <- c("t1_visual_perception", "t2_cubes", "t4_lozenges",
             "t6_paragraph_comprehension")
g1 <- as.matrix(holzinger_swineford[
  holzinger_swineford$school == "Grant-White", battery])
g2 <- as.matrix(holzinger_swineford[
  holzinger_swineford$school == "Pasteur", battery])
ci_mahalanobis(group_1 = g1, group_2 = g2)

# One-sample distance: how far is the Grant-White centroid from a
#     reference vector of (29, 24, 18, 9)?
ci_mahalanobis(group_1 = g1, mu_0 = c(29, 24, 18, 9))

# Pre-computed D^2 (no raw data needed): the two-school distance,
#     reproduced from reported summaries.
ci_mahalanobis(D2 = 0.608, n_1 = 145, n_2 = 156, p = 4)

Confidence Interval for the Noncentrality Parameter of a Noncentral F Distribution

Description

Finds the noncentrality parameters of a noncentral F-distribution that bracket an observed F-value with the requested tail probabilities, giving a confidence interval on the population noncentrality parameter. Together with ci_nc_t and ci_nc_chisq, this is one of the low-level noncentral distribution workhorses on which the ci_* confidence interval functions (e.g., ci_pvaf, ci_snr, ci_R2) are built; most analyses reach it through those functions rather than calling it directly. The function was conf_limits_ncf() in earlier builds of DMAR and is conf.limits.ncf() in MBESS; it is named into the ci_* family because a confidence interval is what it computes.

Usage

ci_nc_F(
  F_value = NULL,
  conf_level = 0.95,
  df_1 = NULL,
  df_2 = NULL,
  alpha_lower = NULL,
  alpha_upper = NULL,
  tol = 1e-09,
  verbose = TRUE,
  ...
)

Arguments

Details

Each confidence limit is the noncentrality parameter \lambda \ge 0 of a noncentral F-distribution with df_1 and df_2 degrees of freedom whose appropriate tail at the observed F_value contains the requested probability:

• the lower limit satisfies P(F \ge \mathtt{F\_value}) = \mathtt{alpha\_lower};

• the upper limit satisfies P(F \le \mathtt{F\_value}) = \mathtt{alpha\_upper}.

The two conditions run in opposite directions in \lambda: the lower-tail probability P(F \le \mathtt{F\_value}) is continuous and strictly decreasing in the noncentrality parameter, so the upper-tail probability P(F \ge \mathtt{F\_value}) is continuous and strictly increasing in it. The lower limit is the \lambda at which the upper tail has grown to alpha_lower, and the upper limit is the \lambda at which the lower tail has shrunk to alpha_upper. Each is therefore the unique non-negative root of a one-dimensional equation, and both are located with uniroot on the decreasing lower-tail scale; extendInt is used to widen the search bracket if needed.

Because the noncentrality parameter is bounded below by zero, the lower limit is set to zero whenever the observed F_value is smaller than the alpha_lower critical value of the central F-distribution (i.e., the data is consistent with \lambda = 0 at the requested confidence level). A warning is issued in that case, and the achieved probabilities reported in the output reflect the actual values at \lambda = 0 rather than the requested alpha_lower. The warning carries the condition class dmar_nc_F_clamp, so a caller that inverts the noncentral F repeatedly can muffle or deduplicate it by class.

Value

A data.frame with one row per confidence limit and the columns:

Author(s)

Ken Kelley kkelley@nd.edu

References

Cumming, G., & Finch, S. (2001). A primer on the understanding, use, and calculation of confidence intervals that are based on central and noncentral distributions. Educational and Psychological Measurement, 61(4), 532–574. doi:10.1177/0013164401614002

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

See Also

ss_aipe_R2, ci_R2, ci_nc_t, ci_nc_chisq, stats::pf(), stats::qf(), uniroot

Other noncentral distribution confidence intervals: ci_nc_chisq(), ci_nc_t()

Examples

ci_nc_F(F_value = 5, conf_level = .95, df_1 = 5, df_2 = 100)

# A one-sided (upper) confidence interval.
ci_nc_F(F_value = 5, conf_level = NULL, df_1 = 5, df_2 = 100,
                alpha_lower = 0, alpha_upper = .05)

Confidence Interval for the Noncentrality Parameter of a Noncentral Chi Square Distribution

Description

Finds the noncentrality parameters of a noncentral chi square distribution that bracket an observed chi square value with the requested tail probabilities, giving a confidence interval on the population noncentrality parameter. Together with ci_nc_t and ci_nc_F, this is one of the low-level noncentral distribution workhorses on which the ci_* confidence interval functions are built; most analyses reach it through those functions rather than calling it directly. The function was conf_limits_nc_chisq() in earlier builds of DMAR and is conf.limits.nc.chisq() in MBESS; it is named into the ci_* family because a confidence interval is what it computes.

Usage

ci_nc_chisq(
  chi_square = NULL,
  conf_level = 0.95,
  df = NULL,
  alpha_lower = NULL,
  alpha_upper = NULL,
  tol = 1e-09,
  verbose = TRUE,
  ...
)

Arguments

Details

Each confidence limit is the noncentrality parameter \lambda \ge 0 of a noncentral chi square distribution with df degrees of freedom whose appropriate tail at the observed chi_square contains the requested probability:

• the lower limit satisfies P(X \ge \mathtt{chi\_square}) = \mathtt{alpha\_lower};

• the upper limit satisfies P(X \le \mathtt{chi\_square}) = \mathtt{alpha\_upper}.

The two conditions run in opposite directions in \lambda: the lower-tail probability P(X \le \mathtt{chi\_square}) is continuous and strictly decreasing in the noncentrality parameter, so the upper-tail probability P(X \ge \mathtt{chi\_square}) is continuous and strictly increasing in it. The lower limit is the \lambda at which the upper tail has grown to alpha_lower, and the upper limit is the \lambda at which the lower tail has shrunk to alpha_upper. Each is therefore the unique non-negative root of a one-dimensional equation, and both are located with uniroot on the decreasing lower-tail scale; extendInt is used to widen the search bracket if needed.

Because the noncentrality parameter is bounded below by zero, the lower limit is set to zero whenever the observed chi_square is smaller than the alpha_lower critical value of the central chi square distribution (i.e., the data is consistent with \lambda = 0 at the requested confidence level). A warning is issued in that case, and the achieved probabilities reported in the output reflect the actual values at \lambda = 0.

Symmetrically, when the observed chi_square is so small that even at \lambda = 0 the lower-tail probability is already at or below alpha_upper, no \lambda \ge 0 places as much as alpha_upper mass at or below chi_square; the upper limit is undefined and is returned as NA, with a warning.

Value

A data.frame with one row per confidence limit and the columns:

Author(s)

Ken Kelley kkelley@nd.edu

References

Cumming, G., & Finch, S. (2001). A primer on the understanding, use, and calculation of confidence intervals that are based on central and noncentral distributions. Educational and Psychological Measurement, 61(4), 532–574. doi:10.1177/0013164401614002

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

See Also

ci_nc_t, ci_nc_F, stats::pchisq(), stats::qchisq(), uniroot

Other noncentral distribution confidence intervals: ci_nc_F(), ci_nc_t()

Examples

# A typical call to the function.
ci_nc_chisq(chi_square = 30, conf_level = .95, df = 15)

# A one-sided (upper) confidence interval.
ci_nc_chisq(chi_square = 30, alpha_lower = 0, alpha_upper = .05,
                     conf_level = NULL, df = 15)

Confidence Interval for the Noncentrality Parameter of a Noncentral t Distribution

Description

Finds the noncentrality parameters of a noncentral t-distribution that bracket an observed t-value with the requested tail probabilities, giving a confidence interval on the population noncentrality parameter. Together with ci_nc_F and ci_nc_chisq, this is one of the low-level noncentral distribution workhorses on which the ci_* confidence interval functions (e.g., ci_smd, ci_smd_c, ci_cv) are built; most analyses reach it through those functions rather than calling it directly. The function was conf_limits_nct() in earlier builds of DMAR and is conf.limits.nct() in MBESS; it is named into the ci_* family because a confidence interval is what it computes.

Usage

ci_nc_t(
  ncp,
  df,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  t_value,
  tol = 1e-09,
  verbose = TRUE,
  ...
)

Arguments

Details

Each confidence limit is the noncentrality parameter of a noncentral t-distribution with df degrees of freedom whose appropriate tail at the observed ncp contains the requested probability:

• the lower limit satisfies P(T \ge \mathrm{ncp}) = \mathtt{alpha\_lower};

• the upper limit satisfies P(T \le \mathrm{ncp}) = \mathtt{alpha\_upper}.

Each tail probability is continuous and strictly monotone in the noncentrality parameter, so each limit is the unique root of a one-dimensional equation. The roots are located with uniroot starting from a bracket centered on ncp with half-width scaled by the asymptotic standard error of the noncentrality estimator; extendInt is used to widen the bracket if needed.

This function is especially useful for forming confidence intervals around standardized mean differences (Cohen's d, Glass's g, Hedges' g), standardized regression coefficients, and coefficients of variation.

Value

A data.frame with one row per confidence limit and the columns:

Warning

As of R 4.0.0, the largest ncp that R can accurately handle is 37.62.

Author(s)

Ken Kelley kkelley@nd.edu

References

Cumming, G., & Finch, S. (2001). A primer on the understanding, use, and calculation of confidence intervals that are based on central and noncentral distributions. Educational and Psychological Measurement, 61(4), 532–574. doi:10.1177/0013164401614002

Kelley, K. (2005). The effects of nonnormal distributions on confidence intervals around the standardized mean difference: Bootstrap and parametric confidence intervals, Educational and Psychological Measurement, 65, 51–69. doi:10.1177/0013164404264850

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Steiger, J. H., & Fouladi, R. T. (1997). Noncentrality interval estimation and the evaluation of statistical methods. In L. L. Harlow, S. A. Mulaik, & J. H. Steiger (Eds.), What if there were no significance tests? (pp. 221–257). Mahwah, NJ: Lawrence Erlbaum.

See Also

stats::pt(), stats::qt(), uniroot, ci_smd, ci_smd_c, ci_nc_F, ci_nc_chisq

Other noncentral distribution confidence intervals: ci_nc_F(), ci_nc_chisq()

Examples

# Suppose observed t-value based on 'df'=126 is 2.83. Finding the lower
# and upper critical values for the population noncentrality parameter
# with a symmetric confidence interval with 95\% confidence is given as:
ci_nc_t(ncp = 2.83, df = 126, conf_level = .95)

# Modifying the above example so that a nonsymmetric 95% confidence interval
# can be formed:
ci_nc_t(ncp = 2.83, df = 126, alpha_lower = .01, alpha_upper = .04, conf_level = NULL)

# Modifying the above example so that a single-sided 95% confidence interval
# can be formed:
ci_nc_t(ncp = 2.83, df = 126, alpha_lower = 0, alpha_upper = .05, conf_level = NULL)

Confidence Interval for Omega Squared (Effect Size for ANOVA)

Description

Computes the point estimate and an exact, noncentrality-based confidence interval for the population omega squared (\omega^2), the proportion of variance in the dependent variable accounted for by a fixed effect. Accepts either the raw ANOVA summary (F, effect df, error df, total N) or a fitted aov / lm object, in which case the function returns a row per effect (partial \omega^2 in factorial designs).

Usage

ci_omega_squared(
  object = NULL,
  F_value = NULL,
  df_effect = NULL,
  df_error = NULL,
  N = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL
)

Arguments

Details

Point estimate. The function reports the usual sample omega squared, which for a one-way design can be written as

\hat{\omega}^2 = \frac{\mathit{SS}_{\text{effect}} - df_{\text{effect}} \cdot \mathit{MS}_{\text{error}}}{\mathit{SS}_{\text{total}} + \mathit{MS}_{\text{error}}} = \frac{df_{\text{effect}} (F - 1)}{df_{\text{effect}} (F - 1) + N}

(Hays, 1994; Keppel, 1991). For factorial designs the same formula applied per effect yields partial omega squared (Olejnik & Algina, 2003); values below zero are truncated to zero.

Confidence interval. The CI is constructed by Steiger's (2004, Proposition 1) confidence interval transformation principle: a CI for the noncentrality parameter \lambda of the F distribution is obtained (via ci_nc_F) and then mapped through

\omega^2_{\text{bound}} = \frac{\lambda_{\text{bound}}}{\lambda_{\text{bound}} + N}.

When the lower CI on \lambda is not identified (i.e., the observed F is below the one-sided critical value), the lower limit on \omega^2 is set to 0, matching the convention used in ci_pvaf. In a one-way design, the interval produced here is identical to the CI for \eta^2 from ci_pvaf; the two estimands coincide in the population and differ only in their sample estimators (an implication of the confidence interval transformation principle of Steiger, 2004).

Sums of squares in factorial designs. When a fitted model is supplied, the function reads the F-values from anova(), which in base R uses Type I (sequential) sums of squares. For balanced designs, Types I, II, and III give identical F-values; for unbalanced designs they differ. If Type II or III F-values are required, compute them with e.g.\ car::Anova(object, type = 3) and pass the relevant F / df into the raw-argument interface.

Value

A data.frame with one row per effect and the columns effect, omega_squared (point estimate), lower_limit, upper_limit, F_value, df_effect, df_error, and N. When the raw-argument interface is used, effect is "overall".

Author(s)

Ken Kelley kkelley@nd.edu

References

Fleishman, A. I. (1980). Confidence intervals for correlation ratios. Educational and Psychological Measurement, 40(3), 659–670.

Hays, W. L. (1994). Statistics (5th ed.). Fort Worth, TX: Harcourt Brace College Publishers.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Preacher, K. J. (2012). On effect size. Psychological Methods, 17, 137–152. doi:10.1037/a0028086

Keppel, G. (1991). Design and analysis: A researcher's handbook (3rd ed.). Prentice Hall.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on \eta^2, Chapter 7 on factorial designs, and Chapter 11 on generalized \eta^2 for within-subjects designs.)

Olejnik, S., & Algina, J. (2003). Generalized eta and omega squared statistics: Measures of effect size for some common research designs. Psychological Methods, 8(4), 434–447. doi:10.1037/1082-989X.8.4.434

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

ci_pvaf, ci_nc_F

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# 1. Raw-argument interface. Bargman's (1970) example, also used in
#        Venables (1975), Fleishman (1980), and Steiger (2004): a 5-group
#        one-way ANOVA with 11 subjects per group, observed F = 11.221.
ci_omega_squared(F_value = 11.221, df_effect = 4, df_error = 50, N = 55)

# Same example with a 90% confidence interval.
ci_omega_squared(
  F_value = 11.221, df_effect = 4, df_error = 50, N = 55,
  conf_level = 0.90
)

# 2. One way ANOVA from a fitted model: mean IQ gain differs across
#        the six grades of the pygmalion data (N = 310).
fit_one <- aov(iq_gain ~ factor(grade), data = pygmalion)
ci_omega_squared(fit_one)

# 3. Two-factor ANOVA: partial omega squared per effect for the
#        manipulated expectancy treatment and the measured grade
#        classification (pygmalion data, N = 310). The treatment by
#        grade interaction is weak here (F = 1.19), so the additive
#        model is used.
fit_additive <- aov(iq_8 ~ treatment + factor(grade), data = pygmalion)
ci_omega_squared(fit_additive)

Confidence Interval for a Single Proportion

Description

The Wilson (1927) score interval for a binomial proportion, the package's default for proportion inference: unlike the textbook Wald interval it cannot escape [0, 1], behaves sensibly at 0 and 1 counts, and holds close to nominal coverage at small n (Brown, Cai, & DasGupta, 2001, recommend it for general use). The Wald interval is available for instruction and comparison.

Usage

ci_proportion(successes, n, conf_level = 0.95, method = c("wilson", "wald"))

Arguments

Value

A data.frame (class dmar_tbl) with rows lower_limit, proportion, upper_limit, successes, and n, so the point estimate sits between its confidence limits.

Author(s)

Ken Kelley kkelley@nd.edu

References

Brown, L. D., Cai, T. T., & DasGupta, A. (2001). Interval estimation for a binomial proportion. Statistical Science, 16(2), 101–133. doi:10.1214/ss/1009213286

Wilson, E. B. (1927). Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association, 22(158), 209–212.

See Also

responder_analysis, which uses this interval for each group's responder proportion.

Examples

ci_proportion(successes = 17, n = 50)

# The Wilson interval stays inside [0, 1] even at the boundary.
ci_proportion(successes = 0, n = 20)

Confidence Interval for the Proportion of Variance Accounted for (in the Dependent Variable by Knowing the Levels of the Factor)

Description

Computes the exact confidence limits for the proportion of variance in the dependent variable accounted for by knowing the levels of the factor (group status in a single factor design) in a fixed effects analysis of variance, so an omnibus F-test is accompanied by an effect size with a statement of its precision.

Usage

ci_pvaf(
  F_value = NULL,
  df_1 = NULL,
  df_2 = NULL,
  N = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  ...
)

Arguments

Details

The confidence level must be specified in one of following two ways: using confidence interval coverage (conf_level), or lower and upper confidence limits (alpha_lower and alpha_upper).

This function uses the confidence interval transformation principle (Steiger, 2004) to transform the confidence limits for the noncentrality parameter to the confidence limits for the population proportion of variance accounted for by knowing the group status. The confidence interval for the noncentral F parameter can be obtained from the function ci_nc_F, which is used within this function.

Value

A 4-row data.frame with columns term, value, prob_less, and prob_greater. The term values are "lower_limit" (the lower confidence limit on the proportion of variance accounted for, on the [0, 1] scale), "pvaf" (the sample proportion of variance accounted for, df_1 * F_value / (df_1 * F_value + df_2), the same value that eta squared reports, so the point estimate sits between its confidence limits), "upper_limit" (the upper confidence limit), and "actual_coverage" (the achieved coverage probability, which equals conf_level when both tail targets are met). The prob_less and prob_greater columns report the achieved tail-error probabilities at the two limits; NA on the "pvaf" and "actual_coverage" rows.

Note

This function can be used for single or factorial ANOVA designs.

Author(s)

Ken Kelley kkelley@nd.edu

References

Fleishman, A. I. (1980). Confidence intervals for correlation ratios. Educational and Psychological Measurement, 40(3), 659–670.

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K. (2008). Sample size planning for the squared multiple correlation coefficient: Accuracy in parameter estimation via narrow confidence intervals. Multivariate Behavioral Research, 43, 524–555. doi:10.1080/00273170802490632

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 3 on R^2 as a model comparison effect size.)

Steiger, J. H. (2004). Beyond the F Test: Effect size confidence intervals and tests of close fit in the Analysis of Variance and Contrast Analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

ci_nc_F

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

## Bargman (1970) gave an example in which a 5-group ANOVA with 11 subjects in each
## group is conducted and the observed F value is 11.221. This example was used
## in Venables (1975),  Fleishman (1980), and Steiger (2004). If one wants to calculate the
## exact confidence interval for the proportion of variance accounted for in that example,
## this function can be used.
ci_pvaf(F_value = 11.221, df_1 = 4, df_2 = 50, N = 55)

ci_pvaf(F_value = 11.221, df_1 = 4, df_2 = 50, N = 55, conf_level = .90)

ci_pvaf(F_value = 11.221, df_1 = 4, df_2 = 50, N = 55, alpha_lower = 0, alpha_upper = .05)

Confidence Interval for an Unstandardized Regression Coefficient

Description

Computes a confidence interval for a population regression coefficient in its raw (unstandardized) metric, by the standard t-based approach or the noncentral t approach. A thin convenience wrapper around ci_reg_coef, which is the general engine; for the standardized coefficient use ci_src.

Usage

ci_rc(
  b_j,
  SE_b_j = NULL,
  s_Y = NULL,
  s_X = NULL,
  N,
  p,
  R2_Y_X = NULL,
  R2_j_X_without_j = NULL,
  conf_level = 0.95,
  R2_Y_X_without_j = NULL,
  t_value = NULL,
  alpha_lower = NULL,
  alpha_upper = NULL,
  noncentral = FALSE,
  ...
)

Arguments

Details

Returns the confidence limits for the regression coefficient of interest from the standard approach to confidence interval formation or from the noncentral approach to confidence interval formation using the noncentral t-distribution.

Value

A 2-row data.frame with columns term, value, prob_less, and prob_greater. The term values are "lower_limit" and "upper_limit", and value holds the confidence limits on the regression coefficient in its raw metric. The prob_less and prob_greater columns report the tail probabilities below and above each limit; when the noncentral t approach is used they are the achieved tail probabilities. Unlike ci_src and ci_reg_coef, which place the point estimate between its limits as a third row, ci_rc returns the two limits only.

Note

Not all of the values need to be specified, only those that contain all of the necessary information in order to compute the confidence interval (options are thus given for the values that need to be specified).

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Maxwell, S. E. (2003). Sample size for multiple regression: Obtaining regression coefficients that are accurate, not simply significant. Psychological Methods, 8(3), 305–321. doi:10.1037/1082-989X.8.3.305

Kelley, K., & Maxwell, S. E. (2008). Sample size planning with applications to multiple regression: Power and accuracy for omnibus and targeted effects. In P. Alasuutari, L. Bickman, & J. Brannen (Eds.), The Sage handbook of social research methods (pp. 166–192). Sage.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 4 on individual comparisons of means and Chapter 6 on trend analysis.)

Smithson, M. (2003). Confidence intervals. Thousand Oaks, CA: Sage Publications.

Steiger, J. H. (2004). Beyond the F Test: Effect size confidence intervals and tests of close fit in the Analysis of Variance and Contrast Analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

ss_aipe_reg_coef, ci_nc_t, ci_reg_coef, ci_src

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

ci_rc(b_j = 0.61319, SE_b_j = 0.16098, N = 30, p = 6, conf_level = 0.95)

Confidence Interval for a Regression Coefficient, Raw or Standardized

Description

The general engine behind ci_rc (unstandardized) and ci_src (standardized): computes a confidence interval for a population regression coefficient by the standard t-based approach or the noncentral t approach, in whichever metric the inputs are supplied.

Usage

ci_reg_coef(
  b_j,
  SE_b_j = NULL,
  s_Y = NULL,
  s_X = NULL,
  N,
  p,
  R2_Y_X = NULL,
  R2_j_X_without_j = NULL,
  conf_level = 0.95,
  R2_Y_X_without_j = NULL,
  t_value = NULL,
  alpha_lower = NULL,
  alpha_upper = NULL,
  noncentral = FALSE,
  ...
)

Arguments

Details

For standardized variables, do not specify the standard deviation of the variables and input the standardized regression coefficient for b_j.

When b_j is reconstructed from squared multiple correlations (that is, from R2_Y_X, R2_Y_X_without_j, and R2_j_X_without_j rather than a supplied b_j, SE_b_j, or t_value), only the magnitude of the coefficient is identifiable; its sign is not. The positive root is returned and a warning is issued. If the coefficient is negative, negate the point estimate and swap and negate the confidence limits, or supply b_j directly.

Value

A 3-row data.frame with columns term, value, prob_less, and prob_greater. The rows are ordered "lower_limit", "reg_coef" (the regression coefficient point estimate), and "upper_limit", so the point estimate sits between its confidence limits. The lower and upper rows give the confidence limits on the regression coefficient. The prob_less and prob_greater columns report the achieved tail probabilities at each limit when the noncentral t method is used (they are NA for the "reg_coef" estimate row).

Note

Not all of the values need to be specified, only those that contain all of the necessary information in order to compute the confidence interval (options are thus given for the values that need to be specified).

The function ci_rc in DMAR also calculates the confidence interval for the population (unstandardized) regression coefficient. The function ci_src also calculates the confidence interval for the population (standardized) regression coefficient. These two functions perform the same tasks as ci_reg_coef does and are preferred to it because of simpler arguments.

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Maxwell, S. E. (2003). Sample size for multiple regression: Obtaining regression coefficients that are accurate, not simply significant. Psychological Methods, 8(3), 305–321. doi:10.1037/1082-989X.8.3.305

Kelley, K., & Maxwell, S. E. (2008). Sample size planning with applications to multiple regression: Power and accuracy for omnibus and targeted effects. In P. Alasuutari, L. Bickman, & J. Brannen (Eds.), The Sage handbook of social research methods (pp. 166–192). Sage.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 4 on individual comparisons of means and Chapter 6 on trend analysis.)

Smithson, M. (2003). Confidence intervals. Thousand Oaks, CA: Sage Publications.

See Also

ss_aipe_reg_coef, ci_nc_t, ci_rc, ci_src

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

ci_reg_coef(b_j = 0.61319, SE_b_j = 0.16098, N = 30, p = 6)

Confidence Interval for the Population Root Mean Square Error of Approximation

Description

Constructs a confidence interval for the population root mean square error of approximation (RMSEA), a population badness-of-fit index for structural equation models. The interval is obtained by inverting the noncentral chi square distribution of the sample fit function T = (N - 1) \hat F_{ML} under the model implied covariance structure, mapping the resulting noncentrality limits to the RMSEA metric (Steiger & Lind, 1980; Browne & Cudeck, 1993).

Usage

ci_rmsea(
  rmsea,
  df,
  N,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL
)

Arguments

Details

The RMSEA expresses the badness of model fit per degree of freedom on the noncentrality scale. Under the noncentral chi square model for the sample fit statistic, the sample T = (N - 1) \hat F_{ML} has approximate noncentral chi square distribution with df degrees of freedom and noncentrality parameter \lambda = (N - 1) df \cdot \mathrm{RMSEA}^2. The CI on \mathrm{RMSEA}^2 is obtained by inverting the noncentral chi square distribution at the requested confidence level (ci_nc_chisq does the inversion); the bounds are then mapped back to the RMSEA scale via the square root. When the lower noncentrality limit hits zero (i.e., the data are compatible with a well-fitting model), the lower RMSEA limit is truncated at zero because RMSEA is non-negative by construction.

The 90 percent CI (rather than the usual 95 percent) is the conventional reporting choice for RMSEA (Browne & Cudeck, 1993) because the upper limit of the 90 percent CI plays a one-sided role in the test of close fit (H_0: \mathrm{RMSEA} \le 0.05). ci_rmsea defaults to conf_level = 0.95 in line with the rest of the package; pass conf_level = 0.90 when the close fit test is the intended use.

Value

A 3-row data.frame with columns term and value. The term values are "lower_limit" (the lower bound of the confidence interval on the population RMSEA, truncated at zero by definition), "rmsea" (the observed point estimate), and "upper_limit" (the upper bound).

Author(s)

Ken Kelley kkelley@nd.edu

References

Browne, M. W., & Cudeck, R. (1993). Alternative ways of assessing model fit. In K. A. Bollen & J. S. Long (Eds.), Testing structural equation models (pp. 136–162). Sage.

Kelley, K., & Lai, K. (2011). Accuracy in parameter estimation for the root mean square error of approximation: Sample size planning for narrow confidence intervals. Multivariate Behavioral Research, 46, 1–32. doi:10.1080/00273171.2011.543027

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge.

Steiger, J. H., & Lind, J. C. (1980). Statistically-based tests for the number of common factors. Paper presented at the annual Spring meeting of the Psychometric Society, Iowa City, IA.

See Also

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_sc(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# 1. A typical 95 percent CI on RMSEA.
ci_rmsea(rmsea = .055, df = 40, N = 425, conf_level = .95)

# 2. The 90 percent CI is the conventional choice when interpretation
#    will follow the Browne and Cudeck (1993) close fit decision rule
#    (the test of H_0: RMSEA <= 0.05 vs. the upper CI limit). Here the
#    upper limit is 0.052, just above the close fit threshold of 0.05
#    that Browne and Cudeck recommend, so close fit is not established
#    even though the point estimate sits comfortably below it.
ci_rmsea(rmsea = .035, df = 40, N = 425, conf_level = .90)

# 3. Wider model with smaller N: more uncertainty, wider CI.
ci_rmsea(rmsea = .055, df = 10, N = 100, conf_level = .90)

Confidence Interval for a Standardized Contrast in a Fixed Effects ANOVA

Description

Computes the exact noncentral t-based confidence interval for a standardized contrast of means in a fixed effects analysis of variance: the contrast of interest divided by the error standard deviation, so groups measured in raw units are compared on a common standardized scale.

Usage

ci_sc(
  means = NULL,
  s_anova = NULL,
  c_weights = NULL,
  n = NULL,
  N = NULL,
  psi = NULL,
  ncp = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  df_error = NULL,
  ...
)

Arguments

Value

A 3-row data.frame with columns term and value. The term values are "lower_limit" (the lower confidence limit on the standardized contrast), "std_contrast" (the standardized contrast), and "upper_limit" (the upper limit).

Note

Be sure to use the standard deviation and not the error variance for s_anova, not the square of this value (the error variance) which would come from the source table (i.e., do not use the variance of the error but rather use its square root, the standard deviation).

Be sure to use fractional c-weights when doing complex contrasts (not integers) to specify c_weights. For example, in an ANCOVA of four groups, if the user wants to compare the mean of group 1 and 2 with the mean of group 3 and 4, c_weights should be specified as c(0.5, 0.5, -0.5, -0.5) rather than c(1, 1, -1, -1). Make sure the sum of the contrast weights are zero.

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Lai, K., & Kelley, K. (2012). Accuracy in parameter estimation for ANCOVA and ANOVA contrasts: Sample size planning via narrow confidence intervals. British Journal of Mathematical and Statistical Psychology, 65, 350–370. doi:10.1111/j.2044-8317.2011.02029.x

Steiger, J. H. (2004). Beyond the F Test: Effect size confidence intervals and tests of close fit in the Analysis of Variance and Contrast Analysis. Psychological Methods, 9(2), 164–182. doi:10.1037/1082-989X.9.2.164

See Also

ci_nc_t, ci_src, ci_smd, ci_smd_c, ci_sm, ci_c ci_c_ancova

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc_ancova(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Here is a four group example. Suppose that the means of groups 1--4 are 2, 4, 9,
# and 13, respectively. Further, let the error variance be .64 and thus the standard
# deviation would be .80 (note we use the standard deviation in the function, not the
# variance). The standardized contrast of interest here is the average of groups 1 and 4
# versus the average of groups 2 and 3.

ci_sc(means = c(2, 4, 9, 13), s_anova = .80, c_weights = c(.5, -.5, -.5, .5),
      n = c(3, 3, 3, 3), N = 12, conf_level = .95)

# Here is an example with two groups.
ci_sc(means = c(1.6, 0), s_anova = .80, c_weights = c(1, -1),
      n = c(10, 10), N = 20, conf_level = .95)

Confidence Interval for a Standardized Contrast in ANCOVA With One Covariate

Description

Calculate the confidence interval for a standardized contrast in ANCOVA with one covariate. The standardizer (i.e., the divisor) can be either the error standard deviation of the ANOVA model (i.e., the model excluding the covariate) or of the ANCOVA model.

Usage

ci_sc_ancova(
  psi = NULL,
  adj_means = NULL,
  s_anova = NULL,
  s_ancova = NULL,
  standardizer = "s_ancova",
  c_weights,
  n,
  cov_means,
  SSwithin_x,
  conf_level = 0.95
)

Arguments

Details

The argument SSwithin_x is the sum of squares within groups for the covariate, taken from the ANOVA source table in which the covariate (not the outcome) is the dependent variable. Published reports do not always print this quantity directly. When a report gives the covariate group means, the group sample sizes, and the F statistic from the one-way ANOVA on the covariate, SSwithin_x can be recovered algebraically. The worked example below follows Lai and Kelley (2012): three groups of sizes 19, 18, and 19 (so N = 56) have covariate means 60.08, 57.08, and 57.97, and the covariate ANOVA reports F = 0.756 with 2 and 53 degrees of freedom. The sum of squares between groups for the covariate, computed from the group means and sample sizes, is approximately 88.5, so the mean square between groups is approximately 88.5 / 2 = 44.3. Because F is the ratio of the mean square between groups to the mean square within groups, the mean square within groups is approximately 44.3 / 0.756 = 58.6, and the sum of squares within groups is that mean square times its degrees of freedom, approximately 58.6 \times 53 = 3103. That recovered value is what you would pass to SSwithin_x. The “Examples” section reproduces this computation in code.

Value

A 3-row data.frame with columns term and value (numeric). The term values are "lower_limit" (the lower confidence limit on the standardized ANCOVA contrast), "psi" (the standardized contrast), and "upper_limit" (the upper limit). The divisor used in standardization (either "s_anova" or "s_ancova") is attached as the "standardizer" attribute of the returned data.frame.

Note

Be sure to use the standard deviations and not the error variances for s_anova and s_ancova, not the squares of these values which would come from the source tables (i.e., do not use the variance of the errors but rather use its square root, the standard deviation).

If n receives a single number, that number is considered as the sample size per group. If n is assigned to a vector, the vector is considered as the sample size of each group.

Be sure to use fractional c-weights when doing complex contrasts (not integers) to specify c_weights. For example, in an ANCOVA of four groups, if the user wants to compare the mean of group 1 and 2 with the mean of group 3 and 4, c_weights should be specified as c(0.5, 0.5, -0.5, -0.5) rather than c(1, 1, -1, -1). Make sure the sum of the contrast weights are zero.

The argument to be assigned to standardizer must be either "s_ancova" or "s_anova".

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Kelley, K., & Rausch, J. R. (2006). Sample size planning for the standardized mean difference: Accuracy in parameter estimation via narrow confidence intervals. Psychological Methods, 11, 363–385. doi:10.1037/1082-989X.11.4.363

Lai, K., & Kelley, K. (2012). Accuracy in parameter estimation for ANCOVA and ANOVA contrasts: Sample size planning via narrow confidence intervals. British Journal of Mathematical and Statistical Psychology, 65, 350–370. doi:10.1111/j.2044-8317.2011.02029.x

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 9.)

Steiger, J. H., & Fouladi, R. T. (1997). Noncentrality interval estimation and the evaluation of statistical methods. In L. L. Harlow, S. A. Mulaik, & J. H. Steiger (Eds.), What if there were no significance tests? (pp. 221–257). Mahwah, NJ: Lawrence Erlbaum.

See Also

ci_c_ancova, ci_sc

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sm(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()

Examples

# Maxwell, Delaney, & Kelley (2027) offer an example that 30 depressive
# individuals are randomly assigned to three groups, 10 in each, and ANCOVA
# is performed on the posttest scores using the participants' pretest
# scores as the covariate. The means of pretest scores of group 1, 2, and 3 are
# 17, 17.7, and 17.4, respectively, whereas the adjusted means of groups 1, 2, and 3
# are 7.5, 12, and 14, respectively. The error variance in ANCOVA is 29 and thus
# 5.385165 is the error standard deviation, with the sum of squares within groups
# from an ANOVA on the covariate is 752.5.

# To obtained the confidence interval for the standardized adjusted mean difference
# between group 1 and 2, using the ANCOVA error standard deviation:
ci_sc_ancova(adj_means = c(7.5, 12, 14), s_ancova = 5.385165, c_weights = c(1, -1, 0),
             n = 10, cov_means = c(17, 17.7, 17.4), SSwithin_x = 752.5)

# Or, with less error in rounding:
ci_sc_ancova(adj_means = c(7.54, 11.98, 13.98), s_ancova = 5.393, c_weights = c(-1, 0, 1),
             n = 10, cov_means = c(17, 17.7, 17.4), SSwithin_x = 752.5)

# Now, using the standard deviation from ANOVA (and not ANCOVA as above), we have:
ci_sc_ancova(adj_means = c(7.54, 11.98, 13.98), s_anova = 6.294, s_ancova = 5.393,
             c_weights = c(-1, 0, 1),n = 10, cov_means = c(17, 17.7, 17.4),
             SSwithin_x = 752.5, standardizer = "s_anova", conf_level = .95)

# Recovering SSwithin_x from a covariate ANOVA F when a report does not print
# it directly (see the Details section). This example follows Lai and Kelley
# (2012): three groups of sizes 19, 18, and 19 have covariate means 60.08,
# 57.08, and 57.97, and the one-way ANOVA on the covariate reports F = 0.756.
cov_means_ex  <- c(60.08, 57.08, 57.97)
n_ex          <- c(19, 18, 19)
grand_x       <- sum(n_ex * cov_means_ex) / sum(n_ex)
ss_between_x  <- sum(n_ex * (cov_means_ex - grand_x)^2)
ms_between_x  <- ss_between_x / (length(n_ex) - 1)
ms_within_x   <- ms_between_x / 0.756
SSwithin_x_ex <- ms_within_x * (sum(n_ex) - length(n_ex))

ci_sc_ancova(adj_means = c(63.88, 62.39, 56.48), s_ancova = 20.267,
             c_weights = c(0.5, 0.5, -1), n = n_ex, cov_means = cov_means_ex,
             SSwithin_x = SSwithin_x_ex)

Scheffe-Adjusted Simultaneous Confidence Intervals for Contrasts

Description

Computes the Scheffe (1953, 1959) simultaneous confidence intervals on user-specified contrasts among the means of a one-way design. The Scheffe procedure controls the family-wise error rate for any set of contrasts, however many and however post-hoc, which makes it more conservative than Tukey-Kramer or Bonferroni for the specific case of all-pairwise comparisons but optimal for arbitrary post-hoc contrasts.

Usage

ci_scheffe(x, group = NULL, contrasts = NULL, conf_level = 0.95)

Arguments

Details

Critical value. For a groups with \nu error degrees of freedom, the Scheffe critical value is

S \;=\; \sqrt{(a - 1) F_{1 - \alpha, a - 1, \nu}},

where F_{1 - \alpha, a - 1, \nu} is the upper \alpha quantile of the central F distribution. The Scheffe simultaneous CI on a contrast \psi = \sum_i c_i \mu_i is

\hat\psi \;\pm\; S \cdot \mathit{SE}(\hat\psi),

where \mathit{SE}(\hat\psi) = \sqrt{\mathit{MS}_E \sum_i c_i^2 / n_i}.

Scope. The Scheffe family-wise coverage holds for any number of contrasts, pairwise, complex, or chosen after looking at the data. The trade-off is conservativeness: for all-pairwise comparisons, Tukey-Kramer is uniformly more powerful.

Value

A data.frame with one row per contrast. Columns: contrast (a printed label), contrast_value, se, F_statistic, lower_limit, upper_limit, p_adjusted.

Author(s)

Ken Kelley kkelley@nd.edu

References

Scheffe, H. (1953). A method for judging all contrasts in the analysis of variance. Biometrika, 40(1/2), 87–104.

Scheffe, H. (1959). The analysis of variance. Wiley.

See Also

cv_scheffe, ci_tukey_kramer, ci_dunnett

Other hypothesis tests: adjusted_means(), ancova(), anova_within(), ci_dunnett(), ci_tukey_kramer(), compare_cov_structures(), contrast_test(), correlations_test(), equivalence_r(), equivalence_smd(), factorial_anova(), manova_split_plot(), mauchly_test(), mixed_anova(), obrien_test(), pairwise_within(), randomization_test(), randomization_test_paired(), regions_of_significance(), simple_effects_AB(), summary_t_test(), welch_t()

Examples

# 1. All pairwise contrasts among the six marketing panels of the
#    test_market data via the default:
fit <- lm(brand_movement ~ panel, data = test_market)
ci_scheffe(fit)

# 2. A contrast chosen after inspecting the means: the two panels with
#    the highest brand movement (5 and 6) against the two with the
#    lowest (1 and 2). The Scheffe coverage holds for a contrast picked
#    this way, and the interval still excludes zero even though none of
#    the pairwise intervals above does.
cmat <- matrix(c(-0.5, -0.5, 0, 0, 0.5, 0.5), nrow = 6,
               dimnames = list(levels(test_market$panel),
                               "panels 5,6 - panels 1,2"))
ci_scheffe(fit, contrasts = cmat)

Confidence Interval for the Standardized Mean

Description

Computes the exact confidence interval for the standardized mean, the mean divided by the standard deviation, by inverting the noncentral t distribution. The standardized mean is the one-sample analog of the standardized mean difference and shares its noncentral interval theory.

Usage

ci_sm(
  sm = NULL,
  mean = NULL,
  sd = NULL,
  ncp = NULL,
  N = NULL,
  conf_level = 0.95,
  alpha_lower = NULL,
  alpha_upper = NULL,
  ...
)

Arguments

Details

The user must specify the standardized mean in one and only one of the three ways: a) mean and standard deviation (mean and sd), b) standardized mean (sm), and c) noncentral parameter (ncp). The confidence level must be specified in one of following two ways: using confidence interval coverage (conf_level), or lower and upper confidence limits (alpha_lower and alpha_upper). This function uses the exact confidence interval method based on noncentral t-distributions. The confidence interval for noncentral t-parameter can be obtained from the ci_nc_t function in DMAR.

Value

A 3-row data.frame with columns term and value. The term values are "lower_limit" (the lower confidence limit on the standardized mean), "std_mean" (the standardized mean), and "upper_limit" (the upper confidence limit).

Note

The standardized mean is the mean divided by the standard deviation.

Author(s)

Ken Kelley kkelley@nd.edu

References

Kelley, K. (2007). Confidence intervals for standardized effect sizes: Theory, application, and implementation. Journal of Statistical Software, 20(8), 1–24. doi:10.18637/jss.v020.i08

Steiger, J. H., & Fouladi, R. T. (1997). Noncentrality interval estimation and the evaluation of statistical methods. In L. L. Harlow, S. A. Mulaik, & J. H. Steiger (Eds.), What if there were no significance tests? (pp. 221–257). Mahwah, NJ: Lawrence Erlbaum.

See Also

ci_nc_t

Other confidence intervals for effect sizes: ci_R2(), ci_c(), ci_c_ancova(), ci_c_ancova_bp(), ci_correlation, ci_cv(), ci_eta_squared(), ci_eta_squared_generalized(), ci_eta_squared_partial(), ci_mahalanobis(), ci_omega_squared(), ci_pvaf(), ci_rc(), ci_reg_coef(), ci_rmsea(), ci_sc(), ci_sc_ancova(), ci_smd(), ci_smd_c(), ci_snr(), ci_src(), ci_srsnr(), contrast_adjusted(), plot_smd()