Package 'visStatistics'

Breusch-Pagan Test for Heteroscedasticity

Description

Performs the Breusch-Pagan test for heteroscedasticity in linear regression models. Tests the null hypothesis that the error variance is constant (homoscedasticity) against the alternative that the error variance depends on the fitted values.

Usage

bp.test(model)

Arguments

model

A fitted linear model object (from lm()).

Details

Implements the Koenker variant of the Breusch-Pagan test, which regresses the squared raw residuals on the fitted values. The test statistic is:

$BP = n \cdot R^2$

where:

• $n$ = sample size

• $R^2$ = coefficient of determination from auxiliary regression of $e_i^2$ on $\hat{y}_i$

Under the null hypothesis of homoscedasticity, the test statistic follows a chi-squared distribution: $BP \sim \chi^2(k-1)$ where $k$ is the number of parameters in the model (including intercept).

Large values of the test statistic (small p-values) provide evidence against homoscedasticity.

Value

An object of class "htest" with components:

statistic	the value of the chi-squared test statistic.
parameter	degrees of freedom.
p.value	the p-value of the test.
method	a character string indicating the test performed.
data.name	a character string giving the name of the model.

References

Breusch, T. S., and Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287-1294. DOI: 10.2307/1911963

Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107-112. DOI: 10.1016/0304-4076(81)90062-2

Examples

# Example with homoscedastic errors
set.seed(123)
x <- runif(100)
y <- 2 + 3 * x + rnorm(100, sd = 1)
model1 <- lm(y ~ x)
bp.test(model1) # Should not reject (p > 0.05)

# Example with heteroscedastic errors (variance increases with x)
set.seed(456)
x <- runif(100)
y <- 2 + 3 * x + rnorm(100, sd = 0.5 + 2 * x)
model2 <- lm(y ~ x)
bp.test(model2) # Should reject (p < 0.05)

---

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

---

Dunn's Post-Hoc Test for Kruskal-Wallis

Description

Performs pairwise comparisons on the rank sums of a single, combined ranking of all groups, as proposed by Dunn (1964). It is the post-hoc procedure matched to kruskal.test(): both rank the observations globally, so a pairwise decision here concerns the same quantity the omnibus test rejected.

Usage

dunn.test(samples, groups, conf.level = 0.95)

Arguments

samples	numeric vector; the dependent variable.
groups	factor or vector; the grouping variable.
conf.level	numeric; confidence level (default: 0.95). Used only for the significant column; the p-values do not depend on it.

Details

All $k$ samples are combined and ranked from smallest to largest, ties receiving the average rank. Writing $\bar R_i$ for the mean rank of group $i$ and $N$ for the total number of observations, the statistic for groups $i$ and $j$ is

$z_{ij} = (\bar R_i - \bar R_j) / \sigma_{ij},$

with

$\sigma_{ij}^2 = \left[\frac{N(N+1)}{12} - \frac{\sum_{s=1}^{r}(t_s^3 - t_s)}{12(N-1)}\right] \left(\frac{1}{n_i} + \frac{1}{n_j}\right),$

where the $r$ groups of tied scores contain $t_s$ observations each; the subtracted term is zero without ties. This is Eq. (3) of Dunn (1964). The function returns two-sided p-values adjusted by Holm's step-down procedure over all $p = k(k-1)/2$ pairwise comparisons. Note that all pairwise comparisons are performed, so $p$ is not chosen in advance as Dunn's formulation assumes.

Value

A data frame with columns:

group1

First group in comparison

group2

Second group in comparison

mean_rank_diff

Difference in mean ranks (group1 - group2)

se

Standard error of the difference in mean ranks

z

Standard normal test statistic

p_value

Unadjusted two-sided p-value

p_adj

Holm-adjusted p-value for multiple comparisons

significant

Logical; TRUE if p_adj < (1 - conf.level)

References

Dunn, O. J. (1964). Multiple Comparisons Using Rank Sums. Technometrics, 6(3), 241-252. doi:10.1080/00401706.1964.10490181.

Examples

# Convert dose to factor
ToothGrowth$dose <- as.factor(ToothGrowth$dose)

# Perform Dunn's test
result <- dunn.test(ToothGrowth$len, ToothGrowth$dose)
print(result)

