Package 'visStatistics'

Title: Automated Selection and Visualisation of Statistical Hypothesis Tests
Description: Automated test selection, visualised. 'visStatistics' automatically selects and visualises statistical hypothesis tests comparing two vectors, based on their class, distribution, and sample size. Visual outputs, including box plots, bar charts, regression lines with confidence bands, mosaic plots, residual plots, and Q-Q plots, are annotated with relevant test statistics, assumption checks, and post-hoc analyses where applicable. The algorithmic workflow shifts attention from ad-hoc test selection to visual diagnostic assessment and statistical interpretation. It is particularly suited for server-side R applications, where end users interact solely through a web interface to select data groups and receive a complete visual statistical analysis automatically. The same automation makes it useful in time-constrained contexts such as statistical consulting, where it reduces effort spent on test selection and leaves more room for interpretation. The implemented tests cover the most frequently applied inferential methods in biomedical research (Hayat et al. (2017) <doi:10.1371/journal.pone.0179032>). The test selection algorithm proceeds as follows: Input vectors of class numeric or integer are considered numerical; those of class factor are considered categorical; those of class ordered are considered ordinal. Assumptions of residual normality and homogeneity of variances are considered met if the corresponding test yields a p-value greater than the significance level alpha = 1 - conf.level. (1) When the response is numerical and the predictor is categorical, a test comparing central tendencies is selected. If every group contains more than 50 observations, the sampling distribution of the group means is assumed approximately normal by the central limit theorem (Lumley et al. (2002) <doi:10.1146/annurev.publhealth.23.100901.140546>); otherwise, residual normality is assessed using shapiro.test() applied to the standardised residuals of lm(). If normality is not met, wilcox.test() is used when the predictor has two levels and kruskal.test() followed by pairwise.wilcox.test() otherwise. If normality is met, levene.test() assesses variance homogeneity. For two-level predictors, Student's t.test(var.equal = TRUE) is applied if variances are homogeneous and Welch's t.test() otherwise. For predictors with more than two levels, aov() followed by TukeyHSD() is applied if variances are homogeneous, and oneway.test() followed by games.howell() otherwise. (2) When both vectors are numerical, lm() is fitted by default (correlation = FALSE). If correlation = TRUE, Spearman rank correlation is performed. (3) When the response is ordinal, it is converted to numeric ranks and the non-parametric path from (1) is followed (Wilcoxon or Kruskal-Wallis). When both variables are ordinal and correlation = TRUE, Kendall's tau_b is used instead. (4) When both vectors are categorical, Cochran's rule (Cochran (1954) <doi:10.2307/3001666>) is applied to test independence either by chisq.test() or fisher.test().
Authors: Sabine Schilling [cre, aut, cph] (ORCID: <https://orcid.org/0000-0002-8318-9421>, year: 2026), Peter Kauf [ctb]
Maintainer: Sabine Schilling <[email protected]>
License: MIT + file LICENSE
Version: 0.2.1
Built: 2026-06-02 15:06:39 UTC
Source: https://github.com/shhschilling/visstatistics

Help Index

Convert data frame of counts to dataframe of cases. The data frame must contain a column with frequencies (counts) as generated by as.data.frame from a contingency table

Description

Convert data frame of counts to dataframe of cases. The data frame must contain a column with frequencies (counts) as generated by as.data.frame from a contingency table

Usage

counts_to_cases(x, countcol = "Freq")

Arguments

x

a data.frame of counts generated from a contingency table.
countcol

character string, name of the column of x containing the counts. Default name of the column is 'Freq'.

Value

data frame of cases of dimension (total number of counts as sum of 'Freq' in x) times 2.

Examples

counts_to_cases(as.data.frame(HairEyeColor[, , 1]), countcol = "Freq")

Compute an effect-size estimate for a visstat result

Description

effect_size() returns the effect-size estimate associated with a visstat() result. If result$effect_size is already present, it is returned unchanged. Otherwise, the estimate is computed from the test object stored in result; for some base R stats results, it is extracted directly from the returned object.

Usage

effect_size(result, x = NULL, y = NULL)

Arguments

result

A list returned by visstat() or a compatible test result object.
x

First input vector, matching the first argument of visstat(x, y). Required when the effect size cannot be extracted from result alone.
y

Second input vector, matching the second argument of visstat(x, y). Required when the effect size cannot be extracted from result alone.

Details

Notation used below: xx and yy are the two variables entering the selected analysis, NN is the total number of non-missing observations, njn_j is the sample size in group jj, kk is the number of groups, yˉj\bar{y}_j is the mean of numeric vector yy in group jj, and sj2s_j^2 is the variance of numeric vector yy in group jj.

The following estimates are computed internally:

  • Student's two-sample t.test(..., var.equal = TRUE): Hedges' gsp=J(N−2)(yˉ1−yˉ2)/spg_{s_p} = J(N-2)(\bar{y}_1-\bar{y}_2)/s_p, where sp=((n1−1)s12+(n2−1)s22)/(N−2)s_p = \sqrt{((n_1-1)s_1^2+(n_2-1)s_2^2)/(N-2)} and J(ν)=Γ(ν/2)/(ν/2Γ((ν−1)/2))J(\nu) = \Gamma(\nu/2)/(\sqrt{\nu/2}\Gamma((\nu-1)/2)).

  • Welch's two-sample t.test(..., var.equal = FALSE): Hedges' gs∗=J(ν∗)(yˉ1−yˉ2)/s∗g_{s^*} = J(\nu^*)(\bar{y}_1-\bar{y}_2)/s^*, where s∗=(s12+s22)/2s^* = \sqrt{(s_1^2+s_2^2)/2} and ν∗=((n1−1)(n2−1)(s12+s22)2)/((n2−1)s14+(n1−1)s24)\nu^* = ((n_1-1)(n_2-1)(s_1^2+s_2^2)^2)/ ((n_2-1)s_1^4+(n_1-1)s_2^4).

  • Wilcoxon rank-sum test: signed rank-biserial correlation r=2W/(n1n2)−1r = 2W/(n_1 n_2)-1, where WW is the statistic returned by wilcox.test() for the first group.

  • Fisher's one-way ANOVA: omega-squared ω2=ν1(F−1)/(ν1F+ν2+1)\omega^2 = \nu_1(F-1)/(\nu_1F+\nu_2+1), where FF is the ordinary one-way ANOVA statistic, ν1=k−1\nu_1=k-1, and ν2=N−k\nu_2=N-k. Negative estimates are truncated to zero.

  • Welch's one-way test: approximate omega-squared-type estimate ν1(FW−1)/(ν1FW+ν2+1)\nu_1(F_W-1)/(\nu_1F_W+\nu_2+1), where FWF_W is the Welch ANOVA statistic, ν1=k−1\nu_1=k-1, and ν2\nu_2 is the usually fractional denominator degree of freedom returned by oneway.test(). Negative estimates are truncated to zero.

  • Kruskal-Wallis test: Kelley-adjusted eta-squared based on HH, ηH2=(H−k+1)/(N−k)\eta_H^2=(H-k+1)/(N-k), where HH is the Kruskal-Wallis statistic. Negative estimates are truncated to zero.

  • Pearson's chi-squared test: Cramer's VV for general R×CR\times C tables, V=χ2/(N⋅(min⁡(R,C)−1))V=\sqrt{\chi^2/(N\cdot(\min(R,C)-1))}, where RR and CC are the numbers of rows and columns. For 2×22\times 2 tables this is phi, χ2/N\sqrt{\chi^2/N}. The chi-squared statistic is used as supplied by chisq.test().

The following estimates are extracted from existing result objects:

  • R2R^2 from summary(lm())$r.squared.

  • Spearman's ρ\rho from cor.test(method = "spearman")$estimate.

  • Kendall's τb\tau_b from cor.test(method = "kendall")$estimate.

  • The conditional maximum-likelihood odds ratio from fisher.test()$estimate and its confidence interval from fisher.test()$conf.int for 2×22\times 2 tables.

Value

A list with components name, estimate, effect_size_method, and optionally conf.int.

References

Hedges, L. V. (1981). Distribution theory for Glass's estimator of effect size and related estimators. Journal of Educational Statistics, 6(2), 107–128. doi:10.3102/10769986006002107.

Delacre, M., Lakens, D., Ley, C., Liu, L., & Leys, C. (2021). Why Hedges' g*s based on the non-pooled standard deviation should be reported with Welch's t-test. PsyArXiv. doi:10.31234/osf.io/tu6mp.

Glass, G. V. (1965). A ranking variable analogue of biserial correlation: Implications for short-cut item analysis. Journal of Educational Measurement, 2(1), 91–95. doi:10.1111/j.1745-3984.1965.tb00396.x.

Albers, C., & Lakens, D. (2018). When power analyses based on pilot data are biased: Inaccurate effect size estimators and follow-up bias. Journal of Experimental Social Psychology, 74, 187–195. doi:10.1016/j.jesp.2017.09.004.

Kelley, T. L. (1935). An unbiased correlation ratio measure. Proceedings of the National Academy of Sciences, 21(9), 554–559. doi:10.1073/pnas.21.9.554.

Cohen, J. (2013). Statistical power analysis for the behavioral sciences. Routledge. doi:10.4324/9780203771587.

Examples

x <- ToothGrowth$supp
y <- ToothGrowth$len
tt <- list("t-test-statistics" = t.test(y ~ x, var.equal = TRUE))
effect_size(tt, x = x, y = y)

kw <- list(
  "Kruskal Wallis rank sum test" = kruskal.test(Petal.Width ~ Species,
                                               data = iris)
)
effect_size(kw, x = iris$Species, y = iris$Petal.Width)

tab <- matrix(c(10, 5, 4, 12), nrow = 2)
effect_size(chisq.test(tab))

## Not run: 
## Large-sample example with a statistically significant Student's
## t-test p-value but a small effect size, measured by Hedges' g
## using the pooled standard deviation. A small mean shift is added
## to noisy normal data. Because N is large, the t-test p-value
## becomes small, while Hedges' g remains close to zero.
## The residual Shapiro-Wilk p-value in the diagnostic panel is NA
## because shapiro.test() is limited to n <= 5000.
set.seed(20260525)
n <- 2501
mean_shift <- 0.1
group <- factor(rep(c("control", "treatment"), each = n))
response <- rnorm(2 * n) + rep(c(0, mean_shift), each = n)
res <- visstat(group, response)
res[["t-test-statistics"]]$method
res[["t-test-statistics"]]$p.value
res$effect_size$effect_size_method
res$effect_size$estimate

## End(Not run)

Cairo wrapper function with plot capture capability

Description

Cairo wrapper function returning NULL if not type is specified. Enhanced version that can capture plots for later replay.

Usage

openGraphCairo(
  width = 640,
  height = 480,
  fileName = NULL,
  type = NULL,
  fileDirectory = getwd(),
  pointsize = 12,
  bg = "transparent",
  canvas = "white",
  units = "px",
  dpi = 150
)

Arguments

width

see Cairo()
height

see Cairo()
fileName

name of file to be created. Does not include both file extension '.type' and file filedirectory. Default file name 'visstat_plot'.
type

Supported output types are 'png', 'jpeg', 'pdf', 'svg', 'ps' and 'tiff'. See Cairo()
fileDirectory

path of directory, where plot is stored. Default current working directory.
pointsize

see Cairo()
bg

see Cairo()
canvas

see Cairo()
units

see Cairo()
dpi

DPI used for the conversion of units to pixels. Default value 150.

Details

openGraphCairo() Cairo() wrapper function. Differences to Cairo: a) prematurely ends the function call to Cairo() returning NULL, if no output type of types 'png', 'jpeg', 'pdf', 'svg', 'ps' or 'tiff' is provided. b) The file argument of the underlying Cairo function is generated by file.path(fileDirectory,paste(fileName,'.', type, sep = '')). c) Can set up plot capture when capture_env is provided.

