Package 'STATassist'

Description

⁠$cells⁠ is the canonical form, one row per cell, because that is the shape that survives being written out as JSON with its labels attached. A table is the shape to read it in, so it is built on request rather than stored twice and left to drift.

Usage

## S3 method for class 'sa_categorical'
as.table(x, ...)

Arguments

Details

The rows dropped for a missing value or for a level outside category_lv are already gone, so this is the table the tests were run on rather than a fresh table() of the input.

Value

A two-dimensional table() of counts, its dimnames named after the two axes of the design.

Examples

res <- compare_categorical_groups(
  data.frame(smoker = rep(c("y", "n"), each = 30),
             grade  = c(rep(c("high", "low"), c(22, 8)),
                        rep(c("high", "low"), c(9, 21))))
)
as.table(res)

## Which is the table the tests were run on, so it feeds `stats` directly.
stats::chisq.test(as.table(res))

Description

Divides, or on the log2 scale subtracts, the centre of the control group out of each feature, so that the control lands on 1 (raw) or 0 (log2) and every other observation reads as its distance from the control rather than as a measurement of its own. The data that comes back is the data that went in with the feats columns replaced, which is what lets the same call be handed straight to a comparison or to draw_heatmap().

Usage

center_by_control(
  data,
  feats,
  group,
  group_lv,
  control_label = group_lv[1],
  fc_mean = c("arith", "geom"),
  input_scale = c("raw", "log2")
)

Arguments

Details

The arguments are the ones compare_two_groups() and compare_multiple_groups() already take, in the same order and with the same defaults, because the point of the function is that one set of arguments describes both steps:

centred <- center_by_control(data, feats, group, group_lv,
                             input_scale = "log2")
compare_two_groups(centred, feats, group, group_lv, input_scale = "log2")

fc_mean is the comparison's own argument rather than a centring method of this function's, and it is resolved the same way: the geometric mean by default on the log2 scale, where it is the convention, and the arithmetic mean otherwise. Passing the same value to both is what keeps the baseline removed here identical to the centre the comparison divides in its effect table.

Two or more levels are accepted, so the same call serves a two-group and a multi-group design. Only the control level takes part in the baseline; the rest are centred on it.

Value

The data that was passed in, as a data.frame, with the feats columns replaced by their control-relative values. Row count, row order, row names and every column that is not a feature are left as they arrived. A feature whose baseline could not be taken comes back as an NA column and is named in a warning.

What the comparison reports afterwards

The transformation is one constant per feature applied to every row of it, so most of what a comparison reports does not move:

fold_change and log2fc

Unchanged. Both group centres are divided by the same baseline, so their ratio survives it. This holds for all four combinations of fc_mean and input_scale.

p-values

Unchanged, in every test family. On the log2 scale the baseline is subtracted, and the tests are shift invariant; on the raw scale it is divided, and dividing by a positive constant leaves the t statistic and every rank untouched.

The reference centre

Becomes 1. It is y_center in compare_two_groups() and ref_center in compare_multiple_groups(), and once it is 1 the other centre in the row is the fold change itself.

mean_diff, hl_shift and trim_diff

Unchanged on the log2 scale. On the raw scale they are divided by the baseline along with the data, since they are differences rather than ratios.

The reference centre lands on 1 because the observations used here and the ones the comparison centres are the same. Under paired = TRUE a comparison keeps complete pairs or complete subjects only, so if any are dropped its control centre is taken on fewer rows than the baseline was and comes back near 1 rather than at it. The ratios are unaffected either way.

See Also

compare_two_groups() and compare_multiple_groups(), whose arguments this shares, and draw_heatmap(), which takes the result with scale = "none" to show departure from the control rather than from each feature's own mean.

Examples

feats <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
group_lv <- levels(iris$Species)

## setosa is the first level and so the control. Every measurement is now a
## multiple of the setosa mean, and setosa itself centres on 1.
centred <- center_by_control(iris, feats, iris$Species, group_lv)
colMeans(centred[iris$Species == "setosa", feats])

## The comparison takes the centred data with its arguments unchanged, and
## reports the fold changes it reported before the centring.
before <- compare_multiple_groups(iris, feats, iris$Species, group_lv,
                                  posthoc = FALSE, diagnose = FALSE)
after <- compare_multiple_groups(centred, feats, iris$Species, group_lv,
                                 posthoc = FALSE, diagnose = FALSE)
cbind(before = before$effect$log2fc, after = after$effect$log2fc)

## What did change is the reference centre, which is now 1 by construction,
## so `extreme_center` reads as the fold change on its own.
after$effect[c("ref_center", "extreme_center", "fold_change")]

## On the log2 scale the baseline is subtracted instead, and the control
## centres on 0. Rounded, because the subtraction leaves floating point dust
## rather than a clean zero.
d <- iris[iris$Species != "setosa", ]
lg <- center_by_control(log2(d[feats]), feats, d$Species,
                        c("versicolor", "virginica"),
                        input_scale = "log2")
round(colMeans(lg[d$Species == "versicolor", ]), 12)

## Pointing the baseline at the other level is one argument rather than a
## rewritten `group_lv`.
virg <- center_by_control(iris, feats, iris$Species, group_lv,
                          control_label = "virginica")
colMeans(virg[iris$Species == "virginica", feats])

Description

Grows a cluster out of every point that has at least min_pts neighbours within eps of it, joining the neighbourhoods that overlap, and leaves the points that never fell into one as noise. How many clusters there are is the answer rather than the question, and it can be none.

Usage

cluster_dbscan(
  data,
  feats = NULL,
  cluster_scale = c("samples", "features"),
  center = TRUE,
  scale = TRUE,
  eps = NULL,
  min_pts = NULL
)

Arguments

Details

The input is the wide format the comparison functions take: one row per sample and one column per feature. Which margin of it becomes the thing being clustered is cluster_scale, and design$point_type reports the answer.

Value

An object of class sa_cluster, a plain list. The slots are cluster_hclust()'s, with analysis "dbscan", fit the dbscan::dbscan() object, and design$n_clusters and design$n_noise both answers rather than arguments.

Noise is a result, not a failure

A point that never joined a neighbourhood gets cluster 0, which is dbscan::dbscan()'s convention and this contract's. It has no silhouette, since noise is not a cluster to be near or far from, and clusters has no row for it; design$n_noise is the count. All points being noise is a possible answer and means the density asked for is not present, which is usually eps being too small or min_pts too large for the data.

Choosing eps

eps is a radius in the units of the matrix being clustered, which is why there is no default constant that could be right. Left as NULL it is derived from the distances actually present: every point's distance to its min_pts - 1th neighbour is collected, and eps is the 95th percentile of those. The rule therefore reads as "assume about one point in twenty is noise, and set the radius that leaves that many outside", which is a statement worth disagreeing with rather than a number off a curve.

The quantity being read is the one dbscan::kNNdistplot() plots and the manual procedure reads a knee off by eye. The knee is not what is taken, because locating it arithmetically fails on a curve that rises gradually rather than turning; the probe in ⁠Test/cursor_test/2026_08_15/⁠ has the case. Either way it is a heuristic: parameters$eps_source records whether the value was supplied or derived, and a derived one is said out loud when it is used.

Note that this eps is a distance. cluster_snn() has an argument of the same name that counts shared neighbours, because the algorithm it comes from named it that; the two are not comparable.

Which margin is clustered

cluster_scale = "samples", the default, puts one point per row of data, and "features" puts one point per column. See cluster_hclust(), which documents the same argument and why transposing data by hand is a different analysis.

See Also

cluster_snn(), the other density method here, which measures closeness by how many neighbours two points share rather than by a radius, and cluster_kmeans(), which is told the number of groups and places every point.

Examples

## Two well-separated blobs and a couple of strays. The count is not an
## argument, so finding two is the result rather than the setup.
set.seed(1)
blobs <- rbind(
  matrix(rnorm(60, mean = -3), ncol = 2),
  matrix(rnorm(60, mean =  3), ncol = 2),
  cbind(c(-12, 12), c(12, -12))
)
colnames(blobs) <- c("x", "y")
res <- cluster_dbscan(blobs, scale = FALSE)
res
table(res$assignments$cluster)

plot(blobs, col = res$assignments$cluster + 1L, pch = 16,
     main = "cluster 0, in black, is noise")

## On a planted two-group structure, with the radius left to be derived.
sim <- simulate_two_groups(n_feats = 30, n_up = 5, n_down = 5, seed = 3)
by_samp <- cluster_dbscan(sim$args$data)
by_samp$parameters$eps_source
table(cluster = by_samp$assignments$cluster, group = sim$args$group)

Description

Merges the two closest points, then the two closest groups, and keeps going until everything is one group; then cuts the resulting tree so that n_clust groups fall out. What comes back is one cluster label per point, beside the stats::hclust() tree the labels were cut from.

Usage

cluster_hclust(
  data,
  feats = NULL,
  cluster_scale = c("samples", "features"),
  center = TRUE,
  scale = TRUE,
  n_clust = 2,
  dist_method = c("euclidean", "manhattan", "correlation"),
  hclust_method = c("average", "complete", "ward.D2")
)

Arguments

Details

The input is the wide format the comparison functions take: one row per sample and one column per feature. Which margin of it becomes the thing being clustered is cluster_scale, and design$point_type reports the answer.

Value

An object of class sa_cluster, a plain list.

analysis

"hclust".

points

Labels of the things that were clustered — samples or features, as cluster_scale asked — in the row order assignments follows.

design

What was clustered: point_type, the feats kept and any dropped_feats, the counts n_samples, n_used, n_dropped and n_feats, and the two counts of the answer, n_clusters and n_noise. n_noise is always 0 here: a tree places every point.

parameters

The choices as they were used.

assignments

One row per point: points, cluster and silhouette.

clusters

One row per cluster: cluster, size and the mean silhouette of its members.

engine

What computed the clustering.

fit

The stats::hclust() object. This is the slot that is not portable; dropping it leaves an object that writes out as JSON.

metadata

Package version, R version, platform and timestamp.

Which margin is clustered

cluster_scale = "samples", the default, puts one point per row of data, and the question is which samples resemble each other. "features" puts one point per column, and the question is which features move together. The second is the one this and draw_heatmap() agree on: the features are standardised first and the transpose is clustered as it stands, which is also what perform_pca(embedding_scale = "features") reports on.

Transposing data by hand instead is a third analysis. scale() standardises columns, so on the transpose it standardises samples, and the picture looks right while answering a different question. cluster_scale exists so that the transpose never has to be taken by hand.

Choosing the distance and the linkage

dist_method decides what "close" means and hclust_method decides what "close" means for two groups rather than two points. "correlation" is 1 - cor(), so two points a constant apart are at distance zero: use it when the shape of a profile is the question and its level is not. The default "average" linkage is the one that neither chains clusters into strings the way single linkage does nor insists they come out round the way "ward.D2" does.

The tree is built on exactly the distance draw_heatmap() builds its dendrogram from, so a heatmap drawn with the same dist_method and hclust_method is showing this clustering and not a near relative of it.

What is dropped before anything runs

Rows that are not complete and finite across feats go before the distance is measured, and design$n_dropped reports how many. This is the listwise deletion the rest of the package uses; nothing is imputed. A feature that takes a single value cannot be scaled, so with scale = TRUE it is left out with a message and named in design$dropped_feats.

See Also

cluster_kmeans(), which is told the same n_clust and answers without a tree, and cluster_dbscan(), which is not told it at all. draw_heatmap() draws the tree this cuts.

Examples

## Three species, and the clustering was not told there were three.
res <- cluster_hclust(iris[1:4], n_clust = 3)
res
table(cluster = res$assignments$cluster, species = iris$Species)

## The tree is `$fit`, so another cut costs nothing.
table(stats::cutree(res$fit, k = 2))

## On a planted two-group structure, clustering the samples.
sim <- simulate_two_groups(n_feats = 30, n_up = 5, n_down = 5, seed = 3)
by_samp <- cluster_hclust(sim$args$data, n_clust = 2)
table(cluster = by_samp$assignments$cluster, group = sim$args$group)

## `cluster_scale = "features"` asks the other question: which features move
## together. The correlation distance is the one that ignores their levels.
by_feat <- cluster_hclust(sim$args$data, cluster_scale = "features",
                          n_clust = 3, dist_method = "correlation")
by_feat
split(by_feat$points, by_feat$assignments$cluster)

Description

Places n_clust centres, gives every point to its nearest one, moves each centre to the mean of the points it was given, and repeats until nothing moves. What comes back is one cluster label per point, beside the stats::kmeans() object holding the centres.

Usage

cluster_kmeans(
  data,
  feats = NULL,
  cluster_scale = c("samples", "features"),
  center = TRUE,
  scale = TRUE,
  n_clust = 2,
  n_start = 25,
  iter_max = 100,
  seed = NULL
)