---

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; or, in raw-data mode, the first input vector x.
x	First input vector, matching the first argument of visstat(x, y). Required when the effect size cannot be extracted from result alone. In raw-data mode this is the second input vector y.
y	Second input vector, matching the second argument of visstat(x, y). Required when the effect size cannot be extracted from result alone.
...	Passed to visstat() in raw-data mode (e.g. correlation, conf.level).

Details

Notation used below: $x$ and $y$ are the two variables entering the selected analysis, $N$ is the total number of non-missing observations, $n_j$ is the sample size in group $j$, $k$ is the number of groups, $\bar{y}_j$ is the mean of numeric vector $y$ in group $j$, and $s_j^2$ is the variance of numeric vector $y$ in group $j$.

The following estimates are computed internally:

• Student's two-sample t.test(..., var.equal = TRUE): Hedges' $g_{s_p} = J(N-2)(\bar{y}_1-\bar{y}_2)/s_p$, where $s_p = \sqrt{((n_1-1)s_1^2+(n_2-1)s_2^2)/(N-2)}$ and $J(\nu) = \Gamma(\nu/2)/(\sqrt{\nu/2}\Gamma((\nu-1)/2))$.

• Welch's two-sample t.test(..., var.equal = FALSE): Hedges' $g_{s^*} = J(\nu^*)(\bar{y}_1-\bar{y}_2)/s^*$, where $s^* = \sqrt{(s_1^2+s_2^2)/2}$ and $\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/(n_1 n_2)-1$, where $W$ is the statistic returned by wilcox.test() for the first group.

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

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

• Kruskal-Wallis test: eta-squared based on $H$, $\widehat{\eta}_H^2=(H-k+1)/(N-k)$, where $H$ is the Kruskal-Wallis statistic. Negative estimates are truncated to zero.

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

The following estimates are extracted from existing result objects:

• $R^2$ from summary(lm())$r.squared, the coefficient of determination $R^2 = 1 - SS_{res}/SS_{tot}$, where $SS_{res}$ is the residual and $SS_{tot}$ the total sum of squares.

• Spearman's $\rho$ from cor.test(method = "spearman")$estimate, the Pearson correlation of the ranks, $\rho = r(rank(x), rank(y))$.

• Kendall's $\tau_b$ from cor.test(method = "kendall")$estimate, $\tau_b = (n_c - n_d)/\sqrt{(n_0 - n_1)(n_0 - n_2)}$, where $n_c$ and $n_d$ are the numbers of concordant and discordant pairs, $n_0 = n(n-1)/2$ the total number of pairs, and $n_1$ and $n_2$ the numbers of pairs tied in the response and in the predictor.

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

Value

A list with components name, estimate, effect_size_method, and optionally conf.int. In raw-data mode (effect_size(x, y)) the list additionally contains selected_test, the name of the test visstat() selected.

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.

Wendt, H. W. (1972). Dealing with a common problem in social science: A simplified rank-biserial coefficient of correlation based on the U statistic. European Journal of Social Psychology, 2(4), 463–465. doi:10.1002/ejsp.2420020412.

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.

Tomczak, M., & Tomczak, E. (2014). The need to report effect size estimates revisited. An overview of some recommended measures of effect size. Trends in Sport Sciences, 1(21), 19–25. https://tss.awf.poznan.pl/The-need-to-report-effect-size-estimates-revisited-An-overview-of-some-recommended,188960,0,2.html

Cohen, J. (2013). Statistical power analysis for the behavioural 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))

## Raw-data mode: select the test with visstat() and return the effect
## size together with the name of the selected test.
effect_size(ToothGrowth$supp, ToothGrowth$len)

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

---

Games-Howell Post-Hoc Test

Description

Performs pairwise comparisons using the package's own implementation of a Games-Howell-type unequal-variance post-hoc procedure. It does not assume equal variances or equal sample sizes and is used after Welch's ANOVA.

Usage

games.howell(samples, groups, conf.level = 0.95)

Arguments

samples	Numeric vector; the dependent variable.
groups	Factor or vector; the grouping variable.
conf.level	Numeric; confidence level for confidence intervals (default: 0.95).

Details

The implementation uses the Welch-Satterthwaite approximation for degrees of freedom and does not pool variances. P-values are adjusted using the Holm method to control family-wise error rate.

Value

A data frame with columns:

group1

First group in comparison

group2

Second group in comparison

mean_diff

Difference in means (group1 - group2)

se

Standard error of the difference

t

t-statistic

df

Degrees of freedom (Welch-Satterthwaite)

p_value

Unadjusted p-value