Value

NULL, if no type is specified. Otherwise see Cairo()

Examples

##  adapted from example in \code{Cairo()}
openGraphCairo(fileName = "normal_dist", type = "pdf", fileDirectory = tempdir())
plot(rnorm(4000), rnorm(4000), col = "#ff000018", pch = 19, cex = 2)
dev.off() # creates a file 'normal_dist.pdf' in the directory specified in fileDirectory
# ## remove the plot from fileDirectory
file.remove(file.path(tempdir(), "normal_dist.pdf"))

Saves Graphical Output with plot capture capability

Description

Closes all graphical devices with dev.off() and saves the output only if both fileName and type are provided. Enhanced version that can capture plots before closing devices.

Usage

saveGraphVisstat(
  fileName = NULL,
  type = NULL,
  fileDirectory = getwd(),
  oldfile = NULL,
  capture_env = NULL
)

Arguments

fileName

name of file to be created in directory fileDirectory without file extension '.type'.
type

see Cairo().
fileDirectory

path of directory, where graphic is stored. Default setting current working directory.
oldfile

old file of same name to be overwritten
capture_env

Environment to store captured plots. If NULL, no capture occurs.

Value

NULL, if no type or fileName is provided, file path if graph is created

Examples

# very simple KDE (adapted from example in Cairo())
openGraphCairo(type = "png", fileDirectory = tempdir())
plot(rnorm(4000), rnorm(4000), col = "#ff000018", pch = 19, cex = 2)
# save file 'norm.png' in directory specified in fileDirectory
saveGraphVisstat("norm", type = "png", fileDirectory = tempdir())
file.remove(file.path(tempdir(), "norm.png")) # remove file 'norm.png'