Arguments

Details

The input is the wide format the comparison functions take: one row per sample and one column per feature. Which margin of it becomes the thing being clustered is cluster_scale, and design$point_type reports the answer.

Value

An object of class sa_cluster, a plain list. The slots are cluster_hclust()'s, with analysis "kmeans", fit the stats::kmeans() object and design$n_noise always 0, since every point is given to a centre.

Reproducibility

The starting centres are random, so two calls can disagree. Two defences are on by default and neither replaces the other. n_start = 25 runs the whole thing twenty-five times from different starts and keeps the one with the smallest within-cluster sum of squares, which makes a bad local optimum unlikely rather than impossible. seed makes the run exactly repeatable, and it restores the random stream afterwards, so seeding this call does not quietly reseed whatever the caller does next.

Neither is a guarantee that n_clust is the right number. k-means always returns the number of groups it was asked for, and clusters$silhouette is where to look for whether they are groups: a cluster near zero is one whose members are about as close to another cluster as to their own.

Which margin is clustered

cluster_scale = "samples", the default, puts one point per row of data, and "features" puts one point per column. See cluster_hclust(), which documents the same argument and why transposing data by hand is a different analysis.

What is dropped before anything runs

Rows that are not complete and finite across feats go before the centres are placed, and design$n_dropped reports how many. A feature that takes a single value cannot be scaled, so with scale = TRUE it is left out with a message and named in design$dropped_feats.

See Also

cluster_hclust(), which is told the same n_clust and returns the tree it cut, and cluster_dbscan(), which derives the count from the density instead and may refuse to place a point.

Examples

## Three species, and the clustering was not told there were three.
res <- cluster_kmeans(iris[1:4], n_clust = 3, seed = 1)
res
table(cluster = res$assignments$cluster, species = iris$Species)

## The centres are on `$fit`, in the scaled units the clustering ran in.
res$fit$centers

## On a planted two-group structure. The group was never shown to the engine.
sim <- simulate_two_groups(n_feats = 30, n_up = 5, n_down = 5, seed = 3)
by_samp <- cluster_kmeans(sim$args$data, n_clust = 2, seed = 1)
table(cluster = by_samp$assignments$cluster, group = sim$args$group)

## Asking for more groups than there are shows up in the silhouettes rather
## than in an error: k-means returns whatever number it was asked for.
cluster_kmeans(sim$args$data, n_clust = 5, seed = 1)$clusters

Description

Builds a graph in which every point keeps its k nearest neighbours, links two points that share at least eps of them, and grows a cluster out of every point with at least min_pts such links. Points that never joined one are left as noise. How many clusters there are is the answer rather than the question.

Usage

cluster_snn(
  data,
  feats = NULL,
  cluster_scale = c("samples", "features"),
  center = TRUE,
  scale = TRUE,
  k = NULL,
  eps = NULL,
  min_pts = NULL
)

Arguments

Details

The input is the wide format the comparison functions take: one row per sample and one column per feature. Which margin of it becomes the thing being clustered is cluster_scale, and design$point_type reports the answer.

Value

An object of class sa_cluster, a plain list. The slots are cluster_hclust()'s, with analysis "snn", fit the dbscan::sNNclust() object, and design$n_clusters and design$n_noise both answers rather than arguments.

eps here counts neighbours, and in cluster_dbscan() it is a distance

The two functions have an argument of the same name meaning two different things, which is the algorithms' doing rather than this package's: dbscan::dbscan() and dbscan::sNNclust() both call their threshold eps, and renaming one of them here would make this package's documentation disagree with the engine's. In cluster_dbscan() eps is a radius in the units of the data and any positive number is meaningful. Here it is a count of shared neighbours, so it is a whole number between 1 and k, and a value of k / 2 says "half of what each of you considers close is the same points".

Why this rather than a radius

One radius has to be right everywhere, so DBSCAN cannot find a tight cluster and a loose one in the same call: the eps that reaches across the loose one merges the tight one into its surroundings. Shared neighbours are relative, so a point in a sparse region is still close to its own neighbours, and clusters of different densities come out together.

The price is dimension. As the number of variables grows, distances concentrate and neighbour lists start to overlap for no reason, so this method will report structure that is an artefact of the width of the table. clusters$silhouette is the check worth making, and agreement with cluster_dbscan() on the same points is worth more than either alone.

Noise

A point that never joined a dense part of the graph gets cluster 0, has no silhouette and no row in clusters; design$n_noise is the count. Every point being noise is a possible answer and means the overlap asked for is not there, which is usually eps or min_pts being too large for the k in use.

See Also

cluster_dbscan(), the other density method here, which measures closeness with one radius rather than by shared neighbours.

Examples

## Two blobs of different spread, which is the case a single radius handles
## badly and shared neighbours handle well.
set.seed(1)
blobs <- rbind(
  matrix(rnorm(80, mean = -6, sd = 0.4), ncol = 2),
  matrix(rnorm(80, mean =  6, sd = 2.5), ncol = 2)
)
colnames(blobs) <- c("x", "y")
res <- cluster_snn(blobs, scale = FALSE, k = 10)
res
table(res$assignments$cluster)

plot(blobs, col = res$assignments$cluster + 1L, pch = 16,
     main = "cluster 0, in black, is noise")

## On a planted two-group structure, with everything left to be derived.
sim <- simulate_two_groups(n_feats = 30, n_up = 5, n_down = 5, seed = 3)
by_samp <- cluster_snn(sim$args$data)
table(cluster = by_samp$assignments$cluster, group = sim$args$group)

Description

The fit element of an sa_model is the engine object, kept so that the model can predict(). These methods make the rest of what anyone does with a fitted model work on it as well: coef() and summary() read through to the stats::lm() or stats::glm() fit that caret::train() built, so what comes back is what the same call on that model would give.

Usage

## S3 method for class 'sa_fit'
coef(object, ...)

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

## S3 method for class 'sa_fit'
predict(object, newdata = NULL, type = "raw", ...)

Arguments

Details

The numbers are also in ⁠$coefficients⁠, as a data.frame with the confidence interval beside them and, for a logistic regression, the odds ratio, and that table is what coef.sa_model() gives for the result as a whole. These methods are for when an lm-shaped answer is what is wanted instead: a named vector to index, or the summary layout that is already familiar.

On a fit_elastic_net() fit, coef() is the coefficients at the penalty that was chosen, which is the model the result describes and the one predict() predicts from, rather than the whole lambda path glmnet keeps. summary() has nothing of the same shape to return there, since a penalized fit has no standard errors and no residual degrees of freedom to test against; what it does have is in ⁠$coefficients⁠ and ⁠$fit_stats⁠.

A fit_rf() or fit_svm() fit has neither. coef() on one is an error naming the importance table to read instead, rather than the NULL coef() on a randomForest object gives or the S4 indexing error it gives on a ksvm one, and summary() reaches a model with no summary method of its own. print(fit$fit) is the readable answer for both, and everything the result reports is in ⁠$coefficients⁠ and ⁠$fit_stats⁠.

predict() here takes newdata as caret prepared it rather than as the fit received it, which for fit_elastic_net() and fit_svm() is a design matrix and not a data frame of predictors. predict.sa_model() is the method that takes the rows themselves, since the result is what knows which columns they were and what levels its factors had.

predict() is otherwise caret's, with one word added. caret::train() knows type = "raw" for the prediction itself and type = "prob" for the class probabilities, and refuses anything else; "response", which is what stats::glm() and stats::lm() call the prediction on the scale of the outcome, is therefore an error on an object that is plainly a logistic regression. It is accepted here and means the same thing in both models: the predicted value for a linear regression, and for a logistic one the probability of outcome_lv[2], the class every coefficient and odds ratio in the result already describes. The probability is caret's own, the second column of type = "prob", so nothing is computed twice.

Value

From coef(), a named numeric vector of coefficients, NA for a term the fit could not estimate and exactly 0 for one a penalty dropped. From summary(), the summary.lm or summary.glm object of the underlying fit. The Call it prints belongs to caret rather than to you: .outcome ~ . on the predictor frame for a linear model, and empty for a logistic one, which caret fits without a formula. What was fitted to what is what print() on the result says.

From predict(), whatever caret::predict.train() returns for that type, except that type = "response" on a classification gives an unnamed numeric vector of probabilities, one per row of newdata.

Examples

fit <- fit_linear_regression(mtcars, outcome = "mpg",
                             predictors = c("wt", "hp"), cv = FALSE)
coef(fit$fit)
summary(fit$fit)

## The same numbers the result reports in its own table.
fit$coefficients

## On a classification, `type = "response"` is the probability of the second
## level of `outcome_lv`, the one the coefficients describe.
iris2 <- iris[iris$Species != "setosa", ]
clf <- fit_logistic_regression(iris2, outcome = "Species",
                               predictors = "Petal.Length",
                               outcome_lv = c("versicolor", "virginica"),
                               cv = FALSE)
head(predict(clf$fit, newdata = iris2, type = "response"))

## A penalized fit answers with the coefficients at the chosen penalty, the
## ones `$coefficients` reports, rather than with the whole lambda path.
pen <- fit_elastic_net(mtcars, outcome = "mpg", penalty = "lasso",
                       lambda = 0.5, cv = FALSE)
coef(pen$fit)

Description

The coefficient table, object$coefficients itself: one row per term, in the order of object$terms, carrying everything the model estimated about each of them rather than the estimate alone.

Usage

## S3 method for class 'sa_model'
coef(object, ...)

Arguments

Details

The lm-shaped answer, a named numeric vector, is coef(x$fit). Two objects are reached by two calls, so they answer in their own terms rather than both giving the same thing: the engine object answers as the engine does, and the result object answers with the table it was assembled to hold. Ask ⁠$fit⁠ when a vector is what is wanted to index or multiply.

What the table holds depends on the model, and the columns say which kind it is. An unpenalized fit carries the standard error, the statistic, the p-value and the confidence limits, and a classification the odds ratio; a penalized one has no standard error to report and carries selected instead, so is.null(coef(x)$pval) is the test for a fit that cannot be asked for inference. A forest and a support vector machine have no coefficient of any kind, and their table is the importance of each term rather than an effect per unit of it; see fit_rf() and fit_svm(). Every term keeps its row either way — a term a penalty dropped with an estimate of exactly 0, one the engine could not estimate with NA — since a table shorter than the model would be a table that had lost terms rather than one describing terms that went to zero.

Value

The coefficients element of object, a data.frame of one row per term.

See Also

fit_linear_regression(), fit_elastic_net(), fit_rf(), fit_svm(), and coef.sa_fit() for the same question asked of the engine object in ⁠$fit⁠.

Examples

fit <- fit_linear_regression(mtcars, outcome = "mpg",
                             predictors = c("wt", "hp"), cv = FALSE)
coef(fit)

## The named vector is one object further in.
coef(fit$fit)

## A penalized fit answers with what it has: no inference, but `selected`.
pen <- fit_elastic_net(mtcars, outcome = "mpg", penalty = "lasso",
                       lambda = 0.5, cv = FALSE)
coef(pen)

Description

Crosses two categorical variables into a contingency table and returns an asymptotic and an exact test of it side by side, together with the measures of how strong the association is. Nothing is chosen on the user's behalf: reporting both makes disagreement between them visible, which is the situation where the choice between an approximation and an exact enumeration actually matters.

Usage

compare_categorical_groups(
  data,
  category_lv = NULL,
  control_label = NULL,
  paired = FALSE,
  conf_level = 0.95,
  correct = TRUE,
  exact = NULL,
  simulate_p_value = FALSE,
  n_resamples = 9999,
  max_levels = 20L,
  seed = NULL,
  diagnose = TRUE
)

Arguments

Details

This is the one scenario in the package with no feature axis. Every other comparison asks its question once per numeric column and returns a table with one row per column; a contingency table is asked about as a whole, so the result carries the table itself in ⁠$cells⁠ and one row per test in ⁠$tests⁠. That is also why the result is not an sa_comparison and does not go through estimate_significance(): a volcano plot needs a signed effect per feature and this scenario has no feature axis to carry one. estimate_categorical_significance() is the counterpart that reads it, one verdict per cell. See "Why this is not a comparison result".

A row is dropped for one of two reasons, and the two are counted apart. A row naming a level category_lv leaves out was measured and excluded, and appears as design$n_dropped. A row missing a value in any variable was not measured, and appears as design$n_incomplete; a table needs the whole row, so the deletion is listwise.

The warning stats::chisq.test() raises about small expected counts is not passed on. ⁠$diagnostics⁠ states the same fact as a number, and the exact test that does not need the approximation is already in the same result, so a warning would be a less precise version of what is sitting next to it.