p_adj

Holm-adjusted p-value for multiple comparisons

ci_lower

Lower bound of confidence interval

ci_upper

Upper bound of confidence interval

significant

Logical; TRUE if p_adj < (1 - conf.level)

References

Games, P. A., & Howell, J. F. (1976). Pairwise Multiple Comparison Procedures with Unequal N's and/or Variances: A Monte Carlo Study. Journal of Educational Statistics, 1(2). doi:10.2307/1164979.

Examples

# Convert dose to factor
ToothGrowth$dose <- as.factor(ToothGrowth$dose)

# Perform Games-Howell test
result <- games.howell(ToothGrowth$len, ToothGrowth$dose)
print(result)

---

Mean-centred Levene Test for Homogeneity of Variance

Description

Performs Levene's original mean-centred test. It tests the null hypothesis that all groups have equal variances by testing whether the absolute deviations from group means are equal across groups. The function reproduces the behaviour of leveneTest(y, g, center = mean, ...) of the car package.

Usage

levene.test(y, g, data = NULL)

Arguments

y	A numeric response vector.
g	A grouping factor.
data	Optional data frame containing 'y' and 'g'.

Details

For each observation $y_{ij}$ in group $i$, compute the absolute deviation from the group mean:

$z_{ij} = |y_{ij} - \bar{y}_i|$

where $\bar{y}_i$ is the mean of group $i$.

The test statistic is the F-statistic from a one-way ANOVA on the $z_{ij}$ values:

$F = \frac{(N-k) \sum_{i=1}^{k} n_i (\bar{z}_i - \bar{z})^2}{ (k-1) \sum_{i=1}^{k} \sum_{j=1}^{n_i} (z_{ij} - \bar{z}_i)^2}$

where:

• $k$ = number of groups

• $N$ = total sample size

• $n_i$ = sample size of group $i$

• $\bar{z}_i$ = mean of absolute deviations in group $i$

• $\bar{z}$ = overall mean of all absolute deviations

Under the null hypothesis of equal variances, the test statistic follows an F-distribution: $F \sim F(k-1, N-k)$.

Value

An object of class "htest" with components:

statistic	the value of the F-statistic.
parameter	degrees of freedom: df1=k-1, df3=N-k, where k is the number of groups and N the total sample size
p.value	the p-value of the test.
method	a character string indicating the test performed.
data.name	a character string giving the name(s) of the data.

References

Gastwirth, J. L., Gel, Y. R., & Miao, W. (2009). The Impact of Levene's Test of Equality of Variances on Statistical Theory and Practice. Statistical Science, 24(3), 343–360. doi:10.1214/09-STS301.

Examples

set.seed(123)
y <- c(rnorm(10), rnorm(10, sd = 2), rnorm(10, sd = 0.5))
g <- factor(rep(1:3, each = 10))
levene.test(y, g)

# Usage with data frame
df <- data.frame(response = y, group = g)
levene.test(response, group, data = df)

# Example with unequal variances (should reject null hypothesis)
set.seed(456)
y_unequal <- c(rnorm(15, sd = 1), rnorm(15, sd = 5), rnorm(15, sd = 0.2))
g_unequal <- factor(rep(c("A", "B", "C"), each = 15))
levene.test(y_unequal, g_unequal)

---

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

---

Plot method for visstat objects

Description

Replays captured plots or reports saved plot file paths from a visstat object.

Usage

## S3 method for class 'visstat'
plot(x, which = NULL, ...)

Arguments

x	An object of class "visstat", returned by visstat().
which	Integer selecting a single plot to display (1, 2, ...). If NULL (the default), all available plots are listed without being drawn.
...	Currently unused. Included for S3 method compatibility.

Details

When called without which, the method lists all available plots (either as file paths or as indices of captured plots). When called with which, the selected plot is displayed: for file-based output the stored path is printed, for captured plots the plot is replayed via replayPlot().

Value

The object x, invisibly. Called for its side effect.

See Also

print.visstat, summary.visstat, visstat

Examples

# File-based output: plot() lists stored paths
anova_path <- visstat(
  npk$block,
  npk$yield,
  graphicsoutput = "png",
  plotDirectory = tempdir()
)

plot(anova_path)

# Interactive output: plot() lists available plots,
# plot(obj, which = n) replays a specific one
linreg <- visstat(trees$Height, trees$Girth)

plot(linreg)
plot(linreg, which = 2)

