tests/regtest_libcoin.R

In libcoin: Linear Test Statistics for Permutation Inference

library("libcoin")
library("coin")
set.seed(29)

n <- 100
p <- 4
q <- 2
X <- matrix(runif(p * n), nc = p)
Y <- matrix(runif(q * n), nc = q)
w <- as.integer(floor(runif(n, max = 4)))
s <- sort(sample(1:n, floor(n/2), replace = TRUE))
b <- sample(gl(2, 2, length = n))

tol <- if (.Platform$OS.type == "unix") {
    sqrt(.Machine$double.eps)
} else {
   .Machine$double.eps^(1/5)
}

cmp <- function(t1, t2) {
    if (is.null(t1$Covariance)) {
        var1 <- t1$Variance
        var2 <- diag(covariance(t2))
    } else {
        var1 <- t1$Covariance
        var2 <- covariance(t2)
        var2 <- var2[!upper.tri(var2)]
    }
    stopifnot(all.equal(
        list(t1$LinearStatistic,
             t1$Expectation,
             var1),
        list(as.vector(statistic(t2, type = "linear")),
             as.vector(expectation(t2)),
             var2),
        check.attributes = FALSE, tolerance = tol
    ))
}

cmp2 <- function(t1, t2) {
    nm <- c("LinearStatistic", "Expectation",
            if(t1$varonly == 1) "Variance" else "Covariance")
    stopifnot(all.equal(t1[nm], t2[nm], tolerance = tol))
}

cmp3 <- function(t1, t2, pvalue = FALSE) {
    stopifnot(all.equal(statistic(t1), t2$TestStatistic, tolerance = tol))
    if (pvalue)
        stopifnot(all.equal(unclass(pvalue(t1)), t2$p.value, check.attributes = FALSE, tolerance = tol))
}

cmp4 <- function(t1, t2)
    stopifnot(all.equal(t1$TestStatistic, sqrt(t2$TestStatistic), tolerance = tol))


t1 <- LinStatExpCov(X, Y)
t1v <- LinStatExpCov(X, Y, varonly = TRUE)
t2 <- independence_test(Y ~ X)
cmp(t1, t2)
cmp(t1v, t2)
cmp3(t2, doTest(t1, teststat = "maximum"), pvalue = TRUE)
t3 <- independence_test(Y ~ X, teststat = "maximum", alternative = "less")
cmp3(t3, doTest(t1, teststat = "maximum", alternative = "less"), pvalue = TRUE)
t3 <- independence_test(Y ~ X, teststat = "maximum", alternative = "greater")
cmp3(t3, doTest(t1, teststat = "maximum", alternative = "greater"), pvalue = TRUE)
t3 <- independence_test(Y ~ X, teststat = "quadratic")
cmp3(t3, doTest(t1, teststat = "quadratic"), pvalue = TRUE)

t1 <- LinStatExpCov(X, Y, weights = w)
t1v <- LinStatExpCov(X, Y, weights = w, varonly = TRUE)
t2 <- independence_test(Y ~ X, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, subset = s)
t1v <- LinStatExpCov(X, Y, subset = s, varonly = TRUE)
t2 <- independence_test(Y ~ X, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, weights = w, subset = s)
t1v <- LinStatExpCov(X, Y, weights = w, subset = s, varonly = TRUE)
t2 <- independence_test(Y ~ X, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, block = b)
t1v <- LinStatExpCov(X, Y, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X  | b)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, weights = w, block = b)
t1v <- LinStatExpCov(X, Y, weights = w, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X | b, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X | b, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, weights = w, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, weights = w, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X | b, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

n <- 100
n1 <- 5
n2 <- 4
p <- 4
q <- 2
X <- rbind(0, matrix(runif(p * n1), nc = p))
Y <- rbind(0, matrix(runif(q * n2), nc = q))
ix <- sample(1:n1, n, replace = TRUE)
iy <- sample(1:n2, n, replace = TRUE)
w <- as.integer(floor(runif(n, max = 4)))
s <- sort(sample(1:n, floor(n/2), replace = TRUE))
b <- sample(gl(2, 2, length = n))

t1 <- LinStatExpCov(X, Y, ix, iy)
t1v <- LinStatExpCov(X, Y, ix, iy, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,])
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,], weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, subset = s)
t1v <- LinStatExpCov(X, Y, ix, iy, subset = s, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,], subset = s)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,], weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, block = b, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,]  | b)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, block = b, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,] | b, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,] | b, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(Y[iy + 1,] ~ X[ix + 1,]| b, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)

