---
title: "varequal"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{varequal}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# 1. Overview

The **`varequal`** R package provides statistical methods for assessing
**homoscedasticity**, that is, whether the variances of a response variable are
equal across independent groups.

Homogeneity of variance is an important assumption for many classical
statistical procedures, including one-way ANOVA and related linear-model
analyses. When the group variances differ substantially, the validity of
procedures that rely on a common variance can be affected.

`varequal` provides a unified interface to several classical, robust,
rank-based, and variance-outlier-oriented procedures. The package also provides
high-level functions for automatically combining several tests into a single
variance-homogeneity decision.

The main functionality can be divided into four levels:

- **Individual tests** for explicit selection of a homogeneity-of-variance
  procedure.
- **`check_var_equal()`** as a wrapper for selecting an individual test by a
  short method code.
- **`is_var_equal()`** as a higher-level assessment combining five tests.
- **Supporting datasets and utilities** for demonstrating and investigating
  the behavior of the procedures.

The currently implemented individual tests are:

| Code | Function | Main characteristic |
|:---|:---|:---|
| `AB` | `Ansari_Bradley_test()` | Rank-based dispersion test |
| `BL` | `Bartlett_test()` | Classical parametric test |
| `BF` | `Brown_Forsythe_test()` | Robust Levene-type procedure |
| `FK` | `Fligner_Killeen_test()` | Robust rank-based procedure |
| `LG` | `Lam_G_test()` | Variance-outlier-oriented procedure |
| `LV` | `Levene_test()` | ANOVA-based robust procedure |
| `MBF` | `Brown_Forsythe_test(method = "MBF")` | Mehrotra-modified Brown--Forsythe |
| `OB` | `O.Brien_test()` | O'Brien variance test |
| `OM` | `O.Neill_Mathews_test()` | Weighted least-squares Levene-type test |

The package source currently describes `varequal` as a toolkit for classical and
robust tests of variance homogeneity, together with supporting measures and
utilities for exploratory analysis and hypothesis testing.

# 2. Installation

Install the released version from CRAN:

```{r eval=FALSE}
install.packages("varequal")
```

Or install the development version from GitHub:

```{r eval=FALSE}
if (!require("pak")) install.packages("pak")
pak::pak("P10911004-NPUST/varequal")
```

Then load the package:

```{r}
library(varequal)

# For reproducibility
set.seed(123)
```

# 3. Quick start

For a quick overall assessment, use `is_var_equal()`.

The package includes the `roGFP` dataset, which contains redox-index
measurements from three experimental batches. Each batch contains a response
variable `ro` and a grouping variable `grp`.

```{r}
data("roGFP")

df1 <- roGFP[[1]]

is_var_equal(df1, ro ~ grp)
```

The result is a logical value:

```{r}
is_var_equal(df1, ro ~ grp)
```

If you want to inspect the component tests, request a summary:

```{r}
is_var_equal(df1, ro ~ grp, summary = TRUE)
```

For explicit test selection, use `check_normality()`:

```{r}
check_var_equal(df1, ro ~ grp, method = "LV")
```

**If you only require a quick check for data heteroscedasticity, use `is_var_equal` and feel free 
to skip the rest of this guide.**

**If you only require a high-level decision about variance homogeneity, start with `is_var_equal()` 
and feel free to skip the rest of this guide. If you need to justify a particular statistical
procedure, use `check_var_equal()` or one of the individual test functions.**

# 4. Understanding homogeneity-of-variance testing

Most procedures in this package test the null hypothesis

$$ H_0:\sigma_1^2 = \sigma_2^2 = \cdots = \sigma_k^2 $$

against an alternative in which at least one group variance differs:

$$ H_1:\text{not all }\sigma_i^2\text{ are equal}. $$

The interpretation follows the usual hypothesis-testing framework.

A small p-value provides evidence against homogeneity of variance. A large p-value means that the 
test did not detect sufficient evidence of unequal variances.

