# test_normality() =========================================================
# TODO:
# 1. Use "parse_formula" for formula interface.
# 2. Enable formaula of form: y1 + y2 + y3 ~ group
# 3. Enable formaula of form: ~ y1 + y2 + y3 | group
#
#
# DONE:
# 1. Order of arguments changed. Now first is `data`
# 2. [+] Prevent from failing if one of the groups do not meet function requirements.
#
# An example of failing:
#
# DATA_1 <- structure(list(
# A = structure(c(2L, 1L, 1L, 1L, 1L, 1L, 1L,
# 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
# 1L, 2L, 2L, NA, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L,
# 2L, 2L, 2L, 2L, 2L, 1L, NA, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L,
# 1L, 2L, 2L, NA, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
# 1L, 1L, 2L, 2L, 1L, 2L, 2L, NA, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
# 2L, 2L, 1L), .Label = c("(-)", "(+)"), class = "factor"),
# B = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 1L,
# 3L, 1L, 1L, 2L, 3L, 3L, 1L, 2L, 2L, 3L, 3L, 3L, 2L, 1L, 2L,
# 3L, 3L, 3L, 1L, 2L, 1L, 3L, 2L, 2L, 3L, 2L, 2L, 1L, 1L, 1L,
# 3L, 2L, 1L, 1L, 1L, 3L, NA, 2L, 2L, 3L, 2L, 3L, 1L, 2L, 1L,
# 1L, 1L, 2L, 2L, NA, 2L, 1L, 1L, 1L, 1L, 2L, NA, 2L, 1L, 2L,
# 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L,
# 1L, 2L, 1L, 1L, 1L, NA), .Label = c("A", "B",
# "C"), class = c("ordered", "factor")),
# Y = c(1.41011635102669,
# 2.67996596035883, 2.93082281827539, 3.47057470748745, 3.93013356916371,
# 4.46463889802315, 4.16879679954147, 0.877400732688826, 2.18144836901425,
# 0.43029448325718, 4.62383351283921, 3.81767823167015, 1.12524434294819,
# 0.668998063328768, 0.632398055004286, 5.67433575537348, 2.61454048710063,
# 3.02047582397111, 0.836036933062681, 1.02234637748381, 0.612431298292482,
# 2.41600148898046, 6.25820428291974, 2.72896425226281, 1.04713994687296,
# 1.02055325724726, 0.871851076216876, 2.29227854194425, 2.53371799679094,
# 2.54263453547687, 0.915853083472324, 2.72562238958133, 1.26251201645741,
# 0.639371807195731, 2.21473393518887, 1.05972500710731, 3.0949115598333,
# 3.45374832012055, 4.46001217389848, 1.84533320388126, 2.50946203663001,
# 2.94273769714618, 5.56578156481604, 2.58221811243787, 0.903003765776914,
# 0.61814786719151, 2.01169911590801, 2.48549359503437, 0.707356744235952,
# 1.30002021889422, 0.567361588379272, 5.00332592197709, 5.1262908784939,
# 4.11646207529042, 5.9301390273134, 5.05449658118008, 1.3153282491963,
# 1.62271144298906, 2.99666457392001, 0.975593543621166, 2.4374333607522,
# 2.72788236111502, 5.03425628052577, 3.8940899952299, 1.80900903415956,
# 0.815976697401983, 1.58629992809868, 2.06685279835032, 2.33730373263884,
# 3.51154948522958, 5.5969942369269, 2.05741437514989, 2.22454410505537,
# 1.57830558076523, 3.03264625533027, 7.87293092517559, 3.96414306364804,
# 2.50924620116726, 1.40667280365032, 1.01336405963084, 4.00607224407719,
# 4.1102592110769, 3.71846260372488, 5.73830971585118, 3.86107565341245,
# 3.20207416723376, 4.1743749518175, 2.9932090278708, 4.68758885818583,
# 5.08063093725737)), class = c("tbl_df", "tbl", "data.frame"
# ), row.names = c(NA, -90L))
#
#
# biostat::test_normality(Y ~ A + B, data = DATA_1)
# # debugonce(biostat:::test_)
# ============================================================================
#' Normality tests by groups
#'
#' Perform a test of normality for each group separately.
#' The available tests: Shapiro-Wilk (default),
#' Lilliefors (Kolmogorov-Smirnov), Anderson-Darling and other.
#'
#' @param y (formula|numeric|character)\cr
#' Either a formula, a numeric vector or a name of a column
#' in `data`. \itemize{
#' \item If `y` is a formula (e.g. `variable ~ factor`), left-hand
#' side provides the name of a variable to be tested. In the right-hand side
#' there are names of factor variables to be used to create subsets.
#' If the left-hand side is empty (e.g. `~ factor`), right-hand
#' side is treated as variable name to test.
#' }
#'
#' @param groups (`NULL`|factor|character) \cr
#' An alternative way to provide groups.
#' If `y` is a numeric vector, `groups` must be
#' a factor (or `NULL`). If `y` is a sting,
#' `groups` must also be a string (or `NULL`).
#'
#' @param data (data frame|`NULL`) \cr
#' Either a data frame that contains the
#' variables mentioned in `y` or `NULL` (if the variables
#' are in the function's environment).
#'
#' @param test (string | function) \cr
#' Either a function that carries out a normality test or
#' a string (case insensitive, maybe unambiguously abbreviated)
#' with a name of a normality test. Possible names of tests:
#' \itemize{
#' \item{"SW", "Shapiro-Wilk" — for Shapiro-Wilk test;}
#' \item{"Lilliefors" — for Kolmogorov-Smirnov test with Lilliefor's correction;}
#' \item{"AD", "Anderson-Darling" — for Anderson-Darling test;}
#' \item{"CVM", "CM", "Cramer-von Mises" — for Cramer-von Mises test;}
#' \item{"SF", "Shapiro-Francia" — for Shapiro-Francia test;}
#' \item{"Chi-squared","Pearsons" — for Pearson's chi-squared test of normality.}
#' }
#'
#' @param signif_stars (logical) \cr
#' If `TRUE`, significance stars are printed.
#'
#' @param legend (logical) \cr
#' If `TRUE`, legend for significance stars
#' is printed.
#'
#' @param digits_stat (integer) \cr
#' Either a number of decimal places or number of
#' significant digits to round test statistic to.
#'
#' @param format_stat (character) \cr
#' Number format "f", "g" or "auto" (default) for test
#' statistic. More about number
#' formats "f" and "g" you can find in documentation of
#' function [base::formatC()] section `format`.
#'
#' @param show_col_method (logical) \cr
#' If `FALSE` column "method" is not
#' printed. Default is `FALSE`.
#' @param hide_error_msg (logical) \cr
#' If `FALSE` column "error_msg" is not printed.
#' If `TRUE` the column is printed if error occurs
#' in calcultations.
#' Default is `FALSE`.
#'
#' @param caption (string|`NULL`|`NA`) \cr
#' A caption for the table with
#' results of a normality test. If `NA` — a default caption
#' is printed (default). If `NULL` – no caption is printed.
#'
#' @param p_adjust_method (`NULL`|string) \cr
#' A name of p value adjustment
#' method for multiple comparisons. For available methods check
#' [stats::p.adjust.methods()].
#' If `NULL`, no adjusted p value is calculated (default).
#' If string (e.g., `"holm"`), an additional column
#' `p.adjust` is calculated.
#'
#' @param ... Further parameters to the function of normatity test.
#'
#' @param x `normality_test` object.
#'
#' @inheritParams format_p_values
#' @inheritParams stats::shapiro.test
#' @inheritParams nortest::lillie.test
#' @inheritParams nortest::pearson.test
#'
#'
#' @return \itemize{
#' \item Function `test_normality` returns a data frame with
#' normality test results for each group.
#' \item `print` and `pander` methods format and print the
#' results. By default, methods `print.test_normality` and
#' `pander.test_normality` do not print column called "method".
#' }
#'
#'
#' @export
#'
#' @importFrom stats shapiro.test
#' @import nortest
#'
#' @seealso \itemize{
#' \item [stats::shapiro.test()] in package \pkg{stats};
#' \item [nortest::lillie.test()] in package \pkg{nortest};
#' \item [nortest::pearson.test()] in package \pkg{nortest};
#' \item [nortest::ad.test()] in package \pkg{nortest};
#' \item [nortest::cvm.test()] in package \pkg{nortest};
#' \item [nortest::sf.test()] in package \pkg{nortest}.
#' }
#'
#' @examples
#' library(biostat)
#' library(pander)
#' data(CO2)
#'
#' # For whole dataset
#' test_normality(~uptake, data = CO2)
#'
#' # For each subgroup
#' test_normality(uptake ~ Treatment, data = CO2)
#'
#' # For each subgroup by several factor variables
#' rez <- test_normality(uptake ~ Type + Treatment,
#' data = CO2)
#' rez
#'
#' # Modify printed output
#' print(rez, rm_zero = TRUE)
#'
#' print(rez, digits_p = 2)
#'
#'
#' # Choose another test of normality
#' rez2 <- test_normality(uptake ~ Type + Treatment,
#' data = CO2,
#' test = "chi-sq")
#' rez2
#'
#'
#' # Print as a 'pandoc' table (for R Markdown reports)
#' pander(rez)
#'
#' pander(rez, digits_stat = 2)
#'
#' pander(rez, digits_p = 2)
#'
#' pander(rez, digits_p = 4, signif_stars = FALSE)
#'
#'
#' \donttest{\dontrun{
#' # View unformatted results in a separate window
#' View(rez)
#' }}
#'
#' # Show object's class
#' class(rez)
#'
test_normality <- function(y,
data = NULL,
test = "Shapiro-Wilk",
p_adjust_method = NULL,
...,
groups = NULL,
ss = signif_syms,
hide_error_msg = FALSE) {
# Choose the test ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (is.function(test)) {
use_test <- test
} else if (checkmate::test_character(test, len = 1)) {
# Possible choices
available_tests <- c(
"SW", "Shapiro-Wilk",
"Lilliefors",
"AD", "Anderson-Darling",
"CVM", "CM", "Cramer-von Mises",
"SF", "Shapiro-Francia",
"Chi-squared", "Pearsons"
)
test <- match.arg(tolower(test), tolower(available_tests))
use_test <- switch(
test,
"sw" = ,
"shapiro.test" = ,
"shapiro-wilk" = stats::shapiro.test,
"lillie.test" = ,
"lilliefors" = nortest::lillie.test,
"ad" = ,
"ad.test" = ,
"anderson-darling" = nortest::ad.test,
"cvm" = ,
"cm" = ,
"cvm.test" = ,
"cramer-von mises" = nortest::cvm.test,
"sf" = ,
"sf.test" = ,
"shapiro-francia" = nortest::sf.test,
"chi-squared" = ,
"pearson.test" = ,
"pearsons" = nortest::pearson.test
)
} else {
stop(
"\n`test` must be either a function",
"\n or a name of a test (a string of length 1)."
)
}
use_test_mod <- function(...) {
safe_rez <- purrr::safely(use_test, otherwise = NULL)(...)
if (is.null(safe_rez$error)) {
out <- safe_rez$result
out$error_msg <- NA_character_
} else {
out <- structure(
data.frame(
method = as.character(safe_rez$error),
stringsAsFactors = FALSE
),
class = "htest"
)
}
out
}
# Output
test_(y,
data = data,
p_adjust_method = p_adjust_method,
...,
groups = groups,
test = use_test_mod,
hide_error_msg = hide_error_msg,
ss = ss
)
}
# test_() =====================================================================
test_ <- function(y,
data = NULL,
p_adjust_method = NULL,
...,
groups = NULL,
test = stats::shapiro.test,
ss = NULL,
hide_error_msg = FALSE)
# na.rm = getOption("na.rm", FALSE)
{
# Make formula according to input type
if (is.numeric(y)) {
if (!is.null(groups)) {
data <- data.frame(y = y, groups = groups)
y <- y ~ groups
} else {
data <- data.frame(y = y)
y <- ~y
}
}
if (is.null(data)) {
data <- rlang::f_env(y)
}
if (rlang::is_formula(y)) {
# Select necessary variables only
data <- stats::model.frame(y, data = data)
# To indicate, that there is no grouping the first column constant
# if (ncol(data) == 1)
# data[["Groups"]] <- "<all values>"
var_names <- names(data)
gr_vars <- rlang::syms(var_names[-1])
if (!is.numeric(data[[1]])) {
stop_glue("Variable `{var_names[1]}` must be numeric.")
}
# The main test
rez <-
data %>%
dplyr::group_by(!!!gr_vars) %>%
dplyr::do(broom::tidy(test(.[[1]], ...))) %>%
dplyr::ungroup() %>%
dplyr::select(method, dplyr::everything()) %>%
as.data.frame()
# Add error message corrently
err <- stringr::str_detect(rez$method, "([Ee]rror)|([Ww]arning)")
rez$error_msg[err] <- rez$method[err]
rez$method[err] <- unique(rez$method[!err])
# If adjusted p value is needed
if (!is.null(p_adjust_method)) {
rez$p.adjust <- p.adjust(rez$p.value, method = p_adjust_method)
}
rez <- structure(rez,
class = c("test_normality", "data.frame"),
test = as.character(unique(rez$method)),
p_adjust_method = p_adjust_method,
hide_error_msg = hide_error_msg,
ss = ss
)
return(rez)
} else {
stop("Incorrect input")
}
}
# print() ----------------------------------------------------------------
#' @rdname test_normality
#' @export
print.test_normality <- function(x,
...,
digits_p = 3,
signif_stars = !is.null(ss),
digits_stat = 3,
format_stat = c("auto", "f", "g"),
rm_zero = FALSE,
legend = TRUE,
show_col_method = FALSE,
hide_error_msg = attr(x, "hide_error_msg"),
ss = attr(x, "ss")) {
format_stat <- match.arg(format_stat)
x <- format_object(
x,
digits_p = digits_p,
digits_stat = digits_stat,
format_stat = format_stat,
signif_stars = signif_stars,
signif_as_separate_column = TRUE,
signif_stars_col_name = " ",
rm_zero = rm_zero,
show_col_method = show_col_method,
ss = ss,
hide_error_msg = isTRUE(hide_error_msg)
)
# Print the name of the method
cat("\n", "The results of", which_test(x), "\n\n")
# Print main results
NextMethod(print, x)
# Print signif. stars legend
if (legend == TRUE && signif_stars == TRUE) {
cat(paste0(
"\n",
# "Legend for p-values: \n",
signif_stars_legend(ss = ss, collapse = ", "), "\n"
))
}
# Print p adjust. method:
p_adjust_method <- attr(x, "p_adjust_method")
if (!is.null(p_adjust_method) && nrow(x) > 1) {
cat(
"The method for p-value adjustment:",
first_capital(p_adjust_method),
"\n"
)
}
}
# pander() ----------------------------------------------------------------
#' @export
#' @rdname test_normality
pander.test_normality <- function(
x,
caption = NA,
...,
digits_p = 3,
signif_stars = !is.null(ss),
digits_stat = 3,
format_stat = c("auto", "f", "g"),
rm_zero = FALSE,
legend = TRUE,
show_col_method = FALSE,
hide_error_msg = attr(x, "hide_error_msg"),
ss = attr(x, "ss")) {
format_stat <- match.arg(format_stat)
x <- format_object(
x,
digits_p = digits_p,
digits_stat = digits_stat,
format_stat = format_stat,
signif_stars = signif_stars,
signif_as_separate_column = FALSE,
rm_zero = rm_zero,
show_col_method = show_col_method,
ss = ss,
hide_error_msg = hide_error_msg
)
# String for p adjust. method:
p_adjust_method <- attr(x, "p_adjust_method")
adj_sting <-
if (!is.null(p_adjust_method) & nrow(x) > 1) {
paste0(
" The method for p-value adjustment: ",
first_capital(p_adjust_method), "."
)
} else {
""
}
# Add caption
caption <-
if (is.null(caption)) {
NULL
} else if (is.na(caption)) {
glue::glue("The results of {which_test(x)}.{adj_sting}")
} else {
caption
}
# Print table of results
NextMethod("pander", x, caption = caption, ...)
# Print legend
if (legend == TRUE && signif_stars == TRUE) {
# cat("Legend for p-values: \n`", signif_stars_legend(ss = ss), "` \n")
cat(paste0(
# "\nLegend for p-values: \n",
signif_stars_legend(ss = ss), "\n"
))
}
}
# helpers --------------------------------------------------------------------
which_test <- function(x) {
# as.character(unique(x$method))
attr(x, "test")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.