Saved graphics

When visstat() is called with graphicsoutput specified, plot() lists the generated file paths instead. All generated graphics can be saved in any file format supported by Cairo() (Urbanek and Horner 2025), including “png”, “jpeg”, “pdf”, “svg”, “ps”, and “tiff”. If plotName is provided, the main result plot uses this name. The assumption-diagnostic plot adds the prefix "glm_assumptions_". If plotName is not provided, file names are generated from the selected plot type and the input variable names.

General linear model definition

In the general linear model, a response \(Y\) is modelled as a linear combination of \(k-1\) predictors \(x_j\). The general linear model for observation \(i,\;i = 1, \ldots, N\) is then

\[\begin{equation} (\#eq:glm) Y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_{k-1} x_{i,k-1} + \varepsilon_i, \end{equation}\]

where \(Y_i\) is the response for observation \(i\), \(x_{ij}\) is the value of predictor \(j\) for observation \(i\), \(\beta_0, \beta_1, \ldots, \beta_{k-1}\) are the \(k\) parameters, and \(\varepsilon_i\) is the model error term assumed to be independent and normally distributed with expectation 0 and constant variance \(\sigma^2\), in short \[\varepsilon_i \sim \mathscr{N}(0, \sigma^2), \quad\mathrm{mutually\; independent}. \]

The variance \(\sigma^2\) represents the variation of the data about the regression, \(\operatorname{Var}(Y_i)=\operatorname{Var}(\varepsilon_i)=\sigma^2\), as both the (unknown) model parameters and predictors are not random.

From Eq. @ref(eq:glm), the special cases used by visstat() follow from the predictor structure:

Student’s t-test uses one binary indicator variable \(x_{i1}\), with \(x_{i1}=0\) for group 1 and \(x_{i1}=1\) for group 2. Let \(\mu_1\) and \(\mu_2\) denote the expected mean values in the two population groups. For the expected values of the response, we then obtain \[E(Y_i \mid x_{i1}=0)=\mu_1=\beta_0\] for group 1 and \[E(Y_i \mid x_{i1}=1)=\mu_2=\beta_0+\beta_1\] for group 2. Testing \(H_0: \beta_1 = 0\) is therefore equivalent to testing \(H_0: \mu_1 = \mu_2\).

Fisher’s ANOVA generalises this coding to \(k-1\) binary indicator variables for \(k\) groups; testing \(H_0: \beta_1 = \cdots = \beta_{k-1} = 0\) is equivalent to testing equality of \(k\) group population means.

Simple linear regression uses one continuous predictor; \(H_0: \beta_1 = 0\) examines whether a linear relationship exists.

General linear model assumption tests

Normality tests The normality of the standardised residuals is formally assessed using both the Shapiro–Wilk (SW) test (Shapiro and Wilk 1965; Royston 1982; Royston 1995) (shapiro.test(); Eq. @ref(eq:shapiro-w)) and the Anderson–Darling test (Anderson and Darling 1952) (ad.test(); Eq. @ref(eq:anderson-a2)). These tests offer complementary strengths: Anderson–Darling is highly sensitive to tail deviations in larger samples (Yap and Sim 2011), while Shapiro–Wilk generally exhibits greater power across non-normal distributions in small samples. The Shapiro–Wilk test is used as the normality gate in the automated test selection (Section @ref(sec:route-1)).

Homoscedasticity tests For grouped central-tendency analyses, variance homogeneity of standardised residuals (Cook and Weisberg 1982) is assessed using the package-implemented mean-centred Levene test (Levene 1960) (levene.test(); Eq. @ref(eq:levene-f)) and Bartlett’s test (Bartlett 1937) (bartlett.test(); Eq. @ref(eq:bartlett-k2)).

Bartlett’s test is powerful under normality but sensitive to non-normality; Levene-type tests trade some power for greater robustness when distributions depart from normality (Brown and Forsythe 1974).

Therefore, Levene’s test is used as the variance gate in the automated workflow.

For simple linear regression, group-based variance tests are not applicable. There, visStatistics uses its package implementation bp.test() of the Breusch–Pagan test (Breusch and Pagan 1979) (Eq. @ref(eq:breusch-pagan-bp)) on raw residuals (Schützenmeister et al. 2012).

Visualisation of the assumptions of the general linear model

Since algorithmic logic based on p-values of assumption tests cannot replace expert visual judgment, visstat() visualises the assumptions of the underlying linear model for the selection of tests of central tendency (Route 1) and for simple linear regression (correlation = FALSE) (Route 3).

For numeric responses with categorical predictors (Route 1), the diagnostic panel displays the residual histogram, the normal Q–Q plot with simultaneous tolerance band (STB) and point-wise tolerance band (TB) (Schützenmeister et al. 2012), and the absolute standardised residuals \(|r_i|\) (Eq. @ref(eq:standardised)) by group. The last panel shows whether residual spread is comparable across factor levels, the pattern assessed formally by the Levene (Eq. @ref(eq:levene-f)) and Bartlett (Eq. @ref(eq:bartlett-k2)) variance checks.

For Route 3 (simple linear regression), the first two diagnostic panels follow the layout of Route 1, whereas the third panels displays z-scaled residuals versus fitted values.

The first row of the outer tile of the diagnostic plot reports p-values of residual-normality checks with the Shapiro–Wilk test and Anderson–Darling tests. The second row reports p-values of variance checks: Levene and Bartlett for grouped central-tendency analyses, or Breusch–Pagan for simple regression.

Note that among the displayed assumption tests, only the Shapiro–Wilk and Levene test results enter automated routing, and only in the central-tendency branch (see Section @ref(sec:route-1)). Anderson–Darling, Bartlett, and Breusch–Pagan are diagnostic output only.

The Route 1 and Route 3 diagnostic-panel designs are illustrated in the examples in Figures @ref(fig:welch-anova-example), left, and @ref(fig:regression-example), left.

Type I simulations: equal means, equal ranks, balanced or unbalanced group sizes

The Type I simulations addresses homo- and heteroscedastic settings: Figure @ref(fig:route1-identical-typeI) adapts the three-group zero-effect one-way ANOVA setting of Blanca et al. (Blanca et al. 2017) to four Route 1 groups with identical distributions; both the equal-means null and the Kruskal–Wallis null are true. Figure @ref(fig:route1-unequal-typeI) keeps the equal-means null and varies the group standard deviations and sample-size pairings as in Delacre et al. (Delacre et al. 2019). For F, W, and L, all columns are therefore Type I checks of the equal-means null. For KW, SW, and SW+L, the first two symmetric columns remain Type I checks without a median or mean-rank shift, whereas the skewed columns change the group medians and mean ranks.

The simulations use Fleishman’s polynomial transformation \(Y=a+bX+cX^2+dX^3\), where \(X\sim N(0,1)\) and the coefficients \(a,b,c,d\) are chosen to obtain mean 0, variance 1, and the target skewness and excess kurtosis (Fleishman 1978).

The Type I simulations use four groups and \(B=50{,}000\) Monte Carlo replications per cell. Across rows, the balanced designs vary the common group size \(n \in \{10,20,50,100\}\), whereas he unbalanced designs vary the target mean group size \(\bar n = n\). Within each unbalanced row, the group-size vector is \(\mathbf n=\bar n\cdot(0.5,0.8,1.2,1.5)\), rounded up component-wise. For an estimated rejection probability \(\hat p\), the Monte Carlo standard error is (Koehler et al. 2009).

\[SE_{\mathrm{MC}}=\sqrt{\hat p(1-\hat p)/B}.\] This is about 0.10 percentage points at \(\hat p=0.05\), with a maximum of about 0.22 percentage points at \(\hat p=0.50\).

Equal mean and equal ranks

In the simulations shown in Fig. @ref(fig:route1-identical-typeI) , all four groups are sampled from the same input distribution, so the equal-means null and the Kruskal–Wallis null are both true. The columns show a normal distribution and four Fleishman distributions with prespecified skewness and excess kurtosis. Panels B and C only differ in sample sizes. Across rows, panel B varies the common group size \(n \in \{10,20,50,100\}\), whereas panel C varies the target mean group size \(\bar n \in \{10,20,50,100\}\). Within panel C, the group-size vector is \(\mathbf n=\bar n\cdot(0.5,0.8,1.2,1.5)\), rounded up component-wise. All group SDs are \(1\), so the group distributions remain identical within each column. :::

Route 1 Type I simulation under identical distributions and identical means, with \(\mathrm{mean}_i = 0\) and \(SD_i = 1\) (\(i = 1, \ldots, 4\)). (A) input distributions. (B) balanced design. (C) unbalanced design. The heatmaps in (B) and (C) report final-test rejection rates at alpha = 5%. Green marks Bradley’s 2.5%–7.5% interval. Cells for SW and SW+L additionally report the route split. In (A), dashed lines mark means and dotted lines mark medians. The index i = 1, …, 4 denotes groups; \(\bar n\) denotes the target mean group size for unbalanced designs. All heatmap numbers are percentages; the first value is the final-test rejection rate, and gated rows list route splits after |.

In this setting, routing to Kruskal–Wallis under non-normality does not change the truth status of the tested null because the equal-rank-distribution null is also true. The expected result is therefore close to nominal Type I error for all strategies. Deviations from the nominal level reflect finite-sample and Monte Carlo behaviour, not a conflict between the mean and rank nulls. This is the cleanest Type I check for the full automatic route.

This is the clean Type I check for the full automatic route. Routing to Kruskal–Wallis under non-normality does not change the truth status of the tested null, because the Kruskal–Wallis null is also true. SW and SW+L therefore differ only by the added Levene gate inside the mean branch. In these cells, their final rejection rates are almost identical: SW ranges from 4.4% to 5.6%, and SW+L from 4.5% to 5.7%.

Equal means, unequal ranks

In Fig. @ref(fig:route1-unequal-typeI), all four group means remain zero, but the common standardised input distribution is multiplied by group-specific SD scale factors. The balanced block (panel B) uses \(\mathbf n=(n,n,n,n)\) and \(\mathbf s=(1.0,1.3,1.7,2.2)\). The unbalanced blocks use \(\mathbf n=\bar n(0.5,0.8,1.2,1.5)\); this group-size vector is paired either with \(\mathbf s=(1.0,1.3,1.7,2.2)\), so larger groups have larger SDs, or with \(\mathbf s=(2.2,1.7,1.3,1.0)\), so larger groups have smaller SDs.

The parametric equal-means null is true in all columns. For the first two symmetric columns, the group means and medians are equal, and SD scaling leaves the group mean ranks equal. In the skewed columns, SD scaling also changes the group medians, so Kruskal–Wallis can reject because the group rank distributions are no longer aligned. Therefore the first two columns are Type I checks for both the parametric tests and Kruskal–Wallis, whereas the skewed columns are Type I checks only for the parametric equal-means tests.

Route 1 simulation: identical means and unequal distributions; \(\mathrm{mean}_i = 0\) (\(i = 1, \ldots, 4\)). Route 1 equal-means simulation with varied group standard deviations and sample-size pairings. (A) input distributions. (B) balanced design. (C) unbalanced design with larger \(n_i\) paired with larger SD. (D) unbalanced design with larger \(n_i\) paired with smaller SD. F, W, and L are mean-null tests; SW and SW+L report final-test rejection after possible routing to Kruskal–Wallis. In (A), dashed lines mark means and dotted lines mark medians. The index i = 1, …, 4 denotes groups; \(\bar n\) denotes the target mean group size for unbalanced designs. All heatmap numbers are percentages; the first value is the final-test rejection rate, and gated rows list route splits after |.

Thus, F, W, and L are still evaluated against a true equal-means null, whereas SW and SW+L may switch to a final test whose null is not the same. High rejection after a switch to Kruskal–Wallis is therefore not a Type I error for the equal-means null. It is the rejection rate of the final selected test after the route has changed the hypothesis. This panel is the stress test for interpreting routed procedures, not a simple comparison of mean tests.

This figure is therefore read column by column. For F, W, and L, every cell is a Type I check for the equal-means null. For KW, SW, and SW+L, the first two symmetric columns are also Type I checks, but the skewed columns are not pure Type I checks for the Kruskal–Wallis route because the scaling changes the group medians and mean ranks. The printed route percentages show how often SW and SW+L leave the mean-based branch.

The power simulation uses the same five fixed input distributions but adds ordered location shifts across the four groups. It uses balanced groups with equal SD, so \(n_i=n\), \(\mathrm{SD}_i=1\), and group means are shifted by \(0,0.25,0.50,0.75\). Because the groups differ by location, the mean-based tests and Kruskal–Wallis all have a true alternative, but they target different summaries: means for F, W, and L, and ranks/location for KW. The main visible result is how often SW+L routes to Fisher, Welch, or Kruskal–Wallis under each input distribution.

The power simulation uses the same five fixed input distributions and adds ordered location shifts across the four groups. It uses balanced groups with equal SD, so \(n_i=n\), \(\mathrm{SD}_i=1\), and group means are shifted by \(0,0.25,0.50,0.75\). The figure therefore asks whether the second gate changes the rejection rate relative to SW, and how often SW+L routes to Fisher, Welch, or Kruskal–Wallis under each input distribution.

Route 1 power simulation with Fleishman input distributions. (A) Input distributions with group mean and median reference lines. (B) Simulated rejection rates for the six testing strategies.

In the normal input-distribution panel, the mean-based tests have slightly higher rejection rates than Kruskal–Wallis at small and moderate \(n\). In the heavy-tailed symmetric panel and in the skewed panels, Kruskal–Wallis rejects more often under the imposed location-shift alternative, and SW+L increasingly follows it because the Shapiro–Wilk gate routes to KW. This is a power gain only for the rank/location alternative; it is not evidence that the routed procedure is a better equal-means test. Once \(n\) is large enough, all strategies approach 100% rejection, so the relevant information is mainly the route split and the small-\(n\) behaviour.

If Route 1 is read strictly as an equal-means machine, the simulations support Welch as the clean default. If Route 1 is allowed to use input shape as a method switch, SW is the simpler automatic route. SW+L is a model-display compromise: it preserves the Fisher/Student GLM branch when both residual normality and homoscedasticity are retained, but its Type I error cost must be read from the Type I simulations.

Thus the simulations do not prove that SW+L is the best equal-means procedure. They show the price and benefit of the automatic route. The price is that, under equal means but unequal distributions, SW and SW+L often report the rejection rate of KW after the hypothesis has changed. The benefit is that the switch is visible and diagnostic-driven. The final design decision is therefore not statistical optimality in the abstract, but which default the package wants to make explicit: Welch for fixed equal-means inference, SW for a parsimonious shape-based switch, or SW+L for the same switch while preserving the GLM Fisher/Student branch when its displayed assumptions pass.

The empirical cost of adding Levene is small in these simulations. In the clean Type I setting with identical group distributions, SW ranges from 4.4% to 5.6% rejection and SW+L from 4.5% to 5.7%. In the power simulation, the median difference SW+L minus SW is 0.004 percentage points, and the maximum is 2.25 percentage points. Thus the second gate changes the route display much more than it changes the final rejection rate.

This supports using SW+L as the default assumption gate in visstat(). The point is not that the gate is the most powerful standalone test. The point is that Fisher/Student is used only when the two displayed GLM assumptions pass, Welch is used when only homoscedasticity fails, and Kruskal–Wallis is used when the residual-normality gate fails. The route percentages in the figures document these decisions directly.

Route 1: Numeric response, categorical predictor

A numeric response with a categorical predictor with \(k\) “levels” (in the following “groups”) asks whether the response differs between groups.

Figure @ref(fig:decision-tree) expands the routing logic for tests of central tendencies.

Decision tree for the default Route 1 test selection (group_test = NULL). Shapiro–Wilk on model residuals determines whether the route remains mean-based or switches to rank-based tests; the Levene test then selects equal-variance or Welch-type procedures.

A linear model of Eq. @ref(eq:glm) is fitted between the numeric response and the categorical predictor, and the model residuals of Eq. @ref(eq:raw-residual) are extracted.

In the default setting (group_test = NULL), Route 1 uses the displayed residual diagnostics of the Shapiro–Wilk (SW) and Levene test (L) (Levene 1960) as automatic gates:

The Shapiro–Wilk (SW) normality test functions as the residual-normality gate, because simulation studies report high power for small to moderate sample sizes (Razali and Wah 2011).

If the SW-test rejects residual normality (\(p_\text{SW} \le \alpha\)), robust non-parametric tests are selected: wilcox.test() (Eq. @ref(eq:wilcoxon-w)) for two groups, or kruskal.test() (Eq. @ref(eq:kruskal-h)) followed by Holm-adjusted pairwise.wilcox.test() for more than two groups.

If residual normality is not rejected, variance homogeneity the mean-centred Levene test (L) (Levene 1960) (Eq. @ref(eq:levene-f)) gates the variance assessment:

For homoscedastic data (\(p_\text{L} > \alpha\)), t.test(var.equal = TRUE) (Eq. @ref(eq:student-t)) is applied for two groups, or Fisher’s aov() (Eq. @ref(eq:fisher-f)) for more than two groups. For heteroscedastic data (\(p_\text{L} \le \alpha\)), Welch’s t.test() (Eq. @ref(eq:welch-t)) is applied for two groups, or Welch’s oneway.test() (Eq. @ref(eq:welch-f)) for more than two groups.

A type I error of the Levene test merely routes homoscedastic data to Welch’s t-test or Welch’s one-way ANOVA, which lose only negligible power when variances are in fact equal simple sizes are large. (Rasch et al. 2011; Delacre et al. 2019).

Route 2: Ordered response

An ordered categorical response with a categorical predictor or ordered categorical predictor is treated as a rank-based group comparison. The ordered response is converted to integer level codes and analysed with the Wilcoxon rank-sum test for two groups or the Kruskal–Wallis test for more than two groups.

Route 3: Numeric response, numeric predictor

Two numeric variables ask whether a numeric response changes with a numeric predictor. By default, visstat() fits a simple linear regression (Eq. @ref(eq:simple-regression-fit)). and the diagnostic panel described in Section @ref(sec:graphical) is displayed. If general linear model assumptions are violated, the corresponding p-values trigger warnings and recommendations, but no automatic model replacement. The regression output is shown in Section @ref(sec:lin-reg).

Route 4: Two unordered factors

Two unordered factors ask whether two categorical variables are independent. visstat() uses Pearson’s \(\chi^2\) test or Fisher’s exact test, depending on expected cell counts following Cochran’s rule (Cochran 1954): the \(\chi^2\) approximation is used if no expected cell count is less than 1 and no more than 20% of cells have expected counts below 5. Yates’ continuity correction is applied by default to \(2 \times 2\) tables when the \(\chi^2\) approximation is used.

Optional rank-correlation mode

The four routes above describe the default, automatic test selection behaviour. For ordered–ordered and numeric–numeric input vectors, the user can instead request a rank-correlation analysis by setting correlation = TRUE.

Both optional analyses test monotone association and are computed by cor.test(): Kendall’s \(\tau_b\) (Eq. @ref(eq:kendall-tau-b)) with method = "kendall", exact = FALSE for two ordered variables, and Spearman’s \(\rho\) (Eq. @ref(eq:spearman-rho)) with method = "spearman" for two numeric variables. Kendall’s \(\tau_b\) corrects for ties present with few ordered levels (Agresti 2010; Xu et al. 2013).

Note that for numeric–numeric input, Pearson correlation is not implemented as a separate optional mode as in simple linear regression with an intercept, the two-sided test of zero slope and the two-sided Pearson correlation test return the same \(p\)-value.

Student’s t-test and Fisher’s one-way ANOVA

Student’s t-test

The ToothGrowth dataset records odontoblast length in 60 guinea pigs given vitamin C by orange juice (OJ) or ascorbic acid (VC). With delivery method as predictor and length as response, the assumption-diagnostic panel shows no residual-normality or variance-homogeneity violation. visstat() therefore selects Student’s t-test, and the result panel shows the two-group box plot with the selected test result.

student_ttest <- visstat(ToothGrowth$supp, ToothGrowth$len)

Student’s t-test applied to the ToothGrowth dataset (len vs. supp). Assumption diagnostics (Shapiro–Wilk does not reject residual normality; Levene does not reject residual variance homogeneity) select the equal-variance mean-based path, followed by box plots with the Student t-test result.

Fisher’s one-way ANOVA with Tukey HSD post-hoc comparisons with methods demonstration

The PlantGrowth dataset records yields (as measured by dried weight of plants) for a control group and two treatment groups. This dataset serves a double purpose: it demonstrates both a branching result and the visstat() S3 methods print(), summary(), and plot() (Section @ref(sec:visstat-methods)).

anova_plantgrowth <- visstat(PlantGrowth$group, PlantGrowth$weight)

In this branch visstat() generates two plots: the assumption-diagnostic panel (which = 1) and the result panel with box plots and post-hoc significance letters (which = 2). Calling plot() without which first lists both available plots:

plot(anova_plantgrowth)

## Plot [1] captured. Use plot(obj, which = 1) to display.

## Plot [2] captured. Use plot(obj, which = 2) to display.

Figure @ref(fig:anova-plantgrowth-panels)(a) replays the assumption-diagnostic panel (which = 1). With control and treatment groups as predictor and plant weight as response, Shapiro–Wilk does not reject normality of the model residuals and levene.test() does not reject homoscedasticity, so visstat() takes the equal-variance mean-based path. Figure @ref(fig:anova-plantgrowth-panels)(b) replays the result panel (which = 2).

plot(anova_plantgrowth, which = 1)
plot(anova_plantgrowth, which = 2)

PlantGrowth Fisher’s one-way ANOVA (weight vs. group). (a) Assumption-diagnostic panel. (b) Result panel with Tukey HSD significance letters (\(\alpha = 0.05\)).

The omnibus F-test is significant at \(\alpha = 0.05\), and the Tukey HSD post-hoc comparison finds no significant difference between the control group and either treatment, but the difference between trt1 and trt2 is significant.

To save the graphics, call visstat() with graphicsoutput; the file paths are returned in the "plot_paths" attribute. Here, plotName is set explicitly so that the output names are stable.

anova_plantgrowth_stored <- visstat(
  PlantGrowth$group,
  PlantGrowth$weight,
  graphicsoutput = "png",
  plotName = "anova_plantgrowth",
  plotDirectory = tempdir()
)
paths <- attr(anova_plantgrowth_stored, "plot_paths")
print(basename(paths))

## [1] "glm_assumptions_anova_plantgrowth.png"
## [2] "anova_plantgrowth.png"

print() lists the returned components:

print(anova_plantgrowth)

## Object of class 'visstat'
## 
## Available components:
## [1] "summary statistics of ANOVA" "post-hoc analysis "         
## [3] "conf.level"                  "effect_size"

summary() prints the full object, including assumption tests, post-hoc comparisons, and effect size.

summary(anova_plantgrowth)

## Summary of visstat object
## 
## --- Named components ---
## [1] "summary statistics of ANOVA" "post-hoc analysis "         
## [3] "conf.level"                  "effect_size"                
## 
## --- Contents ---
## 
## $summary statistics of ANOVA:
##             Df Sum Sq Mean Sq F value Pr(>F)  
## fact         2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $post-hoc analysis :
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = samples ~ fact)
## 
## $fact
##             diff        lwr       upr     p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1  0.865  0.1737839 1.5562161 0.0120064
## 
## 
## $conf.level:
## [1] 0.95
## 
## $effect_size:
## $name
## [1] "omega-squared"
## 
## $estimate
## [1] 0.2040788
## 
## $effect_size_method
## [1] "Omega-squared for one-way ANOVA"

Welch’s t-test and Welch’s one-way ANOVA

Welch’s t-test

The Motor Trend Car Road Tests dataset (mtcars) contains 32 observations, where mpg denotes miles per (US) gallon and am represents the transmission type (0 = automatic, 1 = manual). With binary factor am and continuous response mpg, the assumption-diagnostic panel shows that Shapiro–Wilk does not reject normality of the model residuals, while the Levene test detects heteroscedasticity. The routing therefore leads to Welch’s t-test rather than Student’s t-test, and the result panel shows the corresponding two-group comparison.

mtcars$am <- as.factor(mtcars$am)
t_test_stats <- visstat(mtcars$am, mtcars$mpg)

Welch’s t-test applied to the mtcars dataset (mpg vs. am). Assumption diagnostics (Shapiro–Wilk does not reject residual normality; Levene rejects residual variance homogeneity) select the unequal-variance mean-based path, followed by box plots with the Welch t-test result.

Welch’s heteroscedastic one-way ANOVA with Games–Howell post-hoc comparisons

In the iris dataset, using Species as predictor and Sepal.Length as response, the assumption-diagnostic panel shows that Shapiro–Wilk does not reject normality of the model residuals, whereas the Levene test rejects homoscedasticity at the given \(\alpha = 5\%\). visstat() therefore selects Welch’s heteroscedastic one-way ANOVA (oneway.test()) and applies Games–Howell post-hoc comparisons. The result panel shows the box plots and Games–Howell significance letters.

welch_anova_iris <- visstat(iris$Species, iris$Sepal.Length)

Welch’s heteroscedastic one-way ANOVA applied to the iris dataset (Sepal.Length vs. Species). Assumption diagnostics (Shapiro–Wilk does not reject residual normality; Levene rejects residual variance homogeneity) select the unequal-variance mean-based path, followed by box plots with Games–Howell significance letters (\(\alpha = 0.05\)).

Wilcoxon rank-sum test and Kruskal–Wallis test

Wilcoxon rank-sum test

The warpbreaks dataset records thread breaks during weaving. Using wool type (A or B) as predictor and the number of breaks as response, the assumption-diagnostic panel shows that the Shapiro–Wilk test rejects normality of the model residuals. visstat() therefore selects the Wilcoxon rank-sum test, and the result panel shows the rank-based two-group comparison.

wilcoxon_stats <- visstat(warpbreaks$wool, warpbreaks$breaks)

Wilcoxon rank-sum test applied to the warpbreaks dataset (breaks vs. wool). Assumption diagnostics (Shapiro–Wilk rejects residual normality; non-parametric path selected) and box plots with the Wilcoxon test result.

Kruskal–Wallis rank sum test with pairwise Wilcoxon post-hoc comparisons

In the iris data set, Petal.Width by Species follows a different route than Sepal.Length by Species above (Figure @ref(fig:welch-anova-example)), because the assumption diagnostics differ. The assumption-diagnostic panel shows clear departures from normality, and both normality tests return very small \(p\)-values. Since Shapiro–Wilk falls below \(\alpha\), visstat() switches to kruskal.test() followed by Holm-adjusted pairwise.wilcox.test(). The result panel shows the box plots and Holm-adjusted significance letters; all three species differ significantly in petal width, as indicated by distinct letters.

kruskal_iris <- visstat(iris$Species, iris$Petal.Width)

Kruskal-Wallis test applied to the iris dataset (Petal.Width vs. Species). Assumption diagnostics (Shapiro–Wilk rejects residual normality; non-parametric path selected) and box plots with Holm-adjusted pairwise Wilcoxon significance letters (\(\alpha = 0.05\)).

Ordered response, categorical factor

titanic_df <- counts_to_cases(as.data.frame(Titanic))
titanic_df$Class <- ordered(titanic_df$Class,
  levels = c("1st", "2nd", "3rd", "Crew")
)
wilcox_ordered <- visstat(titanic_df$Sex, titanic_df$Class)

## Warning: Ordered response detected. Converting to integer level codes for
## non-parametric analysis.

Kruskal–Wallis test with ordered response

With three predictor groups, visstat() routes to kruskal.test() followed by Holm-adjusted pairwise.wilcox.test(). The result panel shows the Kruskal–Wallis comparison and Holm-adjusted significance letters on the numeric level scores (see Figure @ref(fig:ordinal-wilcoxon-kruskal-example), right). A synthetic survey records perceived car comfort on a five-point scale across three markets.

set.seed(123)
market <- factor(rep(c("Europe", "North America", "Asia"), each = 50))
comfort_numeric <- c(
  sample(1:5, 50, replace = TRUE, prob = c(0.30, 0.30, 0.20, 0.15, 0.05)),
  sample(1:5, 50, replace = TRUE, prob = c(0.10, 0.20, 0.40, 0.20, 0.10)),
  sample(1:5, 50, replace = TRUE, prob = c(0.05, 0.10, 0.20, 0.35, 0.30))
)
survey_data_3 <- data.frame(
  market = market,
  comfort = ordered(comfort_numeric)
)
kruskal_ordered <- visstat(comfort ~ market, data = survey_data_3)

## Warning: Ordered response detected. Converting to integer level codes for
## non-parametric analysis.

Wilcoxon rank-sum test for ordered passenger class by sex in the expanded Titanic data (left) and its multi-group generalisation, the Kruskal-Wallis test for ordered car comfort ratings by market (right). Holm-adjusted pairwise Wilcoxon post-hoc comparisons are shown as significance letters for the Kruskal-Wallis example (\(\alpha = 0.05\)).

Linear regression

The swiss dataset records standardised fertility and socio-economic indicators for 47 French-speaking Swiss provinces in 1888. We examine how the share of draftees achieving the highest army examination score (Examination) predicts the fertility measure (Fertility), with conf.level = 0.99. The diagnostic panel in Figure @ref(fig:regression-example), left, shows that both normality tests pass and the Breusch–Pagan test confirms homoscedasticity, supporting the linear model. The assumption-diagnostic panel is displayed, but its checks do not trigger automatic model replacement. The regression plot shows the fitted line

\[\begin{equation} \hat{y}_i = b_0 + b_1 x_i (\#eq:simple-regression-fit) \end{equation}\] with the point estimates \(b_0\) and \(b_1\) for the unknown parameters \(\beta_0\) and \(\beta_1\) of the linear regression model in Eq. @ref(eq:glm) with one predictor. It is displayed with pointwise confidence and prediction bands at the specified conf.level.

The returned object contains the regression statistics, residual-normality tests, pointwise confidence and prediction bands, and the coefficient of determination \(R^2\) (Eq. @ref(eq:r-squared)) as effect size.

linreg_swiss <- visstat(swiss$Examination, swiss$Fertility, conf.level = 0.99)

Simple linear regression of Fertility on Examination for the swiss dataset (conf.level = 0.99). Left: residual-diagnostic panel with histogram, normal Q-Q plot with simultaneous tolerance band (STB) and point-wise tolerance band (TB), and residuals versus fitted values. Right: scatter plot with fitted regression line, 99% prediction interval for an individual response, and 99% confidence interval for the mean response.

The airquality ozone example shows the limits of the automated approach when the default linear model is not an adequate final model. visstat() identifies assumption violations and points to analyses outside the automated decision tree. A default visstat() call for ozone concentration (Ozone) as a function of wind speed (Wind) fits the simple linear model.

ozone_lm <- visstat(airquality$Wind, airquality$Ozone)

## Warning: Statistical assumptions violated:
## Normality of residuals violated (Shapiro-Wilk p = 0.00522 )
## Homoscedasticity violated (Breusch-Pagan p = 0.00595 )
## Analysis proceeded but interpret results cautiously.

## RECOMMENDATION: Consider exploring alternatives outside visstat() such as data transformations,
## generalised linear models, or robust regression. For a non-causal alternative
## consider rerunning with correlation = TRUE.

Default simple linear regression for Ozone by Wind in the airquality dataset. Assumption diagnostics flag non-normal model residuals and heteroscedasticity before alternative routes are considered.

The diagnostic output flags non-normal model residuals and heteroscedasticity.

In the “Residual vs. fitted” diagnostic panel we observe an increase in spread from left to right, forming a funnel shape that indicates variance increases with fitted values. The optional Spearman analysis for the same dataset is shown in Section @ref(sec:examples-rank-correlation-mode). The following example shows a Gamma generalised linear model outside visstat().

Model exploration outside visstat()

As a model outside of visstat(), we fit a Gamma generalised linear model with log link. The Gamma family is suited here because Ozone is strictly positive and continuous, and its variance grows with the fitted values — the structure detected by the Breusch–Pagan test. The log link guarantees positive fitted values.

# Gamma model with log mapping
model_gamma <- glm(Ozone ~ Wind, data = airquality, family = Gamma(link = "log"))
model_gamma$aic

## [1] 1040.021

# Comparison with AIC of simple linear regression
model_lm <- glm(Ozone ~ Wind, data = airquality)
model_lm$aic

## [1] 1093.187

Gamma GLM with log link fitted to the airquality dataset Ozone vs. Wind. The red curve shows the fitted Gamma GLM; the y-axis is on a log scale.

For a Gamma generalised linear model with log link, standardised deviance residuals are asymptotically standard normal; we use Shapiro–Wilk and Anderson–Darling as approximate checks of the fitted model:

# Extract standardised deviance residuals
std_dev_res <- rstandard(model_gamma, type = "deviance")
# Validate using the Shapiro-Wilk normality test
shapiro.test(std_dev_res)

## 
##  Shapiro-Wilk normality test
## 
## data:  std_dev_res
## W = 0.99245, p-value = 0.7817

# Validate using the Anderson-Darling normality test
nortest::ad.test(std_dev_res)

## 
##  Anderson-Darling normality test
## 
## data:  std_dev_res
## A = 0.198, p-value = 0.8853

The Gamma model improves the model fit according to the Akaike Information Criterion (Akaike 1974), which decreases from 1093.2 to 1040.0. The increase in the Shapiro–Wilk \(p\)-value from \(p_{SW} = 0.0052\) in the simple linear regression to \(p_{SW} = 0.78\) is more consistent with residual normality. This comparison illustrates how assumption warnings from visstat() can motivate model exploration outside the automated decision tree.

Pearson’s \(\chi^2\) test

For a contingency table with \(R\) response levels and \(C\) predictor levels, Pearson’s \(\chi^2\) test (Eq. @ref(eq:pearson-chi)) shows a grouped column plot of row percentages with the \(p\)-value in the title, followed by a mosaic plot from vcd (Meyer et al. 2006, 2024). Each tile corresponds to one cell of the contingency table. The tile colour represents the Pearson residual value (Eq. @ref(eq:pearson-residual)) on a blue–red colour scale; the tile size reflects the cell count.

With Eye and Hair from HairEyeColor, all expected cell counts exceed the Cochran thresholds (Cochran 1954), so the \(4 \times 4\) \(\chi^2\) approximation is used.

hair_eye_df <- counts_to_cases(as.data.frame(HairEyeColor))
visstat(hair_eye_df$Eye, hair_eye_df$Hair)

Pearson’s \(\chi^2\) test applied to the HairEyeColor dataset. Grouped bar chart of eye colour by hair colour and mosaic plot with tiles coloured by Pearson residuals (blue: over-represented, red: under-represented).

Here, cells for black hair and brown hair, as well as blond hair and blue eyes, show counts above the expectation.

Pearson’s \(\chi^2\) test with Yates’ continuity correction

Restricting HairEyeColor to black or brown hair and brown or blue eyes yields a \(2 \times 2\) table. Cochran’s rule is still satisfied, so visstat() applies Pearson’s \(\chi^2\) test with Yates’ continuity correction. The resulting grouped column plot is shown in Figure @ref(fig:yates-fisher-example), left.

hair_bb_eyes_bb <- HairEyeColor[1:2, 1:2, ]
hair_bb_eyes_bb_df <- counts_to_cases(
  as.data.frame(hair_bb_eyes_bb)
)
yates_stats <- visstat(
  hair_bb_eyes_bb_df$Eye,
  hair_bb_eyes_bb_df$Hair
)

yates_stats$effect_size

## $name
## [1] "phi"
## 
## $estimate
## [1] 0.1709571
## 
## $effect_size_method
## [1] "Phi coefficient for 2 x 2 contingency table"

The returned effect size is \(\phi = 0.17\), which, using Cohen’s benchmarks for \(2 \times 2\) tables (Cohen 2013, 227), is a small association. The p-value instead is below \(\alpha = 0.05\) ($p = 0.0035$) and thus significant. This example underlines the importance of effect sizes: a significant p-value can be accompanied by a small effect size measure.

Fisher’s exact test

Restricting HairEyeColor to male participants with black or brown hair and hazel or green eyes yields a \(2 \times 2\) table where one expected frequency is less than 5, violating Cochran’s rule (Cochran 1954). visstat() therefore applies Fisher’s exact test. The graphical output shows absolute counts with count labels above each bar and the \(p\)-value in the title, so the small cell counts that trigger the exact test remain visible (see Figure @ref(fig:yates-fisher-example), right).

hair_eye_male <- HairEyeColor[, , 1]
black_brown_hazel_green <- hair_eye_male[1:2, 3:4]
black_brown_hazel_green_df <- counts_to_cases(
  as.data.frame(black_brown_hazel_green)
)
fisher_stats <- visstat(
  black_brown_hazel_green_df$Eye,
  black_brown_hazel_green_df$Hair
)

Two \(2 \times 2\) categorical routes in HairEyeColor: Yates-corrected Pearson \(\chi^2\) when Cochran’s rule is satisfied (black/brown hair and brown/blue eyes; left), and Fisher’s exact test when expected counts are too small (male participants, black/brown hair, hazel/green eyes; right). The Yates-corrected plot shows row percentages; the Fisher plot shows absolute counts.

Kendall rank correlation with correlation = TRUE

A hypothetical survey of 150 secondary-school students records alcohol consumption frequency and academic performance on five-point ordinal scales. A negative monotone association is induced by construction: students who consume alcohol more frequently tend to have lower academic performance. The Kendall result is shown in Figure @ref(fig:kendall-spearman-example), left.

set.seed(42)
n <- 150
xs <- sample(1:5, n, replace = TRUE)
ys <- pmin(5, pmax(1, (6 - xs) + sample(-1:1, n, replace = TRUE)))
likert_alc <- c("never", "rarely", "sometimes", "often", "always")
likert_perf <- c("poor", "fair", "ok", "good", "great")
alcohol <- ordered(likert_alc[xs], levels = likert_alc)
performance <- ordered(likert_perf[ys], levels = likert_perf)
kendall_result <- visstat(performance, alcohol, correlation = TRUE)
spearman_air <- visstat(airquality$Wind, airquality$Ozone, correlation = TRUE)

Rank-based correlations: Left: Kendall’s \(\tau_b\) for a hypothetical survey (\(n = 150\)): alcohol consumption frequency vs. academic performance. Right: Spearman rank correlation of Wind and Ozone from the airquality dataset (correlation = TRUE; right). Both plots annotate the corresponding effect measure and \(p\)-value.

Spearman rank correlation with correlation = TRUE

For the ozone example introduced in Section @ref(sec:lin-reg), staying within visstat() with the flag correlation = TRUE gives the Spearman analysis shown in Figure @ref(fig:kendall-spearman-example), right.

Shapiro–Wilk test shapiro.test()

The Shapiro–Wilk test evaluates whether a sample \(x_1,\ldots,x_N\) comes from a normal distribution. Let \(x_{(1)}\le \cdots \le x_{(N)}\) be its order statistics. Introduce a reference sample \(Z_1,\ldots,Z_N\) of independent standard normal random variables, i.e. \(Z_i \sim N(0,1)\) for all \(i\), and let \(Z_{(1)}\le \cdots \le Z_{(N)}\) be their order statistics used to construct the Shapiro–Wilk weights.

Let \(m_i = \operatorname{E}(Z_{(i)})\) and \(v_{ij} = \operatorname{Cov}(Z_{(i)}, Z_{(j)})\) for \(i,j = 1,\ldots,N\). Define \(\mathbf{m} = (m_1,\ldots,m_N)^\top\) and \(V = (v_{ij})_{i,j=1}^N\).

The vector \(\mathbf{m}\) contains the expected standard-normal order statistics, and \(V\) is their covariance matrix. Let \(\mathbf{a}=(a_1,\ldots,a_N)^\top\) be the resulting vector of normalised weights for the ordered observed sample values

\[\mathbf{a} =\frac{V^{-1}\mathbf{m}} {\sqrt{\left(\mathbf{m}^\top V^{-1}V^{-1}\mathbf{m}\right)}}.\] Royston (Royston 1982; Royston 1995) describes the algorithmic approximation used for these weights and for the \(p\)-value calculation. The Shapiro–Wilk statistic (Shapiro and Wilk 1965) is

\[\begin{equation} W=\frac{\left(\sum_{i=1}^{N} a_i x_{(i)}\right)^2} {\sum_{i=1}^{N} (x_i-\bar{x})^2} (\#eq:shapiro-w) \end{equation}\]

In R, shapiro.test() calls the compiled C_SWilk implementation; the weights are not returned at R level. \(W\) takes values in \((0, 1]\); values close to 1 indicate normality.

Anderson–Darling test ad.test()

Let \(z_i = (x_{(i)} - \bar{x})/s,\; i=1,2,\ldots,N\) be the standardised order statistics of \(x_i\), where \(s\) is the sample standard deviation, and let \(\Phi\) denote the standard normal cumulative distribution function. The test statistic is

\[\begin{equation} A^2 = -N - \frac{1}{N}\sum_{i=1}^{N}(2i-1) \left[\ln\Phi(z_i) + \ln\!\left(1 - \Phi(z_{N+1-i})\right)\right] (\#eq:anderson-a2) \end{equation}\] visstat() uses ad.test() from nortest (Gross and Ligges 2015).

The mean-centred Levene test levene.test()

The package implementation uses Levene’s original mean-centred proposal (Levene 1960).

The Levene test statistic is the one-way ANOVA \(F\) statistic, computed on the absolute residuals \(|e_{ij}|\) in place of the responses \(x_{ij}\); the corresponding Fisher ANOVA formula is given in Eq. @ref(eq:fisher-f):

\[\begin{equation} F_L = \frac{\displaystyle\sum_{i=1}^{k} n_i (\overline{|e|}_i - \overline{|e|})^2\;/\;(k-1)} {\displaystyle\sum_{i=1}^{k}\sum_{j=1}^{n_i}(|e_{ij}| - \overline{|e|}_i)^2\;/\;(N-k)}, (\#eq:levene-f) \end{equation}\]

where \(\overline{|e|}_i\) is the within-group mean of the absolute residuals and \(\overline{|e|}\) is their overall mean.

Bartlett’s test bartlett.test()

Bartlett’s test statistic (Bartlett 1937) is

\[\begin{equation} K^2 = \frac{(N-k)\ln s_p^2 - \displaystyle\sum_{i=1}^k (n_i-1)\ln s_i^2} {1 + \dfrac{1}{3(k-1)}\!\left(\displaystyle\sum_{i=1}^k \frac{1}{n_i-1} - \frac{1}{N-k}\right)}, (\#eq:bartlett-k2) \end{equation}\]

where \(s_p^2\) is the pooled variance from Eq. @ref(eq:pooled-variance).

Under the null hypothesis the statistic approximately follows \(\chi^2(k-1)\).

Breusch–Pagan test bp.test()

For simple linear regression, group-based variance tests are not applicable. The package implementation bp.test() performs the Koenker variant (Koenker 1981) of the Breusch–Pagan test (Breusch and Pagan 1979), which tests whether the \(N\) squared residuals \(e_i^2\) vary systematically with the fitted values from the regression model \(\hat{y}_i\).

The Breusch–Pagan statistic is defined as:

\[\begin{equation} BP = N R^2_\text{aux} (\#eq:breusch-pagan-bp), \end{equation}\]

where \(R^2_\text{aux}\) denotes the coefficient of determination from regressing \(e_i^2\) on \(\hat{y}_i\):

\[R^2_\text{aux} = 1 - \frac{\sum_{i=1}^{N} (e_i^2 - \widehat{e_i^2})^2} {\sum_{i=1}^{N} (e_i^2 - \overline{e^2})^2}.\]

Here \(\widehat{e_i^2}\) are the fitted values from this auxiliary regression and \(\overline{e^2}\) is the mean of the squared residuals.

Under the null hypothesis of homoscedasticity, \(BP\) is compared asymptotically to a \(\chi^2(k-1)\) distribution.

Student’s t-test t.test(..., var.equal = TRUE)

Student’s t-test tests the null hypothesis that the population means of two unpaired groups are equal. The test statistic for Student’s t-test (t.test(..., var.equal = TRUE)) is

\[\begin{equation} t = \frac{\bar{x}_1 - \bar{x}_2} {s_p \sqrt{\dfrac{1}{n_1} + \dfrac{1}{n_2}}}, (\#eq:student-t) \end{equation}\]

where \(s_p\) is the square root of the pooled variance in Eq. @ref(eq:pooled-variance) for \(k=2\). The statistic follows a \(t\)-distribution with \(\nu=n_1+n_2-2\) degrees of freedom.

Fisher’s one-way ANOVA aov()

Fisher’s one-way ANOVA generalises the comparison to more than two groups and tests the null hypothesis that the population means of \(k\) groups are equal. Using the grouped-test notation defined above, the between-group sum of squares is \[ SS_\text{between} = \sum_{i=1}^{k} n_i(\bar{x}_i-\bar{x})^2, \] and the within-group sum of squares is \[ SS_\text{within} = \sum_{i=1}^{k}\sum_{j=1}^{n_i}(x_{ij}-\bar{x}_i)^2. \] Dividing these sums of squares by their degrees of freedom gives the mean squares \[ MS_\text{between}=\frac{SS_\text{between}}{k-1}, \qquad MS_\text{within}=\frac{SS_\text{within}}{N-k}. \] The Fisher ANOVA statistic is \[\begin{equation} F=\frac{MS_\text{between}}{MS_\text{within}}. (\#eq:fisher-f) \end{equation}\]

For \(k=2\), the between-group sum of squares can be written as \[ SS_\text{between} = \frac{n_1n_2}{N}(\bar x_1-\bar x_2)^2 . \] The within-group mean square is the pooled variance, \(MS_\text{within}=s_p^2\). Thus \[ F = \frac{SS_\text{between}}{MS_\text{within}} = \frac{\frac{n_1n_2}{N}(\bar x_1-\bar x_2)^2}{s_p^2}. \] Because \[ \frac{n_1n_2}{N} = \frac{1}{1/n_1+1/n_2}, \] this becomes \[ F = \frac{(\bar x_1-\bar x_2)^2} {s_p^2(1/n_1+1/n_2)} =t^2. \] Therefore, in the two-sample case, Student’s \(t\)-test with var.equal = TRUE and Fisher’s one-way ANOVA return identical \(p\)-values.

Under \(H_0: \mu_1 = \cdots = \mu_k\), the statistic follows \(F(k-1, N-k)\).

Welch’s t-test t.test()

Welch’s t-test (t.test(..., var.equal = FALSE)) compares the means of two independent groups when homogeneous variances cannot be assumed. Its statistic is

\[\begin{equation} t = \frac{\bar{x}_1 - \bar{x}_2} {\sqrt{s_1^2/n_1 + s_2^2/n_2}} (\#eq:welch-t) \end{equation}\]

with degrees of freedom approximated by the Welch–Satterthwaite equation (Welch 1947; Satterthwaite 1946):

\[\begin{equation} \nu \approx \frac{\left(\dfrac{s_1^2}{n_1} + \dfrac{s_2^2}{n_2}\right)^2} {\dfrac{(s_1^2/n_1)^2}{n_1-1} + \dfrac{(s_2^2/n_2)^2}{n_2-1}}. (\#eq:welch-satterthwaite-df) \end{equation}\]

Welch’s heteroscedastic ANOVA oneway.test()

Welch’s heteroscedastic ANOVA (oneway.test()) generalises Welch’s t-test to more than two groups by down-weighting groups with large variance. It compares group means using weights based on sample sizes and variances when homogeneous variances cannot be assumed. Its test statistic is

\[\begin{equation} F_W = \frac{\displaystyle\sum_{i=1}^{k} w_i (\bar{x}_i - \bar{x}_w)^2\;/\;(k-1)} {1 + \dfrac{2(k-2)}{k^2-1} \displaystyle\sum_{i=1}^{k} \dfrac{(1-w_i/w)^2}{n_i-1}}, (\#eq:welch-f) \end{equation}\]

where \(w_i = n_i/s_i^2\) are the inverse-variance weights, \(w = \sum_{i=1}^{k} w_i\), and \(\bar{x}_w = \sum_{i=1}^{k} w_i \bar{x}_i / w\) is the weighted grand mean. The numerator degree of freedom is \(k-1\); the denominator degree of freedom is the Satterthwaite-type approximation returned by oneway.test().

Wilcoxon rank-sum test wilcox.test()

For two independent groups, after the combined ranking step described above, let \(R(x_{1,j})\) be the rank assigned to observation \(j\) in the first group supplied to wilcox.test(x, y), and let \[ W_1=\sum_{j=1}^{n_1}R(x_{1,j}) \] be the rank sum of that first group. The rank sum of the second group is \(W_2=N(N+1)/2-W_1\). If observations are tied, R assigns the average of the tied rank positions. The smallest possible rank sum for group 1 is \(1+\cdots+n_1=n_1(n_1+1)/2\). Subtracting this minimum from the observed rank sum \(W_1\) gives the number of cross-group wins for group 1, with ties contributing one half. The statistic returned by wilcox.test() is the Mann–Whitney statistic (Mann and Whitney 1947):

\[\begin{equation} W = U_1 = W_1 - \frac{n_1(n_1+1)}{2} (\#eq:wilcoxon-w) \end{equation}\]

Equivalently, the same statistic can be written as a count over the \(n_1n_2\) possible cross-group pairs. Let \(C_{>}\) be the number of pairs in which the group-1 observation is larger than the group-2 observation, and let \(C_{=}\) be the number of tied pairs. Then \(W=C_{>}+0.5C_{=}\), and dividing by \(n_1n_2\) gives the empirical Mann–Whitney probability: \[\frac{W}{n_1n_2} = \frac{C_{>}+0.5C_{=}}{n_1n_2}.\] Under the null hypothesis that the two groups have the same continuous distribution, neither group is more likely to produce the larger value. Thus \(W/(n_1n_2)\) is centred at \(1/2\); there are no ties in a continuous distribution, so the tie term is zero. If the two group distributions are the same distribution up to an additive constant, the test can be read as a location test and, because the shift moves all quantiles by the same amount, also as a median test (Fay and Proschan 2010).

Because wilcox.test(x, y) uses the first supplied group for \(W_1\), swapping the two groups uses \(W_2\) and reports \(W_2-n_2(n_2+1)/2=n_1n_2-W\). For each cross-group pair, the two directional contributions always sum to \(1\): group 1 larger gives \(1+0\), group 2 larger gives \(0+1\), and a tie gives \(0.5+0.5\). The two-sided \(p\)-value is unchanged, but the reported statistic and one-sided direction change.

The \(p\)-value is the tail probability of the observed \(W\) under the null distribution of the rank-sum statistic. With R’s default settings, wilcox.test() obtains this null distribution exactly when both groups have fewer than 50 finite observations, and otherwise uses a normal approximation with continuity correction.

Kruskal–Wallis test kruskal.test()

For \(k\) independent groups, the observations from all groups are ranked together as described above. Let \(\bar R_i\) be the mean rank assigned back to group \(i\). If all groups have the same rank distribution, each group has expected mean rank \[ \bar R=\frac{N+1}{2}. \] The Kruskal–Wallis statistic measures how far the group mean ranks \(\bar R_i\) are from this common expected rank (Kruskal and Wallis 1952):

\[\begin{equation} H = \frac{12}{N(N+1)} \sum_{i=1}^{k} n_i \left(\bar{R}_i - \bar{R}\right)^2, (\#eq:kruskal-h) \end{equation}\]

The prefactor \(12/[N(N+1)]\) rescales the weighted squared deviations of the group mean ranks by the sample variance of the \(N\) pooled ranks.

Large values of \(H\) occur when at least one group has systematically higher or lower ranks than expected under equal rank distributions. kruskal.test() evaluates \(H\) against the asymptotic \(\chi^2(k-1)\) null distribution. If ties are present, R first divides \(H\) by the tie factor \[ 1-\frac{\sum_j(t_j^3-t_j)}{N^3-N}, \] where \(t_j\) is the number of observations in tie block \(j\). This factor is the proportion of the original rank variance that remains after tied observations have been assigned average ranks.

If the group distributions are the same distribution up to group-specific additive constants, the test can be read as a location test and, because such shifts move all quantiles by the same amount, also as a median test (Hollander et al. 2014).

For \(k=2\), Kruskal–Wallis and Wilcoxon are based on the same pooled ranks. The two group mean ranks are \(\bar R_1=W_1/n_1\) and \(\bar R_2=W_2/n_2\), and \(W=U_1\) is the reported Wilcoxon statistic from Eq. @ref(eq:wilcoxon-w). In the large-sample approximation, the two-group Kruskal–Wallis statistic \(H\) corresponds to a squared, centred, and rescaled form of the reported Wilcoxon statistic \(W\). Therefore the two tests give identical two-sided \(p\)-values only when Wilcoxon is forced to use the uncorrected large-sample approximation, wilcox.test(..., exact = FALSE, correct = FALSE). In visstat(), wilcox.test() is used with R’s default settings, while kruskal.test() uses the large-sample \(\chi^2\) approximation to \(H\). Therefore the two routes should not be expected to return identical two-group \(p\)-values under the defaults.