Fisher's exact test enumerates every table with the observed margins, and on a large r x c one there are more of those than the algorithm's workspace holds. That is a limit of the enumeration rather than a fault in the data, so ⁠$tests$fisher_test$pval⁠ comes back NA with a message saying so, instead of the whole call failing and taking the chi-square result with it. simulate_p_value = TRUE is the way to get an answer there.

Value

An sa_categorical object: a plain list, so it survives being written out as JSON and read back in another language, with an S3 class on top that supplies print(), plot() and as.table(). Its elements are

analysis

"categorical_comparison".

variables

The variable names, in the order category_lv fixed.

design

category_lv, null, paired, pairing, dim, row_var, col_var, n_used, n_dropped and n_incomplete. null is "independence", "symmetry" or "marginal_homogeneity", and it is what expected and the residuals are read under.

parameters

The analysis choices as used, so exact says which branch of McNemar's test ran rather than what was passed.

cells

One row per cell of the table: row_level, col_level, observed, expected, residual, std_residual, prop_total, prop_row and prop_col. This is the canonical form of the table, and what draw_mosaic_plot() reads. as.table() folds it back into a table(), which is the shape to read it in rather than a second copy to keep.

tests

One one-row data.frame per test, named chisq_test, fisher_test, mcnemar_test or cochran_q, each carrying n_used, statistic, df, pval, lower_conf and upper_conf. There is no pval_adj.

test_info

The method id, a readable label and whether the test is a matched one, per element of tests.

association

One row per measure: measure, estimate, lower_conf, upper_conf. Which measures are defined depends on the design and on the size of the table.

diagnostics

The approximation rule this design rests on, or NULL.

metadata

package_version, r_version, platform and an ISO-8601 timestamp.

One test, three questions

The chi-square statistic answers what look like three different questions, and the arithmetic is the same in each.

Independence

Both variables were measured on one sample: is there any association between them? This is the reading the default label uses.

Homogeneity

One variable says which sample a row came from and the other is what was measured: do the samples share the same distribution over the categories? The sampling scheme differs, the null hypothesis is stated differently, and the test is the same one.

Goodness of fit

One variable against a set of expected proportions. This one is genuinely different, being about a single margin rather than a cross-classification, and it is not implemented here.

Which of the first two a caller means is a fact about their sampling scheme, not about their data, and the data cannot be inspected to find out. So this function does not ask, and both readings are the same call.

Which tests run

Decided by paired and by how many variables category_lv names, so there is no argument naming a test.

A matched design reads the columns as repeated measurements of one thing on the same row, so pairing is by row and there is no id argument. It also requires the levels to be binary: the tests of symmetry that generalise McNemar's test past two levels, Bowker's and Stuart-Maxwell's, are not implemented, and a three-level matched design is an error naming them rather than a table quietly collapsed to two levels.

What the cells are expected to hold

design$null names the hypothesis the whole result is about, and ⁠$cells$expected⁠ is read under it. This matters because a contingency table can be held against more than one null and the three designs above hold it against three different ones.

Under independence a cell is expected at the product of its margins over the total. Under symmetry it is expected at the average of it and its transpose, so the diagonal is expected at exactly what it holds and carries no residual: a pair that answered the same way under both conditions says nothing about which condition raises the response. Under marginal homogeneity every condition is expected to show the pooled response rate, which on a condition-by-response table is the same arithmetic as independence and a different claim about the world.

The residuals follow the same hypothesis, which is what keeps them and the p-value beside them talking about one thing. In a matched 2 x 2 the squared Pearson residuals of the two discordant cells sum to (b - c)^2 / (b + c), which is McNemar's uncorrected statistic exactly; the cell table and the test are the same arithmetic read two ways. draw_mosaic_plot() shades on those residuals, so the picture is about the hypothesis the result tested rather than about independence in every case.

std_residual is the exception. Its variance correction is derived for a table held against its own margins and has no counterpart under symmetry, so it is NA there rather than a number that looks referable to a standard normal and is not.

Why this is not a comparison result

sa_new_comparison() holds every table in a comparison to one row per feature, and estimate_significance() reads one log2fc per row off that axis. Neither exists here. There is one table, so there is one p-value and no multiplicity to adjust across, which is why no pval_adj column is carried. And an association has no sign: cramers_v says how far the table sits from independence but not in which direction, because past a 2 x 2 table there is no single direction to name. On a 2 x 2 table there is one, and odds_ratio reports it.

The result is therefore an sa_categorical, which keeps the vocabulary of a comparison – design, parameters, tests beside test_info, metadata – without claiming a contract it cannot meet. draw_mosaic_plot() is what reads it, in the place draw_volcano_plot() holds for the numeric scenarios.

Every sentence above is about the table. A cell does have both axes: observed / expected is a signed departure defined on a table of any shape, and std_residual beside it is referred to a standard normal, so the cells of one table are a family to adjust across. estimate_categorical_significance() is what reads that axis, and it is a function of ⁠$cells⁠ alone, so no slot is carried for it.

Direction

Set once, by the order of the levels in category_lv, and control_label is the second way of stating it. Every quantity that has a direction follows it.

odds_ratio is above 1 when the second level of each variable occurs with the second level of the other more often than independence would give, and phi_coefficient is above zero in the same situation. Pointing control_label at the other level of either variable inverts both, and pointing it at the other level of both leaves them where they were. In a matched design the second level is the response the conditions are compared on, so odds_ratio_paired is above 1 when the later condition raises it.

Which measures are defined

An independent table always reports cramers_v and contingency_coefficient, and a 2 x 2 one adds phi_coefficient and odds_ratio. A matched 2 x 2 reports odds_ratio_paired, risk_difference_paired and cohens_g, all three read off the discordant cells. Three or more matched conditions report kendalls_w.

Only odds_ratio and the paired measures carry an interval. The others are functions of the chi-square statistic whose sampling distribution has no closed-form interval, so their lower_conf and upper_conf are NA: the columns exist because every result table in the package carries them, not because every number in them is finite.

Every measure is built from the uncorrected chi-square statistic, whatever correct was set to. Yates' correction is about referring a discrete statistic to a continuous distribution, which is a statement about a p-value; letting it into an effect size would make the reported strength of an association depend on a choice made about its tail probability.

References

Pearson, K. (1900). On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling. Philosophical Magazine, 50(302), 157-175.

Fisher, R. A. (1935). The logic of inductive inference. Journal of the Royal Statistical Society, 98(1), 39-82.

McNemar, Q. (1947). Note on the sampling error of the difference between correlated proportions or percentages. Psychometrika, 12(2), 153-157.

Cochran, W. G. (1950). The comparison of percentages in matched samples. Biometrika, 37(3-4), 256-266.

Agresti, A. (2002). Categorical Data Analysis, 2nd ed. Wiley.

See Also

simulate_categorical_groups() for a table whose association is known, estimate_categorical_significance() to reduce the result to a verdict per cell, draw_mosaic_plot() to draw it, and compare_multiple_groups() for the case where the measurement is numeric and only the grouping is categorical.

Examples

## Both columns are categorical, so `data` is the whole input.
smoking <- data.frame(
  smoker = rep(c("y", "n"), each = 60),
  grade  = c(rep(c("high", "mid", "low"), c(10, 20, 30)),
             rep(c("high", "mid", "low"), c(30, 20, 10)))
)
res <- compare_categorical_groups(smoking)
res

## The table itself, and the cells that made the statistic what it is.
as.table(res)
res$cells[c("row_level", "col_level", "observed", "expected", "std_residual")]

## How strong the association is, which no test reports.
res$association

## `category_lv` picks the levels and their order, and drops the rest. Two
## levels of `grade` make a 2 x 2 table, which is where the odds ratio exists.
two_by_two <- compare_categorical_groups(
  smoking,
  category_lv = list(smoker = c("n", "y"), grade = c("low", "high"))
)
subset(two_by_two$association, measure == "odds_ratio")

## `control_label` restates the reference without the levels being retyped,
## and pointing it at one variable inverts the odds ratio.
flipped <- compare_categorical_groups(
  smoking,
  category_lv   = list(smoker = c("n", "y"), grade = c("low", "high")),
  control_label = c(smoker = "y")
)
subset(flipped$association, measure == "odds_ratio")

## A matched design: the columns are the same question asked twice, so pairing
## is by row, the null is symmetry rather than independence, and McNemar's test
## reads only the discordant pairs.
before_after <- data.frame(
  before = rep(c("pass", "fail"), c(20, 30)),
  after  = c(rep(c("pass", "fail"), c(18, 2)), rep(c("pass", "fail"), c(14, 16)))
)
matched <- compare_categorical_groups(before_after, paired = TRUE)
matched$design$null

## Which is what the cell table is read under: the diagonal is expected at
## exactly what it holds, so only the discordant cells carry a residual.
matched$cells[c("row_level", "col_level", "observed", "expected", "residual")]

## A simulated table hands over exactly the arguments this function takes.
sim <- simulate_categorical_groups(n_samples = 400, assoc = 0.4, seed = 1)
fit <- do.call(compare_categorical_groups, sim$args)
cbind(planted = sim$truth$cramers_v,
      estimated = fit$association$estimate[1])

Description

The factorial counterpart of compare_multiple_groups(), and the one place in the package where a scenario function fits one model rather than running every applicable test side by side. Two crossed factors are a two-way ANOVA, three are a three-way ANOVA and more than three are a factorial ANOVA, and those are three names for the same fully crossed linear model rather than three procedures to choose between. Which name applies follows from length(factor_lv), so there is no argument for it; the answer is reported in design$anova_type and in the readable label of the test.

Usage

compare_factorial_groups(
  data,
  feats,
  factors,
  factor_lv = NULL,
  control_label = NULL,
  within = NULL,
  id = NULL,
  conf_level = 0.95,
  ss_type = c("III", "II", "I"),
  posthoc = TRUE,
  posthoc_alpha = 0.05,
  posthoc_scope = c("both", "marginal", "simple"),
  fc_mean = c("arith", "geom"),
  input_scale = c("raw", "log2"),
  p_adjust = "BH",
  diagnose = TRUE
)

Arguments

Details

What multiplies instead is the number of answers. One factor asks whether its levels are alike; crossing a second one asks that of each factor and asks whether either one's effect depends on the other, and those questions fail separately. A single p-value cannot carry them, so the result has an axis the sibling functions do not: ⁠$terms⁠, one row per feature and model term.

Features that cannot be tested do not abort the run. Their rows are filled with NA and all such features are reported together in a single warning.

Value

A sa_comparison object with the layout compare_multiple_groups() returns, minus ⁠$pairwise⁠ and plus ⁠$terms⁠ and ⁠$cells⁠. Its elements are

analysis

"factorial_comparison".

features

The feature names, in the row order every table uses.

design

factor_lv, anova_type, n_factors, group_lv (the cell labels), cell_n, n_empty_cells, paired, pairing, within, n_dropped and unmatched_ids.

parameters

The analysis choices as used, plus n_posthoc, the number of features that entered the pairwise stage.

effect

One row per feature: n_used, n_cells, ref_center, extreme_cell, extreme_center, fold_change and log2fc. The reference cell is the one where every factor sits at its first level, which factor_lv or control_label states, and extreme_cell is whichever cell sits furthest from it on the log2 scale, the same definition simulate_factorial_groups() records in truth$extreme_cell.

tests

One table, anova_test, holding the whole-model F test with one row per feature.

terms

One row per feature and model term, holding features, terms, term_order, n_used, df, ss, ms, f_stat, df_error, eta_sq, partial_eta_sq, log2_effect, pval and pval_adj.

cells

One row per feature and cell of the grid, holding features, one column per factor named after it and holding the level name, cell (the dot-joined label, as in design$group_lv), n, mean, sd and se. mean is the arithmetic cell mean the model was fitted on and se is sqrt(ms_error / n), pooled over the whole model, so the variance of a marginal mean over a set of cells is recoverable from it. This is the table draw_interaction_plot() reads, and the one part of the grid the rest of the result cannot be read backwards into.

posthoc

anova_test, one row per feature and contrast, carrying the post-hoc contract columns plus factor and stratum. Absent when posthoc = FALSE.

test_info

The method id, a readable label naming the ANOVA that ran, and the post-hoc procedure that followed it.

diagnostics

Assumption checks over the cells, or NULL.

metadata

package_version, r_version, platform and an ISO-8601 timestamp.

Two axes

The whole-model test in ⁠$tests$anova_test⁠ asks whether a feature responds to the design at all. That is one question per feature, and it is exactly the one-way ANOVA that treats the cells as groups: a fully crossed model is the cell means model written in another basis, so the two fit the same values and leave the same residuals. Keeping it there is what lets estimate_significance(), print(), draw_forest_plot() and draw_volcano_plot() read a factorial result without knowing that it is one.

