4.6.1 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.

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

6.1.1.1 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)

Figure 6.1: 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.

6.1.1.2 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 4.6).

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).

Figure 6.2(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 6.2(b) replays the result panel (which = 2).

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

Figure 6.2: PlantGrowthdata set: 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 stored 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"

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

6.1.2.1 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)

Figure 6.3: 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.

6.1.2.2 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)

Figure 6.4: 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\)).

6.1.3 Wilcoxon rank-sum test and Kruskal–Wallis test

6.1.3.1 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)

Figure 6.5: 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.

6.1.3.2 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 6.4), 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)

Figure 6.6: 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\)).

6.2.1 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.

6.2.1.2 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 6.7, 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.

Figure 6.7: 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\)).

6.3.1 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 6.8, 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 \tag{6.1} \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. (9.1) 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. (F.1)) as effect size.

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

Figure 6.8: 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.

Figure 6.9: 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 6.5. The following example shows a Gamma generalised linear model outside visstat().

6.3.2 Model exploration outside visstat()

As a model outside 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

Figure 6.10: 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.

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

For a contingency table with \(R\) response levels and \(C\) predictor levels, Pearson’s \(\chi^2\) test (Eq. (E.2)) 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. (E.1)) 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)

Figure 6.11: 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.

6.4.2 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 6.12, 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.

6.4.3 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 6.12, 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
)

Figure 6.12: 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.

6.5.1 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 6.13, 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)

Figure 6.13: 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.

6.5.2 Spearman rank correlation with correlation = TRUE

For the ozone example introduced in Section 6.3.1, staying within visstat() with the flag correlation = TRUE gives the Spearman analysis shown in Figure 6.13, right.

9.1.1 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} \tag{9.1} 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. (9.1), 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.

9.1.2 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. (A.1)) and the Anderson–Darling test (Anderson and Darling 1952) (ad.test(); Eq. (A.2)). 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. Among the normality tests compared by Razali and Wah (2011) the Shapiro–Wilk test was the most powerful against both symmetric and asymmetric alternatives, although all of them had low power below about 30 observations. Therefore, the Shapiro–Wilk test is used as the normality gate in the automated test selection (Section 9.2.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. (A.3)) and Bartlett’s test (Bartlett 1937) (bartlett.test(); Eq. (A.4)).

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. (A.5)) on raw residuals (Schützenmeister et al. 2012).

9.1.3 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. (9.5)) by group. The last panel shows whether residual spread is comparable across factor levels, the pattern assessed formally by the Levene (Eq. (A.3)) and Bartlett (Eq. (A.4)) 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 9.2.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 6.4, left, and 6.8, left.

9.2.1 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 9.1 expands the default routing logic for tests of central tendencies.

Figure 9.1: 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. (9.1) is fitted between the numeric response and the categorical predictor, and the model residuals of Eq. (9.2) 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) test is an omnibus normality test, sensitive to both skewness and kurtosis (D’Agostino and Pearson 1973).

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

If residual normality is not rejected, the mean-centred Levene test (L) (Levene 1960) (Eq. (A.3)) gates the variance assessment:

For homoscedastic data (\(p_\text{L} > \alpha\)), t.test(var.equal = TRUE) (Eq. (B.1)) is applied for two groups, or Fisher’s aov() (Eq. (B.2)) for more than two groups. For heteroscedastic data (\(p_\text{L} \le \alpha\)), Welch’s t.test() (Eq. (B.4)) is applied for two groups, or Welch’s oneway.test() (Eq. (B.6)) for more than two groups.

Independent of assumption testing, the user can enforce group mean comparisons by the option group_test = welch which defaults to Welch variants of the t-test (t.test()) and ANOVA (oneway.test()), otherwise the option group_test = rankswitches to the non-parametric alternatives wilcox.test() and kruskal.test().

The two overrides differ in what they display: group_test = welch still fits the linear model and shows the assumption-diagnostic panel, with a warning when the Shapiro–Wilk test rejects residual normality, whereas group_test = rank enters the rank branch directly and does not generate the assumption plot with its corresponding test statistics.

The rationale for the automated gating, and the mean- and rank-based alternatives and the limitations of each approach are discussed in Section 10 and Section 7.

9.2.2 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.

9.2.3 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. (6.1)), and the diagnostic panel described in Section 9.1.3 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 6.3.1.

9.2.4 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.

9.2.5 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. (D.1)) with method = "kendall", exact = FALSE for two ordered variables, and Spearman’s \(\rho\) (Eq. (D.2)) 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.

10.2.1 Bradley’s boundaries in heteroscedastic simulations

Albeit Fisher’s Anova (F) assumes variance homogeneity, it remains within Bradley’s boundaries for all simulated balanced, unequal standard deviation four group comparison (panel B of Fig. 10.2). In unbalanced simulations where smaller group sizes were combined with smaller standard deviations (panel C), F is conservative throughout and sits at the lower bound of 2.5%, which it undercuts under the symmetric input with excess kurtosis 6 at \(\bar n = 10\) by less than the Monte Carlo error.

This is not the case any more in the reversed scenario (panel D): Combining larger groups with smaller standard deviations, F rejects the null hypotheses of equal means between 11.7%–13.2%. In this reversed scenario, also W leaves the Bradley bounds only for skewness 2 with excess kurtosis 6.6 for a mean group size \(\bar n \le 20\) (panel D).

In the first two symmetric columns of Fig. 10.2, both SW and SW+L are still type I tests and remain in Bradley’s boundaries in the balanced simulations (panel B) as well as the unbalanced simulations of panel C. But combining smaller sample sizes with larger standard deviations (panel D), SW leaves the Bradley bound in three out of 8 scenarios with \(\bar n \ge 50\), as in larger samples the Shapiro–Wilk gate is overpowered and sends most replications to KW, whose own rate is affected by unequal variances at unequal group sizes (Zimmerman 2004). Adding the Levene gate (SW+L) increases the type I error in smaller samples even further, as the Levene gate at \(\bar n=10\) is underpowered and routes nearly half the replications to F. SW+L leaves the bounds in seven of eight cells of panel D, at worst at about twice the nominal 5%.

A.1.1 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} \tag{A.1} \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.

A.1.2 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] \tag{A.2} \end{equation}\] visstat() uses ad.test() from nortest (Gross and Ligges 2015).

A.2.1 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. (B.2):

\[\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)}, \tag{A.3} \end{equation}\]

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

A.2.2 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)}, \tag{A.4} \end{equation}\]

where \(s_p^2\) is the pooled variance from Eq. (10.1).

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

A.2.3 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} \tag{A.5}, \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.

B.1.1 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}}}, \tag{B.1} \end{equation}\]

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

B.1.2 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}}. \tag{B.2} \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)\).

B.2.1 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}} \tag{B.4} \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}}. \tag{B.5} \end{equation}\]

B.2.2 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}}, \tag{B.6} \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().

C.1.1 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} \tag{C.1} \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.

C.1.2 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, \tag{C.2} \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. (C.1). 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.