A large p-value does **not** prove that all population variances are exactly equal. The result 
should be interpreted together with:

- the sample size in each group;
- the presence or absence of outliers;
- the distributional shape of the response;
- the balance of the experimental design;
- the scientific magnitude of the variance differences; and
- the assumptions of the downstream analysis.

Different tests have different sensitivities. In particular, Bartlett's test can be very 
sensitive to non-normality, whereas robust and rank-based procedures are generally less 
affected by departures from normality.

# 5. Data input

The package accepts either:

1. a data frame together with a formula; or
2. a list of numeric vectors representing groups.

## 5.1 Data frame and formula

The standard interface is:

```text
response ~ group
```

For example:

```{r}
Levene_test(df1, ro ~ grp)
```

The response is the dependent variable and the right-hand side identifies the
independent grouping variable.

A data frame containing several experimental variables can therefore be used
directly:

```{r}
check_var_equal(df1, ro ~ grp, method = "BF")
```

## 5.2 A list of numeric vectors

A list can be used when the observations are already separated into groups.

```{r}
x <- list(A = rnorm(12), B = rnorm(12), C = rnorm(12))
is_var_equal(x)
```

# 6. The individual tests

## 6.1 Ansari--Bradley test

`Ansari_Bradley_test()` is a rank-based procedure for assessing equality of dispersion.

```{r}
ab <- Ansari_Bradley_test(df1, ro ~ grp)
```

The implementation uses an approximate chi-square statistic. The package documentation notes that 
results may be inaccurate for very small total sample sizes or when there are many tied values. 
The test is useful as a nonparametric comparison of dispersion, but it should not automatically be 
treated as interchangeable with a robust ANOVA-based homogeneity test.

## 6.2 Bartlett's test

`Bartlett_test()` implements the classical Bartlett test.

Bartlett's test is the most powerful test when no outliers exist and the observations follow normal
distribution, but it is sensitive to departures from normality and to outliers. 

For data that satisfy the normality assumption reasonably well and are free of important outliers, 
Bartlett's test is an appropriate classical procedure.

```{r}
bl <- Bartlett_test(df1, ro ~ grp, summary = TRUE)
```

## 6.3 Brown--Forsythe test

`Brown_Forsythe_test()` provides both the original Brown--Forsythe procedure and the 
Mehrotra-modified version.

The default method is `"MBF"`:

```{r}
mbf <- Brown_Forsythe_test(df1, ro ~ grp)
```

The original Brown--Forsythe procedure can be selected with `method = "BF"`:

```{r}
bf <- Brown_Forsythe_test(df1, ro ~ grp, method = "BF")
```

The Brown--Forsythe approach is a robust modification of Levene's procedure that uses a robust 
location measure, conventionally the group median.

The package also permits a custom transformation of the response. The default is based on absolute 
deviations from the group median:

$$
y'_{ij} = |y_{ij} - \tilde y_i|.
$$

For example:

```{r}
Brown_Forsythe_test(df1, ro ~ grp, transform = function(x) abs(x - mean(x)))
```

The `method` argument distinguishes the two implementations:

- `"BF"`: original Brown--Forsythe procedure;
- `"MBF"`: Mehrotra modification, which adjusts the degrees of freedom.

The modified procedure is designed to reduce Type I error inflation, with a
potential trade-off in power.

## 6.4 Fligner--Killeen test

`Fligner_Killeen_test()` is a rank-based procedure that is particularly useful
when robustness to non-normality and outliers is important.

```{r}
fk <- Fligner_Killeen_test(df1, ro ~ grp)
```

The procedure works with ranks of absolute deviations from the group median and uses a normal-score 
transformation. `Fligner--Killeen` is an important general-purpose choice when the normality 
assumption is questionable.

## 6.5 't Lam's G test