Which part of the design a feature responds to is the question a crossed model was fitted to answer, and it has one answer per term. Those live in ⁠$terms⁠, whose terms and term_order columns are the ones simulate_factorial_groups() writes into truth_term, so a result and an answer key merge on c("features", "terms") with neither side renamed.

pval_adj is adjusted within each term, across features. A term is one family: asking of five hundred features whether the treatment matters is five hundred instances of one question, while asking whether the treatment matters and whether it depends on sex are two questions and pooling them would correct each for the other's multiplicity.

There is no ⁠$pairwise⁠ slot. Its keys are contrast labels, and a factorial design has a two-dimensional contrast axis, factor by stratum, that a flat list keyed by label cannot hold without a naming convention. Until there is one, estimate_significance(by = "contrast") reports that the result has no pairwise stage, which is the true statement.

The size of a term

eta_sq and partial_eta_sq say how much of the variance a term accounts for, and neither has a sign, so neither can be read as a direction. log2_effect can: it is the largest ANOVA component of the term, with its sign, taken by decomposing log2() of the same cell centres effect is built from. It is what estimate_significance() puts on the effect axis in a term reading, and so what draw_volcano_plot() plots a term panel against when passed the verdict of estimate_significance(fact, by = "term").

A component is a deviation from what the rest of the model already predicts, not a difference between two levels. A two-level factor whose levels differ by one log2 unit has components -0.5 and +0.5, so ⁠|log2_effect|⁠ is 0.5 where the marginal fold change is 2. Exact and near ties in absolute value (within the factorial tolerance) keep the earlier, reference cell's sign; do not read that sign as a treatment-versus-control direction. The definition is kept because it is the one simulate_factorial_groups() records in truth_term$max_abs_delta, which makes the column scorable against an answer key; the consequence is that cutoffs meant for a fold change, log2fc_cutoff = 1 among them, are stricter here than they look.

Components do not depend on which cell is the reference, since subtracting one cell from every cell cancels out of the deviation. That is the one place the term axis is simpler than effect$log2fc, which is read against the reference cell by definition.

Unbalanced cells

With equal cell sizes the three types of sums of squares are identical and ss_type changes nothing. They part company when the cells are unequal, because the factors are then no longer orthogonal and some of the variation can be attributed to more than one of them.

The default is Type III, each term adjusted for every other term. It is the type whose main effects are statements about the levels rather than about how many observations happened to land in each, which is what makes it agree with the unweighted marginal means the post-hoc stage compares and with the decomposition simulate_factorial_groups() plants. Type II adjusts a term for every term that does not contain it, leaving a main effect unadjusted for the interaction it is part of. Type I is sequential and therefore depends on the order factors was written in; it is what stats::aov() reports, which makes it the type to ask for when checking these numbers against an external implementation on unbalanced data.

A cell holding no usable observation leaves the crossed model with nothing to estimate there, so that feature's rows are NA and the reason is reported in one warning. design$n_empty_cells counts the cells that hold no rows at all.

Marginal contrasts and simple effects

A marginal contrast compares two levels of one factor with the other factors averaged away, which is the comparison a main effect is a statement about. A simple effect compares the same two levels inside one combination of the other factors, and it is the only one of the two that means anything when an interaction is real: if the treatment helps males and harms females, the average of the two is a number that describes nobody. stratum tells them apart, NA for a marginal contrast and the levels of the other factors joined by a dot for a simple one, and factor names the factor being compared. Both columns are the ones truth_contrast carries.

The marginal mean is the unweighted mean of the cell means rather than the mean of the observations, so a level is not pulled towards whichever combination of the other factors was sampled most heavily. Every contrast is scaled by the mean square error of the whole model, which is what makes the pairwise stage consistent with the term tests instead of a second analysis of the same data, and judged against the studentised range over the number of levels of its own factor. Those p-values are already family-wise within a block, so there is no posthoc_p_adjust argument to apply and parameters$posthoc_p_adjust is NA.

Which contrasts are computed is decided term by term. A marginal contrast of a factor runs when that factor's main effect cleared posthoc_alpha, and a simple effect runs when the interaction of the factor with the factors held fixed cleared it. Gating everything on the whole-model test instead would compare the levels of a factor that the model says has no effect, on the strength of a different factor that has one.

References

Fisher, R. A. (1935). The Design of Experiments.

Yates, F. (1934). The analysis of multiple classifications with unequal numbers in the different classes. Journal of the American Statistical Association, 29(185), 51-66.

Tukey, J. W. (1949). Comparing individual means in the analysis of variance. Biometrics, 5(2), 99-114.

Kramer, C. Y. (1956). Extension of multiple range tests to group means with unequal numbers of replications. Biometrics, 12(3), 307-310.

See Also

compare_multiple_groups() for a single factor, simulate_factorial_groups() for data whose answer per term is known, and draw_forest_plot() to draw the result.

Examples

## Two crossed factors, so a two-way ANOVA. The factors name columns of `data`.
res <- compare_factorial_groups(
  data    = warpbreaks,
  feats   = "breaks",
  factors = list(wool = "wool", tension = "tension")
)
res

## The whole-model test, one row per feature: does this feature respond at all.
res$tests$anova_test

## The answer per term, which is what the design was crossed to get.
res$terms

## Every fold change is read against the cell where each factor sits at its
## first level. `control_label` moves one of those levels without the rest
## having to be listed, so the reference cell becomes A.M rather than A.L.
medium <- compare_factorial_groups(
  data          = warpbreaks,
  feats         = "breaks",
  factors       = list(wool = "wool", tension = "tension"),
  control_label = list(tension = "M"),
  posthoc       = FALSE
)
medium$design$group_lv[1]
medium$effect[c("ref_center", "extreme_cell", "log2fc")]

## The whole-model volcano names the reference cell on the x axis.
draw_volcano_plot(estimate_significance(res, log2fc_cutoff = 0.1))

## Marginal contrasts and simple effects, told apart by `stratum`.
subset(res$posthoc$anova_test, factor == "tension" & is.na(stratum))

## A simulated design hands over exactly the arguments this function takes,
## and `truth_term` is the answer key for `$terms`.
sim <- simulate_factorial_groups(n_feats = 12, n_per_cell = 8, seed = 1)
fac <- do.call(compare_factorial_groups, sim$args)
scored <- merge(fac$terms, sim$truth_term, by = c("features", "terms"))
with(scored, table(called = pval_adj <= 0.05, planted = is_effect, terms))

Description

Compares three or more group levels across any number of numeric features and returns the omnibus tests side by side, each followed by the post-hoc procedure that shares its assumptions. As with compare_two_groups(), nothing is chosen on the user's behalf: reporting the parametric, the rank-based and the robust result together makes disagreement between them visible, which is the situation where the choice of test actually matters.

Usage

compare_multiple_groups(
  data,
  feats,
  group,
  group_lv,
  control_label = group_lv[1],
  id = NULL,
  paired = FALSE,
  conf_level = 0.95,
  tr = 0.2,
  posthoc = TRUE,
  posthoc_alpha = 0.05,
  fc_mean = c("arith", "geom"),
  input_scale = c("raw", "log2"),
  p_adjust = "BH",
  posthoc_p_adjust = "holm",
  diagnose = TRUE
)

Arguments

Details

Which family runs depends on paired:

Each one is followed by its own post-hoc procedure, never by a borrowed one: one-way ANOVA by Tukey's HSD, Welch's ANOVA by Games-Howell, the trimmed mean ANOVA by pairwise Yuen tests, Kruskal-Wallis by Dunn's test, repeated measures ANOVA by pairwise paired t-tests and Friedman by Conover's test. A rank-based omnibus test is therefore never followed by a parametric comparison.

Direction is set once, by the order of group_lv, and control_label is the second way of stating that order: the level it names moves to the front and the rest keep the order they were given. The move carries further here than it does in compare_two_groups(), since the reference is the denominator of the fold change and the subtracted side of every post-hoc contrast at the same time, and the remaining levels keep contrasting each other in display order. Naming a level that group_lv already holds is a correction rather than a contradiction, which is why it is accepted here and refused by the models, whose outcome_lv holds the two classes and nothing else.

Features that cannot be tested do not abort the run. Their row is filled with NA and all such features are reported together in a single warning.

Value

A sa_comparison object with the same layout compare_two_groups() returns, plus the posthoc and pairwise slots that only a comparison of three or more levels has. Its elements are

analysis

"multi_group_comparison".

features

The feature names, in the row order every table uses.

design

group_lv, paired, pairing, n_dropped and unmatched_ids.

parameters

The analysis choices as used, plus n_posthoc, the number of features that entered the pairwise stage per test.

effect

One row per feature, described below.

tests

The omnibus tables, one row per feature.

posthoc

One table per omnibus test, one row per feature and pair of levels, carrying features, contrast, group1, group2, n1, n2, estimate, stderr, statistic, df, pval, pval_adj, lower_conf and upper_conf. Absent when posthoc = FALSE.

pairwise

The same numbers one contrast at a time, described below. Absent when posthoc = FALSE.

test_info

Per test, the method id, a readable label and the post-hoc procedure that followed it.

diagnostics

Assumption checks, or NULL.

metadata

package_version, r_version, platform and an ISO-8601 timestamp.

The effect table holds n_used, n_groups, ref_center (the centre of group_lv[1]), extreme_level, extreme_center, fold_change and log2fc. extreme_level is whichever level sits furthest from the reference on the log2 scale, and the ratio puts it over the reference, so a positive log2fc means that level is the higher one. group_lv[1] is the denominator here exactly as it is in compare_two_groups(): the first level is the reference in both. The post-hoc contrasts subtract the reference for the same reason, so a feature the treatment raised is positive in effect$log2fc and in its post-hoc estimate alike. Keeping the column named log2fc means estimate_significance() and draw_volcano_plot() work on a multi-group result unchanged.

Omnibus intervals

The omnibus tables carry lower_conf and upper_conf as required by the result contract, and both are NA. An omnibus test states that the levels are not all alike; it does not state by how much, and there is no single quantity for an interval to be about. The intervals of a multi-group comparison live in ⁠$posthoc⁠, where each row is one contrast and does have a scale of its own.

Post-hoc stage

Only features whose omnibus pval_adj clears posthoc_alpha are compared pairwise, and a feature that did not qualify is absent from the post-hoc table rather than present with NA. An absent row means the question was never asked; an NA row means it was asked and could not be answered, and the two should not look the same. parameters$n_posthoc records how many features entered each stage.

Tukey's HSD and Games-Howell control the error rate over the whole set of contrasts through the studentised range, so pval_adj equals pval for those two and posthoc_p_adjust is not applied. Dunn, Conover, pairwise Yuen and pairwise paired t-tests are adjusted across the contrasts within each feature.

Reading one contrast at a time

⁠$posthoc⁠ stacks every contrast of every feature into one long table, which is the honest record of what was asked but an awkward shape for a reader who came for a single comparison. ⁠$pairwise⁠ holds the same numbers keyed first by test and then by contrast, so res$pairwise$anova_test[["virginica - setosa"]] is that one comparison across all features.

Those tables are rectangular where ⁠$posthoc⁠ is ragged: each holds every feature, in the order features fixes, so contrasts can be lined up against each other and against the omnibus tables by position. A feature that did not qualify is present with its inference columns NA.

They add fold_change and log2fc, which no post-hoc procedure reports, being the ratio of the two group centres rather than anything a test produced. The ratio puts group1 over group2, matching the direction estimate reads in, so the two always agree in sign. Note that this is a different quantity from effect$log2fc, which compares the most extreme level rather than a named pair, though the two now point the same way: both divide by the reference. A feature that never entered the pairwise stage still has its ratio: dividing two centres does not require a test to have been run.

estimate_significance() reads these tables with by = "contrast", whose significance element is then one verdict table per contrast rather than one table.

Repeated conditions

A within-subject omnibus test needs a complete rectangle, so id is required and subjects missing any condition are dropped whole rather than partially used. The number dropped is reported in design$unmatched_ids. Missing values are then handled per feature: a subject with NA on one feature is left out of that feature only, which is why n_used can differ between features.

Row order pairing, which compare_two_groups() allows, is deliberately not offered here. With two groups it is at least well defined; with three or more it would also have to assume every condition is stored in the same subject order, and there is no way to notice when it is not.

References

Welch, B. L. (1951). On the comparison of several mean values: an alternative approach. Biometrika, 38(3-4), 330-336.

Kruskal, W. H. and Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. Journal of the American Statistical Association, 47(260), 583-621.

Tukey, J. W. (1949). Comparing individual means in the analysis of variance. Biometrics, 5(2), 99-114.

Games, P. A. and 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), 113-125.