#### X factor
n <- 10000
p <- 40
q <- 20
X <- diag(p)[fx <- unclass(factor(sample(1:p, n, replace = TRUE))),]
Y <- matrix(runif(q * n), nc = q)
w <- as.integer(floor(runif(n, max = 4)))
s <- sort(sample(1:n, floor(n/2), replace = TRUE))
b <- sample(gl(2, 2, length = n))

t1 <- LinStatExpCov(X, Y)
t1v <- LinStatExpCov(X, Y, varonly = TRUE)
t2 <- independence_test(Y ~ X)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y)
t1vf <- LinStatExpCov(fx, Y, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, weights = w)
t1v <- LinStatExpCov(X, Y, weights = w, varonly = TRUE)
t2 <- independence_test(Y ~ X, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, weights = w)
t1vf <- LinStatExpCov(fx, Y, weights = w, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, subset = s)
t1v <- LinStatExpCov(X, Y, subset = s, varonly = TRUE)
t2 <- independence_test(Y ~ X, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, subset = s)
t1vf <- LinStatExpCov(fx, Y, subset = s, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, weights = w, subset = s)
t1v <- LinStatExpCov(X, Y, weights = w, subset = s, varonly = TRUE)
t2 <- independence_test(Y ~ X, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, weights = w, subset = s)
t1vf <- LinStatExpCov(fx, Y, weights = w, subset = s, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, block = b)
t1v <- LinStatExpCov(X, Y, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X  | b)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, block = b)
t1vf <- LinStatExpCov(fx, Y, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)


t1 <- LinStatExpCov(X, Y, weights = w, block = b)
t1v <- LinStatExpCov(X, Y, weights = w, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X | b, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, weights = w, block = b)
t1vf <- LinStatExpCov(fx, Y, weights = w, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X | b, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, subset = s, block = b)
t1vf <- LinStatExpCov(fx, Y, subset = s, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, weights = w, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, weights = w, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(Y ~ X | b, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(fx, Y, weights = w, subset = s, block = b)
t1vf <- LinStatExpCov(fx, Y, weights = w, subset = s, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

n1 <- 5
n2 <- 7
X <- rbind(0, diag(n1))
Y <- rbind(0, matrix(runif(q * n2), nc = q))
ix <- sample(1:n1, n, replace = TRUE)
iy <- sample(1:n2, n, replace = TRUE)
w <- as.integer(floor(runif(n, max = 4)))
s <- sort(sample(1:n, floor(n/2), replace = TRUE))
b <- sample(gl(2, 2, length = n))

YY <- Y[iy + 1,]
XX <- X[ix + 1,]

t1 <- LinStatExpCov(X, Y, ix, iy)
t1v <- LinStatExpCov(X, Y, ix, iy, varonly = TRUE)
t2 <- independence_test(YY ~ XX)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, varonly = TRUE)
t2 <- independence_test(YY ~ XX, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, ix, iy, subset = s)
t1v <- LinStatExpCov(X, Y, ix, iy, subset = s, varonly = TRUE)
t2 <- independence_test(YY ~ XX, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, subset = s)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, subset = s, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)


t1 <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s, varonly = TRUE)
t2 <- independence_test(YY ~ XX, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, subset = s)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, subset = s, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, ix, iy, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, block = b, varonly = TRUE)
t2 <- independence_test(YY ~ XX  | b)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, block = b)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, block = b, varonly = TRUE)
t2 <- independence_test(YY ~ XX | b, weights = ~ w)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, block = b)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, ix, iy, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(YY ~ XX | b, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, subset = s, block = b)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, subset = s, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

t1 <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s, block = b)
t1v <- LinStatExpCov(X, Y, ix, iy, weights = w, subset = s, block = b, varonly = TRUE)
t2 <- independence_test(YY ~ XX| b, weights = ~w, subset = s)
cmp(t1, t2)
cmp(t1v, t2)
t1f <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, subset = s, block = b)
t1vf <- LinStatExpCov(numeric(0), Y, ix, iy, weights = w, subset = s, block = b, varonly = TRUE)
cmp2(t1, t1f)
cmp2(t1v, t1vf)