---

Print method for visstat objects

Description

Displays a brief summary of the statistical test results and, if available, assumption tests and post hoc comparisons.

Usage

## S3 method for class 'visstat'
print(x, ...)

Arguments

x	An object of class "visstat", returned by visstat().
...	Currently unused. Included for S3 method compatibility.

Details

Quick overview of the statistical analysis results.

Value

The object x, invisibly.

See Also

summary.visstat, plot.visstat, visstat

Examples

anova <- visstat(npk$block, npk$yield)
print(anova)

---

Simulated Q-Q envelopes for linear model residuals

Description

Computes pointwise and simultaneous Q-Q envelopes for internally studentised residuals from an unweighted lm or aov object by repeatedly simulating responses from the fitted normal-error model and refitting the model. The simultaneous envelope follows the Monte Carlo tolerance-band idea of Schützenmeister et al. (2012); see vis_lm_assumptions and visstat for where these bands are drawn.

Usage

qq_lm_envelope(
  model,
  conf.level = 0.95,
  nsim = getOption("visStatistics.qq_nsim", 5000L),
  q.type = 2L,
  tol = 1e-04,
  max.iter = 100L
)

Arguments

model	An unweighted lm or aov object.
conf.level	Numeric confidence level for the envelopes.
nsim	Integer number of simulated refits.
q.type	Integer quantile type passed to stats::quantile().
tol	Numeric tolerance for the simultaneous-band coverage search.
max.iter	Integer maximum number of bisection iterations.

Value

A list of class qq_lm_envelope with the observed sorted residuals, theoretical quantiles, pointwise bounds, simultaneous bounds, the simulated sorted residual matrix, and the achieved simultaneous coverage.

References

Schützenmeister, A., Jensen, U., & Piepho, H.-P. (2012). Checking Normality and Homoscedasticity in the General Linear Model Using Diagnostic Plots. Communications in Statistics - Simulation and Computation, 41(2). doi:10.1080/03610918.2011.582560.

Examples

fit <- lm(mpg ~ wt, data = mtcars)
env <- qq_lm_envelope(fit, nsim = 100)
str(env)

---

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 or if no graphics file was created, 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'

---

Summary method for visstat objects

Description

Displays the full statistical test results and, if available, assumption tests and post hoc comparisons.

Usage

## S3 method for class 'visstat'
summary(object, ...)

Arguments

object	An object of class "visstat", returned by visstat().
...	Currently unused. Included for S3 method compatibility.

Details

This method provides a full textual report of the statistical test results returned by visstat(), and prints the contents of posthoc_summary if present.

Value

The object object, invisibly.

See Also

print.visstat, visstat

Examples

anova <- visstat(npk$block, npk$yield)
summary(anova)

---

Visualisation of linear model assumption diagnostics

Description

Checks the residual diagnostics in the general linear model Student's t-test (t.test,var=EQUAL) Fisher oneway ANOVA (aov) or simple linear regression. Performs the Shapiro-Wilk and Anderson-Darling tests for normality, and for grouped data also levene.test and Bartlett's test for homogeneity of variances. For simple linear regression, heteroscedasticity is assessed with bp.test, the studentised Breusch-Pagan test, which regresses squared raw residuals on fitted values. The normality tests, the grouped variance tests, and the histogram and Q-Q panels are computed from the internally studentised residuals r_i = e_i / (SE_res sqrt(1 - h_i)), which remove the leverage-dependent variance of the raw residuals (Var(e_i) = sigma^2 (1 - h_i)). The residuals-vs-fitted panel (regression mode) uses the z-residuals z_i = e_i / SE_res, which retain the leverage-dependent spread.

Usage

vis_lm_assumptions(
  samples,
  fact,
  cex = 1,
  correlation = FALSE,
  conf.level = 0.95,
  qq_nsim = getOption("visStatistics.qq_nsim", 5000L),
  plot_args = list()
)

Arguments

samples	Numeric vector; the dependent variable.
fact	Factor; the independent variable.
cex	Numeric; scaling factor for plot text and symbols (default: 1).
correlation	Logical. If FALSE and fact is numeric, regression diagnostics are shown. If TRUE, no regression diagnostics are shown. Default is FALSE.
conf.level	Numeric confidence level for the simulated Q-Q envelopes.
qq_nsim	Integer number of simulated refits for the Q-Q envelopes.
plot_args	Optional named list of base graphics parameters.

Details