`Lam_G_test()` implements 't Lam's G procedure. The procedure is related to Cochran-type 
variance-outlier assessment rather than being a conventional homogeneity test in the same 
sense as Levene's or Bartlett's test.

```{r}
lg <- Lam_G_test(df1, ro ~ grp)
```


It supports three alternatives:

```{r}
Lam_G_test(df1, ro ~ grp, alternative = "two.sided")
Lam_G_test(df1, ro ~ grp, alternative = "greater")
Lam_G_test(df1, ro ~ grp, alternative = "less")
```

The interpretation is directional:

- `"greater"` focuses on unusually large variances;
- `"less"` focuses on unusually small variances;
- `"two.sided"` considers both directions.

The method is highly sensitive to outliers, so it should be used with an explicit understanding of 
that property.

## 6.6 Levene's test

`Levene_test()` provides an ANOVA-based test of homogeneity.

```{r}
lv <- Levene_test(df1, ro ~ grp)
```

The implementation uses a transformation of each observation into a measure of within-group 
deviation. Its default transformation is the absolute deviation from the group median:

$$ y'_{ij} = |y_{ij} - \tilde y_i|. $$

A custom transformation can be supplied:

```{r}
Levene_test(df1, ro ~ grp, transform = function(x) (x - median(x)) ^ 2)
```

Other useful transformations include:

```{r}
# Absolute deviations from the mean
Levene_test(df1, ro ~ grp, transform = function(x) abs(x - mean(x)))

# Square-root absolute deviations
Levene_test(df1, ro ~ grp, transform = function(x) sqrt(abs(x - median(x))))
```

The transformation is subsequently analyzed using an ANOVA-style decomposition.

## 6.7 O'Brien's test

`O.Brien_test()` implements O'Brien's variance-homogeneity procedure.

```{r}
ob <- O.Brien_test(df1, ro ~ grp)
```

The default transformation is based on squared deviations from the group median:

$$ y'_{ij} = (y_{ij} - \tilde y_i)^2. $$

O'Brien's test can be conservative and may have relatively low power for some forms of 
heteroscedasticity. The package documentation therefore places Levene-type and Brown--Forsythe 
procedures ahead of O'Brien's test for general homogeneity assessment.

## 6.8 O'Neill--Mathews test

`O.Neill_Mathews_test()` implements the weighted least-squares approach to Levene's test.

```{r}
om <- O.Neill_Mathews_test(df1, ro ~ grp)
```

Its default transformation is again based on absolute deviations from the group median:

$$ y'_{ij} = |y_{ij} - \tilde y_i|. $$

This procedure is useful when a weighted least-squares formulation is desired, particularly when 
group sizes are unequal.

# 7. Selecting a test with `check_var_equal()`

`check_var_equal()` provides a common wrapper around the individual tests. The available method 
codes are:

| Code  | Procedure                 |
| :---  | :-------                  |
| `AB`  | Ansari--Bradley           |
| `BL`  | Bartlett                  |
| `FK`  | Fligner--Killeen          |
| `LG`  | 't Lam's G                |
| `LV`  | Levene                    |
| `MBF` | Mehrotra--Brown--Forsythe |
| `BF`  | Brown--Forsythe           |
| `OB`  | O'Brien                   |
| `OM`  | O'Neill--Mathews          |

This is particularly useful when the method is selected programmatically. For example:

```{r}
method <- "FK"
result <- check_var_equal(df1, ro ~ grp, method = method, silent = TRUE)
result
```

# 8. Automatic assessment with `is_var_equal()`

`is_var_equal()` combines five procedures:

- Brown--Forsythe;
- Fligner--Killeen;
- 't Lam's G;
- Levene;
- O'Neill--Mathews.

The `sensitivity` argument ranges from 1 to 5. A larger value requires more component tests to 
agree that the variances are equal:

```{r}
for (s in 1:5) {
  cat("sensitivity =", s, "\n")
  print(is_var_equal(df1, ro ~ grp, sensitivity = s))
}
```

