https://github.com/cran/coin
Tip revision: 2b221e4e7ade59d5c33d4f9f277e802c90a7656c authored by Torsten Hothorn on 26 April 2013, 00:00:00 UTC
version 1.0-22
version 1.0-22
Tip revision: 2b221e4
MTP.R
### single step maxT multiple testing procedure
singlestep <- function(object, ...) {
if (object@statistic@alternative == "two.sided") {
ts <- abs(statistic(object, "standardized"))
} else {
ts <- statistic(object, "standardized")
}
### <FIXME> iterate over unique(ts) only and remap
ret <- 1 - sapply(ts, pperm, object = object, ...)
### </FIXME>
ret <- matrix(ret, nrow = nrow(ts), ncol = ncol(ts))
rownames(ret) <- rownames(ts)
colnames(ret) <- colnames(ts)
ret
}
### algorithm 2.8 (Free Step-Down Resampling Method) in
### Westfall & Young (1993), page 66 _using standardized
### statistics instead of p-values_!
### <FIXME>
sdmaxT <- function(pls, ts) {
### order of original statistics
rts <- order(ts)
### algorithm 2.8 (Free Step-Down Resampling Method) in
### Westfall & Young (1993), page 66 _using standardized
### statistics instead of p-values_!
q <- pls[,rts, drop = FALSE]
if (ncol(q) > 1) {
for (j in 2:ncol(q))
q[,j] <- pmax(q[,j], q[,j-1])
}
ret <- matrix(rowMeans(GE(t(q), ts[rts]))[rank(ts)],
nrow = nrow(ts), ncol = ncol(ts))
rownames(ret) <- rownames(ts)
colnames(ret) <- colnames(ts)
ret
}
### stepdown maxT multiple testing procedure
stepdown <- function(object, ...) {
if (!(extends(class(object), "MaxTypeIndependenceTest") &&
extends(class(object@distribution), "ApproxNullDistribution")))
stop(sQuote("object"), " is not of class ",
sQuote("MaxTypeIndependenceTest"),
" or distribution was not approximated via Monte-Carlo")
### raw simulation results, scores have been handled already
pls <- support(object, raw = TRUE)
### standardize
dcov <- sqrt(variance(object))
expect <- expectation(object)
pls <- t((pls - expect) / dcov)
ts <- statistic(object, "standardized")
if (object@statistic@alternative == "two.sided") {
pls <- abs(pls)
ts <- abs(ts)
}
if (object@statistic@alternative == "less") {
pls <- -pls
ts <- -ts
}
sdmaxT(pls, ts)
}
### Bonferroni permutation method (Westfall & Wolfinger, 1997, AmStat 51, 3-8)
dbonf <- function(object, ...) {
### <FIXME> this should be possible when the _exact_ marginal
### distributions are available
### </FIXME>
if (!(extends(class(object), "MaxTypeIndependenceTest") &&
extends(class(object@distribution), "ApproxNullDistribution")))
stop(sQuote("object"), " is not of class ",
sQuote("MaxTypeIndependenceTest"),
" or distribution was not approximated via Monte-Carlo")
alternative <- object@statistic@alternative
### standardize
dcov <- sqrt(variance(object))
expect <- expectation(object)
### raw simulation results, scores have been handled already
pls <- support(object, raw = TRUE)
pls <- (pls - expect) / dcov
ts <- (statistic(object, "standardized"))
pvals <- switch(alternative,
"less" = rowMeans(LE(pls, as.vector(drop(ts)))),
"greater" = rowMeans(GE(pls, as.vector(drop(ts)))),
"two.sided" = rowMeans(GE(abs(pls), as.vector(abs(drop(ts))))))
foo <- function(x, t)
switch(alternative,
"less" = mean(LE(x, t)),
"greater" = mean(GE(x, t)),
"two.sided" = mean(GE(abs(x), abs(t))))
p <- vector(mode = "list", length = nrow(pls))
for (i in 1:nrow(pls)) {
ux <- unique(pls[i,])
p[[i]] <- sapply(ux, foo, x = pls[i,])
}
### Bonferroni adjustment (Westfall & Wolfinger, 1997)
adjp <- rep(1, length(ts))
for (i in 1:length(pvals)) {
for (q in 1:length(p)) {
x <- p[[q]][p[[q]] <= pvals[i]]
if (length(x) > 0)
adjp[i] <- adjp[i] * (1 - max(x))
}
}
ret <- matrix(1 - pmin(adjp, 1), nrow = nrow(ts), ncol = ncol(ts))
rownames(ret) <- rownames(ts)
colnames(ret) <- colnames(ts)
ret
}