The Q-Q panel shows the simultaneous and point-wise tolerance bands computed by qq_lm_envelope, which documents their construction and gives the reference for it.

Value

A list with elements:

summary_anova

Summary of the ANOVA model.

shapiro_test

Result from shapiro.test().

ad_test

Result from nortest::ad.test() or a character message if n < 7.

levene_test

Result from levene.test() (grouped diagnostics only).

bartlett_test

Result from bartlett.test() (grouped diagnostics only).

bp_test

Result from bp.test() (regression diagnostics only).

References

Schützenmeister, A., Jensen, U., & Piepho, H.-P. (2012). Checking Normality and Homoscedasticity in the General Linear Model Using Diagnostic Plots. Communications in Statistics - Simulation and Computation, 41(2). doi:10.1080/03610918.2011.582560. (Q-Q simultaneous tolerance band, see qq_lm_envelope.)

See levene.test and bp.test for the references of those two tests.

Examples

ToothGrowth$dose <- as.factor(ToothGrowth$dose)
vis_lm_assumptions(ToothGrowth$len, ToothGrowth$dose, qq_nsim = 100L)

---

Wrapper for visstat_core Allowing Three Different Input Styles

Description

visstat() is a wrapper around visstat_core that provides three input styles and dispatches to one of four analysis routes. Route 1 handles a numeric response with a categorical predictor. By default, Route 1 uses residual-based assumption diagnostics: the Shapiro–Wilk test gates mean-based versus rank-based analysis, and the Levene test then gates equal-variance versus Welch-type mean tests. Above 5000 observations, where shapiro.test() is undefined, the Anderson–Darling test takes over as the residual-normality gate. If the group sizes are unbalanced, the largest group standard deviations occur in the smallest groups, and the selected route is the equal-variance test or a rank-based test, the default route warns and points at group_test = "welch", because those two routes can exceed the nominal significance level in that configuration. Alternatively, group_test = "welch" keeps Route 1 on the mean scale and forces Welch-type tests, while group_test = "rank" forces Wilcoxon/Kruskal–Wallis rank tests. Route 2 handles ordered responses with Wilcoxon/Kruskal–Wallis tests. Route 3 handles two numeric variables with linear regression by default, or Spearman correlation when correlation = TRUE. Route 4 handles two unordered factors with Pearson's $\chi^2$ test or Fisher's exact test.

Usage

visstat(
  x,
  y,
  ...,
  data = NULL,
  conf.level = 0.95,
  correlation = FALSE,
  numbers = TRUE,
  minpercent = 0.05,
  group_test = "gated",
  graphicsoutput = NULL,
  plotName = NULL,
  plotDirectory = getwd(),
  qq_nsim = getOption("visStatistics.qq_nsim", 5000L)
)

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.
...	Named graphical arguments are forwarded to the plots selected by the analysis route for all three input forms. Only in the backward-compatible form, the first unnamed extra argument must specify the predictor or grouping column name in x. Ignored otherwise.
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.
group_test	Character. For a numeric response and factor predictor, "gated" (the default) keeps the assumption gates, "welch" forces Welch-type mean tests, but still displays the assumption-diagnostic plot and warns when residual normality is rejected, whereas "rank" forces Wilcoxon/Kruskal-Wallis rank tests without assessing the assumptions.
graphicsoutput	Saves plot(s) of type "png", "jpeg", "pdf", "svg", "ps" or "tiff" in directory specified in plotDirectory. If NULL, no plots are saved. Any other value is not supported by Cairo(): it triggers a warning, no file is written, and no plot_paths attribute is returned.
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.
qq_nsim	Integer number of simulated refits for the Q-Q envelopes in the assumption diagnostics. Defaults to the option visStatistics.qq_nsim, or 5000 if that is unset.

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 is Route 1: residual normality and residual variance diagnostics select between Student/Welch mean tests and Wilcoxon/Kruskal–Wallis rank tests, unless group_test explicitly forces "welch" or "rank".

If both variables are numeric, a linear model is fitted with y as the response and x as the predictor, unless correlation = TRUE, which requests Spearman rank correlation.

If both variables are factors, unordered factors enter the Pearson-$\chi^2$/Fisher exact association route. Ordered responses with categorical predictors are analysed as rank-based group comparisons; if both variables are ordered and correlation = TRUE, Kendall's $\tau_b$ is used.

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

