R Under development (unstable) (2019-01-14 r75992) -- "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 multiple adjustments
>
> set.seed(290875)
> library("coin")
Loading required package: survival
> isequal <- coin:::isequal
> options(useFancyQuotes = FALSE)
>
> ### example from Westfall & Wolfinger (1997), Table 4
> tab <- as.table(matrix(c(12, 1, 3, 8, 17, 12, 9, 9, 16, 24), nrow = 2,
+ dimnames = list(group = c("Placebo", "Active"),
+ response = c("Very Poor", "Poor", "Fair",
+ "Good", "Excellent"))))
> df <- coin:::table2df(tab)
>
> it <- independence_test(response ~ group, data = df,
+ distribution = approximate(nresample = 100000))
>
> ### Table 5, first column: OK
> pvalue(it, method = "unadjusted")
Very Poor Poor Fair Good Excellent
Placebo 0.00229 0.11624 0.39398 1 0.07899
>
> ### Table 5, next-to-last column: OK
> pvalue(it, method = "step-down", distribution = "marginal", type = "Sidak")
Very Poor Poor Fair Good Excellent
Placebo 0.006991189 0.2500849 0.5782464 1 0.2309167
>
> ### Table 5, last column: OK
> pvalue(it, method = "step-down")
Very Poor Poor Fair Good Excellent
Placebo 0.00689 0.24334 0.56908 1 0.21351
>
> ### example from Westfall & Wolfinger (1997), Table 1
> df <- data.frame(group = factor(c(rep("Control", 50), rep("Treatment", 48))),
+ V1 = c(rep(0, 50), rep(1, 0), rep(0, 48 - 5), rep(1, 5)),
+ V2 = c(rep(0, 50 - 4), rep(1, 4), rep(0, 48 - 3), rep(1, 3)),
+ V3 = c(rep(0, 50), rep(1, 0), rep(0, 48 - 4), rep(1, 4)),
+ V4 = c(rep(0, 50 - 6), rep(1, 6), rep(0, 48 - 4), rep(1, 4)))
>
> ### alternative: less
> it <- independence_test(V1 + V2 + V3 + V4 ~ group, data = df,
+ distribution = approximate(nresample = 100000),
+ alternative = "less")
>
> ### page 4, 2nd column: adjusted p-value = 0.03665 for V1
> pvalue(it, method = "single-step", distribution = "marginal", type = "Sidak")
V1 V2 V3 V4
Control 0.03666069 0.9842865 0.08844018 0.997373
>
> ### page 4, 2nd column: adjusted p-value = 0.03698 for V1
> ### Note: 0.02521 + 0.00532 + 0 + 0.00645 = 0.03698
> pvalue(it, method = "single-step", distribution = "marginal")
V1 V2 V3 V4
Control 0.03699 1 0.09074 1
>
> ### alternative: two.sided
> it <- independence_test(V1 + V2 + V3 + V4 ~ group, data = df,
+ distribution = approximate(nresample = 100000))
>
> ### page 5, 1st column: adjusted p-value = 0.05261 for V1
> pvalue(it, method = "single-step", distribution = "marginal", type = "Sidak")
V1 V2 V3 V4
Control 0.05407491 1 0.1350998 0.964179
>
> ### page 5, 1st column: adjusted p-value = 0.05352 for V1
> ### Note: 0.02521 + 0.01254 + 0 + 0.01577 = 0.05352
> pvalue(it, method = "single-step", distribution = "marginal")
V1 V2 V3 V4
Control 0.05503 1 0.14162 1
>
> ### artificial example, checked against 'multtest:mt.maxT'
>
> set.seed(290875)
>
> gr <- gl(2, 50)
> x1 <- rnorm(100) + (as.numeric(gr) - 1) * 0.5
> x2 <- rnorm(100) - (as.numeric(gr) - 1) * 0.5
>
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "two.sided"),
+ method = "single-step")
x1 x2
1 0.009466134 0.01080496
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "less"),
+ method = "single-step")
x1 x2
1 0.004742767 0.9999959
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "greater"),
+ method = "single-step")
x1 x2
1 0.9999969 0.005414957
>
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "two.sided"),
+ method = "step-down")
x1 x2
1 0.009466134 0.009466134
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "less"),
+ method = "step-down")
x1 x2
1 0.004742767 0.9972905
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "greater"),
+ method = "step-down")
x1 x2
1 0.9976271 0.005414957
>
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "two.sided",
+ distribution = approximate(nresample = 10000)),
+ method = "single-step")
x1 x2
1 0.0088 0.0097
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "less",
+ distribution = approximate(nresample = 10000)),
+ method = "single-step")
x1 x2
1 0.0048 1
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "greater",
+ distribution = approximate(nresample = 10000)),
+ method = "single-step")
x1 x2
1 1 0.0052
>
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "two.sided",
+ distribution = approximate(nresample = 10000)),
+ method = "step-down")
x1 x2
1 0.0089 0.0089
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "less",
+ distribution = approximate(nresample = 10000)),
+ method = "step-down")
x1 x2
1 0.0053 0.9971
> pvalue(independence_test(x1 + x2 ~ gr, alternative = "greater",
+ distribution = approximate(nresample = 10000)),
+ method = "step-down")
x1 x2
1 0.9981 0.0057
>
> if (FALSE) {
+ #library("multtest")
+ #a <- mt.maxT(t(cbind(x1, x2)), as.numeric(gr) - 1)
+ #a[order(a$index),]
+ #a <- mt.maxT(t(cbind(x1, x2)), as.numeric(gr) - 1, side = "upper")
+ #a[order(a$index),]
+ #a <- mt.maxT(t(cbind(x1, x2)), as.numeric(gr) - 1, side = "lower")
+ #a[order(a$index),]
+ }
>
> ### Monte Carlo distribution
>
> y <- rnorm(20)
> x <- runif(20)
>
> mt <- maxstat_test(y ~ x, distribution = approximate())
> pvalue(mt)
[1] 0.026
99 percent confidence interval:
0.02207744 0.03038090
> pperm(mt, 1)
[1] 0.0452
> qperm(mt, 0.9)
[1] 2.414349
> dperm(mt, qperm(mt, 0.9))
[1] 0.0078
> #support(mt)
>
> mt <- maxstat_test(y ~ x, distribution = approximate(), alternative = "greater")
> pvalue(mt)
[1] 0.0148
99 percent confidence interval:
0.01187122 0.01820033
> pperm(mt, 1)
[1] 0.4143
> qperm(mt, 0.9)
[1] 2.146543
> dperm(mt, qperm(mt, 0.9))
[1] 1e-04
> #support(mt)
>
> mt <- maxstat_test(y ~ x, distribution = approximate(), alternative = "less")
> pvalue(mt)
[1] 0.9739
99 percent confidence interval:
0.9695115 0.9778303
> pperm(mt, 1)
[1] 0.0027
> qperm(mt, 0.9)
[1] -0.1277333
> dperm(mt, qperm(mt, 0.9))
[1] 1e-04
> #support(mt)
>
> ### unadjusted
>
> set.seed(290875)
>
> gr <- gl(3, 50)
> x1 <- rnorm(150) + (as.numeric(gr) - 1) * 0.5
> x2 <- rnorm(150) - (as.numeric(gr) - 1) * 0.5
>
> pvalue(it1 <- independence_test(x1 + x2 ~ gr),
+ method = "unadjusted")
x1 x2
1 1.847081e-05 4.416624e-06
2 5.281476e-01 4.753754e-01
3 2.604284e-04 1.057581e-04
> pvalue(it2 <- independence_test(x1 + x2 ~ gr, alternative = "less"),
+ method = "unadjusted")
x1 x2
1 9.235407e-06 9.999978e-01
2 7.359262e-01 2.376877e-01
3 9.998698e-01 5.287905e-05
> pvalue(it3 <- independence_test(x1 + x2 ~ gr, alternative = "greater"),
+ method = "unadjusted")
x1 x2
1 0.9999907646 2.208312e-06
2 0.2640737994 7.623123e-01
3 0.0001302142 9.999471e-01
>
> pvalue(it4 <- independence_test(x1 + x2 ~ gr,
+ distribution = approximate(nresample = 10000)),
+ method = "unadjusted")
x1 x2
1 <0.0001 <0.0001
2 0.5428 0.4729
3 0.0001 0.0001
> pvalue(it5 <- independence_test(x1 + x2 ~ gr, alternative = "less",
+ distribution = approximate(nresample = 10000)),
+ method = "unadjusted")
x1 x2
1 <0.0001 1.0000
2 0.7370 0.2385
3 0.9998 <0.0001
> pvalue(it6 <- independence_test(x1 + x2 ~ gr, alternative = "greater",
+ distribution = approximate(nresample = 10000)),
+ method = "unadjusted")
x1 x2
1 1.0000 <0.0001
2 0.2653 0.7536
3 0.0003 1.0000
>
> ### consistency of minimum p-value for "global"/"single-step"/"step-down"
>
> set.seed(290875); pg1 <- pvalue(it1)[1]
> set.seed(290875); pss1 <- pvalue(it1, method = "single-step")
> set.seed(290875); psd1 <- pvalue(it1, method = "step-down")
> identical(pg1, min(pss1))
[1] TRUE
> identical(pg1, min(psd1))
[1] TRUE
>
> set.seed(290875); pg2 <- pvalue(it2)[1]
> set.seed(290875); pss2 <- pvalue(it2, method = "single-step")
> set.seed(290875); psd2 <- pvalue(it2, method = "step-down")
> identical(pg2, min(pss2))
[1] TRUE
> identical(pg2, min(psd2))
[1] TRUE
>
> set.seed(290875); pg3 <- pvalue(it3)[1]
> set.seed(290875); pss3 <- pvalue(it3, method = "single-step")
> set.seed(290875); psd3 <- pvalue(it3, method = "step-down")
> identical(pg3, min(pss3))
[1] TRUE
> identical(pg3, min(psd3))
[1] TRUE
>
> pg4 <- pvalue(it4)[1]
> pss4 <- pvalue(it4, method = "single-step")
> psd4 <- pvalue(it4, method = "step-down")
> identical(pg4, min(pss4))
[1] TRUE
> identical(pg4, min(psd4))
[1] TRUE
>
> pg5 <- pvalue(it5)[1]
> pss5 <- pvalue(it5, method = "single-step")
> psd5 <- pvalue(it5, method = "step-down")
> identical(pg5, min(pss5))
[1] TRUE
> identical(pg5, min(psd5))
[1] TRUE
>
> pg6 <- pvalue(it6)[1]
> pss6 <- pvalue(it6, method = "single-step")
> psd6 <- pvalue(it6, method = "step-down")
> identical(pg6, min(pss6))
[1] TRUE
> identical(pg6, min(psd6))
[1] TRUE
>
> ### adjusted marginal asymptotic p-values
>
> pvalue(it1, method = "single-step", distribution = "marginal")
x1 x2
1 0.0001108249 2.649974e-05
2 1.0000000000 1.000000e+00
3 0.0015625705 6.345486e-04
> pvalue(it1, method = "single-step", distribution = "marginal", type = "Sidak")
x1 x2
1 0.0001108198 2.649945e-05
2 0.9889633564 9.791507e-01
3 0.0015615535 6.343809e-04
> pvalue(it1, method = "step-down", distribution = "marginal")
x1 x2
1 9.235407e-05 2.649974e-05
2 9.507509e-01 9.507509e-01
3 7.812853e-04 4.230324e-04
> pvalue(it1, method = "step-down", distribution = "marginal", type = "Sidak")
x1 x2
1 9.235066e-05 2.649945e-05
2 7.247691e-01 7.247691e-01
3 7.810818e-04 4.229653e-04
>
> pvalue(it2, method = "single-step", distribution = "marginal")
x1 x2
1 5.541244e-05 1.0000000000
2 1.000000e+00 1.0000000000
3 1.000000e+00 0.0003172743
> pvalue(it2, method = "single-step", distribution = "marginal", type = "Sidak")
x1 x2
1 5.541116e-05 1.0000000000
2 9.996609e-01 0.8037555042
3 1.000000e+00 0.0003172324
> pvalue(it2, method = "step-down", distribution = "marginal")
x1 x2
1 5.541244e-05 1.0000000000
2 1.000000e+00 0.9507508919
3 1.000000e+00 0.0002643953
> pvalue(it2, method = "step-down", distribution = "marginal", type = "Sidak")
x1 x2
1 5.541116e-05 0.9999999830
2 9.815848e-01 0.6622995367
3 1.000000e+00 0.0002643673
>
> pvalue(it3, method = "single-step", distribution = "marginal")
x1 x2
1 1.0000000000 1.324987e-05
2 1.0000000000 1.000000e+00
3 0.0007812853 1.000000e+00
> pvalue(it3, method = "single-step", distribution = "marginal", type = "Sidak")
x1 x2
1 1.000000000 0.0000132498
2 0.841143280 0.9998196814
3 0.000781031 1.0000000000
> pvalue(it3, method = "step-down", distribution = "marginal")
x1 x2
1 1.000000000 1.324987e-05
2 1.000000000 1.000000e+00
3 0.000651071 1.000000e+00
> pvalue(it3, method = "step-down", distribution = "marginal", type = "Sidak")
x1 x2
1 0.9999999972 0.0000132498
2 0.7066831177 0.9865717243
3 0.0006509015 0.9999999972
>
> ## ### mcp
>
> ## YOY <- data.frame(length = c(46, 28, 46, 37, 32, 41, 42, 45, 38, 44,
> ## 42, 60, 32, 42, 45, 58, 27, 51, 42, 52,
> ## 38, 33, 26, 25, 28, 28, 26, 27, 27, 27,
> ## 31, 30, 27, 29, 30, 25, 25, 24, 27, 30),
> ## site = factor(c(rep("I", 10), rep("II", 10),
> ## rep("III", 10), rep("IV", 10))))
>
> ## ### permutation based Dunnett
> ## it <- independence_test(length ~ site, data = YOY,
> ## xtrafo = mcp_trafo(site = "Dunnett"),
> ## distribution = approximate(10000),
> ## alternative = "two.sided")
> ## pvalue(it, method = "npmcp")
>
> ## ### asymptotic Dunnett
> ## it <- independence_test(length ~ site, data = YOY,
> ## xtrafo = mcp_trafo(site = "Dunnett"),
> ## alternative = "two.sided")
> ## pvalue(it, method = "npmcp")
>
> ## ### asymptotic Dunnett, user-defined w/o column names
> ## cm <- rbind("II vs I" = c(-1, 1, 0, 0),
> ## "III vs I" = c(-1, 0, 1, 0),
> ## "IV vs I" = c(-1, 0, 0, 1))
> ## it <- independence_test(length ~ site, data = YOY,
> ## xtrafo = mcp_trafo(site = cm),
> ## alternative = "two.sided")
> ## pvalue(it, method = "npmcp")
>
> proc.time()
user system elapsed
1.76 0.15 2.07
Generated by dwww version 1.15 on Sun Jun 16 16:34:18 CEST 2024.