R Under development (unstable) (2019-02-27 r76167) -- "Unsuffered Consequences" Copyright (C) 2019 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. > ### Regression tests for the r x c x K problem, i.e., > ### testing the independence of a factor > ### 'y' and a factor factor 'x' (possibly blocked) > > suppressWarnings(RNGversion("3.5.2")) > set.seed(290875) > library("coin") Loading required package: survival > isequal <- coin:::isequal > options(useFancyQuotes = FALSE) > > > ### generate data: 2 x 2 x K > dat <- data.frame(x = gl(2, 50), y = gl(2, 50)[sample(1:100)], + block = gl(10, 10)[sample(1:100)])[sample(1:100, 75),] > > ### Pearsons Chisq Test, asymptotic distribution > ptwo <- chisq.test(table(dat$x, dat$y), correct = FALSE)$p.value > > stopifnot(isequal(pvalue(chisq_test(y ~ x, data = dat)), ptwo)) > stopifnot(isequal(pvalue(chisq_test(table(dat$y, dat$x))), ptwo)) > > ### Cochran-Mantel-Haenzel Test, asymptotic distribution > ptwo <- drop(mantelhaen.test(table(dat$x, dat$y, dat$block), + correct = FALSE)$p.value) > > stopifnot(isequal(pvalue(cmh_test(y ~ x | block, data = dat)), ptwo)) > stopifnot(isequal(pvalue(cmh_test(table(dat$y, dat$x, dat$block))), ptwo)) > > > ### generate data: r x c x K > dat <- data.frame(x = gl(4, 25), y = gl(4, 25)[sample(1:100)], + block = gl(2, 50)[sample(1:100)]) > > ### Cochran-Mantel-Haenzel Test, asymptotic distribution > ### (was wrong in R < 2.1.0) > ptwo <- drop(mantelhaen.test(table(dat$y, dat$x, dat$block), + correct = FALSE)$p.value) > > stopifnot(isequal(pvalue(cmh_test(y ~ x | block, data = dat)), ptwo)) > stopifnot(isequal(pvalue(cmh_test(table(dat$y, dat$x, dat$block))), ptwo)) > > ### Linear-by-linear Test, asymptotic distribution > lbl_test(y ~ x | block, data = dat) Asymptotic Linear-by-Linear Association Test data: y (ordered) by x (1 < 2 < 3 < 4) stratified by block Z = -0.35009, p-value = 0.7263 alternative hypothesis: two.sided > lbl_test(table(dat$y, dat$x, dat$block)) Asymptotic Linear-by-Linear Association Test data: Var2 (ordered) by Var1 (1 < 2 < 3 < 4) stratified by Var3 Z = -0.35009, p-value = 0.7263 alternative hypothesis: two.sided > > > ### generate data: r x c x K > dat <- data.frame(x = gl(4, 25), y = gl(5, 20)[sample(1:100)], + block = gl(2, 50)[sample(1:100)]) > > ### Cochran-Mantel-Haenzel Test, asymptotic distribution > ### (was wrong in R < 2.1.0) > ptwo <- drop(mantelhaen.test(table(dat$y, dat$x, dat$block), + correct = FALSE)$p.value) > > stopifnot(isequal(pvalue(cmh_test(y ~ x | block, data = dat)), ptwo)) > stopifnot(isequal(pvalue(cmh_test(table(dat$y, dat$x, dat$block))), ptwo)) > > ### Linear-by-linear Test, asymptotic distribution > lbl_test(y ~ x | block, data = dat) Asymptotic Linear-by-Linear Association Test data: y (ordered) by x (1 < 2 < 3 < 4) stratified by block Z = 1.1436, p-value = 0.2528 alternative hypothesis: two.sided > lbl_test(table(dat$y, dat$x, dat$block)) Asymptotic Linear-by-Linear Association Test data: Var2 (ordered) by Var1 (1 < 2 < 3 < 4 < 5) stratified by Var3 Z = 1.1436, p-value = 0.2528 alternative hypothesis: two.sided > > > ### 2x2 table and maxstat > x <- c(rep(1,51), rep(2,49)) > y <- factor(c(rep(0,49), rep(1,51)))[sample(1:100)] > sit <- as.vector(statistic(independence_test(table(x, y)))) > stopifnot(isequal(as.vector(statistic(maxstat_test(y ~ x))), sit)) > stopifnot(isequal(as.vector(statistic(maxstat_test(table(x, y)))), sit)) > > > ### maxstat for multiple, ordered and unordered covariates > dat <- data.frame(w = rnorm(100), x = runif(100), y = gl(4, 25)[sample(1:100)], + z = ordered(gl(4, 25)[sample(1:100)])) > > mt <- maxstat_test(w ~ x, data = dat) > mt Asymptotic Generalized Maximally Selected Statistics data: w by x maxT = 2.3738, p-value = 0.184 alternative hypothesis: two.sided sample estimates: "best" cutpoint: <= 0.9060508 > est <- mt@estimates$estimate$cutpoint > stopifnot(isequal(statistic(mt), + abs(statistic(independence_test(w ~ (x <= est), data = dat))))) > > mt <- maxstat_test(w ~ y, data = dat) > mt Asymptotic Generalized Maximally Selected Statistics data: w by y (1, 2, 3, 4) maxT = 0.7294, p-value = 0.8905 alternative hypothesis: two.sided sample estimates: "best" cutpoint: {1, 2, 3} vs. {4} > est <- mt@estimates$estimate$cutpoint > xx <- dat$y %in% est > stopifnot(isequal(statistic(mt), + abs(statistic(independence_test(w ~ xx, data = dat))))) > > mt <- maxstat_test(w ~ z, data = dat) > mt Asymptotic Generalized Maximally Selected Statistics data: w by z (1 < 2 < 3 < 4) maxT = 0.52218, p-value = 0.9125 alternative hypothesis: two.sided sample estimates: "best" cutpoint: {1} vs. {2, 3, 4} > est <- mt@estimates$estimate$cutpoint > xx <- dat$z <= est > stopifnot(isequal(statistic(mt), + abs(statistic(independence_test(w ~ xx, data = dat))))) > > mt <- maxstat_test(w ~ x + y + z, data = dat) > mt Asymptotic Generalized Maximally Selected Statistics data: w by x, y, z(ordered) maxT = 2.3738, p-value = 0.286 alternative hypothesis: two.sided sample estimates: "best" cutpoint: <= 0.9060508 covariable: x > est <- mt@estimates$estimate > xsel <- dat[[est[[1]]]] > if (is.factor(xsel) && !is.ordered(xsel)) { + xx <- xsel %in% est[[2]] + } else { + xx <- xsel <= est[[2]] + } > stopifnot(isequal(statistic(mt), + abs(statistic(independence_test(w ~ xx, data = dat))))) > > > ### marginal homogeneity > rating <- c("low", "moderate", "high") > x <- as.table(matrix(c(20, 10, 5, + 3, 30, 15, + 0, 5, 40), + ncol = 3, byrow = TRUE, + dimnames = list(Rater1 = rating, Rater2 = rating))) > ### test statistic W_0 = 13.76 > ### see http://ourworld.compuserve.com/homepages/jsuebersax/mcnemar.htm > stopifnot(all.equal(round(statistic(mh_test(x)), 2), 13.76)) > > proc.time() user system elapsed 2.20 0.12 2.31
Generated by dwww version 1.15 on Sun Jun 16 16:57:50 CEST 2024.