R Under development (unstable) (2020-12-20 r79658) -- "Unsuffered Consequences"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
>
> library("libcoin")
> library("coin")
Loading required package: survival
> 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)))
Approximative Generalized Maximally Selected Statistics
data: y by x
maxT = 2.1537, p-value = 0.26
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: <= 0.53
> lev <- LinStatExpCov(ix, y, nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum"))
$TestStatistic
[1] 2.153664
$p.value
[1] 0.245
$index
[1] 36
> cmp3(mt, tst)
> ux[tst$index]
[1] 0.53
> (tst2 <- doTest(lev, teststat = "quadratic"))
$TestStatistic
[1] 4.63827
$p.value
[1] 0.245
$index
[1] 36
> cmp4(tst, tst2)
> ux[tst2$index]
[1] 0.53
> lev <- LinStatExpCov(ix, y, nresample = 1000, varonly = TRUE)
> (tst <- doTest(lev, teststat = "maximum"))
$TestStatistic
[1] 2.153664
$p.value
[1] 0.253
$index
[1] 36
> cmp3(mt, tst)
> ux[tst$index]
[1] 0.53
> (tst2 <- doTest(lev, teststat = "quadratic"))
$TestStatistic
[1] 4.63827
$p.value
[1] 0.253
$index
[1] 36
> cmp4(tst, tst2)
> ux[tst2$index]
[1] 0.53
>
> (mt <- maxstat_test(y ~ x | blk, distrib = approximate(nresample = 1000)))
Approximative Generalized Maximally Selected Statistics
data: y by x
stratified by blk
maxT = 2.2764, p-value = 0.221
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: <= 0.53
> lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum"))
$TestStatistic
[1] 2.276446
$p.value
[1] 0.192
$index
[1] 36
> cmp3(mt, tst)
> ux[tst$index]
[1] 0.53
> (tst2 <- doTest(lev, teststat = "quadratic"))
$TestStatistic
[1] 5.182204
$p.value
[1] 0.192
$index
[1] 36
> cmp4(tst, tst2)
> ux[tst$index]
[1] 0.53
> lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000, varonly = TRUE)
> try(tst <- doTest(lev, teststat = "maximum"))
Error in doTest(lev, teststat = "maximum") :
need covariance for maximally statistics with blocks
>
> (mt <- maxstat_test(y + y2 ~ x , distrib = approximate(nresample = 1000)))
Approximative Generalized Maximally Selected Statistics
data: y, y2 by x
maxT = 2.1537, p-value = 0.445
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: <= 0.53
> lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum"))
$TestStatistic
[1] 2.153664
$p.value
[1] 0.471
$index
[1] 36
> cmp3(mt, tst)
> ux[tst$index]
[1] 0.53
> (tst <- doTest(lev, teststat = "quadratic"))
$TestStatistic
[1] 5.720153
$p.value
[1] 0.412
$index
[1] 40
> ux[tst$index]
[1] 0.57
> lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000, varonly = TRUE)
> (tst <- doTest(lev, teststat = "maximum"))
$TestStatistic
[1] 2.153664
$p.value
[1] 0.476
$index
[1] 36
> cmp3(mt, tst)
> ux[tst$index]
[1] 0.53
> (tst <- doTest(lev, teststat = "quadratic"))
$TestStatistic
[1] 5.720153
$p.value
[1] 0.429
$index
[1] 40
> ux[tst$index]
[1] 0.57
>
> (mt <- maxstat_test(y + y2 ~ x | blk, distrib = approximate(nresample = 1000)))
Approximative Generalized Maximally Selected Statistics
data: y, y2 by x
stratified by blk
maxT = 2.2764, p-value = 0.38
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: <= 0.53
> lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum"))
$TestStatistic
[1] 2.276446
$p.value
[1] 0.394
$index
[1] 36
> cmp3(mt, tst)
> ux[tst$index]
[1] 0.53
> (tst <- doTest(lev, teststat = "quadratic"))
$TestStatistic
[1] 5.406605
$p.value
[1] 0.504
$index
[1] 40
> ux[tst$index]
[1] 0.57
> lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 1000, varonly = TRUE)
> try(tst <- doTest(lev, teststat = "maximum"))
Error in doTest(lev, teststat = "maximum") :
need covariance for maximally statistics with blocks
>
> 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)))
Approximative Generalized Maximally Selected Statistics
data: y by x (1, 2, 3, 4, 5)
maxT = 1.7088, p-value = 0.443
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: {1, 3} vs. {2, 4, 5}
> lev <- LinStatExpCov(ix, y, nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
$TestStatistic
[1] 1.708775
$p.value
[1] 0.485
$index
[1] 0 1 0 1 1
> cmp3(mt, tst)
> (tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
$TestStatistic
[1] 2.919912
$p.value
[1] 0.485
$index
[1] 0 1 0 1 1
> cmp4(tst, tst2)
> lev <- LinStatExpCov(ix, y, nresample = 1000, varonly = TRUE)
> (tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
$TestStatistic
[1] 1.708775
$p.value
[1] 0.476
$index
[1] 0 1 0 1 1
> cmp3(mt, tst)
> (tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
$TestStatistic
[1] 2.919912
$p.value
[1] 0.476
$index
[1] 0 1 0 1 1
> cmp4(tst, tst2)
>
> (mt <- maxstat_test(y ~ x | blk, distrib = approximate(nresample = 1000)))
Approximative Generalized Maximally Selected Statistics
data: y by x (1, 2, 3, 4, 5)
stratified by blk
maxT = 1.6489, p-value = 0.543
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: {1, 3} vs. {2, 4, 5}
> lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
$TestStatistic
[1] 1.648947
$p.value
[1] 0.464
$index
[1] 0 1 0 1 1
> cmp3(mt, tst)
> (tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
$TestStatistic
[1] 2.719028
$p.value
[1] 0.464
$index
[1] 0 1 0 1 1
> cmp4(tst, tst2)
> lev <- LinStatExpCov(ix, y, block = blk, nresample = 1000, varonly = TRUE)
> try(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
Error in doTest(lev, teststat = "maximum", ordered = FALSE) :
need covariance for maximally statistics with blocks
>
> (mt <- maxstat_test(y + y2 ~ x , distrib = approximate(nresample = 1000)))
Approximative Generalized Maximally Selected Statistics
data: y, y2 by x (1, 2, 3, 4, 5)
maxT = 2.358, p-value = 0.285
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: {1, 3, 5} vs. {2, 4}
> lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000)
> (tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
$TestStatistic
[1] 2.357973
$p.value
[1] 0.301
$index
[1] 0 1 0 1 0
> cmp3(mt, tst)
> (tst <- doTest(lev, teststat = "quadratic", ordered = FALSE))
$TestStatistic
[1] 8.373114
$p.value
[1] 0.143
$index
[1] 0 1 0 1 1
> lev <- LinStatExpCov(ix, cbind(y, y2), nresample = 1000, varonly = TRUE)
> (tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
$TestStatistic
[1] 2.357973
$p.value
[1] 0.275
$index
[1] 0 1 0 1 0
> cmp3(mt, tst)
> (tst <- doTest(lev, teststat = "quadratic", ordered = FALSE))
$TestStatistic
[1] 8.373114
$p.value
[1] 0.134
$index
[1] 0 1 0 1 1
>
> (mt <- maxstat_test(y + y2 ~ x | blk, distrib = approximate(nresample = 1000)))
Approximative Generalized Maximally Selected Statistics
data: y, y2 by x (1, 2, 3, 4, 5)
stratified by blk
maxT = 2.4473, p-value = 0.218
alternative hypothesis: two.sided
sample estimates:
"best" cutpoint: {1, 3, 5} vs. {2, 4}
> lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 50)
> (tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
$TestStatistic
[1] 2.447272
$p.value
[1] 0.2
$index
[1] 0 1 0 1 0
> cmp3(mt, tst)
> (tst2 <- doTest(lev, teststat = "quadratic", ordered = FALSE))
$TestStatistic
[1] 8.115438
$p.value
[1] 0.12
$index
[1] 0 1 0 1 0
>
> xx <- factor(x %in% levels(x)[tst2$index == 1])
> (it <- independence_test(y + y2 ~ xx | blk, teststat = "quadratic"))
Asymptotic General Independence Test
data: y, y2 by xx (FALSE, TRUE)
stratified by blk
chi-squared = 8.1154, df = 2, p-value = 0.01729
> cmp3(it, tst2)
>
> lev <- LinStatExpCov(ix, cbind(y, y2), block = blk, nresample = 1000, varonly = TRUE)
> try(tst <- doTest(lev, teststat = "maximum", ordered = FALSE))
Error in doTest(lev, teststat = "maximum", ordered = FALSE) :
need covariance for maximally statistics with blocks
>
> proc.time()
user system elapsed
4.17 0.34 4.56
Generated by dwww version 1.15 on Sat May 18 08:46:41 CEST 2024.