diff --git a/NAMESPACE b/NAMESPACE index 233ea5b6..3f138851 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -238,6 +238,7 @@ export(t_to_w) export(tschuprows_t) export(v_to_t) export(v_to_w) +export(variance_ratio) export(vd_a) export(w_to_c) export(w_to_fei) diff --git a/NEWS.md b/NEWS.md index bfd8f3f0..0875f261 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,7 +6,7 @@ - `repeated_measures_d()` to compute standardized mean differences (SMD) for repeated measures data. - Also supported in `effectsize()` - +- `variance_ratio()` to compute variance ratios for independent or paired samples. ## Bug fixes - `nnt()` now properly accepts the `y` argument. diff --git a/R/sysdata.rda b/R/sysdata.rda index 50806d57..9c9957fc 100644 Binary files a/R/sysdata.rda and b/R/sysdata.rda differ diff --git a/R/vars_ratio.R b/R/vars_ratio.R new file mode 100644 index 00000000..d3f7ffc9 --- /dev/null +++ b/R/vars_ratio.R @@ -0,0 +1,146 @@ +#' Ratio of Variances +#' +#' Computes the ratio of two variances of independent or paired samples. For +#' paired data, this can also be thought of as the ratio of marginal (squared) +#' scales of a bivariate normal distribution. For independent samples, this is a +#' convenient wrapper around [stats::var.test()] (for a paired version of this +#' test - the Pitman-Morgan test, see `PairedData::Var.test()`). +#' +#' @inheritParams cohens_d +#' @inheritParams means_ratio +#' +#' @details +#' +#' # Confidence (Compatibility) Intervals (CIs) +#' For independent samples, confidence intervals are estimated using +#' [stats::var.test()]. For paired samples, confidence intervals are estimated +#' using the standard parametric method (see Pitman, 1939; Morgan, 1929). +#' +#' @inheritSection effectsize_CIs CIs and Significance Tests +#' @inheritSection print.effectsize_table Plotting with `see` +#' +#' @return A data frame with the effect size (`Vars_ratio`) and its CIs +#' (`CI_low` and `CI_high`). +#' +#' @seealso [means_ratio()], [stats::var.test()] and `PairedData::Var.test()` +#' for the _Pitman-Morgan Test_ for paired samples. +#' +#' @references +#' - Morgan, W. A. (1939). A test for the significance of the difference between +#' the two variances in a sample from a normal bivariate population. Biometrika, +#' 31(1/2), 13-19. +#' - Pitman, E. J. G. (1939). A note on normal correlation. Biometrika, 31(1/2), 9-12. +#' +#' @examples +#' # Two Independent Samples ---------- +#' +#' x <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29) +#' y <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30) +#' +#' variance_ratio(x, y) +#' +#' # More options: +#' variance_ratio(x, y, alternative = "less") +#' variance_ratio(x, y, log = TRUE) +#' +#' # The ratio is scale invariant +#' variance_ratio(x * 3, y * 3) +#' +#' +#' +#' # Paired Samples ---------- +#' +#' sleep2 <- reshape(sleep, +#' direction = "wide", +#' idvar = "ID", timevar = "group" +#' ) +#' variance_ratio(Pair(extra.1, extra.2) ~ 1, data = sleep2) +#' +#' @export +variance_ratio <- function(x, y = NULL, data = NULL, + paired = FALSE, log = FALSE, + ci = 0.95, alternative = "two.sided", + verbose = TRUE, ...) { + alternative <- .match.alt(alternative) + out <- .get_data_2_samples(x, y, data, paired = paired, verbose = verbose, ...) + x <- out[["x"]] + y <- out[["y"]] + paired <- out[["paired"]] + + if (is.null(y)) { + insight::format_error("'y' is missing.") + } + + if (paired) { + out <- .var_ratio_dep(x, y, ci, alternative) + } else { + out <- .var_ratio_ind(x, y, ci, alternative) + } + + if (log) { + i <- intersect(colnames(out), c("Vars_ratio", "CI_low", "CI_high")) + out[i] <- sapply(out[i], log) + colnames(out)[1] <- "log_Vars_ratio" + } + + out +} + +#' @keywords internal +.var_ratio_ind <- function(x, y, ci, alternative) { + conf.level <- if (is.null(ci)) 0.95 else ci + vt <- stats::var.test(x, y, alternative = alternative, conf.level = conf.level) + pars <- parameters::model_parameters(vt) + pars$CI <- ci + out <- pars[intersect(colnames(pars), c("Estimate", "CI", "CI_low", "CI_high"))] + if (is.null(ci)) { + out$CI_low <- out$CI_high <- NULL + } + colnames(out)[1] <- "Vars_ratio" + ci_method <- list(method = "F") + + class(out) <- c("effectsize_table", "see_effectsize_table", "data.frame") + .someattributes(out) <- .nlist( + paired = FALSE, ci, ci_method, alternative, + approximate = FALSE + ) + return(out) +} + +#' @keywords internal +.var_ratio_dep <- function(x, y, ci, alternative) { + v1 <- stats::var(x) + v2 <- stats::var(y) + w <- v1 / v2 + out <- data.frame(Vars_ratio = w) + + if (.test_ci(ci)) { + # Add cis + out$CI <- ci + ci.level <- .adjust_ci(ci, alternative) + alpha <- 1 - ci.level + + # Pitman 1939, pp 11 + n <- length(x) + df <- n - 2 + r <- stats::cor(x, y) + qs <- qt(alpha / 2, df = n - 2) + K <- 1 + (2 * (1 - r^2) * qs^2) / (n - 2) + CONF <- w * (K + c(-1, 1) * sqrt(K^2 - 1)) + + out$CI_low <- CONF[1] + out$CI_high <- CONF[2] + ci_method <- list(method = "t") + out <- .limit_ci(out, alternative, 0, Inf) + } else { + ci_method <- alternative <- NULL + } + + + class(out) <- c("effectsize_table", "see_effectsize_table", "data.frame") + .someattributes(out) <- .nlist( + paired = TRUE, ci, ci_method, alternative, + approximate = FALSE + ) + return(out) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 4defa812..077e81fc 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -19,6 +19,7 @@ reference: - rank_biserial - p_superiority - means_ratio + - variance_ratio - subtitle: "For contingency tables" desc: > diff --git a/data-raw/es_info.R b/data-raw/es_info.R index 88044de6..4c5c413b 100644 --- a/data-raw/es_info.R +++ b/data-raw/es_info.R @@ -12,8 +12,10 @@ es_info <- tibble::tribble( "Mahalanobis_D", "Mahalanobis' D", NA, "onetail", 0, Inf, 0, "Means_ratio", "Means Ratio", NA, "twotail", 0, Inf, 1, "Means_ratio_adjusted", "Means Ratio (adj.)", NA, "twotail", 0, Inf, 1, - "log_Means_ratio", "log(Means Ratio)", NA, "twotail", 0, Inf, 1, - "log_Means_ratio_adjusted", "log(Means Ratio, adj.)", NA, "twotail", 0, Inf, 1, + "log_Means_ratio", "log(Means Ratio)", NA, "twotail", -Inf, Inf, 0, + "log_Means_ratio_adjusted", "log(Means Ratio, adj.)", NA, "twotail", -Inf, Inf, 0, + "Vars_ratio", "Variance Ratio", NA, "twotail", 0, Inf, 1, + "log_Vars_ratio", "log(Variance Ratio)", NA, "twotail", -Inf, Inf, 0, ## xtab cor "Cohens_w", "Cohen's w", NA, "onetail", 0, Inf, 0, diff --git a/man/variance_ratio.Rd b/man/variance_ratio.Rd new file mode 100644 index 00000000..66d1023f --- /dev/null +++ b/man/variance_ratio.Rd @@ -0,0 +1,126 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/vars_ratio.R +\name{variance_ratio} +\alias{variance_ratio} +\title{Ratio of Variances} +\usage{ +variance_ratio( + x, + y = NULL, + data = NULL, + paired = FALSE, + log = FALSE, + ci = 0.95, + alternative = "two.sided", + verbose = TRUE, + ... +) +} +\arguments{ +\item{x, y}{A numeric vector, or a character name of one in \code{data}. +Any missing values (\code{NA}s) are dropped from the resulting vector. +\code{x} can also be a formula (see \code{\link[stats:t.test]{stats::t.test()}}), in which case \code{y} is +ignored.} + +\item{data}{An optional data frame containing the variables.} + +\item{paired}{If \code{TRUE}, the values of \code{x} and \code{y} are considered as paired. +This produces an effect size that is equivalent to the one-sample effect +size on \code{x - y}. See also \code{\link[=repeated_measures_d]{repeated_measures_d()}} for more options.} + +\item{log}{Should the log-ratio be returned? Defaults to \code{FALSE}. +Normally distributed and useful for meta-analysis.} + +\item{ci}{Confidence Interval (CI) level} + +\item{alternative}{a character string specifying the alternative hypothesis; +Controls the type of CI returned: \code{"two.sided"} (default, two-sided CI), +\code{"greater"} or \code{"less"} (one-sided CI). Partial matching is allowed (e.g., +\code{"g"}, \code{"l"}, \code{"two"}...). See \emph{One-Sided CIs} in \link{effectsize_CIs}.} + +\item{verbose}{Toggle warnings and messages on or off.} + +\item{...}{Arguments passed to or from other methods. When \code{x} is a formula, +these can be \code{subset} and \code{na.action}.} +} +\value{ +A data frame with the effect size (\code{Vars_ratio}) and its CIs +(\code{CI_low} and \code{CI_high}). +} +\description{ +Computes the ratio of two variances of independent or paired samples. For +paired data, this can also be thought of as the ratio of marginal (squared) +scales of a bivariate normal distribution. For independent samples, this is a +convenient wrapper around \code{\link[stats:var.test]{stats::var.test()}} (for a paired version of this +test - the Pitman-Morgan test, see \code{PairedData::Var.test()}). +} +\section{Confidence (Compatibility) Intervals (CIs)}{ +For independent samples, confidence intervals are estimated using +\code{\link[stats:var.test]{stats::var.test()}}. For paired samples, confidence intervals are estimated +using the standard parametric method (see Pitman, 1939; Morgan, 1929). +} + +\section{CIs and Significance Tests}{ +"Confidence intervals on measures of effect size convey all the information +in a hypothesis test, and more." (Steiger, 2004). Confidence (compatibility) +intervals and p values are complementary summaries of parameter uncertainty +given the observed data. A dichotomous hypothesis test could be performed +with either a CI or a p value. The 100 (1 - \eqn{\alpha})\% confidence +interval contains all of the parameter values for which \emph{p} > \eqn{\alpha} +for the current data and model. For example, a 95\% confidence interval +contains all of the values for which p > .05. +\cr\cr +Note that a confidence interval including 0 \emph{does not} indicate that the null +(no effect) is true. Rather, it suggests that the observed data together with +the model and its assumptions combined do not provided clear evidence against +a parameter value of 0 (same as with any other value in the interval), with +the level of this evidence defined by the chosen \eqn{\alpha} level (Rafi & +Greenland, 2020; Schweder & Hjort, 2016; Xie & Singh, 2013). To infer no +effect, additional judgments about what parameter values are "close enough" +to 0 to be negligible are needed ("equivalence testing"; Bauer & Kiesser, +1996). +} + +\section{Plotting with \code{see}}{ + +The \code{see} package contains relevant plotting functions. See the \href{https://easystats.github.io/see/articles/effectsize.html}{plotting vignette in the \code{see} package}. +} + +\examples{ +# Two Independent Samples ---------- + +x <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29) +y <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30) + +variance_ratio(x, y) + +# More options: +variance_ratio(x, y, alternative = "less") +variance_ratio(x, y, log = TRUE) + +# The ratio is scale invariant +variance_ratio(x * 3, y * 3) + + + +# Paired Samples ---------- + +sleep2 <- reshape(sleep, + direction = "wide", + idvar = "ID", timevar = "group" +) +variance_ratio(Pair(extra.1, extra.2) ~ 1, data = sleep2) + +} +\references{ +\itemize{ +\item Morgan, W. A. (1939). A test for the significance of the difference between +the two variances in a sample from a normal bivariate population. Biometrika, +31(1/2), 13-19. +\item Pitman, E. J. G. (1939). A note on normal correlation. Biometrika, 31(1/2), 9-12. +} +} +\seealso{ +\code{\link[=means_ratio]{means_ratio()}}, \code{\link[stats:var.test]{stats::var.test()}} and \code{PairedData::Var.test()} +for the \emph{Pitman-Morgan Test} for paired samples. +} diff --git a/tests/testthat/test-rov.R b/tests/testthat/test-rov.R new file mode 100644 index 00000000..27b758e0 --- /dev/null +++ b/tests/testthat/test-rov.R @@ -0,0 +1,37 @@ +test_that("variance_ratio | independent", { + x <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29) + y <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30) + + rov1 <- variance_ratio(x, y) + rov2 <- variance_ratio(y, x) + + expect_equal(1 / rov2[[1]], rov1[[1]]) + expect_equal(1 / rov2$CI_low, rov1$CI_high) + expect_equal(1 / rov2$CI_high, rov1$CI_low) + + lrov1 <- variance_ratio(x, y, log = TRUE) + lrov2 <- variance_ratio(y, x, log = TRUE) + + expect_equal(-lrov2[[1]], lrov1[[1]]) + expect_equal(-lrov2$CI_low, lrov1$CI_high) + expect_equal(-lrov2$CI_high, lrov1$CI_low) + + expect_error(variance_ratio(x), regexp = "missing") +}) + + +test_that("variance_ratio | paired", { + sleep2 <- reshape(sleep, + direction = "wide", + idvar = "ID", timevar = "group" + ) + pdat1 <- Pair(sleep2$extra.1, sleep2$extra.2) + pdat2 <- Pair(sleep2$extra.2, sleep2$extra.1) + + rov1 <- variance_ratio(pdat1) + rov2 <- variance_ratio(pdat2) + + expect_equal(1 / rov2[[1]], rov1[[1]]) + expect_equal(1 / rov2$CI_low, rov1$CI_high) + expect_equal(1 / rov2$CI_high, rov1$CI_low) +})