Dunn, O. J. (1964). Multiple comparisons using rank sums. Technometrics, 6(3), 241-252.

Conover, W. J. (1999). Practical Nonparametric Statistics, 3rd edition.

See Also

compare_two_groups() for exactly two levels, diagnose_distribution() for the assumption checks on their own, and draw_forest_plot() to draw the result.

Examples

## Independent samples: all three iris species
feats <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
res <- compare_multiple_groups(
  data     = iris,
  feats    = feats,
  group    = iris$Species,
  group_lv = c("setosa", "versicolor", "virginica")
)
res
res$tests$anova_test
res$tests$kruskal_test

## The pairwise stage, one row per feature and pair
head(res$posthoc$anova_test)

## The same numbers one contrast at a time, every feature present. Setosa is
## the reference, so it is the level every contrast subtracts.
names(res$pairwise$anova_test)
res$pairwise$anova_test[["virginica - setosa"]]

## Setosa is the reference, so a positive log2fc means the most extreme
## other species is the larger one.
res$effect

## Repeated conditions: each ChickWeight chick weighed at several times
chicks <- ChickWeight[ChickWeight$Time %in% c(0, 6, 12, 18), ]
rep_res <- compare_multiple_groups(
  data     = data.frame(weight = chicks$weight),
  feats    = "weight",
  group    = paste0("day", chicks$Time),
  group_lv = c("day0", "day6", "day12", "day18"),
  id       = chicks$Chick,
  paired   = TRUE
)
rep_res$tests$anova_test   # Mauchly, Greenhouse-Geisser and Huynh-Feldt too
rep_res$posthoc$kruskal_test

Description

Tests each feature against mu and returns a parametric, a rank-based and a proportion result side by side, in the same shape compare_two_groups() uses. There is no second group here, so the reference is a number the user supplies rather than a set of observations.

Usage

compare_one_sample(
  data,
  feats,
  mu = 0,
  p = 0.5,
  success = 1,
  alternative = c("two.sided", "less", "greater"),
  conf_level = 0.95,
  fc_mean = c("arith", "geom"),
  input_scale = c("raw", "log2"),
  p_adjust = "BH",
  diagnose = TRUE
)

Arguments

Details

t_test

One-sample t-test on the mean.

wilcox_test

One-sample Wilcoxon signed-rank test, with the Hodges-Lehmann pseudo-median. It tests location without assuming normality, but it does assume the distribution of x - mu is symmetric.

prop_test

Score test of the proportion of successes against p. It only applies to features that are binary, and features that are not come back as NA rather than being silently coerced.

Missing values are dropped per feature, so n_used varies between features when the data are incomplete. Features that cannot be tested do not abort the run: their row is filled with NA and all such features are reported together in a single warning.

Value

A sa_comparison object with the layout described in compare_two_groups(), with two differences. design carries mu, p and success instead of group_lv, since there are no groups, and effect holds n_used, center (the centre fc_mean selected), mu, diff, fold_change and log2fc, the ratio being that centre over mu. Every column of effect is on the original measurement scale, so with input_scale = "log2" its center, mu and diff are back-transformed and differ from the same-named columns of tests$t_test, which stay on the scale the tests ran on. There is no posthoc or pairwise slot: a single sample has no pair of levels to contrast.

Every test table starts with n_used and carries pval, pval_adj, lower_conf and upper_conf. The remaining columns are:

t_test

center, mu, diff, stderr, t_stat, df and cohens_d, the mean difference over the sample standard deviation.

wilcox_test

hl_shift, the Hodges-Lehmann pseudo-median, and v_stat.

prop_test

n_success, proportion, p, diff, chi_sq, df and cohens_h. The interval is a Wilson score interval, which stays inside ⁠[0, 1]⁠ where a Wald interval does not.

Fold change against a hypothesised value

fold_change is center / mu and log2fc its base-2 logarithm, so both are undefined when mu is zero, which is also its most common value. Both columns are NA in that case rather than infinite, a message says so, and estimate_significance() will call every feature undecided. Reporting Inf would read as an infinitely large increase when what actually happened is that the question has no answer.

This cannot arise under input_scale = "log2". There the reference is 2^mu, which is positive whatever mu is, and mu = 0 simply means a reference of 1. Under the default fc_mean = "geom" the whole row reduces to log2fc = mean(x) - mu.

References

Student (1908). The probable error of a mean. Biometrika, 6(1), 1-25.

Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1(6), 80-83.

Wilson, E. B. (1927). Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association, 22(158), 209-212.

See Also

compare_two_groups() for two groups and compare_multiple_groups() for three or more.

Examples

## Are the iris measurements different from 3 cm?
res <- compare_one_sample(
  data  = iris,
  feats = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"),
  mu    = 3
)
res
res$tests$t_test
res$tests$wilcox_test
res$effect

## A binary feature reaches the proportion test; a continuous one does not
## and comes back NA with a warning rather than being coerced.
cars <- mtcars
cars$am <- as.numeric(cars$am)
suppressWarnings(
  compare_one_sample(cars, c("am", "vs", "mpg"), mu = 0.5, p = 0.5)
)$tests$prop_test

## One-sided: is mileage above 20 mpg?
compare_one_sample(mtcars, "mpg", mu = 20,
                   alternative = "greater")$tests$t_test

Description

Compares exactly two group levels across any number of numeric features and returns a parametric, a rank-based and a robust test side by side, together with the fold change between the two groups. Nothing is chosen on the user's behalf: reporting all of them makes disagreement between them visible, which is the situation where the choice of test actually matters.

Usage

compare_two_groups(
  data,
  feats,
  group,
  group_lv,
  control_label = group_lv[1],
  id = NULL,
  alternative = c("two.sided", "less", "greater"),
  paired = FALSE,
  conf_level = 0.95,
  tr = 0.2,
  fc_mean = c("arith", "geom"),
  input_scale = c("raw", "log2"),
  p_adjust = "BH",
  diagnose = TRUE
)

Arguments

Details

Which member of each family is used depends on paired:

Direction is set once, by the order of group_lv, and every quantity in the result follows it. The first level is the reference, as it is in compare_multiple_groups(), so alternative = "greater" tests whether group_lv[2] exceeds group_lv[1] in all three families, and mean_diff, hl_shift, trim_diff, fold_change and relative_effect are all above their null value when group_lv[2] is the larger group.

control_label is the second way of stating that order. It names the level to hold as the reference and moves it to the front, so the direction can be flipped without the two levels being retyped in the other order. Because group_lv says which levels take part as well as which one is the reference, naming the second one here is a correction rather than a contradiction, which is where the comparisons differ from the models: an outcome_lv holds nothing but the two classes, so fit_logistic_regression() rejects a control_label that disagrees with it.

Features that cannot be tested do not abort the run. Their row is filled with NA and all such features are reported together in a single warning. Informational engine warnings, such as an exact p-value being unavailable because of ties, are grouped into one message().

Missing values are handled per feature: independent samples drop NA within each group, paired samples keep only complete pairs. n_x, n_y and n_used therefore vary between features when the data are incomplete. The fold change is computed from those same observations, so it never rests on a different subset of the data than the p-value beside it.

Value

A sa_comparison object: a plain list, so it survives being written out as JSON and read back in another language, with an S3 class on top that only supplies print(). Its elements are

analysis

"two_group_comparison".

features

The feature names, in the row order every table uses.

design

group_lv, paired, pairing ("order", "id" or NA when not paired), n_dropped (rows removed for belonging to a level outside group_lv) and unmatched_ids.

parameters

alternative, conf_level, tr, fc_mean, input_scale and p_adjust, as used.

effect

One row per feature: x_center, y_center (the centres fc_mean selected for group_lv[2] and for the reference group_lv[1]), fold_change and log2fc.

tests

t_test, wilcox_test and robust_test, described below.

test_info

Per test, the method id, a readable label and whether it was the paired variant.

diagnostics

Assumption checks, or NULL when diagnose is FALSE.

metadata

package_version, r_version, platform and an ISO-8601 timestamp.

There is no posthoc or pairwise slot. With two groups the omnibus comparison is already the only contrast there is, so the question a post-hoc stage answers never arises. compare_multiple_groups() is where those two slots appear.

Every table in tests has one row per feature and starts with features, the per-group sample sizes n_x / n_y and n_used (total observations for independent samples, complete pairs for paired samples), and carries pval, pval_adj, lower_conf and upper_conf. Every column named for x or y reads the same way throughout: x is group_lv[2] and y is the reference group_lv[1]. The remaining columns are:

t_test

x_mean, y_mean, mean_diff, stderr, t_stat, df. Group means are computed directly, so the columns are identical for paired and independent designs.

wilcox_test

hl_shift (Hodges-Lehmann location shift, the pseudo-median of differences when paired) and w_stat.

robust_test, independent

relative_effect (P(X > Y) + 0.5 * P(X = Y), above 0.5 when group_lv[2] is the larger group), bm_stat and df. The interval is on the probability scale, so a one-sided alternative leaves it open at 0 or at 1 rather than at infinity.

robust_test, paired

x_trim_mean, y_trim_mean, trim_diff, stderr, yuen_stat, df and robust_dz, a robust counterpart of Cohen's dz.

Pairing

With paired = TRUE and no id, the only available information is row order: the first row of group_lv[1] is matched with the first row of group_lv[2], and so on. Both groups must then have the same number of rows. Note that this cannot detect rows that have been reordered, so a data set sorted differently in each group would produce wrong pairs and no complaint. Passing id removes that failure mode: pairs are matched on the key, ids appearing in only one group are dropped with a message, and an id repeated within a group is an error.

Log-transformed input

A ratio only means something on the scale the measurement was made on. Dividing two means of already logged values answers a different question and can even come out with the wrong sign: two log2 centres of -1 and -2 are a two-fold increase, but their ratio, 0.5, reads as a two-fold decrease. With input_scale = "log2" each observation is raised back through 2^x before the centres are taken, so fold_change and log2fc mean the same thing they do for raw input and fold_change still equals x_center / y_center.

Under the default fc_mean = "geom" this reduces exactly to log2fc = mean(x) - mean(y), the usual definition for log-scale data. fc_mean = "arith" is accepted but averages the back-transformed values, which the largest observations dominate, and is not the conventional log2 fold change.

Only the effect table is converted. The tests run on the values as supplied, which is the reason for logging them in the first place, so with input_scale = "log2" the centres in effect are on the original scale while x_mean and y_mean in tests$t_test are on the log2 scale. For raw input with fc_mean = "arith" those columns hold the same numbers; here they do not.

References

Welch, B. L. (1947). The generalization of Student's problem when several different population variances are involved. Biometrika, 34(1-2), 28-35.

Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1(6), 80-83.

Mann, H. B. and Whitney, D. R. (1947). On a test of whether one of two random variables is stochastically larger than the other. Annals of Mathematical Statistics, 18(1), 50-60.

Brunner, E. and Munzel, U. (2000). The nonparametric Behrens-Fisher problem: asymptotic theory and a small-sample approximation. Biometrical Journal, 42(1), 17-25.

Yuen, K. K. (1974). The two-sample trimmed t for unequal population variances. Biometrika, 61(1), 165-170.

See Also

estimate_significance() to reduce the result to one significance verdict per feature, draw_grouped_boxplot() to visualise the same input, and simulate_two_groups() for input whose answer is known in advance.

Examples

## Independent samples: two of the three iris species
iris2 <- iris[iris$Species != "setosa", ]
res <- compare_two_groups(
  data     = iris2,
  feats    = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"),
  group    = iris2$Species,
  group_lv = c("versicolor", "virginica")
)
res
res$tests$t_test
res$tests$robust_test
res$effect

## setosa is perfectly separated from the others on petal size, which leaves
## the Brunner-Munzel variance at zero. Those rows come back NA with a
## warning rather than aborting the other features.
suppressWarnings(
  compare_two_groups(
    data     = iris[iris$Species != "virginica", ],
    feats    = c("Sepal.Length", "Petal.Length"),
    group    = iris$Species[iris$Species != "virginica"],
    group_lv = c("setosa", "versicolor")
  )$tests$robust_test
)

## Paired samples: sleep holds 10 subjects under 2 drugs, listed in the same
## subject order within each group.
paired_res <- compare_two_groups(
  data     = sleep["extra"],
  feats    = "extra",
  group    = sleep$group,
  group_lv = c("1", "2"),
  paired   = TRUE,
  alternative = "less"
)
paired_res$tests$robust_test

