#' One-sample permutation test for the center of a distribution
#'
#' @param x Vector of observations
#' @param reps Number of replications to approximate distribution (default 1000)
#' @param center Optional; specify null center of distribution (default 0)
#'
#' @return A vector of length `reps` containing the permutation distribution
#' @examples
#' x <- -5:5 + rnorm(11)
#' distr <- one_sample(x)
#' pvalue <- t2p(mean(x), distr, alternative = 'two-sided')
#' pvalue
#'
#' distr2 <- one_sample(x, center = 5)
#' pvalue2 <- t2p(mean(x), distr2, alternative = 'greater')
#' pvalue2
one_sample <- function(x, reps = 1000, center = 0) {
if (!is.vector(x)) {
stop("x must be a vector")
}
if (!is.numeric(x)) {
stop("x must be numeric")
}
x_center <- x - center
n <- length(x)
tst <- function(z) {
mean(x_center * z + center)
}
distr <- replicate(reps, {
perm <- 1 - 2 * rbinom(n, 1, 0.5)
tst(perm)
})
return(distr)
}
#' Two-sample permutation test for equality of means
#'
#' @param x Vector of observations in group 1
#' @param y Vector of observations in group 2
#' @param reps Number of replications to approximate distribution (default 1000)
#'
#' @return A vector of length `reps` containing the permutation distribution
two_sample <- function(x, y, reps = 1000) {
if (!is.vector(x) | !is.vector(y)) {
stop("x and y must be vectors")
}
if (!is.numeric(x) | !is.numeric(y)) {
stop("x and y must be numeric")
}
tst <- function(x, y) {
mean(x) - mean(y)
}
z <- c(x, y)
n <- length(x)
distr <- replicate(reps, {
z.star <- sample(z)
x.star <- z.star[1:n]
y.star <- z.star[-(1:n)]
tst(x.star, y.star)
})
return(distr)
}
#' Difference in means within groups
#'
#' Compute difference in mean residual between treated and control within each group
#' @param group Vector of group memberships or treatment conditions
#' @param response Vector of measured outcomes, same length as group
#' @param stratum Vector of stratum assignments, same length as group
#' @param groups Vector of unique group labels
#' @param strata Vector of unique stratum labels
#'
#' @return a vector of differences
within_group_mean <- function(group, response, stratum, groups, strata) {
tt <- (group == groups[1])
sapply(strata, function(s) {
ind <- stratum == s
treated <- response[tt == 1 & ind]
ctrl <- response[tt == 0 & ind]
return(mean(treated, na.rm = TRUE) - mean(ctrl, na.rm = TRUE))
})
}
#' Stratified two-sample permutation test for equality of means
#'
#' @param group Vector of group memberships or treatment conditions
#' @param response Vector of measured outcomes, same length as group
#' @param stratum Vector of stratum assignments, same length as group
#' @param stat The test statistic. Default is 'mean'. See details for other options.
#' @param reps Number of replications to approximate distribution (default 1000)
#' @details
#' \enumerate{
#' \item If stat == 'mean', the test statistic is (mean(x) - mean(y))
#' (equivalently, sum(x), since those are monotonically related), omitting
#' NaNs, which therefore can be used to code non-responders
#' \item If stat == 't', the test statistic is the two-sample t-statistic--
#' but the p-value is still estimated by the randomization,
#' approximating the permutation distribution.
#' The t-statistic is computed using t.test(...,var.equal=TRUE)
#' \item If stat == 'mean_within_strata', the test statistic is the difference
#' in means within each stratum, added across strata.
#' \item If stat is a function (a callable object), the test statistic is
#' that function. The function should take a permutation of the pooled
#' data and compute the test function from it.
#' }
#' @return A vector of length `reps` containing the permutation distribution
stratified_two_sample <- function(group, response, stratum, stat = "mean", reps = 1000) {
if (!is.vector(group) | !is.vector(response) | !is.vector(stratum)) {
stop("inputs must be vectors")
}
if (!is.numeric(response)) {
stop("response must be numeric")
}
if (length(unique(group)) > 2) {
stop("two samples only")
}
groups <- unique(group)
strata <- unique(stratum)
ordering <- order(group)
response <- response[ordering]
stratum <- stratum[ordering]
group <- group[ordering]
ntreat <- table(group)[1]
N <- length(group)
# If stat is callable, use it as the test function. Otherwise, look in the
# dictionary
stats <- list(mean = function(u) {
mean(u[1:ntreat], na.rm = TRUE) - mean(u[(ntreat + 1):N], na.rm = TRUE)
}, t = function(u) {
t.test(u[1:ntreat], u[(ntreat + 1):N], var.equal = TRUE)$statistic
}, mean_within_strata = function(u) {
sum(abs(within_group_mean(group, u, stratum, groups, strata)))
})
if (is.function(stat)) {
tst_fun <- stat
} else {
if (stat %in% names(stats)) {
tst_fun <- stats[[stat]]
} else {
stop("stat must be in the dictionary of stats or a function")
}
}
distr <- replicate(reps, {
tst_fun(permute_within_groups(response, strata))
})
return(distr)
}
#' Confidence interval for the additive shift in means for a one- or two-sample problem
#'
#' Invert the permutation test to get a \eqn{1-\alpha} confidence interval for either the shift parameter in a two-sample problem or the center of the distribution in a one-sample problem
#' @inheritParams two_sample
#' @param side Type of interval, either 'both', 'upper', or 'lower'. Default is 'both'.
#' @param alpha Significance level
#' @return a vector containing the desired limits of the confidence set
CI_mean <- function(x, y = NULL, reps = 1000, side = "both", alpha = 0.05) {
if (!(side %in% c("upper", "lower", "both"))) {
stop("side must be one of 'upper', 'lower', or 'both'")
}
if (side == "both") {
d1 <- CI_mean(x = x, y = y, reps = reps, side = "lower", alpha/2)
d2 <- CI_mean(x = x, y = y, reps = reps, side = "upper", alpha/2)
return(c(d1, d2))
}
alternative <- ifelse(side == "upper", "less", "greater")
# initialize
shift_pval <- 1
if (is.null(y)) {
observed <- mean(x)
shift_permtest <- function(ss) {
shift_fun <- function(shift) one_sample(x, reps = reps, center = shift)
distr <- shift_fun(ss)
return(t2p(observed, shift_fun(ss), alternative = alternative))
}
} else {
observed <- mean(x) - mean(y)
shift_permtest <- function(ss) {
shift_fun <- function(shift) two_sample(x - shift, y, rep = reps)
distr <- shift_fun(ss)
return(t2p(observed - ss, shift_fun(ss), alternative = alternative))
}
}
d_prev <- observed
d_next <- observed
incr <- diff(range(c(x, y)))/2
# Conduct permutation test for H0: shift = d until we reject H0
while (shift_pval > alpha) {
d_prev <- d_next
d_next <- ifelse(side == "upper", d_prev + incr, d_prev - incr)
shift_pval <- shift_permtest(d_next)
}
# Bisection
if (side == "upper") {
l_int <- d_prev
u_int <- d_next
} else {
l_int <- d_next
u_int <- d_prev
}
d <- uniroot(function(z) {
shift_permtest(z) - alpha
}, lower = l_int, upper = u_int, tol = alpha/reps)
return(d$root)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.