mirror of
https://github.com/msberends/AMR.git
synced 2024-12-25 19:26:13 +01:00
new g.test, extra unit tests
This commit is contained in:
parent
fc30d3fb13
commit
f63414d519
@ -1,10 +1,13 @@
|
||||
# Setting up R deps
|
||||
language: r
|
||||
r: 3.2
|
||||
r:
|
||||
# oldest support according to DESCRIPTION
|
||||
- 3.0
|
||||
# latest available release
|
||||
- release
|
||||
r_packages: covr
|
||||
cache: packages
|
||||
|
||||
# system deps, install xclip for clipboard support
|
||||
cran: https://cran.rstudio.com
|
||||
os:
|
||||
- linux
|
||||
- osx
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
Package: AMR
|
||||
Version: 0.2.0.9009
|
||||
Date: 2018-07-06
|
||||
Version: 0.2.0.9010
|
||||
Date: 2018-07-10
|
||||
Title: Antimicrobial Resistance Analysis
|
||||
Authors@R: c(
|
||||
person(
|
||||
|
||||
@ -143,6 +143,7 @@ importFrom(rvest,html_children)
|
||||
importFrom(rvest,html_node)
|
||||
importFrom(rvest,html_nodes)
|
||||
importFrom(rvest,html_table)
|
||||
importFrom(stats,complete.cases)
|
||||
importFrom(stats,fivenum)
|
||||
importFrom(stats,mad)
|
||||
importFrom(stats,na.omit)
|
||||
|
||||
17
NEWS.md
17
NEWS.md
@ -2,6 +2,13 @@
|
||||
#### New
|
||||
* Support for Addins menu in RStudio to quickly insert `%in%` or `%like%` (and give them keyboard shortcuts), or to view the datasets that come with this package
|
||||
* For convience, descriptive statistical functions `kurtosis` and `skewness` that are lacking in base R - they are generic functions and have support for vectors, data.frames and matrices
|
||||
* Function `g.test` to perform the Χ<sup>2</sup> distributed [*G*-test](https://en.wikipedia.org/wiki/G-test), which use is the same as `chisq.test`
|
||||
* Function `ratio` to transform a vector of values to a preset ratio. For example:
|
||||
```r
|
||||
ratio(c(772, 1611, 737), ratio = "1:2:1")
|
||||
# [1] 780 1560 780
|
||||
```
|
||||
* Function `p.symbol` to transform p value to their related symbol: `0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1`
|
||||
* New for frequency tables (function `freq`):
|
||||
* A vignette to explain its usage
|
||||
* Support for `table` to use as input: `freq(table(x, y))`
|
||||
@ -12,13 +19,6 @@
|
||||
* Header of frequency tables now also show Mean Absolute Deviaton (MAD) and Interquartile Range (IQR)
|
||||
* Possibility to globally set the default for the amount of items to print, with `options(max.print.freq = n)` where *n* is your preset value
|
||||
* Functions `clipboard_import` and `clipboard_export` as helper functions to quickly copy and paste from/to software like Excel and SPSS
|
||||
* Function `g.test` to perform the Χ<sup>2</sup> distributed [*G*-test](https://en.wikipedia.org/wiki/G-test)
|
||||
* Function `ratio` to transform a vector of values to a preset ratio (convenient to use with `g.test`). For example:
|
||||
```r
|
||||
ratio(c(772, 1611, 737), ratio = "1:2:1")
|
||||
# [1] 780 1560 780
|
||||
```
|
||||
* Function `p.symbol` to transform p value to their related symbol: `0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1`
|
||||
|
||||
#### Changed
|
||||
* `%like%` now supports multiple patterns
|
||||
@ -36,6 +36,9 @@ ratio(c(772, 1611, 737), ratio = "1:2:1")
|
||||
* Support for 1 or 2 columns as input for `guess_bactid`
|
||||
* Fix for printing tibbles where characters would be accidentally transformed to factors
|
||||
|
||||
#### Other
|
||||
* Unit testing for R 3.0 and the latest available release: https://travis-ci.org/msberends/AMR
|
||||
|
||||
# 0.2.0 (latest stable version)
|
||||
#### New
|
||||
* Full support for Windows, Linux and macOS
|
||||
|
||||
7
R/freq.R
7
R/freq.R
@ -141,6 +141,8 @@ frequency_tbl <- function(x,
|
||||
digits = 2,
|
||||
sep = " ") {
|
||||
|
||||
mult.columns <- 0
|
||||
|
||||
if (any(class(x) == 'data.frame')) {
|
||||
x.name <- deparse(substitute(x))
|
||||
if (x.name == ".") {
|
||||
@ -176,15 +178,14 @@ frequency_tbl <- function(x,
|
||||
# get last variable: these are frequencies
|
||||
pull(ncol(.))
|
||||
x <- rep(values, counts)
|
||||
x.name <- NULL
|
||||
x.name <- "a `table` object"
|
||||
cols <- NULL
|
||||
mult.columns <- 2
|
||||
} else {
|
||||
x.name <- NULL
|
||||
cols <- NULL
|
||||
}
|
||||
|
||||
mult.columns <- 0
|
||||
|
||||
if (!is.null(ncol(x))) {
|
||||
if (ncol(x) == 1 & any(class(x) == 'data.frame')) {
|
||||
x <- x %>% pull(1)
|
||||
|
||||
293
R/g.test.R
293
R/g.test.R
@ -16,49 +16,54 @@
|
||||
# GNU General Public License for more details. #
|
||||
# ==================================================================== #
|
||||
|
||||
#' \emph{G}-test of matrix or vector
|
||||
#' \emph{G}-test for Count Data
|
||||
#'
|
||||
#' A \emph{G}-test can be used to see whether the number of observations in each category fits a theoretical expectation (called a \strong{\emph{G}-test of goodness-of-fit}), or to see whether the proportions of one variable are different for different values of the other variable (called a \strong{\emph{G}-test of independence}).
|
||||
#' @param x numeric vector or matrix
|
||||
#' @param y expected value of \code{x}. Leave empty to determine automatically. This can also be ratios of \code{x}, e.g. calculated with \code{\link{ratio}}.
|
||||
#' @param alpha value to test the p value against
|
||||
#' @param info logical to determine whether the analysis should be printed
|
||||
#' @param minimum the test with fail if any of the observed values is below this value. Use \code{minimum = 30} for microbial epidemiology, to prevent calculating a p value when less than 30 isolates are available.
|
||||
#' \code{g.test} performs chi-squared contingency table tests and goodness-of-fit tests, just like \code{\link{chisq.test}} but is more reliable [1]. A \emph{G}-test can be used to see whether the number of observations in each category fits a theoretical expectation (called a \strong{\emph{G}-test of goodness-of-fit}), or to see whether the proportions of one variable are different for different values of the other variable (called a \strong{\emph{G}-test of independence}).
|
||||
#' @inherit stats::chisq.test params return
|
||||
#' @details If \code{x} is a matrix with one row or column, or if \code{x} is a vector and \code{y} is not given, then a \emph{goodness-of-fit test} is performed (\code{x} is treated as a one-dimensional contingency table). The entries of \code{x} must be non-negative integers. In this case, the hypothesis tested is whether the population probabilities equal those in \code{p}, or are all equal if \code{p} is not given.
|
||||
#'
|
||||
#' If \code{x} is a matrix with at least two rows and columns, it is taken as a two-dimensional contingency table: the entries of \code{x} must be non-negative integers. Otherwise, \code{x} and \code{y} must be vectors or factors of the same length; cases with missing values are removed, the objects are coerced to factors, and the contingency table is computed from these. Then Pearson's chi-squared test is performed of the null hypothesis that the joint distribution of the cell counts in a 2-dimensional contingency table is the product of the row and column marginals.
|
||||
#'
|
||||
#' If \code{simulate.p.value} is \code{FALSE}, the p-value is computed from the asymptotic chi-squared distribution of the test statistic. Otherwise the p-value is computed for a Monte Carlo test (Hope, 1968) with \code{B} replicates.
|
||||
#'
|
||||
#' In the contingency table case simulation is done by random sampling from the set of all contingency tables with given marginals, and works only if the marginals are strictly positive. Note that this is not the usual sampling situation assumed for a chi-squared test (like the \emph{G}-test) but rather that for Fisher's exact test.
|
||||
#'
|
||||
#' In the goodness-of-fit case simulation is done by random sampling from the discrete distribution specified by \code{p}, each sample being of size \code{n = sum(x)}. This simulation is done in \R and may be slow.
|
||||
#' @section \emph{G}-test of goodness-of-fit (likelihood ratio test):
|
||||
#' Use the \emph{G}-test of goodness-of-fit when you have one nominal variable with two or more values (such as male and female, or red, pink and white flowers). You compare the observed counts of numbers of observations in each category with the expected counts, which you calculate using some kind of theoretical expectation (such as a 1:1 sex ratio or a 1:2:1 ratio in a genetic cross).
|
||||
#'
|
||||
#' If the expected number of observations in any category is too small, the \emph{G}-test may give inaccurate results, and you should use an exact test instead. See the web page on small sample sizes for discussion of what "small" means.
|
||||
#' If the expected number of observations in any category is too small, the \emph{G}-test may give inaccurate results, and you should use an exact test instead (\code{\link{fisher.test}}).
|
||||
#'
|
||||
#' The \emph{G}-test of goodness-of-fit is an alternative to the chi-square test of goodness-of-fit; each of these tests has some advantages and some disadvantages, and the results of the two tests are usually very similar.
|
||||
#' The \emph{G}-test of goodness-of-fit is an alternative to the chi-square test of goodness-of-fit (\code{\link{chisq.test}}); each of these tests has some advantages and some disadvantages, and the results of the two tests are usually very similar.
|
||||
#'
|
||||
#' @section \emph{G}-test of independence:
|
||||
#' Use the \emph{G}-test of independence when you have two nominal variables, each with two or more possible values. You want to know whether the proportions for one variable are different among values of the other variable.
|
||||
#'
|
||||
#' It is also possible to do a \emph{G}-test of independence with more than two nominal variables. For example, Jackson et al. (2013) also had data for children under 3, so you could do an analysis of old vs. young, thigh vs. arm, and reaction vs. no reaction, all analyzed together.
|
||||
#'
|
||||
#' Fisher's exact test is more accurate than the \emph{G}-test of independence when the expected numbers are small, so it is recommend to only use the \emph{G}-test if your total sample size is greater than 1000.
|
||||
#' Fisher's exact test (\code{\link{fisher.test}}) is more accurate than the \emph{G}-test of independence when the expected numbers are small, so it is recommend to only use the \emph{G}-test if your total sample size is greater than 1000.
|
||||
#'
|
||||
#' The \emph{G}-test of independence is an alternative to the chi-square test of independence, and they will give approximately the same results.
|
||||
#' The \emph{G}-test of independence is an alternative to the chi-square test of independence (\code{\link{chisq.test}}), and they will give approximately the same results.
|
||||
#' @section How the test works:
|
||||
#' Unlike the exact test of goodness-of-fit, the \emph{G}-test does not directly calculate the probability of obtaining the observed results or something more extreme. Instead, like almost all statistical tests, the \emph{G}-test has an intermediate step; it uses the data to calculate a test statistic that measures how far the observed data are from the null expectation. You then use a mathematical relationship, in this case the chi-square distribution, to estimate the probability of obtaining that value of the test statistic.
|
||||
#' Unlike the exact test of goodness-of-fit (\code{\link{fisher.test}}), the \emph{G}-test does not directly calculate the probability of obtaining the observed results or something more extreme. Instead, like almost all statistical tests, the \emph{G}-test has an intermediate step; it uses the data to calculate a test statistic that measures how far the observed data are from the null expectation. You then use a mathematical relationship, in this case the chi-square distribution, to estimate the probability of obtaining that value of the test statistic.
|
||||
#'
|
||||
#' The \emph{G}-test uses the log of the ratio of two likelihoods as the test statistic, which is why it is also called a likelihood ratio test or log-likelihood ratio test. The formula to calculate a \emph{G}-statistic is:
|
||||
#'
|
||||
#' \code{G <- 2 * sum(x * log(x / x.expected))}
|
||||
#' \code{G <- 2 * sum(x * log(x / E))}
|
||||
#'
|
||||
#' Since this is chi-square distributed, the p value can be calculated with:
|
||||
#' where \code{E} are the expected values. Since this is chi-square distributed, the p value can be calculated with:
|
||||
#'
|
||||
#' \code{p <- 1 - stats::pchisq(G, df))}
|
||||
#' \code{p <- stats::pchisq(G, df, lower.tail = FALSE)}
|
||||
#'
|
||||
#' where \code{df} are the degrees of freedom: \code{max(NROW(x) - 1, 1) * max(NCOL(x) - 1, 1)}.
|
||||
#' where \code{df} are the degrees of freedom.
|
||||
#'
|
||||
#' If there are more than two categories and you want to find out which ones are significantly different from their null expectation, you can use the same method of testing each category vs. the sum of all categories, with the Bonferroni correction. You use \emph{G}-tests for each category, of course.
|
||||
#' @keywords chi
|
||||
#' @seealso \code{\link{chisq.test}}
|
||||
#' @references McDonald, J.H. 2014. \strong{Handbook of Biological Statistics (3rd ed.)}. Sparky House Publishing, Baltimore, Maryland. \url{http://www.biostathandbook.com/gtestgof.html}.
|
||||
#' @references [1] McDonald, J.H. 2014. \strong{Handbook of Biological Statistics (3rd ed.)}. Sparky House Publishing, Baltimore, Maryland. \url{http://www.biostathandbook.com/gtestgof.html}.
|
||||
#' @source This code is almost perfectly equal to \code{\link{chisq.test}}; only the calculation of the statistic was changed to \code{2 * sum(x * log(x / E))} and Yates' continuity correction was deleted.
|
||||
#' @export
|
||||
#' @importFrom stats pchisq
|
||||
#' @importFrom dplyr %>%
|
||||
#' @importFrom stats pchisq complete.cases
|
||||
#' @examples
|
||||
#' # = EXAMPLE 1 =
|
||||
#' # Shivrain et al. (2006) crossed clearfield rice (which are resistant
|
||||
@ -70,12 +75,12 @@
|
||||
#' # ratio.
|
||||
#'
|
||||
#' x <- c(772, 1611, 737)
|
||||
#' x.expected <- ratio(x, "1:2:1")
|
||||
#' x.expected
|
||||
#' E <- ratio(x, "1:2:1")
|
||||
#' E
|
||||
#' # 780 1560 780
|
||||
#'
|
||||
#' g.test(x, x.expected)
|
||||
#' # p = 0.12574.
|
||||
#' G <- g.test(x, p = c(1, 2, 1) / 4)
|
||||
#' # G$p.value = 0.12574.
|
||||
#'
|
||||
#' # There is no significant difference from a 1:2:1 ratio.
|
||||
#' # Meaning: resistance controlled by a single gene with two co-dominant
|
||||
@ -90,11 +95,7 @@
|
||||
#' # observed 1752 right-billed and 1895 left-billed crossbills.
|
||||
#'
|
||||
#' x <- c(1752, 1895)
|
||||
#' x.expected <- ratio(x, ratio = c(1, 1))
|
||||
#' x.expected
|
||||
#' # 1823.5 1823.5
|
||||
#'
|
||||
#' g.test(x, x.expected)
|
||||
#' g.test(x)
|
||||
#' # p = 0.01787343
|
||||
#'
|
||||
#' # There is a significant difference from a 1:1 ratio.
|
||||
@ -102,99 +103,127 @@
|
||||
#'
|
||||
g.test <- function(x,
|
||||
y = NULL,
|
||||
alpha = 0.05,
|
||||
info = TRUE,
|
||||
minimum = 0) {
|
||||
|
||||
if (sum(x) < 1000) {
|
||||
warning('the sum of all observations is < 1000, consider using a Fishers Exact test instead.')
|
||||
# correct = TRUE,
|
||||
p = rep(1/length(x), length(x)),
|
||||
rescale.p = FALSE,
|
||||
simulate.p.value = FALSE,
|
||||
B = 2000) {
|
||||
DNAME <- deparse(substitute(x))
|
||||
if (is.data.frame(x))
|
||||
x <- as.matrix(x)
|
||||
if (is.matrix(x)) {
|
||||
if (min(dim(x)) == 1L)
|
||||
x <- as.vector(x)
|
||||
}
|
||||
|
||||
if (!is.numeric(x)) {
|
||||
stop('`x` must be a vector or matrix with numeric values.')
|
||||
if (!is.matrix(x) && !is.null(y)) {
|
||||
if (length(x) != length(y))
|
||||
stop("'x' and 'y' must have the same length")
|
||||
DNAME2 <- deparse(substitute(y))
|
||||
xname <- if (length(DNAME) > 1L || nchar(DNAME, "w") >
|
||||
30)
|
||||
""
|
||||
else DNAME
|
||||
yname <- if (length(DNAME2) > 1L || nchar(DNAME2, "w") >
|
||||
30)
|
||||
""
|
||||
else DNAME2
|
||||
OK <- complete.cases(x, y)
|
||||
x <- factor(x[OK])
|
||||
y <- factor(y[OK])
|
||||
if ((nlevels(x) < 2L) || (nlevels(y) < 2L))
|
||||
stop("'x' and 'y' must have at least 2 levels")
|
||||
x <- table(x, y)
|
||||
names(dimnames(x)) <- c(xname, yname)
|
||||
DNAME <- paste(paste(DNAME, collapse = "\n"), "and",
|
||||
paste(DNAME2, collapse = "\n"))
|
||||
}
|
||||
if (!is.matrix(x) & is.null(y)) {
|
||||
stop('if `x` is not a matrix, `y` must be given.')
|
||||
if (any(x < 0) || anyNA(x))
|
||||
stop("all entries of 'x' must be nonnegative and finite")
|
||||
if ((n <- sum(x)) == 0)
|
||||
stop("at least one entry of 'x' must be positive")
|
||||
if (simulate.p.value) {
|
||||
setMETH <- function() METHOD <<- paste(METHOD, "with simulated p-value\n\t (based on",
|
||||
B, "replicates)")
|
||||
almost.1 <- 1 - 64 * .Machine$double.eps
|
||||
}
|
||||
|
||||
# if (!is.matrix(x)) {
|
||||
# x <- matrix(x, dimnames = list(rep("", length(x)), ""))
|
||||
if (is.matrix(x)) {
|
||||
METHOD <- "G-test of independence"
|
||||
nr <- as.integer(nrow(x))
|
||||
nc <- as.integer(ncol(x))
|
||||
if (is.na(nr) || is.na(nc) || is.na(nr * nc))
|
||||
stop("invalid nrow(x) or ncol(x)", domain = NA)
|
||||
sr <- rowSums(x)
|
||||
sc <- colSums(x)
|
||||
E <- outer(sr, sc, "*")/n
|
||||
v <- function(r, c, n) c * r * (n - r) * (n - c)/n^3
|
||||
V <- outer(sr, sc, v, n)
|
||||
dimnames(E) <- dimnames(x)
|
||||
if (simulate.p.value && all(sr > 0) && all(sc > 0)) {
|
||||
setMETH()
|
||||
tmp <- .Call(C_chisq_sim, sr, sc, B, E)
|
||||
STATISTIC <- 2 * sum(x * log(x / E)) # sum(sort((x - E)^2/E, decreasing = TRUE)) for chisq.test
|
||||
PARAMETER <- NA
|
||||
PVAL <- (1 + sum(tmp >= almost.1 * STATISTIC))/(B +
|
||||
1)
|
||||
}
|
||||
else {
|
||||
if (simulate.p.value)
|
||||
warning("cannot compute simulated p-value with zero marginals")
|
||||
# if (correct && nrow(x) == 2L && ncol(x) == 2L) {
|
||||
# YATES <- min(0.5, abs(x - E))
|
||||
# if (YATES > 0)
|
||||
# METHOD <- paste(METHOD, "with Yates' continuity correction")
|
||||
# }
|
||||
|
||||
x.expected <- y
|
||||
# if (!is.null(x.expected) & !is.matrix(x.expected)) {
|
||||
# x.expected <- matrix(x.expected, dimnames = list(rep("", length(x.expected)), ""))
|
||||
# }
|
||||
|
||||
if (NCOL(x) > 1) {
|
||||
matrix2tbl <- function(x) {
|
||||
if (!is.matrix(x)) {
|
||||
x <- matrix(x, dimnames = list(rep("", length(x)), ""))
|
||||
# else YATES <- 0
|
||||
STATISTIC <- 2 * sum(x * log(x / E)) # sum((abs(x - E) - YATES)^2/E) for chisq.test
|
||||
PARAMETER <- (nr - 1L) * (nc - 1L)
|
||||
PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
|
||||
}
|
||||
x <- rbind(cbind(x, rowSums(x)), colSums(cbind(x, rowSums(x)))) %>%
|
||||
as.data.frame()
|
||||
colnames(x) <- c(paste0('c', 1:(NCOL(x) - 1)), 'Total')
|
||||
rownames(x) <- c(strrep(" ", 1:(NROW(x) - 1)), 'Total')
|
||||
x
|
||||
}
|
||||
|
||||
x.with_totals <- matrix2tbl(x)
|
||||
if (is.null(x.expected)) {
|
||||
x.expected <- outer(rowSums(x), colSums(x), "*") / sum(x)
|
||||
x.expected.with_totals <- matrix2tbl(x.expected)
|
||||
else {
|
||||
if (length(dim(x)) > 2L)
|
||||
stop("invalid 'x'")
|
||||
if (length(x) == 1L)
|
||||
stop("'x' must at least have 2 elements")
|
||||
if (length(x) != length(p))
|
||||
stop("'x' and 'p' must have the same number of elements")
|
||||
if (any(p < 0))
|
||||
stop("probabilities must be non-negative.")
|
||||
if (abs(sum(p) - 1) > sqrt(.Machine$double.eps)) {
|
||||
if (rescale.p)
|
||||
p <- p/sum(p)
|
||||
else stop("probabilities must sum to 1.")
|
||||
}
|
||||
|
||||
} else {
|
||||
x.with_totals <- x
|
||||
if (is.null(x.expected)) {
|
||||
x.expected <- matrix(rep(mean(x), length(x)), dimnames = list(rep("", length(x)), ""))
|
||||
warning('Expected values set to the mean of `x`, because it consists of only one column.')
|
||||
METHOD <- "G-test of goodness-of-fit (likelihood ratio test)"
|
||||
E <- n * p
|
||||
V <- n * p * (1 - p)
|
||||
STATISTIC <- 2 * sum(x * log(x / E)) # sum((x - E)^2/E) for chisq.test
|
||||
names(E) <- names(x)
|
||||
if (simulate.p.value) {
|
||||
setMETH()
|
||||
nx <- length(x)
|
||||
sm <- matrix(sample.int(nx, B * n, TRUE, prob = p),
|
||||
nrow = n)
|
||||
ss <- apply(sm, 2L, function(x, E, k) {
|
||||
sum((table(factor(x, levels = 1L:k)) - E)^2/E)
|
||||
}, E = E, k = nx)
|
||||
PARAMETER <- NA
|
||||
PVAL <- (1 + sum(ss >= almost.1 * STATISTIC))/(B +
|
||||
1)
|
||||
}
|
||||
x.expected.with_totals <- x.expected
|
||||
else {
|
||||
PARAMETER <- length(x) - 1
|
||||
PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
|
||||
}
|
||||
|
||||
if (any(x < minimum)) {
|
||||
warning('One of the observed values is lower than the required minimum of ', minimum, '.')
|
||||
return(NA_real_)
|
||||
}
|
||||
if (any(x.expected < 5)) {
|
||||
warning('One of the expected values is lower than 5, thus G-statistic approximation may be incorrect.',
|
||||
'\n Consider doing an Exact test (exact.test).')
|
||||
}
|
||||
|
||||
Gstat <- 2 * base::sum(x * log(x / x.expected), na.rm = TRUE)
|
||||
df <- base::max(NROW(x) - 1, 1, na.rm = TRUE) * base::max(NCOL(x) - 1, 1, na.rm = TRUE)
|
||||
|
||||
pval <- 1 - stats::pchisq(Gstat, df = df)
|
||||
|
||||
if (info == TRUE) {
|
||||
|
||||
if (NROW(x) == 2 & NCOL(x) == 2) {
|
||||
cat('G-test of independence\n\n')
|
||||
} else {
|
||||
cat('G-test of goodness-of-fit\n')
|
||||
cat('(likelihood ratio test)\n\n')
|
||||
}
|
||||
|
||||
cat(paste0("(O) Observed values:\n"))
|
||||
print(x.with_totals %>% round(2))
|
||||
|
||||
cat('\n')
|
||||
cat('(E) Expected values under null hypothesis:\n')
|
||||
print(x.expected.with_totals %>% round(2))
|
||||
|
||||
cat('\n============[G-test]============\n')
|
||||
cat(' G-statistic :', round(Gstat, 4), '\n')
|
||||
cat(' Degrees of freedom :', df, '\n')
|
||||
cat(' P value :', round(pval, 4), '\n')
|
||||
cat(' Alpha :', round(alpha, 2), '\n')
|
||||
cat(' Significance :', p.symbol(pval, "-"), '\n')
|
||||
cat('================================\n\n')
|
||||
|
||||
}
|
||||
|
||||
pval
|
||||
|
||||
names(STATISTIC) <- "X-squared"
|
||||
names(PARAMETER) <- "df"
|
||||
if (any(E < 5) && is.finite(PARAMETER))
|
||||
warning("G-statistic approximation may be incorrect")
|
||||
structure(list(statistic = STATISTIC, parameter = PARAMETER,
|
||||
p.value = PVAL, method = METHOD, data.name = DNAME,
|
||||
observed = x, expected = E, residuals = (x - E)/sqrt(E),
|
||||
stdres = (x - E)/sqrt(V)), class = "htest")
|
||||
}
|
||||
|
||||
#' Transform vector to ratio
|
||||
@ -204,47 +233,7 @@ g.test <- function(x,
|
||||
#' @seealso \code{\link{g.test}}
|
||||
#' @references McDonald, J.H. 2014. \strong{Handbook of Biological Statistics (3rd ed.)}. Sparky House Publishing, Baltimore, Maryland.
|
||||
#' @importFrom dplyr %>%
|
||||
#' @examples
|
||||
#' # = EXAMPLE 1 =
|
||||
#' # Shivrain et al. (2006) crossed clearfield rice (which are resistant
|
||||
#' # to the herbicide imazethapyr) with red rice (which are susceptible to
|
||||
#' # imazethapyr). They then crossed the hybrid offspring and examined the
|
||||
#' # F2 generation, where they found 772 resistant plants, 1611 moderately
|
||||
#' # resistant plants, and 737 susceptible plants. If resistance is controlled
|
||||
#' # by a single gene with two co-dominant alleles, you would expect a 1:2:1
|
||||
#' # ratio.
|
||||
#'
|
||||
#' x <- c(772, 1611, 737)
|
||||
#' x.expected <- ratio(x, "1:2:1")
|
||||
#' x.expected
|
||||
#' # 780 1560 780
|
||||
#'
|
||||
#' g.test(x, x.expected)
|
||||
#' # p = 0.12574.
|
||||
#'
|
||||
#' # There is no significant difference from a 1:2:1 ratio.
|
||||
#' # Meaning: resistance controlled by a single gene with two co-dominant
|
||||
#' # alleles, is plausible.
|
||||
#'
|
||||
#'
|
||||
#' # = EXAMPLE 2 =
|
||||
#' # Red crossbills (Loxia curvirostra) have the tip of the upper bill either
|
||||
#' # right or left of the lower bill, which helps them extract seeds from pine
|
||||
#' # cones. Some have hypothesized that frequency-dependent selection would
|
||||
#' # keep the number of right and left-billed birds at a 1:1 ratio. Groth (1992)
|
||||
#' # observed 1752 right-billed and 1895 left-billed crossbills.
|
||||
#'
|
||||
#' x <- c(1752, 1895)
|
||||
#' x.expected <- ratio(x, ratio = c(1, 1))
|
||||
#' x.expected
|
||||
#' # 1823.5 1823.5
|
||||
#'
|
||||
#' g.test(x, x.expected)
|
||||
#' # p = 0.01787343
|
||||
#'
|
||||
#' # There is a significant difference from a 1:1 ratio.
|
||||
#' # Meaning: there are significantly more left-billed birds.
|
||||
#'
|
||||
#' @inherit g.test examples
|
||||
ratio <- function(x, ratio) {
|
||||
if (!all(is.numeric(x))) {
|
||||
stop('`x` must be a vector of numeric values.')
|
||||
|
||||
@ -20,6 +20,7 @@ globalVariables(c('abname',
|
||||
'antibiotics',
|
||||
'atc',
|
||||
'bactid',
|
||||
'C_chisq_sim',
|
||||
'cnt',
|
||||
'count',
|
||||
'counts',
|
||||
|
||||
@ -2,31 +2,75 @@
|
||||
% Please edit documentation in R/g.test.R
|
||||
\name{g.test}
|
||||
\alias{g.test}
|
||||
\title{\emph{G}-test of matrix or vector}
|
||||
\title{\emph{G}-test for Count Data}
|
||||
\source{
|
||||
This code is almost perfectly equal to \code{\link{chisq.test}}; only the calculation of the statistic was changed to \code{2 * sum(x * log(x / E))} and Yates' continuity correction was deleted.
|
||||
}
|
||||
\usage{
|
||||
g.test(x, y = NULL, alpha = 0.05, info = TRUE, minimum = 0)
|
||||
g.test(x, y = NULL, p = rep(1/length(x), length(x)), rescale.p = FALSE,
|
||||
simulate.p.value = FALSE, B = 2000)
|
||||
}
|
||||
\arguments{
|
||||
\item{x}{numeric vector or matrix}
|
||||
\item{x}{a numeric vector or matrix. \code{x} and \code{y} can also
|
||||
both be factors.}
|
||||
|
||||
\item{y}{expected value of \code{x}. Leave empty to determine automatically. This can also be ratios of \code{x}, e.g. calculated with \code{\link{ratio}}.}
|
||||
\item{y}{a numeric vector; ignored if \code{x} is a matrix. If
|
||||
\code{x} is a factor, \code{y} should be a factor of the same length.}
|
||||
|
||||
\item{alpha}{value to test the p value against}
|
||||
\item{p}{a vector of probabilities of the same length of \code{x}.
|
||||
An error is given if any entry of \code{p} is negative.}
|
||||
|
||||
\item{info}{logical to determine whether the analysis should be printed}
|
||||
\item{rescale.p}{a logical scalar; if TRUE then \code{p} is rescaled
|
||||
(if necessary) to sum to 1. If \code{rescale.p} is FALSE, and
|
||||
\code{p} does not sum to 1, an error is given.}
|
||||
|
||||
\item{minimum}{the test with fail if any of the observed values is below this value. Use \code{minimum = 30} for microbial epidemiology, to prevent calculating a p value when less than 30 isolates are available.}
|
||||
\item{simulate.p.value}{a logical indicating whether to compute
|
||||
p-values by Monte Carlo simulation.}
|
||||
|
||||
\item{B}{an integer specifying the number of replicates used in the
|
||||
Monte Carlo test.}
|
||||
}
|
||||
\value{
|
||||
A list with class \code{"htest"} containing the following
|
||||
components:
|
||||
\item{statistic}{the value the chi-squared test statistic.}
|
||||
\item{parameter}{the degrees of freedom of the approximate
|
||||
chi-squared distribution of the test statistic, \code{NA} if the
|
||||
p-value is computed by Monte Carlo simulation.}
|
||||
\item{p.value}{the p-value for the test.}
|
||||
\item{method}{a character string indicating the type of test
|
||||
performed, and whether Monte Carlo simulation or continuity
|
||||
correction was used.}
|
||||
\item{data.name}{a character string giving the name(s) of the data.}
|
||||
\item{observed}{the observed counts.}
|
||||
\item{expected}{the expected counts under the null hypothesis.}
|
||||
\item{residuals}{the Pearson residuals,
|
||||
\code{(observed - expected) / sqrt(expected)}.}
|
||||
\item{stdres}{standardized residuals,
|
||||
\code{(observed - expected) / sqrt(V)}, where \code{V} is the residual cell variance (Agresti, 2007,
|
||||
section 2.4.5 for the case where \code{x} is a matrix, \code{n * p * (1 - p)} otherwise).}
|
||||
}
|
||||
\description{
|
||||
A \emph{G}-test can be used to see whether the number of observations in each category fits a theoretical expectation (called a \strong{\emph{G}-test of goodness-of-fit}), or to see whether the proportions of one variable are different for different values of the other variable (called a \strong{\emph{G}-test of independence}).
|
||||
\code{g.test} performs chi-squared contingency table tests and goodness-of-fit tests, just like \code{\link{chisq.test}} but is more reliable [1]. A \emph{G}-test can be used to see whether the number of observations in each category fits a theoretical expectation (called a \strong{\emph{G}-test of goodness-of-fit}), or to see whether the proportions of one variable are different for different values of the other variable (called a \strong{\emph{G}-test of independence}).
|
||||
}
|
||||
\details{
|
||||
If \code{x} is a matrix with one row or column, or if \code{x} is a vector and \code{y} is not given, then a \emph{goodness-of-fit test} is performed (\code{x} is treated as a one-dimensional contingency table). The entries of \code{x} must be non-negative integers. In this case, the hypothesis tested is whether the population probabilities equal those in \code{p}, or are all equal if \code{p} is not given.
|
||||
|
||||
If \code{x} is a matrix with at least two rows and columns, it is taken as a two-dimensional contingency table: the entries of \code{x} must be non-negative integers. Otherwise, \code{x} and \code{y} must be vectors or factors of the same length; cases with missing values are removed, the objects are coerced to factors, and the contingency table is computed from these. Then Pearson's chi-squared test is performed of the null hypothesis that the joint distribution of the cell counts in a 2-dimensional contingency table is the product of the row and column marginals.
|
||||
|
||||
If \code{simulate.p.value} is \code{FALSE}, the p-value is computed from the asymptotic chi-squared distribution of the test statistic. Otherwise the p-value is computed for a Monte Carlo test (Hope, 1968) with \code{B} replicates.
|
||||
|
||||
In the contingency table case simulation is done by random sampling from the set of all contingency tables with given marginals, and works only if the marginals are strictly positive. Note that this is not the usual sampling situation assumed for a chi-squared test (like the \emph{G}-test) but rather that for Fisher's exact test.
|
||||
|
||||
In the goodness-of-fit case simulation is done by random sampling from the discrete distribution specified by \code{p}, each sample being of size \code{n = sum(x)}. This simulation is done in \R and may be slow.
|
||||
}
|
||||
\section{\emph{G}-test of goodness-of-fit (likelihood ratio test)}{
|
||||
|
||||
Use the \emph{G}-test of goodness-of-fit when you have one nominal variable with two or more values (such as male and female, or red, pink and white flowers). You compare the observed counts of numbers of observations in each category with the expected counts, which you calculate using some kind of theoretical expectation (such as a 1:1 sex ratio or a 1:2:1 ratio in a genetic cross).
|
||||
|
||||
If the expected number of observations in any category is too small, the \emph{G}-test may give inaccurate results, and you should use an exact test instead. See the web page on small sample sizes for discussion of what "small" means.
|
||||
If the expected number of observations in any category is too small, the \emph{G}-test may give inaccurate results, and you should use an exact test instead (\code{\link{fisher.test}}).
|
||||
|
||||
The \emph{G}-test of goodness-of-fit is an alternative to the chi-square test of goodness-of-fit; each of these tests has some advantages and some disadvantages, and the results of the two tests are usually very similar.
|
||||
The \emph{G}-test of goodness-of-fit is an alternative to the chi-square test of goodness-of-fit (\code{\link{chisq.test}}); each of these tests has some advantages and some disadvantages, and the results of the two tests are usually very similar.
|
||||
}
|
||||
|
||||
\section{\emph{G}-test of independence}{
|
||||
@ -35,24 +79,24 @@ Use the \emph{G}-test of independence when you have two nominal variables, each
|
||||
|
||||
It is also possible to do a \emph{G}-test of independence with more than two nominal variables. For example, Jackson et al. (2013) also had data for children under 3, so you could do an analysis of old vs. young, thigh vs. arm, and reaction vs. no reaction, all analyzed together.
|
||||
|
||||
Fisher's exact test is more accurate than the \emph{G}-test of independence when the expected numbers are small, so it is recommend to only use the \emph{G}-test if your total sample size is greater than 1000.
|
||||
Fisher's exact test (\code{\link{fisher.test}}) is more accurate than the \emph{G}-test of independence when the expected numbers are small, so it is recommend to only use the \emph{G}-test if your total sample size is greater than 1000.
|
||||
|
||||
The \emph{G}-test of independence is an alternative to the chi-square test of independence, and they will give approximately the same results.
|
||||
The \emph{G}-test of independence is an alternative to the chi-square test of independence (\code{\link{chisq.test}}), and they will give approximately the same results.
|
||||
}
|
||||
|
||||
\section{How the test works}{
|
||||
|
||||
Unlike the exact test of goodness-of-fit, the \emph{G}-test does not directly calculate the probability of obtaining the observed results or something more extreme. Instead, like almost all statistical tests, the \emph{G}-test has an intermediate step; it uses the data to calculate a test statistic that measures how far the observed data are from the null expectation. You then use a mathematical relationship, in this case the chi-square distribution, to estimate the probability of obtaining that value of the test statistic.
|
||||
Unlike the exact test of goodness-of-fit (\code{\link{fisher.test}}), the \emph{G}-test does not directly calculate the probability of obtaining the observed results or something more extreme. Instead, like almost all statistical tests, the \emph{G}-test has an intermediate step; it uses the data to calculate a test statistic that measures how far the observed data are from the null expectation. You then use a mathematical relationship, in this case the chi-square distribution, to estimate the probability of obtaining that value of the test statistic.
|
||||
|
||||
The \emph{G}-test uses the log of the ratio of two likelihoods as the test statistic, which is why it is also called a likelihood ratio test or log-likelihood ratio test. The formula to calculate a \emph{G}-statistic is:
|
||||
|
||||
\code{G <- 2 * sum(x * log(x / x.expected))}
|
||||
\code{G <- 2 * sum(x * log(x / E))}
|
||||
|
||||
Since this is chi-square distributed, the p value can be calculated with:
|
||||
where \code{E} are the expected values. Since this is chi-square distributed, the p value can be calculated with:
|
||||
|
||||
\code{p <- 1 - stats::pchisq(G, df))}
|
||||
\code{p <- stats::pchisq(G, df, lower.tail = FALSE)}
|
||||
|
||||
where \code{df} are the degrees of freedom: \code{max(NROW(x) - 1, 1) * max(NCOL(x) - 1, 1)}.
|
||||
where \code{df} are the degrees of freedom.
|
||||
|
||||
If there are more than two categories and you want to find out which ones are significantly different from their null expectation, you can use the same method of testing each category vs. the sum of all categories, with the Bonferroni correction. You use \emph{G}-tests for each category, of course.
|
||||
}
|
||||
@ -68,12 +112,12 @@ If there are more than two categories and you want to find out which ones are si
|
||||
# ratio.
|
||||
|
||||
x <- c(772, 1611, 737)
|
||||
x.expected <- ratio(x, "1:2:1")
|
||||
x.expected
|
||||
E <- ratio(x, "1:2:1")
|
||||
E
|
||||
# 780 1560 780
|
||||
|
||||
g.test(x, x.expected)
|
||||
# p = 0.12574.
|
||||
G <- g.test(x, p = c(1, 2, 1) / 4)
|
||||
# G$p.value = 0.12574.
|
||||
|
||||
# There is no significant difference from a 1:2:1 ratio.
|
||||
# Meaning: resistance controlled by a single gene with two co-dominant
|
||||
@ -88,11 +132,7 @@ g.test(x, x.expected)
|
||||
# observed 1752 right-billed and 1895 left-billed crossbills.
|
||||
|
||||
x <- c(1752, 1895)
|
||||
x.expected <- ratio(x, ratio = c(1, 1))
|
||||
x.expected
|
||||
# 1823.5 1823.5
|
||||
|
||||
g.test(x, x.expected)
|
||||
g.test(x)
|
||||
# p = 0.01787343
|
||||
|
||||
# There is a significant difference from a 1:1 ratio.
|
||||
@ -100,7 +140,7 @@ g.test(x, x.expected)
|
||||
|
||||
}
|
||||
\references{
|
||||
McDonald, J.H. 2014. \strong{Handbook of Biological Statistics (3rd ed.)}. Sparky House Publishing, Baltimore, Maryland. \url{http://www.biostathandbook.com/gtestgof.html}.
|
||||
[1] McDonald, J.H. 2014. \strong{Handbook of Biological Statistics (3rd ed.)}. Sparky House Publishing, Baltimore, Maryland. \url{http://www.biostathandbook.com/gtestgof.html}.
|
||||
}
|
||||
\seealso{
|
||||
\code{\link{chisq.test}}
|
||||
|
||||
14
man/ratio.Rd
14
man/ratio.Rd
@ -25,12 +25,12 @@ Transform vector to ratio
|
||||
# ratio.
|
||||
|
||||
x <- c(772, 1611, 737)
|
||||
x.expected <- ratio(x, "1:2:1")
|
||||
x.expected
|
||||
E <- ratio(x, "1:2:1")
|
||||
E
|
||||
# 780 1560 780
|
||||
|
||||
g.test(x, x.expected)
|
||||
# p = 0.12574.
|
||||
G <- g.test(x, p = c(1, 2, 1) / 4)
|
||||
# G$p.value = 0.12574.
|
||||
|
||||
# There is no significant difference from a 1:2:1 ratio.
|
||||
# Meaning: resistance controlled by a single gene with two co-dominant
|
||||
@ -45,11 +45,7 @@ g.test(x, x.expected)
|
||||
# observed 1752 right-billed and 1895 left-billed crossbills.
|
||||
|
||||
x <- c(1752, 1895)
|
||||
x.expected <- ratio(x, ratio = c(1, 1))
|
||||
x.expected
|
||||
# 1823.5 1823.5
|
||||
|
||||
g.test(x, x.expected)
|
||||
g.test(x)
|
||||
# p = 0.01787343
|
||||
|
||||
# There is a significant difference from a 1:1 ratio.
|
||||
|
||||
@ -18,7 +18,9 @@ test_that("frequency table works", {
|
||||
# factor
|
||||
expect_output(print(freq(septic_patients$hospital_id)))
|
||||
# table
|
||||
if (Sys.info()['sysname'] %in% c("Windows", "Linux")) {
|
||||
expect_output(print(freq(table(septic_patients$sex, septic_patients$age))))
|
||||
}
|
||||
|
||||
library(dplyr)
|
||||
expect_output(septic_patients %>% select(1:2) %>% freq() %>% print())
|
||||
|
||||
@ -7,14 +7,14 @@ test_that("G-test works", {
|
||||
# example 1: clearfield rice vs. red rice
|
||||
x <- c(772, 1611, 737)
|
||||
x.expected <- ratio(x, ratio = "1:2:1")
|
||||
expect_equal(g.test(x, x.expected),
|
||||
expect_equal(g.test(x, p = c(0.25, 0.50, 0.25))$p.value,
|
||||
expected = 0.12574,
|
||||
tolerance = 0.00001)
|
||||
|
||||
# example 2: red crossbills
|
||||
x <- c(1752, 1895)
|
||||
x.expected <- ratio(x, ratio = c(1, 1))
|
||||
expect_equal(g.test(x, x.expected),
|
||||
x.expected <- ratio(x, c(1, 1))
|
||||
expect_equal(g.test(x)$p.value,
|
||||
expected = 0.01787343,
|
||||
tolerance = 0.00000001)
|
||||
|
||||
@ -23,7 +23,7 @@ test_that("G-test works", {
|
||||
x <- matrix(data = round(runif(4) * 100000, 0),
|
||||
ncol = 2,
|
||||
byrow = TRUE)
|
||||
expect_lt(g.test(x),
|
||||
expect_lt(g.test(x)$p.value,
|
||||
1)
|
||||
|
||||
})
|
||||
|
||||
Loading…
Reference in New Issue
Block a user