Wrapper for visstat_core Allowing Three Different Input Styles

Description

visstat() is a wrapper around visstat_core that provides three alternative input styles: a formula interface, a standardised vector interface, and a backward-compatible data frame interface. visstat_core defines the decision logic for statistical hypothesis testing and visualisation between two variables of class "numeric", "integer", or "factor".

Usage

visstat(
  x,
  y,
  ...,
  data = NULL,
  conf.level = 0.95,
  correlation = FALSE,
  numbers = TRUE,
  minpercent = 0.05,
  graphicsoutput = NULL,
  plotName = NULL,
  plotDirectory = getwd()
)

Arguments

x

For the formula interface: a formula of the form y ~ x, where y is the response variable and x is the predictor or grouping variable (requires data argument). For the standardised form: a vector of class "numeric", "integer", or "factor" representing the predictor or grouping variable. For the backward-compatible form: a data.frame containing the relevant columns.
y

For the formula interface: not used (variables are extracted from the formula). For the standardised form: a vector of class "numeric", "integer", or "factor" representing the response variable. For the backward-compatible form: a character string specifying the name of the response variable column in x.
...

For the backward-compatible form only: a character string specifying the name of the predictor or grouping variable column in x. Ignored for formula and standardised input styles.
data

A data.frame containing the variables specified in the formula. Required when using the formula interface. Ignored for other input styles.
conf.level