Conceptually:

- `sensitivity = 1` is permissive;
- `sensitivity = 3` is the default;
- `sensitivity = 5` is conservative.

Set `summary = TRUE` to obtain the component-test results:

```{r}
result <- is_var_equal(df1, ro ~ grp, summary = TRUE)
result$is_var_equal
result$summary
```

The summary contains the individual test decisions, p-values, test statistics, and critical values.

# 9. Comparing the tests

No single homogeneity-of-variance test is uniformly optimal.

A practical comparison is:

| Situation | Useful starting point |
|:---|:---|
| Approximately normal data, no important outliers | Bartlett |
| Normality is uncertain | Fligner--Killeen |
| Outliers may be present | Brown--Forsythe or Fligner--Killeen |
| A conventional Levene framework is desired | Levene |
| Unequal group sizes with weighted least squares | O'Neill--Mathews |
| Variance-outlier behavior is of interest | 't Lam's G |
| Explicit comparison with O'Brien's transformation | O'Brien |
| Rank-based dispersion comparison | Ansari--Bradley |

These are guidelines rather than universal rules. The final choice should depend on the data and 
the assumptions of the downstream analysis.

# 10. Why multiple tests can disagree

Different procedures weight distributional features differently.

For example, a dataset can contain:

- approximately equal variances but strong skewness;
- a single extreme observation in one group;
- unequal variances that differ primarily in one group;
- several small groups with low statistical power; or
- unequal sample sizes.

Under such conditions, two valid tests can produce different decisions.

Therefore, the question should not simply be:

> Which test gives a significant *p*-value?

A better question is:

> Which variance-homogeneity procedure is appropriate for the distribution,
> sample size, outlier structure, and experimental design?

`is_var_equal()` is intended to provide a practical aggregate assessment, while the individual 
functions allow the analyst to inspect the behavior of specific procedures.

# 11. Inspecting the data before testing

Formal tests should not replace graphical inspection.

For example, boxplots can reveal differences in spread and possible outliers:

```{r, fig.width=6}
text = sprintf("Variance equal: %s", is_var_equal(df1, ro ~ grp))
boxplot(ro ~ grp, 
        data = df1, 
        main = text, 
        horizontal = TRUE, 
        xlab = "Redox index", 
        ylab = "Group")
points(x = df1$ro, 
       y = jitter(as.numeric(df1$grp), amount = 0.15))
```

The same graphical strategy can be used with the `CYCB1` dataset.

```{r, fig.width=6}
data("CYCB1")
df2 <- CYCB1[[2]]
text = sprintf("Variance equal: %s", is_var_equal(df2, cells ~ grp))
boxplot(cells ~ grp,
        data = df2,
        main = text,
        horizontal = TRUE,
        xlab = "Cell number",
        ylab = "Group"
)
points(x = df2$cells, 
       y = jitter(as.numeric(df2$grp), amount = 0.15))
```

# 12. Built-in datasets

## 12.1 `roGFP`

`roGFP` is a list containing three experimental batches. Each data frame
contains:

- `TEMP`: air temperature;
- `RGF1`: RGF1 peptide concentration;
- `treatment`: combined treatment;
- `grp`: group label;
- `ro`: reduced--oxidized redox index.

The redox index ranges from -1 (reduced) to 1 (oxidized).

```{r}
data("roGFP")
str(roGFP[[1]])
```

The dataset is based on measurements from *Arabidopsis thaliana* roots.

## 12.2 `CYCB1`

`CYCB1` is a list containing three experimental batches of meristematic root cell counts. 
Each data frame contains:

- `TEMP`: air temperature;
- `RGF1`: RGF1 peptide concentration;
- `treatment`: combined treatment;
- `grp`: group label;
- `cells`: number of meristematic root cells.

```{r}
data("CYCB1")
str(CYCB1[[1]])
```