## Same data with the second group shuffled. Row order pairing silently uses
## the wrong partners, while `id` recovers the correct result.
shuffled <- rbind(sleep[1:10, ], sleep[10 + c(4, 9, 1, 7, 2, 10, 3, 6, 8, 5), ])
by_order <- compare_two_groups(
  data = shuffled["extra"], feats = "extra",
  group = shuffled$group, group_lv = c("1", "2"), paired = TRUE
)
by_id <- compare_two_groups(
  data = shuffled["extra"], feats = "extra",
  group = shuffled$group, group_lv = c("1", "2"),
  id = shuffled$ID, paired = TRUE
)
## The means agree but the paired standard error, and so the p-value, do not.
rbind(order = by_order$tests$t_test,
      id    = by_id$tests$t_test)

Description

Runs the normality tests, the homogeneity of variance tests and the outlier screen together, and reports them as one object. The same checks are attached automatically to compare_two_groups() and compare_multiple_groups(), so an assumption is never silently ignored; this function is for looking at them on their own, before a test has been chosen.

Usage

diagnose_distribution(
  data,
  feats,
  group = NULL,
  group_lv = NULL,
  alpha = 0.05,
  criterion = c("iqr", "robust_z", "grubbs"),
  iqr_multiplier = 1.5,
  z_threshold = 3.5,
  center = c("median", "mean", "trimmed"),
  trim = 0.1
)

Arguments

Details

A failed assumption never blocks an analysis and never causes a test to be swapped for another one. It changes which member of the reported test family deserves the most weight, and that judgement stays with the user. Skewed groups make the rank-based and robust members more trustworthy than the parametric one; unequal variances make Welch's and Brunner-Munzel's treatments of the same data more trustworthy than the pooled ones.

Each assumption is checked twice on purpose, by tests that fail differently:

Normality

Shapiro-Wilk is the more powerful of the two and is the one to read first. The Kolmogorov-Smirnov test is fitted against a normal with the sample's own mean and standard deviation, which makes its p-value anti-conservative, so it disagreeing with Shapiro-Wilk usually means the departure is in the tails.

Homogeneity of variance

The Levene test is centred on the median, the Brown-Forsythe variant, and tolerates skew. Bartlett's test is more powerful when the groups really are normal and cannot tell unequal variances apart from heavy tails when they are not. The two disagreeing is itself evidence about normality.

Value

A sa_diagnosis object: a plain list carrying analysis, features, design, parameters, metadata and four tables.

normality

One row per feature and group level: features, group, n_used, shapiro_stat, shapiro_pval, ks_stat, ks_pval, skewness and excess_kurtosis. The level column is named group to match summarize_descriptive_stats(), and is NA when no grouping was supplied.

variance

One row per feature: features, n_used, n_groups, levene_stat, levene_df1, levene_df2, levene_pval, bartlett_stat, bartlett_df and bartlett_pval. Zero rows when no group was supplied.

outliers

The screen_outliers() table, one row per flagged observation.

summary

One row per feature: features, n_levels, n_outliers, min_shapiro_pval, normal_ok and variance_ok. The two flags are NA when the corresponding test could not be run.

The result is deliberately not an sa_comparison. Normality is a property of one sample and homogeneity a property of a set of them, so the two tables have different numbers of rows and would not fit a contract built around one row per feature.

References

Shapiro, S. S. and Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3-4), 591-611.

Massey, F. J. (1951). The Kolmogorov-Smirnov test for goodness of fit. Journal of the American Statistical Association, 46(253), 68-78.

Brown, M. B. and Forsythe, A. B. (1974). Robust tests for the equality of variances. Journal of the American Statistical Association, 69(346), 364-367.

Bartlett, M. S. (1937). Properties of sufficiency and statistical tests. Proceedings of the Royal Society A, 160(901), 268-282.

See Also

screen_outliers() for the outlier stage on its own, and compare_multiple_groups(), whose ⁠$diagnostics⁠ slot holds the same tables for the data it tested.

Examples

feats <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")

## Within species, so a genuine species difference is not read as skew
d <- diagnose_distribution(iris, feats, iris$Species)
d
d$normality
d$variance
d$summary

## Petal length pooled over species is strongly bimodal, which every
## normality test picks up once the grouping is taken away.
diagnose_distribution(iris, "Petal.Length")$normality

Description

Draws the distribution of a single feature for two group levels back to back, on a shared set of breaks. The first level runs left from the centre line and the second one right, so the two shapes can be compared bin by bin instead of read off two separate panels.

Usage

draw_butterfly_hist(
  data,
  feat,
  group,
  group_lv,
  breaks = "Sturges",
  scale = c("count", "proportion", "density"),
  type = c("freq", "dens", "both"),
  dens_adjust = 1,
  dens_lwd = 2,
  dens_col = NULL,
  dens_alpha = 0.45,
  col = c("#4575B4", "#D73027"),
  border = "white",
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  cex.lab = 1.3,
  cex.axis = 1.2,
  cex.main = 1.3,
  cex.legend = 1.1,
  legend.position = "topright",
  xlim = NULL,
  ylim = NULL,
  margin = c(5, 5, 4, 3),
  out_statistics = TRUE,
  ...
)

Arguments

Details

The left group is drawn at negative coordinates so that both distributions share one axis, but the tick labels are the absolute values, so a bar is read the same way on either side.

A density curve and a bar can only be read against the same axis when the bar is a density too: a count or a proportion per bin scales with the bin width, which the curve knows nothing about. So type = "dens" and type = "both" move scale to "density", and reject a scale that was asked for explicitly and says otherwise, rather than drawing two incomparable shapes.

The plot starts from the device background rather than imposing a theme of its own, as draw_volcano_plot() does. Only two colours are needed here, so there is no dark variant; that is reserved for draw_grouped_boxplot(), where three or more group colours have to stay apart.

Only the margins are restored on exit, not every graphical parameter, so the coordinate system survives the call and graphics::abline(), graphics::points() and friends can still add to the finished plot.

Value

If out_statistics = FALSE, NULL invisibly. Otherwise a list, invisibly:

bin_summary_stats

One row per bin, with bin_start, bin_end, bin_mid and one column per group level holding the plotted bar length. What that length means follows scale.

group_summary_stats

One column per group level with rows n (finite values used), n_dropped (missing or non-finite values left out), min and max.

group_hists

One "histogram" object per group level, named by the level, exactly as graphics::hist() returns them, so plot(res$group_hists[[lv]]) redraws that group on its own. Both groups are binned on the shared breaks, so these are not what calling graphics::hist() on one group alone would give: that would pick its own breaks from that group only.

group_densities

Present only when type is not "freq". One "density" object per group level, named by the level, as stats::density() returns them. The bins are still reported in bin_summary_stats even when no bars were drawn.

See Also

draw_grouped_boxplot() for several features at once, and compare_two_groups() to test the difference the plot shows.

Examples

iris2 <- iris[iris$Species != "setosa", ]

draw_butterfly_hist(
  data     = iris2,
  feat     = "Petal.Length",
  group    = iris2$Species,
  group_lv = c("versicolor", "virginica"),
  breaks   = 12,
  scale    = "proportion",
  main     = "virginica vs versicolor"
)

## Counts, plus the numbers behind the bars
res <- draw_butterfly_hist(iris2, "Sepal.Width", iris2$Species,
                           c("versicolor", "virginica"))
res$bin_summary_stats
res$group_summary_stats

## Each group also comes back as a plain histogram object
plot(res$group_hists$virginica)

## Bars with a kernel density estimate on top, both on the density scale
res <- draw_butterfly_hist(
  data     = iris2,
  feat     = "Petal.Length",
  group    = iris2$Species,
  group_lv = c("versicolor", "virginica"),
  breaks   = 12,
  type     = "both",
  main     = "Petal length, bars and density"
)
res$group_densities$virginica

## The smoothed shapes on their own
draw_butterfly_hist(iris2, "Petal.Length", iris2$Species,
                    c("versicolor", "virginica"), type = "dens")

Description

A corrplot: one cell per pair of features, coloured by the coefficient on a fixed -1 to 1 scale, with both axes in one order so that the diagonal runs corner to corner and blocks of features that move together sit next to each other. Given the p-values as well, the pairs that did not clear sig_level are drawn as blank cells, so that what is left coloured is what there is evidence for.

Usage

draw_corrplot(
  cor_matrix,
  method = NULL,
  pvalue = NULL,
  use_adjusted = TRUE,
  sig_level = 0.05,
  cluster = TRUE,
  hclust_method = c("average", "complete", "ward.D2"),
  zlim = c(-1, 1),
  anno = TRUE,
  main = NULL,
  cex.anno = 1,
  cex.axis = 0.9,
  cex.main = 1.5,
  cex.legend = 1.2,
  ...
)

Arguments

Details

draw_heatmap() draws it. This function decides what it is handed: the matrix unscaled, the colour range fixed at the range a correlation can take, one clustering shared by the two axes, and the blanking applied afterwards.

The distance the clustering runs on is 1 - cor(), the same one cluster_hclust() and draw_heatmap() mean by dist_method = "correlation", so a corrplot and a heatmap of the same features group them the same way. It is computed once and both axes are permuted by it. A matrix holding an NA, which is what a feature with no variance leaves behind, has no distance for that feature, and rather than fail the features are left in their input order with a message saying so.

Blanking happens after the clustering rather than before it. A cell removed for its p-value would otherwise change the tree, and the picture would no longer be the matrix the reader is being shown.

A cell whose p-value is NA, a pair that could not be tested, is left as it arrived rather than blanked: there is no evidence against it either, and the coefficient beside it is usually already NA.

The matrix is symmetric, so the transpose draw_heatmap() takes on the way in leaves it unchanged and the features come out on both axes.

Value

A list, invisibly: everything draw_heatmap() returns, and beside it

corr

The matrix as it was drawn, in the drawn order and with the blanked cells NA.

pvalue

The p-values in that same order, or NULL.

order

The permutation of the input the clustering chose.

hclust

The stats::hclust() object behind it, or NULL when the features were not clustered.

n_masked

How many cells were blanked.

See Also

summarize_association_stats(), which produces what this draws, and draw_heatmap(), which draws it.

Examples

feats <- c("mpg", "cyl", "disp", "hp", "drat", "wt")
res <- summarize_association_stats(mtcars, feats, methods = "pearson")

## Every pair, with the features clustered so that the blocks sit together
drawn <- draw_corrplot(res, main = "mtcars")
drawn$order

## Only the pairs that cleared the adjustment at 1%
draw_corrplot(res, sig_level = 0.01, cex.anno = 0.8, cex.axis = 0.7)

## A bare matrix, in the order it arrives and with nothing to blank
draw_corrplot(stats::cor(mtcars[feats]), cluster = FALSE)

Description

Plots two coordinates of a perform_pca(), perform_tsne() or perform_umap() result against each other. A clustering of the same frame colours the points and a known grouping shapes them, so what the data was found to say and what the caller already knew are read off one picture.

Usage

draw_dim_reduction_plot(
  reduction_result,
  group = NULL,
  group_lv = NULL,
  cluster_result = NULL,
  cluster_lv = NULL,
  dims = c(1L, 2L),
  anno_points = FALSE,
  dark = FALSE,
  asp = NULL,
  col = NULL,
  pch = NULL,
  cex = 1.2,
  xlim = NULL,
  ylim = NULL,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  cex.axis = 1.2,
  cex.lab = 1.3,
  cex.main = 1.3,
  cex.legend = 1.1,
  cex.anno = NULL
)

## S3 method for class 'sa_reduction'
plot(x, ...)

Arguments

Value

A data.frame of the points as they were drawn, invisibly: points, x, y, col, pch, and cluster and group when those were given. The resolved colouring is carried as a "view" attribute, one of "both", "cluster", "group" or "plain".

Colour and shape say two different things

cluster_result takes the colours and group takes the shapes when both are given, and giving both is the point of the arrangement rather than a conflict to be resolved. With one channel alone, that channel takes the colours when col is named; otherwise a lone clustering is coloured and a lone grouping is shaped. The question a clustering usually raises is whether it recovered a grouping that was known all along, and on two channels that is read directly: one colour per shape is a clustering that found the groups, and a shape split across colours is a group the data does not see as one thing. Give neither and the points are drawn in the foreground colour.

The clustering has to be of the same points, which the two contracts already promise: sa_cluster and sa_reduction both read their input through the same function, so a mismatch is refused rather than lined up by position.

Noise

cluster_dbscan() and cluster_snn() can leave a point in no cluster at all. Those points are grey rather than a palette colour, since a point left out is the absence of a cluster and not a cluster of its own, and the legend counts them on a line of their own.

The axes

A principal component analysis reports what share of the variance each of its components carries, so its axis labels carry it too: ⁠PC1 (30.2%)⁠ is read from ⁠$variance⁠ rather than recomputed. An embedding has no such number and its axes are labelled with their names alone. asp = 1 is what makes one unit of the vertical axis the same length as one unit of the horizontal, which is worth setting when the distance between two points is what is being read.