The Q-Q envelopes in the assumption diagnostics are simulated (see qq_lm_envelope). The number of simulated refits is taken from the qq_nsim argument, which itself defaults to the option visStatistics.qq_nsim and to 5000 if that is unset. Lower it to trade precision for speed, either per call with qq_nsim = 1000L or session wide with options(visStatistics.qq_nsim = 1000L).

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, or using dev.new(width=10, height=6) for a larger plot window.

References

Welch, B. L. (1951). On the Comparison of Several Mean Values: An Alternative Approach. Biometrika, 38(3/4). doi:10.2307/2332579.

Kruskal, W. H., & Wallis, W. A. (1952). Use of Ranks in One-Criterion Variance Analysis. Journal of the American Statistical Association, 47(260). doi:10.2307/2280779.

Bradley, J. V. (1978). Robustness? British Journal of Mathematical and Statistical Psychology, 31(2). doi:10.1111/j.2044-8317.1978.tb00581.x. (Liberal robustness bounds of 2.5%-7.5%, used to judge Type I error control in the size and variance-pairing warnings below.)

Brunner, E., Konietschke, F., Pauly, M., & Puri, M. L. (2017). Rank-Based Procedures in Factorial Designs: Hypotheses About Non-Parametric Treatment Effects. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(5). doi:10.1111/rssb.12222. (Sensitivity of rank-based tests to unequal variances at unequal group sizes, underlying the adverse-pairing warning.)

Delacre, M., Leys, C., Mora, Y. L., & Lakens, D. (2019). Taking Parametric Assumptions Seriously: Arguments for the Use of Welch's F-test instead of the Classical F-test in One-Way ANOVA. International Review of Social Psychology, 32(1). doi:10.5334/irsp.198. (Sample-size thresholds below which Welch's Type I error control is not guaranteed under skewness, used in the small-group warning under group_test = "welch".)

See qq_lm_envelope and vis_lm_assumptions for the Q-Q tolerance-band construction, and levene.test and bp.test for the two variance tests implemented by this package.

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

old_qq_nsim <- getOption("visStatistics.qq_nsim")
options(visStatistics.qq_nsim = 100L)

# 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
# dunn.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, cex = 0.7)

## 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())

options(visStatistics.qq_nsim = old_qq_nsim)

---

Automated Visualization of Statistical Hypothesis Testing

Description

visstat_core() implements the decision tree used by visstat. It receives a data.frame and two column names, determines the corresponding analysis route, creates the diagnostic and result plots, and returns the selected test results as a visstat object.

Usage

visstat_core(
  dataframe,
  varsample,
  varfactor,
  conf.level = 0.95,
  correlation = FALSE,
  numbers = TRUE,
  minpercent = 0.05,
  group_test = "gated",
  graphicsoutput = NULL,
  plotName = NULL,
  plotDirectory = getwd(),
  qq_nsim = getOption("visStatistics.qq_nsim", 5000L),
  plot_args = list()
)

Arguments

dataframe	data.frame with at least two columns.
varsample	character string matching a column name in dataframe. Interpreted as the response if the referenced column is of class numeric or integer and the column named by varfactor is of class factor.
varfactor	character string matching a column name in dataframe. Interpreted as the grouping variable if the referenced column is of class factor and the column named by varsample is of class numeric or integer.
conf.level	Confidence level
correlation	Logical. If FALSE (default), performs simple linear regression analysis with confidence and prediction bands. If TRUE, performs Spearman correlation analysis with trend line only (no regression interpretation).
numbers	a logical indicating whether to show numbers in mosaic count plots.
minpercent	number between 0 and 1 indicating minimal fraction of total count data of a category to be displayed in mosaic count plots.
group_test	Character. For Route 1 only, "gated" (the default) keeps the assumption gates, "welch" forces Welch-type mean tests, but still displays the assumption-diagnostic plot and warns when residual normality is rejected, whereas "rank" forces Wilcoxon/Kruskal-Wallis rank tests without assessing the assumptions.
graphicsoutput	saves plot(s) of type "png", "jpeg", "pdf", "svg", "ps" or "tiff" in directory specified in plotDirectory. If graphicsoutput=NULL, no plots are saved. Any other value is not supported by Cairo(): it triggers a warning and no file is written.
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.
qq_nsim	Integer number of simulated refits for the Q-Q envelopes in the assumption diagnostics. Defaults to the option visStatistics.qq_nsim, or 5000 if that is unset.
plot_args	Optional named list of base graphics parameters.