These datasets are included as practical examples for applying variance homogeneity tests to 
experimental biological data.

# 13. Benchmarking and method performance

For comparing variance-homogeneity procedures, the benchmark uses four possible outcomes:

- **O_O**: variance is homogeneous and the method predicts homogeneous;
- **O_X**: variance is homogeneous but the method predicts heterogeneous;
- **X_X**: variance is heterogeneous and the method predicts heterogeneous;
- **X_O**: variance is heterogeneous but the method predicts homogeneous.

From these quantities:

$$ \text{Type I error} = \frac{O\_X}{O\_O + O\_X}\times100\%, $$

$$ \text{Type II error} = \frac{X\_O}{X\_X + X\_O}\times100\%, $$

$$ \text{Accuracy} = \frac{O\_O + X\_X} {O\_O + O\_X + X\_X + X\_O}\times100\%. $$

## Scenario 1: normal data, moderate sample size

For normally distributed, outlier-free data with group sizes from 8 to 20:

| Test | Type I error | Type II error | Accuracy |
|---|---:|---:|---:|
| Ansari–Bradley | 2.6% | 25.1% | 86.15% |
| Bartlett | 6.1% | 1.5% | 96.20% |
| Modified Brown–Forsythe | 2.3% | 23.3% | 87.20% |
| Brown–Forsythe | 2.7% | 17.7% | 89.80% |
| Fligner–Killeen | 2.4% | 17.6% | 90.00% |
| 't Lam's G | 8.1% | 1.8% | 95.05% |
| Levene | 2.7% | 13.6% | 91.85% |
| O'Brien | 1.7% | 37.5% | 80.40% |
| O'Neill–Mathews | 1.7% | 18.0% | 90.15% |

This scenario illustrates the advantage of Bartlett's test when its normality assumptions are 
satisfied. Its low Type II error produces high accuracy in this benchmark.

## Scenario 2: normal data, small sample size

For group sizes from 3 to 7, the reported Type II errors become much larger for many procedures. 
For example:

| Test | Type I error | Type II error | Accuracy |
|---|---:|---:|---:|
| Ansari–Bradley | 3.3% | 89.6% | 53.55% |
| Bartlett | 4.3% | 56.0% | 69.85% |
| Modified Brown–Forsythe | 0.9% | 93.9% | 52.60% |
| Brown–Forsythe | 1.4% | 92.5% | 53.05% |
| Fligner–Killeen | 2.2% | 90.3% | 53.75% |
| 't Lam's G | 16.0% | 39.3% | 72.35% |
| Levene | 2.7% | 87.9% | 54.70% |
| O'Brien | 0.4% | 97.0% | 51.30% |
| O'Neill–Mathews | 0.6% | 94.2% | 52.60% |

The central lesson is not that one test is universally best. Rather, small samples can make 
variance testing intrinsically difficult. A non-significant result should therefore be interpreted 
cautiously when group sizes are very small.

## Scenario 3: normal data with outliers

The benchmark also considers moderate sample sizes with one or two outliers. In that setting, the 
Bartlett and 't Lam's G can become extremely sensitive to the introduced outliers, while robust 
procedures can maintain substantially better Type I error control.

For example, the reported results include:

| Test | Type I error | Type II error | Accuracy |
|---|---:|---:|---:|
| Bartlett | 100.0% | 0.0% | 50.00% |
| Modified Brown–Forsythe | 0.0% | 41.0% | 79.50% |
| Brown–Forsythe | 0.0% | 25.5% | 87.25% |
| Fligner–Killeen | 2.6% | 9.0% | 94.20% |
| Levene | 0.0% | 21.0% | 89.50% |

## Practical interpretation

The benchmark supports three general principles:

1. **Normal, clean data:** Bartlett can be highly effective.
2. **Non-normal or contaminated data:** robust procedures such as Brown–Forsythe and Fligner–Killeen 
are often preferable.
3. **Very small samples:** all tests can have poor ability to detect genuine heteroscedasticity, so 
graphical and substantive assessment becomes especially important.