See Also

perform_pca() for the coordinates and cluster_kmeans() for the labels, and draw_heatmap(), which shows the same wide input a cell at a time rather than a point at a time.

Examples

res <- perform_pca(iris[1:4])

## What was known all along, as shapes. Name `col` to colour the levels
## instead.
draw_dim_reduction_plot(res, group = iris$Species)

draw_dim_reduction_plot(res, group = iris$Species,
                      col = c("#E69F00", "#56B4E9", "#009E73"))

## What the data was found to say, as colours. The clustering was never shown
## the species, so one colour per shape is its own finding.
cl <- cluster_kmeans(iris[1:4], n_clust = 3, seed = 1)
drawn <- draw_dim_reduction_plot(res, group = iris$Species,
                                 cluster_result = cl)
table(cluster = drawn$cluster, species = drawn$group)

## Names for the clusters the legend would otherwise number.
draw_dim_reduction_plot(res, cluster_result = cl,
                      cluster_lv = c("A", "B", "C"))

## A density method can place no point at all, and those are grey.
db <- cluster_dbscan(iris[1:4])
draw_dim_reduction_plot(res, cluster_result = db, asp = 1)

## The features as the points, labelled, on the third and fourth components.
by_feat <- perform_pca(iris[1:4], embedding_scale = "features")
draw_dim_reduction_plot(by_feat, dims = c(1, 3), anno_points = TRUE)

Description

Draws the estimate of each feature beside its confidence interval, falling back to a bar of -log10(pval_adj) for a table that has no interval to draw, which is every omnibus test.

Usage

draw_forest_plot(
  comparison_result,
  test = names(comparison_result$tests)[1],
  type = c("auto", "estimate", "posthoc", "pvalue"),
  feats = NULL,
  use_adjusted = TRUE,
  alpha = 0.05,
  sort_by = c("none", "pvalue"),
  dark = FALSE,
  xlim = NULL,
  xlab = NULL,
  main = NULL,
  col_signif = "#D1495B",
  col_plain = "#7F8C8D",
  cex.axis = 0.9,
  cex.lab = 1.1,
  cex.main = 1.2,
  cex.legend = 1,
  ...
)

## S3 method for class 'sa_comparison'
plot(x, ...)

Arguments

Details

Only the columns the result contract guarantees are read, which is why one function covers every scenario: it never asks whether the object came from two groups, three groups or a single sample, only whether the table it was handed has intervals or p-values, and every table has one or the other. plot() on an sa_comparison is the same function under the name R users reach for first.

Three views are available:

"estimate"

A forest plot of the effect estimate and its confidence interval, one row per feature. Available whenever the chosen table has finite intervals, which the two-group and one-sample scenarios always do.

"posthoc"

The same forest plot for the pairwise contrasts of a multi-group comparison, one row per contrast. Only the first feature of the post-hoc table is drawn unless feats names the ones to draw. An estimate reads as group1 - group2, the direction the row label spells out, and the reference level is the one subtracted, so a point to the right of the guide agrees in sign with the log2fc a volcano plot of the same comparison draws.

"pvalue"

-log10() of the p-value per feature with the alpha threshold marked. The fallback when a table has no interval to draw, which is the case for every omnibus test.

type = "auto", the default, picks the first of those three that the chosen table can actually support.

The function changes graphical parameters and the panel layout, and restores both on exit, so the caller's device is left as it was found. The legend lives in a narrow panel of its own on the right, where it cannot cover the rows it describes.

The estimate view marks the null value with a vertical line, at zero for a difference and at one for the Brunner-Munzel relative effect, which is the one quantity in the package whose null is not zero.

Value

The plotted data.frame, invisibly, in the row order it was drawn, with the view that type = "auto" resolved to attached as the attribute "view".

See Also

draw_volcano_plot(), which plots effect size against significance rather than estimate against interval, and draw_grouped_boxplot() for the input data.

Examples

iris2 <- iris[iris$Species != "setosa", ]
feats <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
res <- compare_two_groups(iris2, feats, iris2$Species,
                          c("virginica", "versicolor"))

## Mean differences with their intervals
draw_forest_plot(res)

## plot() reaches the same function
plot(res)

## The same comparison read through the rank-based test
draw_forest_plot(res, test = "wilcox_test", sort_by = "pvalue")

## Two features only, drawn in the order they are named, against the
## unadjusted p-value
draw_forest_plot(res, feats = c("Petal.Width", "Petal.Length"),
                 use_adjusted = FALSE)

## An axis range of your own, so that two plots can be read against each other
draw_forest_plot(res, xlim = c(0, 2))

## A multi-group omnibus table has no interval, so "auto" falls through to
## the pairwise contrasts.
multi <- compare_multiple_groups(iris, feats, iris$Species,
                                 levels(iris$Species))
draw_forest_plot(multi, feats = "Petal.Length")

## Several features at once label each contrast with the feature it belongs to
draw_forest_plot(multi, feats = c("Petal.Length", "Sepal.Width"))
draw_forest_plot(multi, type = "pvalue")

Description

Draws one bar per feature and group level, the levels side by side inside each feature's cluster and a legend in a narrow panel on the right. The bar heights are one column of summarize_descriptive_stats(), read from the same wide input the comparison functions take, so a bar and a row of that table are the same number and neither has to be recomputed to check the other.

Usage

draw_grouped_barplot(
  data,
  feats,
  group = NULL,
  group_lv = NULL,
  control_label = NULL,
  mainbar = c("mean", "median", "n", "n_missing", "sd", "var", "se", "cv", "mad",
    "skewness", "excess_kurtosis"),
  errorbar = c("none", "se", "sd", "ci"),
  conf_level = 0.95,
  gap = 1,
  lwd = 1.5,
  col = NULL,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  ylim = NULL,
  dark = FALSE,
  grid_lty = 1,
  grid_lwd = 0.25,
  cex.lab = 1.3,
  cex.axis = 1.2,
  cex.main = 1.3,
  cex.legend = 1.1,
  out_statistics = TRUE
)

Arguments

Details

This is the summary counterpart of draw_grouped_boxplot(). A box shows the distribution a group's observations have; a bar shows one number standing for them, which is less of the data and is what a figure wants when the point being made is about a location rather than a spread.

A bar is read against the zero it stands on, so a derived ylim always includes zero and puts its headroom on the side the bars run to. Heights that go both ways, which "skewness" and a "mean" of centred features do, get the baseline drawn as well.

A bar whose height is NA leaves a blank rather than shifting the ones beside it, which is what makes a group too small for a shape estimate visible instead of silently absent.

The function changes graphical parameters and the panel layout, and restores both on exit, so the caller's device is left as it was found.

Value

If out_statistics = FALSE, NULL invisibly. Otherwise, invisibly, a data.frame of one row per bar in the order they were drawn, which is the row order of summarize_descriptive_stats(): a feature's levels stay together, in group_lv order.

features,group

Which bar the row is.

n

Finite observations the bar was computed from.

value

The bar height, the mainbar column.

lower,upper

The ends of the interval, NA under errorbar = "none" and for a bar whose interval was not defined.

Which height and which interval were drawn are attached as the attributes "mainbar" and "errorbar".

What a bar may carry an interval for

mainbar names any of the summary columns, but only two of them are locations that an interval either side says something about, so errorbar is read under mainbar rather than independently of it.

"mean"

Takes every bar. "se" and "sd" are one standard error and one standard deviation either side, and "ci" is Student's interval at conf_level: stats::qt() on n - 1 degrees of freedom times the standard error.

"median"

Takes "ci" only, the notch interval ⁠median +/- 1.58 * IQR / sqrt(n)⁠. That is the same interval draw_grouped_boxplot() returns as median_confidence_stats, so a bar and the notch of the box beside it are the same width on the same data. "se" and "sd" describe the observations' spread about their mean and are refused here rather than drawn around a median they are not about.

Everything else

Takes "none". A height that is itself a spread, a count or a shape has no second quantity for an interval to be about.

A refused combination is an error rather than a silently dropped interval, because the alternative is a figure that answers a question other than the one it was asked.

References

McGill, R., Tukey, J. W. and Larsen, W. A. (1978). Variations of box plots. The American Statistician, 32(1), 12-16.

See Also

summarize_descriptive_stats() for the table the heights come from, draw_grouped_boxplot() for the observations behind them, and compare_two_groups() or compare_multiple_groups() to test the difference a pair of bars shows.

Examples

feats <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")

## Mean per species, one standard error either side
bars <- draw_grouped_barplot(iris, feats, iris$Species, errorbar = "se")
head(bars)

## The heights are the table's own column, so neither has to be checked
## against the other by eye.
summ <- summarize_descriptive_stats(iris, feats, iris$Species)
all.equal(bars$value, summ$mean)

## Student's interval, at a level of your own
draw_grouped_barplot(iris, feats, iris$Species, errorbar = "ci",
                     conf_level = 0.99)

## A median takes the notch interval draw_grouped_boxplot() notches with
draw_grouped_barplot(iris, feats, iris$Species, mainbar = "median",
                     errorbar = "ci")

## A count carries no interval, and `group_lv` picks the levels and their
## order
draw_grouped_barplot(iris, feats, iris$Species,
                     group_lv = c("virginica", "versicolor"),
                     mainbar = "n")

## `control_label` draws the named level first without rewriting `group_lv`
draw_grouped_barplot(iris, feats, iris$Species,
                     control_label = "setosa")

Description

Draws clusters of boxes with the levels of one factor side by side inside each cluster and a legend in a narrow panel on the right. A single factor gives one cluster per feature. A crossed design gives one panel per feature, with the factors after the first along the x axis, so that the crossing sits inside a panel where an interaction can be seen. Optionally returns the summary statistics behind the boxes.

Usage

draw_grouped_boxplot(
  data,
  feats,
  group = NULL,
  group_lv = NULL,
  factors = NULL,
  factor_lv = NULL,
  control_label = NULL,
  panel_by = c("feature", "factor"),
  panel_nrow = NULL,
  gap = 1,
  lwd = 1.5,
  xlab = NULL,
  ylab = NULL,
  cex.lab = 1.3,
  cex.axis = 1.2,
  cex.main = 1.3,
  ylim = NULL,
  main = NULL,
  dark = FALSE,
  grid_lty = 1,
  grid_lwd = 0.25,
  cex.legend = 1.1,
  out_statistics = TRUE
)

Arguments

Details

The function changes graphical parameters and the panel layout, and restores both on exit, so the caller's device is left as it was found.

Value

If out_statistics = FALSE, NULL invisibly. Otherwise a list of two elements, invisibly:

box_summary_stats

One data.frame per feature with rows min, lower_bound, Q1, median, Q3, upper_bound, max and one column per box: the group levels under group, and the cell labels of the design under factors. The bounds are the Tukey whisker fences Q1 - 1.5 * IQR and Q3 + 1.5 * IQR, not the drawn whisker ends.

median_confidence_stats

One data.frame per feature with rows n, lower_conf, upper_conf, the notch interval ⁠median +/- 1.58 * IQR / sqrt(n)⁠. n counts non-missing values.

What is drawn where

A crossed design has three categorical axes to place and two dimensions to place them in, so one of them has to become the panels. Which one decides whether the picture shows an interaction.

panel_by = "feature", the default, gives one panel per feature. Inside it the remaining factors run along the x axis and the primary factor is the coloured boxes within each of those clusters, so the two factors are side by side in one panel: an effect of the treatment that reverses between the levels of the other factor is a pattern of colours that visibly flips a couple of centimetres away, which is what an interaction is.

panel_by = "factor" gives the transpose: one panel per combination of the remaining factors, with the features along the x axis. It puts every feature of one cell together, which is what to ask for when the question is about the features rather than about the crossing. The cost is that the two factors are then split between the legend and the panels, and reading an interaction means comparing a colour profile in one panel against the same in another.

Either way the boxes are the same boxes and the returned statistics are identical; only the grouping of them into panels and clusters differs. The cells, their order and their labels come from the same helpers compare_factorial_groups() uses, so a box of this plot and a row of that result are the same observations, and the columns of the returned statistics are the cell labels the comparison and simulate_factorial_groups() both key on. A cell holding no observation leaves its box blank rather than shifting the ones beside it, and is reported in a message.

References

McGill, R., Tukey, J. W. and Larsen, W. A. (1978). Variations of box plots. The American Statistician, 32(1), 12-16.

See Also

compare_two_groups() and compare_factorial_groups() to test the same input.

Examples

## All three iris species across the four measurements
stats <- draw_grouped_boxplot(
  data     = iris,
  feats    = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"),
  group    = iris$Species,
  group_lv = c("setosa", "versicolor", "virginica"),
  ylab     = "cm",
  dark     = FALSE
)
stats$box_summary_stats$Sepal.Length
stats$median_confidence_stats$Sepal.Length

