R语言如何查看函数源码[例:wilcox.test]
程序员文章站
2022-04-16 15:43:55
...
R语言如何查看函数源码[例:wilcox.test]
命令
## 查看源码
>stats:::待查看的对象
以wilcox.test为例
## 查看wilcox.test的源码
>stats:::wilcox.test#命令行
function (x, ...)
UseMethod("wilcox.test")
<bytecode: 0x0000024da5b4a4f0>
<environment: namespace:stats>
> stats:::wilcox.test.default#命令行
function (x, y = NULL, alternative = c("two.sided", "less",
"greater"), mu = 0, paired = FALSE, exact = NULL, correct = TRUE,
conf.int = FALSE, conf.level = 0.95, ...)
{
alternative <- match.arg(alternative)
if (!missing(mu) && ((length(mu) > 1L) || !is.finite(mu)))
stop("'mu' must be a single number")
if (conf.int) {
if (!((length(conf.level) == 1L) && is.finite(conf.level) &&
(conf.level > 0) && (conf.level < 1)))
stop("'conf.level' must be a single number between 0 and 1")
}
if (!is.numeric(x))
stop("'x' must be numeric")
if (!is.null(y)) {
if (!is.numeric(y))
stop("'y' must be numeric")
DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
if (paired) {
if (length(x) != length(y))
stop("'x' and 'y' must have the same length")
OK <- complete.cases(x, y)
x <- x[OK] - y[OK]
y <- NULL
}
else {
x <- x[is.finite(x)]
y <- y[is.finite(y)]
}
}
else {
DNAME <- deparse(substitute(x))
if (paired)
stop("'y' is missing for paired test")
x <- x[is.finite(x)]
}
if (length(x) < 1L)
stop("not enough (finite) 'x' observations")
CORRECTION <- 0
if (is.null(y)) {
METHOD <- "Wilcoxon signed rank test"
x <- x - mu
ZEROES <- any(x == 0)
if (ZEROES)
x <- x[x != 0]
n <- as.double(length(x))
if (is.null(exact))
exact <- (n < 50)
r <- rank(abs(x))
STATISTIC <- setNames(sum(r[x > 0]), "V")
TIES <- length(r) != length(unique(r))
if (exact && !TIES && !ZEROES) {
PVAL <- switch(alternative, two.sided = {
p <- if (STATISTIC > (n * (n + 1)/4)) psignrank(STATISTIC -
1, n, lower.tail = FALSE) else psignrank(STATISTIC,
n)
min(2 * p, 1)
}, greater = psignrank(STATISTIC - 1, n, lower.tail = FALSE),
less = psignrank(STATISTIC, n))
if (conf.int) {
x <- x + mu
alpha <- 1 - conf.level
diffs <- outer(x, x, "+")
diffs <- sort(diffs[!lower.tri(diffs)])/2
cint <- switch(alternative, two.sided = {
qu <- qsignrank(alpha/2, n)
if (qu == 0) qu <- 1
ql <- n * (n + 1)/2 - qu
achieved.alpha <- 2 * psignrank(trunc(qu) -
1, n)
c(diffs[qu], diffs[ql + 1])
}, greater = {
qu <- qsignrank(alpha, n)
if (qu == 0) qu <- 1
achieved.alpha <- psignrank(trunc(qu) - 1,
n)
c(diffs[qu], +Inf)
}, less = {
qu <- qsignrank(alpha, n)
if (qu == 0) qu <- 1
ql <- n * (n + 1)/2 - qu
achieved.alpha <- psignrank(trunc(qu) - 1,
n)
c(-Inf, diffs[ql + 1])
})
if (achieved.alpha - alpha > alpha/2) {
warning("requested conf.level not achievable")
conf.level <- 1 - signif(achieved.alpha, 2)
}
attr(cint, "conf.level") <- conf.level
ESTIMATE <- c(`(pseudo)median` = median(diffs))
}
}
else {
NTIES <- table(r)
z <- STATISTIC - n * (n + 1)/4
SIGMA <- sqrt(n * (n + 1) * (2 * n + 1)/24 - sum(NTIES^3 -
NTIES)/48)
if (correct) {
CORRECTION <- switch(alternative, two.sided = sign(z) *
0.5, greater = 0.5, less = -0.5)
METHOD <- paste(METHOD, "with continuity correction")
}
z <- (z - CORRECTION)/SIGMA
PVAL <- switch(alternative, less = pnorm(z), greater = pnorm(z,
lower.tail = FALSE), two.sided = 2 * min(pnorm(z),
pnorm(z, lower.tail = FALSE)))
if (conf.int) {
x <- x + mu
alpha <- 1 - conf.level
if (n > 0) {
mumin <- min(x)
mumax <- max(x)
W <- function(d) {
xd <- x - d
xd <- xd[xd != 0]
nx <- length(xd)
dr <- rank(abs(xd))
zd <- sum(dr[xd > 0]) - nx * (nx + 1)/4
NTIES.CI <- table(dr)
SIGMA.CI <- sqrt(nx * (nx + 1) * (2 * nx +
1)/24 - sum(NTIES.CI^3 - NTIES.CI)/48)
if (SIGMA.CI == 0)
warning("cannot compute confidence interval when all observations are zero or tied",
call. = FALSE)
CORRECTION.CI <- if (correct) {
switch(alternative, two.sided = sign(zd) *
0.5, greater = 0.5, less = -0.5)
}
else 0
(zd - CORRECTION.CI)/SIGMA.CI
}
Wmumin <- W(mumin)
Wmumax <- if (!is.finite(Wmumin))
NA
else W(mumax)
}
if (n == 0 || !is.finite(Wmumax)) {
cint <- structure(c(if (alternative == "less") -Inf else NaN,
if (alternative == "greater") +Inf else NaN),
conf.level = 0)
ESTIMATE <- if (n > 0)
c(midrange = (mumin + mumax)/2)
else NaN
}
else {
wdiff <- function(d, zq) W(d) - zq
root <- function(zq) {
uniroot(wdiff, lower = mumin, upper = mumax,
f.lower = Wmumin - zq, f.upper = Wmumax -
zq, tol = 1e-04, zq = zq)$root
}
cint <- switch(alternative, two.sided = {
repeat {
mindiff <- Wmumin - qnorm(alpha/2, lower.tail = FALSE)
maxdiff <- Wmumax - qnorm(alpha/2)
if (mindiff < 0 || maxdiff > 0) alpha <- alpha *
2 else break
}
if (alpha >= 1 || 1 - conf.level < alpha *
0.75) {
conf.level <- 1 - pmin(1, alpha)
warning("requested conf.level not achievable")
}
if (alpha < 1) {
l <- root(zq = qnorm(alpha/2, lower.tail = FALSE))
u <- root(zq = qnorm(alpha/2))
c(l, u)
} else {
rep(median(x), 2)
}
}, greater = {
repeat {
mindiff <- Wmumin - qnorm(alpha, lower.tail = FALSE)
if (mindiff < 0) alpha <- alpha * 2 else break
}
if (alpha >= 1 || 1 - conf.level < alpha *
0.75) {
conf.level <- 1 - pmin(1, alpha)
warning("requested conf.level not achievable")
}
l <- if (alpha < 1) root(zq = qnorm(alpha,
lower.tail = FALSE)) else median(x)
c(l, +Inf)
}, less = {
repeat {
maxdiff <- Wmumax - qnorm(alpha/2)
if (maxdiff > 0) alpha <- alpha * 2 else break
}
if (alpha >= 1 || 1 - conf.level < alpha *
0.75) {
conf.level <- 1 - pmin(1, alpha)
warning("requested conf.level not achievable")
}
u <- if (alpha < 1) root(zq = qnorm(alpha)) else median(x)
c(-Inf, u)
})
attr(cint, "conf.level") <- conf.level
correct <- FALSE
ESTIMATE <- c(`(pseudo)median` = uniroot(W,
lower = mumin, upper = mumax, tol = 1e-04)$root)
}
}
if (exact && TIES) {
warning("cannot compute exact p-value with ties")
if (conf.int)
warning("cannot compute exact confidence interval with ties")
}
if (exact && ZEROES) {
warning("cannot compute exact p-value with zeroes")
if (conf.int)
warning("cannot compute exact confidence interval with zeroes")
}
}
}
else {
if (length(y) < 1L)
stop("not enough 'y' observations")
METHOD <- "Wilcoxon rank sum test"
r <- rank(c(x - mu, y))
n.x <- as.double(length(x))
n.y <- as.double(length(y))
if (is.null(exact))
exact <- (n.x < 50) && (n.y < 50)
STATISTIC <- c(W = sum(r[seq_along(x)]) - n.x * (n.x +
1)/2)
TIES <- (length(r) != length(unique(r)))
if (exact && !TIES) {
PVAL <- switch(alternative, two.sided = {
p <- if (STATISTIC > (n.x * n.y/2)) pwilcox(STATISTIC -
1, n.x, n.y, lower.tail = FALSE) else pwilcox(STATISTIC,
n.x, n.y)
min(2 * p, 1)
}, greater = {
pwilcox(STATISTIC - 1, n.x, n.y, lower.tail = FALSE)
}, less = pwilcox(STATISTIC, n.x, n.y))
if (conf.int) {
alpha <- 1 - conf.level
diffs <- sort(outer(x, y, "-"))
cint <- switch(alternative, two.sided = {
qu <- qwilcox(alpha/2, n.x, n.y)
if (qu == 0) qu <- 1
ql <- n.x * n.y - qu
achieved.alpha <- 2 * pwilcox(trunc(qu) - 1,
n.x, n.y)
c(diffs[qu], diffs[ql + 1])
}, greater = {
qu <- qwilcox(alpha, n.x, n.y)
if (qu == 0) qu <- 1
achieved.alpha <- pwilcox(trunc(qu) - 1, n.x,
n.y)
c(diffs[qu], +Inf)
}, less = {
qu <- qwilcox(alpha, n.x, n.y)
if (qu == 0) qu <- 1
ql <- n.x * n.y - qu
achieved.alpha <- pwilcox(trunc(qu) - 1, n.x,
n.y)
c(-Inf, diffs[ql + 1])
})
if (achieved.alpha - alpha > alpha/2) {
warning("Requested conf.level not achievable")
conf.level <- 1 - achieved.alpha
}
attr(cint, "conf.level") <- conf.level
ESTIMATE <- c(`difference in location` = median(diffs))
}
}
else {
NTIES <- table(r)
z <- STATISTIC - n.x * n.y/2
SIGMA <- sqrt((n.x * n.y/12) * ((n.x + n.y + 1) -
sum(NTIES^3 - NTIES)/((n.x + n.y) * (n.x + n.y -
1))))
if (correct) {
CORRECTION <- switch(alternative, two.sided = sign(z) *
0.5, greater = 0.5, less = -0.5)
METHOD <- paste(METHOD, "with continuity correction")
}
z <- (z - CORRECTION)/SIGMA
PVAL <- switch(alternative, less = pnorm(z), greater = pnorm(z,
lower.tail = FALSE), two.sided = 2 * min(pnorm(z),
pnorm(z, lower.tail = FALSE)))
if (conf.int) {
alpha <- 1 - conf.level
mumin <- min(x) - max(y)
mumax <- max(x) - min(y)
W <- function(d) {
dr <- rank(c(x - d, y))
NTIES.CI <- table(dr)
dz <- sum(dr[seq_along(x)]) - n.x * (n.x +
1)/2 - n.x * n.y/2
CORRECTION.CI <- if (correct) {
switch(alternative, two.sided = sign(dz) *
0.5, greater = 0.5, less = -0.5)
}
else 0
SIGMA.CI <- sqrt((n.x * n.y/12) * ((n.x + n.y +
1) - sum(NTIES.CI^3 - NTIES.CI)/((n.x + n.y) *
(n.x + n.y - 1))))
if (SIGMA.CI == 0)
warning("cannot compute confidence interval when all observations are tied",
call. = FALSE)
(dz - CORRECTION.CI)/SIGMA.CI
}
wdiff <- function(d, zq) W(d) - zq
Wmumin <- W(mumin)
Wmumax <- W(mumax)
root <- function(zq) {
f.lower <- Wmumin - zq
if (f.lower <= 0)
return(mumin)
f.upper <- Wmumax - zq
if (f.upper >= 0)
return(mumax)
uniroot(wdiff, lower = mumin, upper = mumax,
f.lower = f.lower, f.upper = f.upper, tol = 1e-04,
zq = zq)$root
}
cint <- switch(alternative, two.sided = {
l <- root(zq = qnorm(alpha/2, lower.tail = FALSE))
u <- root(zq = qnorm(alpha/2))
c(l, u)
}, greater = {
l <- root(zq = qnorm(alpha, lower.tail = FALSE))
c(l, +Inf)
}, less = {
u <- root(zq = qnorm(alpha))
c(-Inf, u)
})
attr(cint, "conf.level") <- conf.level
correct <- FALSE
ESTIMATE <- c(`difference in location` = uniroot(W,
lower = mumin, upper = mumax, tol = 1e-04)$root)
}
if (exact && TIES) {
warning("cannot compute exact p-value with ties")
if (conf.int)
warning("cannot compute exact confidence intervals with ties")
}
}
}
names(mu) <- if (paired || !is.null(y))
"location shift"
else "location"
RVAL <- list(statistic = STATISTIC, parameter = NULL, p.value = as.numeric(PVAL),
null.value = mu, alternative = alternative, method = METHOD,
data.name = DNAME)
if (conf.int)
RVAL <- c(RVAL, list(conf.int = cint, estimate = ESTIMATE))
class(RVAL) <- "htest"
RVAL
}
<bytecode: 0x0000024da7f81db8>
<environment: namespace:stats>
> stats:::wilcox.test.formula#命令行
function (formula, data, subset, na.action, ...)
{
if (missing(formula) || (length(formula) != 3L) || (length(attr(terms(formula[-2L]),
"term.labels")) != 1L))
stop("'formula' missing or incorrect")
m <- match.call(expand.dots = FALSE)
if (is.matrix(eval(m$data, parent.frame())))
m$data <- as.data.frame(data)
m[[1L]] <- quote(stats::model.frame)
m$... <- NULL
mf <- eval(m, parent.frame())
DNAME <- paste(names(mf), collapse = " by ")
names(mf) <- NULL
response <- attr(attr(mf, "terms"), "response")
g <- factor(mf[[-response]])
if (nlevels(g) != 2L)
stop("grouping factor must have exactly 2 levels")
DATA <- setNames(split(mf[[response]], g), c("x", "y"))
y <- do.call("wilcox.test", c(DATA, list(...)))
y$data.name <- DNAME
y
}
<bytecode: 0x0000024da7fb4a40>
<environment: namespace:stats>