The benchmark is a simulation study rather than a universal ranking of methods. Results depend on the exact simulation design, and performance under a user's data-generating process can differ.

<br>

# Summary

The **`varequal`** package provides a unified R interface for assessing homogeneity of variance 
across independent groups.

The principal functions are:

- `is_var_equal()` for an aggregate assessment based on multiple tests;
- `check_var_equal()` for selecting a specific procedure by method code;
- `Brown_Forsythe_test()` for a robust Levene-type approach;
- `Fligner_Killeen_test()` for a robust rank-based procedure;
- `Levene_test()` for an ANOVA-based approach;
- `Bartlett_test()` for classical normal-theory testing;
- `Ansari_Bradley_test()` for rank-based dispersion assessment;
- `Lam_G_test()` for variance-outlier-oriented assessment;
- `O.Brien_test()` for O'Brien's variance test; and
- `O.Neill_Mathews_test()` for a weighted least-squares Levene-type approach.

No single test should be treated as universally best. The appropriate method
depends on distributional assumptions, outliers, sample size, group balance,
and the statistical model that follows.

<br>

# References

Ansari, A. R., & Bradley, R. A. (1960). Rank-sum tests for dispersions.
*The Annals of Mathematical Statistics*, 31, 1174--1189.

Bartlett, M. S. (1937). Properties of sufficiency and statistical tests.
*Proceedings of the Royal Society of London. Series A*, 160(901), 268--282.
<https://doi.org/10.1098/rspa.1937.0109>

Brown, M. B., & Forsythe, A. B. (1974). Robust tests for the equality of variances. 
*Journal of the American Statistical Association*, 69(346), 364--367.
<https://doi.org/10.1080/01621459.1974.10482955>

Conover, W. J., Johnson, M. E., & Johnson, M. M. (1981). A comparative study of tests for 
homogeneity of variances, with applications to the outer continental shelf bidding data. 
*Technometrics*, 23, 351--361.
<https://doi.org/10.1080/00401706.1981.10487680>

Fligner, M. A., & Killeen, T. J. (1976). Distribution-free two-sample tests for scale. 
*Journal of the American Statistical Association*, 71(353), 210--213. <https://doi.org/10.1080/01621459.1976.10481517>

Mehrotra, D. V. (1997). Improving the Brown--Forsythe solution to the generalized Behrens--Fisher 
problem. *Communications in Statistics--Simulation and Computation*, 26, 1139--1145.
<https://doi.org/10.1080/03610919708813431>

O'Brien, R. G. (1981). A simple test for variance effects in experimental designs. 
*Psychological Bulletin*, 89(3), 570--574.
<https://doi.org/10.1037/0033-2909.89.3.570>

O'Neill, M. E., & Mathews, K. (2000). Theory & Methods: A Weighted Least Squares Approach to 
Levene's Test of Homogeneity of Variance.
*Australian & New Zealand Journal of Statistics*, 42(1), 81--100.
<https://doi.org/10.1111/1467-842X.00109>

Sharma, D., & Kibria, B. M. G. (2013). On some test statistics for testing homogeneity of variances: 
A comparative study.  *Journal of Statistical Computation and Simulation*, 83, 1944--1963.
<https://doi.org/10.1080/00949655.2012.675336>

't Lam, R. U. E. (2010). Scrutiny of variance results for outliers: Cochran's test optimized. 
*Analytica Chimica Acta*, 659(1--2), 68--84.
<https://doi.org/10.1016/j.aca.2009.11.032>

Zhou, Y., Zhu, Y., & Wong, W. K. (2023). Statistical tests for homogeneity of variance for clinical 
trials and recommendations. *Contemporary Clinical Trials Communications*, 33, 101119.
<https://doi.org/10.1016/j.conctc.2023.101119>