## A crossed design: one panel per feature, the sexes along the x axis and
## the treatment the colours, so an effect that flips with sex is local.
sim <- simulate_factorial_groups(n_feats = 6, n_per_cell = 8, seed = 1)
cells <- draw_grouped_boxplot(
  data      = sim$args$data,
  feats     = sim$args$feats,
  factors   = sim$args$factors,
  factor_lv = sim$args$factor_lv,
  ylab      = "log2 abundance"
)
## One column per cell, labelled the way the answer key labels them.
cells$box_summary_stats$prot_1

## The transpose: one panel per sex, every feature together inside it.
draw_grouped_boxplot(
  data      = sim$args$data,
  feats     = sim$args$feats,
  factors   = sim$args$factors,
  factor_lv = sim$args$factor_lv,
  panel_by  = "factor",
  ylab      = "log2 abundance"
)

## `control_label` draws the named level first in its own factor, and takes
## the same argument the comparison does, so the two agree on where the
## reference cell sits.
draw_grouped_boxplot(
  data          = sim$args$data,
  feats         = sim$args$feats[1:2],
  factors       = sim$args$factors,
  control_label = list(sex = "female"),
  ylab          = "log2 abundance"
)

Description

Draws one cell per feature and sample, with the sample groups as a coloured annotation strip above the columns and a dendrogram on each axis that was clustered. The input is the wide format the comparison functions take, one row per observation and one column per feature, and it is transposed so that features run down the rows the way an expression heatmap is usually read. stats::heatmap() draws the cells, the strip and the trees; this function decides what goes into them and adds the colour key beside them.

Usage

draw_heatmap(
  data,
  group = NULL,
  group_lv = NULL,
  feats = NULL,
  scale = c("feature", "sample", "none"),
  zlim = NULL,
  dist_method = c("euclidean", "correlation", "manhattan"),
  hclust_method = c("average", "complete", "ward.D2"),
  cluster_feats = TRUE,
  cluster_samples = TRUE,
  feat_labels = NULL,
  sample_labels = NULL,
  show_feat_names = TRUE,
  show_sample_names = TRUE,
  anno = FALSE,
  cex.anno = 1,
  n_colors = 101,
  main = NULL,
  cex.axis = 0.9,
  cex.main = 1.5,
  cex.legend = 1.2
)

Arguments

Details

Features are z-scored across the samples by default. Without it a single high-abundance feature takes the whole colour range and everything else is left white, since the colour scale is shared by every cell in the plot and features are not measured on a common scale. The plot does not say which scale ran: naming it would take a line of text wider than the key it titles, and the numbers beside the key are the scaled values either way. A feature with no variance would divide by zero, so it is only centred and ends up flat at the middle of the scale; how many were is reported in a message.

A diverging palette needs a meaningful midpoint. When zlim is not given, zero is that midpoint if the values being drawn have both signs, which is always the case after z-scoring, and the range is made symmetric around it. Values of one sign only have no such point, so the range is the range of the values and white falls in the middle of it.

Missing values are drawn as grey cells rather than dropped, so a gap is visible as a gap. Clustering leaves them to stats::dist(), which measures a pair on the values they share. A pair that shares none has no distance at all, and rather than fail, that axis is left in its input order with a message saying so.

The clustering is done here rather than left to stats::heatmap(), which would standardise after clustering and divide a flat feature by zero, and the trees are handed over already built. What comes back is therefore the tree the plot shows. The panel proportions are stats::heatmap()'s, including the square cells its respect = TRUE layout asks for. That layout keeps the block of panels square and centres it, so a device much wider than it is tall is left with empty space to the right of the plot and its key.

stats::heatmap() draws no colour key, so this function keeps a strip beside the plot for one and overlays it after: an upright bar with its numbers to the right of it, and the group legend to the right of those, both starting level with the top row of cells. The strip is as wide as those need rather than a fixed fraction of the device, and it is placed against the feature labels rather than at the edge, since the space the centred layout leaves over would otherwise read as a gap between the plot and its key. The panel layout, the margins and the plot region are all restored on exit, so the caller's device is left as it was found.

Value

A list, invisibly:

matrix

The scaled matrix as it was drawn: features in rows, samples in columns, both in the order they appear on the plot, and the row and column names as they were labelled. Values are not clamped to zlim, which is a decision about colour rather than about the data.

feat_order,sample_order

The permutations the clustering chose, as indices into feats and into the rows of data that were kept.

feat_hclust,sample_hclust

The stats::hclust() objects behind the dendrograms, or NULL for an axis that was not clustered.

zlim

The colour range, derived or as supplied.

group_colors

The colour drawn for each group level, or NULL when no group was supplied.

See Also

draw_grouped_boxplot(), which takes the same wide input.

Examples

## Features are z-scored, so the two planted directions separate by colour
## and the clustering of the columns finds the two groups on its own.
sim <- simulate_two_groups(n_feats = 30, n_up = 5, n_down = 5, seed = 3)
out <- draw_heatmap(
  data     = as.matrix(sim$args$data),
  group    = sim$args$group,
  group_lv = sim$args$group_lv,
  show_sample_names = FALSE,
  main     = "30 features, 10 of them planted"
)

## The clustering is on the result, so what the picture shows can be checked
## rather than eyeballed.
head(rownames(out$matrix))
out$sample_hclust$method

## Only the planted features, in the order they were named
draw_heatmap(
  data          = sim$args$data,
  group         = sim$args$group,
  group_lv      = sim$args$group_lv,
  feats         = sim$truth$features[sim$truth$direction != "none"],
  cluster_feats = FALSE,
  show_sample_names = FALSE
)

## Unscaled values, with the four measurements on their own scale
draw_heatmap(
  data     = iris[1:4],
  group    = iris$Species,
  group_lv = levels(iris$Species),
  scale    = "none",
  show_sample_names = FALSE,
  main     = "iris, cm"
)

## No group strip or group legend when both are omitted
draw_heatmap(
  data = as.matrix(sim$args$data),
  show_sample_names = FALSE
)

Description

Joins the cell means of one factor across the levels of another, one line per level of the tracing factor, which is the picture an interaction term is a test of: parallel lines are additive factors and crossing or fanning lines are the interaction the term reports a p-value for.

Usage

draw_interaction_plot(
  comparison_result,
  x = NULL,
  trace = NULL,
  facet = NULL,
  type = c("auto", "pairwise", "matrix", "facet"),
  feats = NULL,
  errorbar = c("none", "se", "ci"),
  panel_nrow = NULL,
  dark = FALSE,
  ylim = NULL,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  col = NULL,
  lwd = 2,
  cex.axis = 1.2,
  cex.lab = 1.3,
  cex.main = 1.3,
  cex.legend = 1.1
)

Arguments

Details

The means come from comparison_result$cells, so the plot and the F tests beside it describe the same fit rather than two passes over the data. Where the design holds more factors than the two being drawn, the remaining ones are averaged away, unweighted, which is the same marginal mean the post-hoc stage contrasts: the gap between two points of one line is the estimate its contrast reports in comparison_result$posthoc.

Three views are available:

"pairwise"

One pair of factors, trace against x, with every other factor averaged away. One panel per feature, so several features can be read side by side.

"matrix"

Every pair of factors at once, in an upper triangle: the row is the tracing factor and the column the one on the x axis. One feature only, since the panels are spent on the factors.

"facet"

trace against x again, but with the levels of facet kept apart in panels of their own rather than averaged away, which is what shows a two-factor interaction that itself depends on a third factor. One feature only, for the same reason.

type = "auto", the default, reads the arguments: naming facet asks for the facet view, naming x or trace for the pairwise one, and naming none of them gives the pairwise view for two factors and the matrix for three or more.

The function changes graphical parameters and the panel layout, and restores both on exit, so the caller's device is left as it was found.

Each level of the tracing factor gets a colour, a plotting symbol and a line type of its own, so the lines stay apart in greyscale and for a reader who cannot tell the two colours apart.

Value

The plotted means, invisibly, one row per point in the order they were drawn, with columns panel, features, x_factor, x_level, trace_factor, trace_level, n_cells, mean, se, lower_conf and upper_conf. n_cells is how many cells were averaged into the point, so 1 marks a cell mean and anything more a marginal one. The view that type = "auto" resolved to is attached as the attribute "view".

What an error bar here is not

The bars are built from cells$se, which is sqrt(ms_error / n) pooled over the whole model, so a bar on a marginal mean over a set of cells has half-width sqrt(sum(se^2) / length(se)^2) times one or the qt() multiplier. That is the standard error of the mean being drawn, not of the difference between two of them, and the difference is what the interaction is about. Two bars that overlap do not settle the term test, which is in ⁠$terms⁠, and two points whose gap you want tested are in ⁠$posthoc⁠.

See Also

compare_factorial_groups() for the result this reads, draw_grouped_boxplot() for the observations behind the means, and draw_volcano_plot() with estimate_significance(by = "term") for which features have an interaction at all.

Examples

sim <- simulate_factorial_groups(n_feats = 4, n_per_cell = 8, seed = 1)
res <- do.call(compare_factorial_groups, c(sim$args, list(diagnose = FALSE)))

## One panel per feature, the first factor tracing the second
draw_interaction_plot(res)

## The same two factors the other way round, with standard error bars
draw_interaction_plot(res, x = "treatment", trace = "sex", errorbar = "se")

## Two features only, stacked
draw_interaction_plot(res, feats = res$features[1:2], panel_nrow = 2)

## Every pair of factors of a three-factor design, for one feature
sim3 <- simulate_factorial_groups(
  n_feats = 2, n_per_cell = 8, seed = 2,
  factor_lv = list(treatment = c("control", "treat_A"),
                   sex       = c("male", "female"),
                   site      = c("north", "south", "east"))
)
res3 <- do.call(compare_factorial_groups,
                c(sim3$args, list(diagnose = FALSE)))
draw_interaction_plot(res3, type = "matrix")

## The third factor kept apart in panels instead of averaged away
draw_interaction_plot(res3, x = "site", trace = "treatment", facet = "sex",
                      errorbar = "ci")

Description

Splits the x axis by the first variable's marginal shares and each strip by the second variable's conditional shares, so the area of a tile is the cell's share of the table. Three things are drawn on top of that geometry, and each of them answers a question the geometry alone leaves open.

Usage

draw_mosaic_plot(
  categorical_comparison_result,
  shade = TRUE,
  residual = c("pearson", "standardized"),
  expected_line = TRUE,
  anno_cells = c("auto", "count", "percent", "both", "none"),
  gap = 0.015,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  cex.lab = 1.3,
  cex.axis = 1.2,
  cex.main = 1.3,
  cex.legend = 1.1,
  cex.anno = 1,
  dark = FALSE
)

## S3 method for class 'sa_categorical'
plot(x, ...)

Arguments

Details

The shading says which cells made the statistic what it is. It reads the Pearson residual, the quantity that squares and sums to that statistic, at the conventional cuts of 2 and 4. The two colours are the ones draw_volcano_plot() already uses for a feature that moved up and one that moved down, so "more than expected" and "less than expected" are the same pair of colours the rest of the package reads as a direction.

The expected line says what "expected" was. A dashed segment sits at each boundary the tiles of a strip would have had under the null hypothesis, so the departure is the distance between a tile edge and the line beside it rather than something to be inferred by comparing strips by eye. Under independence the lines fall at the same heights in every strip, which is what makes an association visible at a glance; under symmetry they do not, because there the expectation is a cell against its own transpose.

The annotation says how many observations a tile stands for, which area cannot: a wide short tile and a narrow tall one can hold the same count.

The level names on the y axis are read off the reference strip, the first one, which the control_label and category_lv of compare_categorical_groups() decide. No single set of positions can label every strip, since the whole content of a mosaic is that the strips are cut at different heights, so labelling one of them and saying which is the honest version of the choice. The first strip is the one the rest of the package already treats as the reference.

Only the par() values this function sets are put back. A blanket par(no.readonly = TRUE) snapshot also carries fin, pin and mai, which are absolute sizes, and restoring them pins the next plot to the size this one happened to be drawn at.

Value

Invisibly, a list of the picture as it was drawn.

cells

The cell table with x1, x2, y1, y2 and fill added, in the order the tiles were painted.

widths

Marginal share of each strip, named by row level. These are the widths before the gaps are taken out of them.

heights

Conditional share of each tile within its strip, as a matrix of row level by column level.

expected_prop

The same shares the null hypothesis expects, which is what the dashed segments were drawn from.

empty_levels

row and col levels that hold no observation, so drew no tile and took no axis label.

null

Which hypothesis the shading and the segments are about.

residual,residual_breaks,colors

The scale the shading read.