Confidence level for statistical inference; default is 0.95.
correlation

Logical. If FALSE (default), performs simple linear regression analysis with confidence and prediction bands when both variables are numeric. If TRUE, performs Spearman correlation analysis with trend line only (no regression interpretation).
numbers

Logical. Whether to annotate plots with numeric values.
minpercent

Number between 0 and 1 indicating minimal fraction of total count data of a category to be displayed in mosaic count plots.
graphicsoutput

Saves plot(s) of type "png", "jpg", "tiff" or "bmp" in directory specified in plotDirectory. If NULL, no plots are saved.
plotName

Graphical output is stored following the naming convention "plotName.graphicsoutput" in plotDirectory. Without specifying this parameter, plotName is automatically generated following the convention "statisticalTestName_varsample_varfactor".
plotDirectory

Specifies directory where generated plots are stored. Default is current working directory.

Details

This wrapper supports three input formats:

(1) Formula interface: visstat(y ~ x, data = df), where the formula specifies the response (y) and predictor (x) variables, and data is a data frame containing these variables.

(2) Standardised form: visstat(x, y), where both x and y are vectors of class "numeric", "integer", or "factor". Here x is the predictor or grouping variable and y is the response variable.

(3) Backward-compatible form: visstat(dataframe, "name_of_y", "name_of_x"), where the first character string refers to the response variable and the second to the predictor or grouping variable. Both must be column names in dataframe. This form gives a warning and may be removed in a future version.

The interpretation of x and y depends on the variable classes. Throughout, data of class numeric or integer are referred to as numeric, while data of class factor are referred to as categorical:

If one variable is numeric and the other a factor, the numeric vector is the response (y) and the factor is the grouping variable (x). This supports tests of central tendencies (e.g., t-test, Welch's ANOVA, Wilcoxon, Kruskal-Wallis).

If both variables are numeric, a linear model is fitted with y as the response and x as the predictor.

If both variables are factors, an association test (Chi-squared or Fisher's exact) is used. The test result is invariant to variable order, but visualisations (e.g., axis layout, bar orientation) depend on the roles of x and y.

This wrapper standardises the input and calls visstat_core, which selects and executes the appropriate test with visual output and assumption diagnostics.

Value

An object of class "visstat" containing the results of the automatically selected statistical test. The specific contents depend on which test was performed. Additionally, the returned object includes two attributes:

  • plot_paths: Character vector of file paths where plots were saved (if graphicsoutput was specified)

  • captured_plots: List of captured plot objects for programmatic access

In case of insufficient data, returns a list with an error element and basic input summary information.

Note

For best visualization, ensure that the RStudio Plots pane is adequately sized. If you get "figure margins too large" errors, try expanding the Plots pane in RStudio, using dev.new(width=10, height=6) for a larger plot window, or reducing the cex parameter.

See Also

visstat_core defining the decision logic, the package's vignette vignette("visStatistics") explaining the decision logic accompanied by illustrative examples, and the accompanying webpage https://shhschilling.github.io/visStatistics/.

Examples

# Formula interface
mtcars$am <- as.factor(mtcars$am)
visstat(mpg ~ am, data = mtcars)

# Standardised usage
visstat(mtcars$am, mtcars$mpg)

## Student's t-test (equal variances, two groups)
# When residuals are normally distributed and Levene's test indicates
# homoscedasticity, the classic Student's t-test with pooled variance is used
df <- droplevels(subset(PlantGrowth, group %in% c("ctrl", "trt1")))
visstat(df$group,df$weight)

## Welch's t-test (unequal variances, two groups)
# When residuals are normally distributed but Levene's test indicates
# heteroscedasticity, Welch's t-test is used
visstat(mtcars$am, mtcars$mpg)

## Wilcoxon rank sum test (non-normal, two groups)
# When residuals are not normally distributed
grades_gender <- data.frame(
  Sex = as.factor(c(rep("Girl", 20), rep("Boy", 20))),
  Grade = c(
    19.3, 18.1, 15.2, 18.3, 7.9, 6.2, 19.4, 20.3, 9.3, 11.3,
    18.2, 17.5, 10.2, 20.1, 13.3, 17.2, 15.1, 16.2, 17.3, 16.5,
    5.1, 15.3, 17.1, 14.8, 15.4, 14.4, 7.5, 15.5, 6.0, 17.4,
    7.3, 14.3, 13.5, 8.0, 19.5, 13.4, 17.9, 17.7, 16.4, 15.6
  )
)
visstat(grades_gender$Sex, grades_gender$Grade)

## Fisher's ANOVA (equal variances, >2 groups)
# When residuals are normally distributed and Levene's test indicates
# homoscedasticity, classic Fisher's ANOVA with TukeyHSD post-hoc is used.
# Different green letters indicate significant differences between groups. 
visstat(PlantGrowth$group, PlantGrowth$weight)

## Welch's one-way ANOVA (unequal variances, >2 groups)
set.seed(123)
values <- c(rnorm(20, 10, 1),rnorm(20, 15, 5),rnorm(20, 12, 2))
groups <- factor(rep(c("A", "B", "C"), each = 20))
visstat(groups, values)
## Kruskal-Wallis (non-normal, >2 groups)
# When residuals are not normally distributed, kruskal.test() is followed by 
# pairwise.wilcox.test. 
visstat(iris$Species, iris$Petal.Width)

## Simple linear regression (both numeric)
visstat(trees$Height, trees$Girth, conf.level = 0.99)

## Pearson's Chi-squared test (both factors, large expected counts)
HairEyeColorDataFrame <- counts_to_cases(as.data.frame(HairEyeColor))
visstat(HairEyeColorDataFrame$Eye, HairEyeColorDataFrame$Hair)

## Fisher's exact test (both factors, small expected counts)
HairEyeColorMaleFisher <- HairEyeColor[, , 1]
blackBrownHazelGreen <- HairEyeColorMaleFisher[1:2, 3:4]
blackBrownHazelGreen <- counts_to_cases(as.data.frame(blackBrownHazelGreen))
visstat(blackBrownHazelGreen$Eye, blackBrownHazelGreen$Hair)

## Save PNG
visstat(blackBrownHazelGreen$Hair, blackBrownHazelGreen$Eye,
        graphicsoutput = "png", plotDirectory = tempdir())

## Custom plot name
visstat(iris$Species, iris$Petal.Width,
        graphicsoutput = "pdf", plotName = "kruskal_iris", plotDirectory = tempdir())