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.