Details

The decision logic is organised into four routes. Route 1 handles a numeric response with a categorical predictor. By default, Route 1 uses residual-based test selection: Shapiro–Wilk on model residuals gates mean-based versus rank-based analysis, and Levene gates equal-variance versus Welch-type mean tests inside the mean branch. Above 5000 observations, where shapiro.test() is undefined, the Anderson–Darling test takes over as the residual-normality gate. Alternatively, group_test = "welch" forces Welch-type mean tests, and group_test = "rank" forces Wilcoxon/Kruskal–Wallis tests.

Route 2 handles ordered responses with categorical predictors by converting the ordered response to integer level codes and applying Wilcoxon or Kruskal–Wallis tests. Route 3 handles two numeric variables by fitting lm() by default, or Spearman rank correlation when correlation = TRUE. Route 4 handles two unordered factors with Pearson's $\chi^2$ test or Fisher's exact test, depending on expected counts. If both variables are ordered and correlation = TRUE, Kendall's $\tau_b$ is used.

The significance level alpha is defined as 1 - conf.level. Assumption tests are interpreted relative to this threshold.

Under the default group_test = "gated", Route 1 issues a warning when the group sizes are unbalanced, the largest group standard deviations occur in the smallest groups, and the route selected is either the equal-variance test or a rank-based test. In that configuration those two routes can exceed the nominal significance level, whereas a selected Welch test does not and is therefore not flagged; see the Route 1 simulations in vignette("visStatistics").

Implemented main tests:

t.test(), wilcox.test(), aov(), oneway.test(), lm(), kruskal.test(), fisher.test(), chisq.test().

Implemented tests for assumptions:

• Normality: shapiro.test() and ad.test()

• Heteroscedasticity: bartlett.test() and levene.test() and bp.test()

For the general linear model the Shapiro-Wilk, Anderson-Darling, Levene and Bartlett tests are applied to the internally studentised residuals r_i = e_i / (SE_res sqrt(1 - h_i)), which remove the leverage-dependent variance of the raw residuals (Var(e_i) = sigma^2 (1 - h_i)).

Implemented post hoc tests:

• TukeyHSD() for aov()

• games.howell() for oneway.test()

• dunn.test() for kruskal.test()

The Q-Q envelopes in the assumption diagnostics are simulated (see qq_lm_envelope). The number of simulated refits is taken from the qq_nsim argument, which itself defaults to the option visStatistics.qq_nsim and to 5000 if that is unset. Lower it to trade precision for speed, either per call with qq_nsim = 1000L or session wide with options(visStatistics.qq_nsim = 1000L).

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

See Also

The package's vignette vignette("visStatistics") for a description of the decision logic, illustrated with numerous examples. The package is accompanied by its webpage https://shhschilling.github.io/visStatistics/. The main function visstat provides a detailed description of the return value.

Examples

old_qq_nsim <- getOption("visStatistics.qq_nsim")
options(visStatistics.qq_nsim = 100L)

# Welch Two Sample t-test (t.test())
visstat_core(mtcars, "mpg", "am")

## Wilcoxon rank sum test (wilcox.test())
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_core(grades_gender, "Grade", "Sex")

## Welch's oneway ANOVA not assuming equal variances (oneway.test())
anova_npk <- visstat_core(npk, "yield", "block")
anova_npk # prints summary of tests

## Kruskal-Wallis rank sum test (kruskal.test())
visstat_core(iris, "Petal.Width", "Species")
visstat_core(InsectSprays, "count", "spray")

## Simple linear regression  (lm())
visstat_core(trees, "Girth", "Height", conf.level = 0.99)

## Pearson's Chi-squared test (chisq.test())
### Transform array to data.frame
HairEyeColorDataFrame <- counts_to_cases(as.data.frame(HairEyeColor))
visstat_core(HairEyeColorDataFrame, "Hair", "Eye")

## Fisher's exact test (fisher.test())
HairEyeColorMaleFisher <- HairEyeColor[, , 1]
### slicing out a 2 x2 contingency table
blackBrownHazelGreen <- HairEyeColorMaleFisher[1:2, 3:4]
blackBrownHazelGreen <- counts_to_cases(as.data.frame(blackBrownHazelGreen))
fisher_stats <- visstat_core(blackBrownHazelGreen, "Hair", "Eye")

options(visStatistics.qq_nsim = old_qq_nsim)

---