### and now maximally selected statistics
n <- 100
B <- 1000
x <- round(runif(n), 2)
y <- rnorm(n, mean = x < .5, sd = 2.6)
y2 <- runif(n)
blk <- gl(4, n/4)
ux <- sort(unique(x))
ix <- unclass(cut(x, breaks = c(-Inf, ux[-length(ux)] + diff(ux) / 2, Inf)))

cmp3 <- function(t1, t2)
    stopifnot(all.equal(statistic(t1), t2$TestStatistic, tolerance = tol))

cmp4 <- function(t1, t2)
    stopifnot(all.equal(t1$TestStatistic, sqrt(t2$TestStatistic), tolerance = tol))

(mt <- maxstat_test(y ~ x , distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, y, nresample = 1000)
(tst <- doTest(lev, teststat = "maximum"))
cmp3(mt, tst)
ux[tst$index]
(tst2 <- doTest(lev, teststat = "quadratic"))
cmp4(tst, tst2)
ux[tst2$index]
lev <- LinStatExpCov(ix, y, nresample = 1000, varonly = TRUE)
(tst <- doTest(lev, teststat = "maximum"))
cmp3(mt, tst)
ux[tst$index]
(tst2 <- doTest(lev, teststat = "quadratic"))
cmp4(tst, tst2)
ux[tst2$index]

(mt <- maxstat_test(y ~ x | blk, distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000)
(tst <- doTest(lev, teststat = "maximum"))
cmp3(mt, tst)
ux[tst$index]
(tst2 <- doTest(lev, teststat = "quadratic"))
cmp4(tst, tst2)
ux[tst$index]
lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000, varonly = TRUE)
try(tst <- doTest(lev, teststat = "maximum"))

(mt <- maxstat_test(y + y2 ~ x , distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000)
(tst <- doTest(lev, teststat = "maximum"))
cmp3(mt, tst)
ux[tst$index]
(tst <- doTest(lev, teststat = "quadratic"))
ux[tst$index]
lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000, varonly = TRUE)
(tst <- doTest(lev, teststat = "maximum"))
cmp3(mt, tst)
ux[tst$index]
(tst <- doTest(lev, teststat = "quadratic"))
ux[tst$index]

(mt <- maxstat_test(y + y2 ~ x | blk, distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 1000)
(tst <- doTest(lev, teststat = "maximum"))
cmp3(mt, tst)
ux[tst$index]
(tst <- doTest(lev, teststat = "quadratic"))
ux[tst$index]
lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 1000, varonly = TRUE)
try(tst <- doTest(lev, teststat = "maximum"))

x <- sample(gl(5, n))
y <- rnorm(length(x), mean = x %in% levels(x)[c(1, 3, 5)], sd = 4.5)
y2 <- runif(length(x))
ix <- unclass(x)
blk <- gl(5, n)

(mt <- maxstat_test(y ~ x , distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, y, nresample = 1000)
(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
cmp3(mt, tst)
(tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
cmp4(tst, tst2)
lev <- LinStatExpCov(ix, y, nresample = 1000, varonly = TRUE)
(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
cmp3(mt, tst)
(tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
cmp4(tst, tst2)

(mt <- maxstat_test(y ~ x | blk, distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000)
(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
cmp3(mt, tst)
(tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
cmp4(tst, tst2)
lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000, varonly = TRUE)
try(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))

(mt <- maxstat_test(y + y2 ~ x , distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000)
(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
cmp3(mt, tst)
(tst <- doTest(lev, teststat = "quadratic", ordered = FALSE))
lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000, varonly = TRUE)
(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
cmp3(mt, tst)
(tst <- doTest(lev, teststat = "quadratic", ordered = FALSE))

(mt <- maxstat_test(y + y2 ~ x | blk, distrib = approximate(nresample = 1000)))
lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 50)
(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
cmp3(mt, tst)
(tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))

xx <- factor(x %in% levels(x)[tst2$index == 1])
(it <- independence_test(y + y2 ~ xx | blk, teststat = "quadratic"))
cmp3(it, tst2)

lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 1000, varonly = TRUE)
try(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))