Nothing
# Copyright 2014 Google Inc. All rights reserved.
#
# Use of this source code is governed by a BSD-style
# license that can be found in the LICENSE file or at
# http://opensource.org/licenses/BSD-3-Clause
bootstrap2 <- function(data, statistic,
treatment, data2 = NULL,
R = 10000, ratio = FALSE, args.stat = NULL,
seed = NULL,
sampler = samp.bootstrap,
label = NULL, statisticNames = NULL,
block.size = 100, trace = FALSE)
{
# Two-sample nonparametric bootstrap.
# Specify either (data, data2) or (data, treatment).
#
# Args:
# data: vector, matrix, or data frame.
# Let 'nObs' = length of vector, or nrow otherwise.
# statistic: a function, or expression (e.g. mean(data, trim = .2)
# If data is a data frame, can refer to variables in it.
# This may be a vector; let 'd' be its length.
# treatment: a vector of length nObs, with two unique values.
# data2: like data.
# R: number of replications.
# ratio: logical, if FALSE then bootstrap the difference between
# data and data2 (or first and second treatment).
# If TRUE then bootstrap the ratio.
# args.stat: additional arguments to pass to the function.
# seed: old value of .Random.seed, or argument to set.seed.
# sampler: a function for resampling, see help(samp.bootstrap).
# label: used for labeling plots.
# statisticNames: names used for printing, character vector of length 'd'.
# block.size: replicates are done 'block.size' at a time.
# trace: logical, if TRUE an indication of progress is printed.
Call <- match.call()
dimData <- dim(data)
n <- IfElse(is.null(dimData), length(data), dimData[1])
stopifnot(length(dimData) <= 2)
resultsBoth <- vector("list", 2)
if(is.null(data2)) { # use treatment
if(is.data.frame(data))
treatment <- eval(substitute(treatment),
envir = data, enclos = parent.frame())
treatmentInds <- split(1:n, treatment)
treatmentNames <- names(treatmentInds)
if(length(treatmentInds) != 2)
stop("treatment must have 2 unique values, observed: ",
paste(treatmentNames, collapse = " "))
for(k in 1:2) {
# Create a smaller data set corresponding to this treatment
dataK <- IfElse(length(dimData) == 2,
data[treatmentInds[[k]], , drop = FALSE],
data[treatmentInds[[k]]])
resampleFun <-
.resampleMakeFunction(dataK, statistic,
substitute(data), substitute(statistic),
args.stat)
resultsBoth[[k]] <-
resample(dataK, resampleFun, sampler = sampler, R = R,
seed = IfElse(k == 1, seed, .Random.seed),
statisticNames = statisticNames,
block.size = block.size, trace = trace, call = Call)
}
} else { # data and data2
if(length(dim(data2)) != length(dimData) ||
class(data2) != class(data))
stop("data and data2 must have similar structure")
treatmentNames <- c(IfElse(is.name(substitute(data)),
substitute(data), "data"),
IfElse(is.name(substitute(data2)),
substitute(data2), "data2"))
for(k in 1:2) {
resampleFun <-
.resampleMakeFunction(IfElse(k == 1, data, data2), statistic,
as.name("data"), substitute(statistic), args.stat)
resultsBoth[[k]] <-
resample(IfElse(k == 1, data, data2),
resampleFun, sampler = sampler, R = R,
seed = IfElse(k == 1, seed, .Random.seed),
statisticNames = statisticNames,
block.size = block.size, trace = trace, call = Call)
}
}
for(k in 1:2) {
resultsBoth[[k]]$stats <- .BootStats(resultsBoth[[k]])
resultsBoth[[k]]$call <- paste("bootstrap for data set", k)
class(resultsBoth[[k]]) <- c("bootstrap", "resample")
}
names(resultsBoth) <- treatmentNames
if(!.CheckCompatibleObserved(resultsBoth))
return(resultsBoth)
# combine results
op <- get(IfElse(ratio, "/", "-"))
result <- list(observed =
op(resultsBoth[[1]]$observed, resultsBoth[[2]]$observed),
replicates =
op(resultsBoth[[1]]$replicates, resultsBoth[[2]]$replicates),
n = c(resultsBoth[[1]]$n, resultsBoth[[2]]$n),
p = resultsBoth[[1]]$p,
R = R,
seed = resultsBoth[[1]]$seed,
call = Call,
resultsBoth = resultsBoth,
ratio = ratio)
result$failures <- unique(c(resultsBoth[[1]]$failures,
resultsBoth[[2]]$failures))
statisticNames <- .StatisticNames2(statisticNames, treatmentNames, ratio,
resultsBoth[[1]]$observed, Call)
names(result$observed) <- statisticNames
colnames(result$replicates) <- statisticNames
result$stats <- .BootStats(result)
class(result) <- c("bootstrap2", "bootstrap", "resample")
result
}
# Note difference:
# bootstrap(x, mean(x))
# bootstrap2(x1, data2 = x2, mean(data))
.CheckCompatibleObserved <- function(x){
# x is a pair of objects from bootstrap2 or permutationTest2.
# Check that observed statistics have the same names
# return(TRUE) if they do.
if(length(x[[1]]$observed) == length(x[[2]]$observed) &&
identical(names(x[[1]]$observed), names(x[[2]]$observed)))
return(TRUE)
warning("Statistics returned from the two samples are not compatible,",
"will return results so far, without combining them.")
print(all.equal(x[[1]]$observed, x[[2]]$observed))
return(FALSE)
}
.StatisticNames2 <- function(statisticNames, treatmentNames, ratio,
observed1, call){
# Create statisticNames like "mean: data1 - data2"
# statisticNames correspond to one sample
if(is.null(statisticNames))
statisticNames <- names(observed1)
p <- length(observed1)
if(is.null(statisticNames) && is.call(call) && is.name(call$statistic))
statisticNames <- IfElse(p == 1, as.character(call$statistic),
paste0(as.character(call$statistic), 1:p))
if(is.null(statisticNames))
statisticNames <- paste0("stat", 1:p)
temp <- which(statisticNames == "")
if(length(temp))
statisticNames[temp] <- paste0("stat", temp)
statisticNames <- paste0(statisticNames, ": ",
paste(treatmentNames,
collapse = IfElse(ratio, "/", "-")))
}
# print.bootstrap2 <- function(x, ...) {
# cat0n("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"))
# catn("Replications:", x$R)
# catn("Two groups, sample sizes are", x$n[1], x$n[2])
# catn("\nSummary Statistics for the", IfElse(x$ratio, "ratio", "difference"),
# "between samples 1 and 2:")
# print(x$stats, ...)
# invisible(x)
# }
# same as print.permutatonTest2
if(FALSE) {
x9 <- 1:9
xDF <- data.frame(X = x9, Y = 2 * x9)
t9 <- letters[c(1, 1, 1, 2, 1, 2, 1, 2, 1)]
xDF1 <- xDF[t9 == "a", ]
xDF2 <- xDF[t9 == "b", ]
x91 <- x9[t9 == "a"]
x92 <- x9[t9 == "b"]
source("~/Rlang/resample/resample/R/bootstrap2.R")
##### treatment
### statistic by name
# base case: data by name, statistic is function by name
temp <- bootstrap2(treatment = t9, x9, mean, R=100)
temp
sapply(temp$resultsBoth, class)
temp$resultsBoth
# data expression
bootstrap2(treatment = t9, (x9), mean, R=100)
# args.stat
bootstrap2(treatment = t9, x9, mean, args.stat = list(trim = .25), R=100)
# inline function
bootstrap2(treatment = t9, x9, function(z) mean(z), R=100)
# data frame
bootstrap2(treatment = t9, xDF, colMeans, R=100)
# data expression, data frame
bootstrap2(treatment = t9, (xDF), colMeans, R=100)
# data expression, matrix
bootstrap2(treatment = t9, as.matrix(xDF), colMeans, R=100)
### statistic expression
# data by name
bootstrap2(treatment = t9, x9, mean(x9), R=100)
# data as expression, refer to 'data'
bootstrap2(treatment = t9, (x9), mean(data), R=100)
# data frame
bootstrap2(treatment = t9, xDF, mean(X), R=100)
# data frame expression
bootstrap2(treatment = t9, (xDF), mean(X), R=100)
# See if results reproduce
temp <- .Last.value
.Random.seed <- temp$seed
all.equal(temp, eval(temp$call))
##### data and data2
### statistic by name
# base case: data by name, statistic is function by name
temp <- bootstrap2(x91, data2 = x92, mean, R=100)
temp
sapply(temp$resultsBoth, class)
temp$resultsBoth
# data expression
bootstrap2((x91), data2 = (x92), mean, R=100)
# args.stat
bootstrap2(x91, data2 = x92, mean, args.stat = list(trim = .25), R=100)
# inline function
bootstrap2(x91, data2 = x92, function(z) mean(z), R=100)
# data frame
bootstrap2(xDF1, data2 = xDF2, colMeans, R=100)
# data expression, data frame
bootstrap2((xDF1), data2 = (xDF2), colMeans, R=100)
# data expression, matrix
bootstrap2(as.matrix(xDF1), data2 = as.matrix(xDF2), colMeans, R=100)
### statistic expression
# data by name
bootstrap2(x91, data2 = x92, mean(data)) # use 'data' in this case
# data as expression, refer to 'data'
bootstrap2((x91), data2 = (x92), mean(data), R=100)
# data frame
bootstrap2(xDF1, data2 = xDF2, mean(X), R=100)
# data frame expression
bootstrap2((xDF1), data2 = (xDF2), mean(X), R=100)
source("~/Rlang/resample/resample/R/bootstrap2.R")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.