R Under development (unstable) (2018-11-02 r75535) -- "Unsuffered Consequences" Copyright (C) 2018 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. > ### compare output for examples from StatXact 6 manual with > ### StatXact 6 output as given in the manual > > library("coin") Loading required package: survival > set.seed(290875) > isequal <- coin:::isequal > options(useFancyQuotes = FALSE) > > ### marginal homogenity > > ### StatXact 6 manual, page 315 > presidents <- as.table(matrix(c(28, 7, 13, 27), nrow = 2, + dimnames = list(BeforeTVDebate = c("Carter", "Reagan"), + AfterTVDebate = c("Carter", "Reagan")))) > > bta <- mh_test(presidents) > > # test statistic, page 306 > stopifnot(isequal(round(sqrt(statistic(bta)), 3), 1.342)) > > # two-sided asymptotic p-value, page 306 > stopifnot(isequal(round(pvalue(bta), 4), 0.1797)) > > # exact p-value, page 306 > btMC <- mh_test(presidents, distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(btMC), "conf.int") > stopifnot(pci[1] < 0.2632 & pci[2] > 0.2632) > > > ### StatXact 6 manual, page 308 > endometrial_cancer <- as.table(matrix(c(6, 9, 9, 12, + 2, 4, 2, 1, + 3, 2, 3, 2, + 1, 1, 1, 1), nrow = 4, + dimnames = list(Cases = c(0, 0.2, 0.5125, 0.7), + Controls = c(0, 0.2, 0.5125, 0.7)))) > > bta <- mh_test(endometrial_cancer, + scores = list(response = c(0, 0.2, 0.5125, 0.7))) > > # test statistic, page 311 > stopifnot(isequal(round(statistic(bta), 3), 3.735)) > > # two-sided asymptotic p-value, page 311 > stopifnot(isequal(round(pvalue(bta), 4), 0.0002)) > > ### StatXact 6 manual, page 312 > pathologists <- as.table(matrix(c(22, 5, 0, 0, 0, + 2, 7, 2, 1, 0, + 2, 14, 36, 14, 3, + 0, 0, 0, 7, 0, + 0, 0, 0, 0, 3), nrow = 5, + dimnames = list(Pathologist1 = paste("Level", 1:5, sep = "-"), + Pathologist2 = paste("Level", 1:5, sep = "-")))) > > bta <- mh_test(pathologists, + scores = list(response = 1:5)) > > # test statistic, page 313 > stopifnot(isequal(round(statistic(bta), 3), 1.152)) > > # two-sided asymptotic p-value, page 313 > stopifnot(isequal(round(pvalue(bta), 4), 0.2492)) > > # exact p-value, page 313 > btMC <- mh_test(pathologists, scores = list(response = 1:5), + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(btMC), "conf.int") > stopifnot(pci[1] < 0.3073 & pci[2] > 0.3073) > > > ### Two independent samples > > ### StatXact 6 manual, page 340 > load("bloodp.rda") > > wta <- wilcox_test(bp ~ group, data = bloodp) > > # test statistic, page 342 > stopifnot(isequal(round(statistic(wta), 3), 1.720)) > > # two-sided asymptotic p-value, page 342 > stopifnot(isequal(round(pvalue(wta), 4), 0.0853)) > > wte <- wilcox_test(bp ~ group, data = bloodp, distribution = "exact") > > # two-sided exact p-value, page 342 > stopifnot(isequal(round(pvalue(wte), 4), 0.0989)) > > wtel <- wilcox_test(bp ~ group, data = bloodp, distribution = "exact", + alternative = "greater") > > # one-sided exact p-value, page 342 > stopifnot(isequal(round(pvalue(wtel), 4), 0.0542)) > > # two-sided approximated p-value > wtMC <- wilcox_test(bp ~ group, data = bloodp, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(wtMC), "conf.int") > > stopifnot(pci[1] < pvalue(wte) & pci[2] > pvalue(wte)) > > # one-sided approximated p-value > wtMC <- wilcox_test(bp ~ group, data = bloodp, alternative = "greater", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(wtMC), "conf.int") > > stopifnot(pci[1] < pvalue(wtel) & pci[2] > pvalue(wtel)) > > # deal with ties as described in StatXact 6 manual, page 329ff > nta <- normal_test(bp ~ group, data = bloodp, ties.method = "average") > > # test statistic, page 354 > stopifnot(isequal(round(statistic(nta), 3), 1.789)) > > # two-sided asymptotic p-value, page 354 > stopifnot(isequal(round(pvalue(nta), 4), 0.0737)) > > nte <- normal_test(bp ~ group, data = bloodp, distribution = "exact", + ties.method = "average") > > # two-sided exact p-value, page 354 > stopifnot(isequal(round(pvalue(nte), 4), 0.0799)) > > ntel <- normal_test(bp ~ group, data = bloodp, distribution = "exact", + ties.method = "average", alternative = "greater") > > # one-sided exact p-value, page 354 > stopifnot(isequal(round(pvalue(ntel), 4), 0.0462)) > > # two-sided approximated p-value > ntMC <- normal_test(bp ~ group, data = bloodp, ties.method = "average", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ntMC), "conf.int") > > stopifnot(pci[1] < pvalue(nte) & pci[2] > pvalue(nte)) > > # one-sided approximated p-value > ntMC <- normal_test(bp ~ group, data = bloodp, alternative = "greater", + distribution = approximate(nresample = 10000), + ties.method = "average") > pci <- attr(pvalue(ntMC), "conf.int") > > stopifnot(pci[1] < pvalue(ntel) & pci[2] > pvalue(ntel)) > > pta <- oneway_test(bp ~ group, data = bloodp) > > # test statistic, page 394 > stopifnot(isequal(round(statistic(pta), 3), 1.612)) > > # two-sided asymptotic p-value, page 394 > stopifnot(isequal(round(pvalue(pta), 4), 0.1070)) > > pte <- oneway_test(bp ~ group, data = bloodp, distribution = "exact") > > # two-sided exact p-value, page 394 > stopifnot(isequal(round(pvalue(pte), 4), 0.1040)) > > ptel <- oneway_test(bp ~ group, data = bloodp, distribution = "exact", + alternative = "greater") > > # one-sided exact p-value, page 394 > stopifnot(isequal(round(pvalue(ptel), 4), 0.0564)) > > # two-sided approximated p-value > ptMC <- oneway_test(bp ~ group, data = bloodp, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ptMC), "conf.int") > > stopifnot(pci[1] < pvalue(pte) & pci[2] > pvalue(pte)) > > # one-sided approximated p-value > ptMC <- oneway_test(bp ~ group, data = bloodp, alternative = "greater", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ptMC), "conf.int") > > stopifnot(pci[1] < pvalue(ptel) & pci[2] > pvalue(ptel)) > > > ### StatXact 6 manual, page 345 > load("employment.rda") > > wta <- wilcox_test(Salary ~ Gender | Year, data = employment) > > # test statistic, page 346 > stopifnot(isequal(round(statistic(wta), 3), -1.897)) > > # two-sided asymptotic p-value, page 346 > stopifnot(isequal(round(pvalue(wta), 4), 0.0578)) > > wte <- wilcox_test(Salary ~ Gender | Year, data = employment, + distribution = "exact") > > # two-sided exact p-value, page 346 > stopifnot(isequal(round(pvalue(wte), 4), 0.04)) > > wtel <- wilcox_test(Salary ~ Gender | Year, data = employment, + distribution = "exact", alternative = "less") > > # one-sided exact p-value, page 346 > stopifnot(isequal(round(pvalue(wtel), 4), 0.04)) > > # two-sided approximated p-value > wtMC <- wilcox_test(Salary ~ Gender | Year, data = employment, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(wtMC), "conf.int") > > stopifnot(pci[1] < pvalue(wte) & pci[2] > pvalue(wte)) > > # one-sided approximated p-value > wtMC <- wilcox_test(Salary ~ Gender | Year, data = employment, + alternative = "less", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(wtMC), "conf.int") > > stopifnot(pci[1] < pvalue(wtel) & pci[2] > pvalue(wtel)) > > nta <- normal_test(Salary ~ Gender | Year, data = employment, + ties.method = "average") > > # test statistic, page 358 > stopifnot(isequal(round(statistic(nta), 3), -1.802)) > > # two-sided asymptotic p-value, page 358 > stopifnot(isequal(round(pvalue(nta), 4), 0.0716)) > > # blocks not yet implemented > # nte <- normal_test(Salary ~ Gender | Year, data = employment, > # distribution = "exact") > # > # two-sided exact p-value, page 358 > # stopifnot(isequal(round(pvalue(nte), 4), 0.04)) > # > # ntel <- normal_test(Salary ~ Gender | Year, data = employment, > # distribution = "exact", alternative = "less") > # > # one-sided exact p-value, page 358 > # stopifnot(isequal(round(pvalue(ntel), 4), 0.04)) > > # two-sided approximated p-value > ntMC <- normal_test(Salary ~ Gender | Year, data = employment, + ties.method = "average", distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ntMC), "conf.int") > > stopifnot(pci[1] < 0.04 & pci[2] > 0.04) > > # one-sided approximated p-value > ntMC <- normal_test(Salary ~ Gender | Year, data = employment, + alternative = "less", distribution = approximate(nresample = 10000), + ties.method = "average") > pci <- attr(pvalue(ntMC), "conf.int") > > stopifnot(pci[1] < 0.04 & pci[2] > 0.04) > > pta <- oneway_test(Salary ~ Gender | Year, data = employment) > > # test statistic, page 396 > stopifnot(isequal(round(statistic(pta), 3), -1.873)) > > # two-sided asymptotic p-value, page 396 > stopifnot(isequal(round(pvalue(pta), 4), 0.0610)) > > pte <- oneway_test(Salary ~ Gender | Year, data = employment, + distribution = "exact") > > # two-sided exact p-value, page 396 > stopifnot(isequal(round(pvalue(pte), 4), 0.04)) > > ptel <- oneway_test(Salary ~ Gender | Year, data = employment, + distribution = "exact", alternative = "less") > > # one-sided exact p-value, page 396 > stopifnot(isequal(round(pvalue(ptel), 4), 0.04)) > > # two-sided approximated p-value > ptMC <- oneway_test(Salary ~ Gender | Year, data = employment, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ptMC), "conf.int") > > stopifnot(pci[1] < pvalue(pte) & pci[2] > pvalue(pte)) > > # one-sided approximated p-value > ptMC <- oneway_test(Salary ~ Gender | Year, data = employment, + alternative = "less", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ptMC), "conf.int") > > stopifnot(pci[1] < pvalue(ptel) & pci[2] > pvalue(ptel)) > > > ### StatXact 9 manual, page 182 > machines <- data.frame(cereal = c(10.8, 11.1, 10.4, 10.1, 11.3, + 10.8, 10.5, 11.0, 10.9, 10.8, + 10.7, 10.8), + machine = factor(rep(1:2, c(5, 7)), + labels = c("Present", "New"))) > > ata <- ansari_test(cereal ~ machine, data = machines, + ties.method = "average") > > # test statistic, page 182 > stopifnot(isequal(round(statistic(ata), 4), -1.9983)) > > # two-sided asymptotic p-value, page 182 > stopifnot(isequal(round(pvalue(ata), 4), 0.0457)) > > ate <- ansari_test(cereal ~ machine, data = machines, + ties.method = "average", distribution = "exact") > > # two-sided exact p-value, page 182 > stopifnot(isequal(round(pvalue(ate), 4), 0.0581)) > > atMC <- ansari_test(cereal ~ machine, data = machines, + ties.method = "average", distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(atMC), "conf.int") > > # two-sided approximated p-value, page 182 > stopifnot(pci[1] < pvalue(ate) & pci[2] > pvalue(ate)) > > # Note: StatXact has '.LE.' here since *small* test statistics furnish evidence > # for the alternative that sample 1 is *more* variable than sample 2, but > # 'coin' relieves the user from thinking about that detail > > atag <- ansari_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater") > > # one-sided asymptotic p-value, page 182 > stopifnot(isequal(round(pvalue(atag), 4), 0.0228)) > > ateg <- ansari_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater", + distribution = "exact") > > # one-sided exact p-value, page 182 > stopifnot(isequal(round(pvalue(ateg), 4), 0.0253)) > > atMC <- ansari_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(atMC), "conf.int") > > # one-sided approximated p-value, page 182 > stopifnot(pci[1] < pvalue(ateg) & pci[2] > pvalue(ateg)) > > > ### StatXact 9 manual, page 186 > > kta <- klotz_test(cereal ~ machine, data = machines, + ties.method = "average") > > # test statistic, page 186 > stopifnot(isequal(round(statistic(kta), 4), 2.3082)) > > # two-sided asymptotic p-value, page 186 > stopifnot(isequal(round(pvalue(kta), 4), 0.0210)) > > kte <- klotz_test(cereal ~ machine, data = machines, + ties.method = "average", distribution = "exact") > > # two-sided exact p-value, page 186 > stopifnot(isequal(round(pvalue(kte), 4), 0.0101)) > > ktMC <- klotz_test(cereal ~ machine, data = machines, + ties.method = "average", distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ktMC), "conf.int") > > # two-sided approximated p-value, page 186 > stopifnot(pci[1] < pvalue(kte) & pci[2] > pvalue(kte)) > > ktag <- klotz_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater") > > # one-sided asymptotic p-value, page 186 > stopifnot(isequal(round(pvalue(ktag), 4), 0.0105)) > > kteg <- klotz_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater", + distribution = "exact") > > # one-sided exact p-value, page 186 > stopifnot(isequal(round(pvalue(kteg), 4), 0.0101)) > > ktMC <- klotz_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ktMC), "conf.int") > > # one-sided approximated p-value, page 186 > stopifnot(pci[1] < pvalue(kteg) & pci[2] > pvalue(kteg)) > > > ### StatXact 9 manual, page 190 > > mta <- mood_test(cereal ~ machine, data = machines, + ties.method = "average") > > # test statistic, page 190 > stopifnot(isequal(round(statistic(mta), 4), 2.2715)) > > # two-sided asymptotic p-value, page 190 > stopifnot(isequal(round(pvalue(mta), 4), 0.0231)) > > mte <- mood_test(cereal ~ machine, data = machines, + ties.method = "average", distribution = "exact") > > # two-sided exact p-value, page 190 > stopifnot(isequal(round(pvalue(mte), 4), 0.0202)) > > mtMC <- mood_test(cereal ~ machine, data = machines, + ties.method = "average", distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(mtMC), "conf.int") > > # two-sided approximated p-value, page 190 > stopifnot(pci[1] < pvalue(mte) & pci[2] > pvalue(mte)) > > mtag <- mood_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater") > > # one-sided asymptotic p-value, page 190 > stopifnot(isequal(round(pvalue(mtag), 4), 0.0116)) > > mteg <- mood_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater", + distribution = "exact") > > # one-sided exact p-value, page 190 > stopifnot(isequal(round(pvalue(mteg), 4), 0.0126)) > > mtMC <- mood_test(cereal ~ machine, data = machines, + ties.method = "average", alternative = "greater", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(mtMC), "conf.int") > > # one-sided approximated p-value, page 190 > stopifnot(pci[1] < pvalue(mteg) & pci[2] > pvalue(mteg)) > > > ### StatXact 9 manual, page 195 > load("FAILURE.rda") > > cta <- conover_test(ftime ~ tyretype, data = failure) > > # test statistic, page 195 > stopifnot(isequal(round(statistic(cta), 4), -1.5274)) > > # two-sided asymptotic p-value, page 195 > stopifnot(isequal(round(pvalue(cta), 4), 0.1267)) > > cte <- conover_test(ftime ~ tyretype, data = failure, + distribution = "exact") > > # two-sided exact p-value, page 195 > stopifnot(isequal(round(pvalue(cte), 4), 0.1300)) > > ctMC <- conover_test(ftime ~ tyretype, data = failure, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ctMC), "conf.int") > > # two-sided approximated p-value, page 195 > stopifnot(pci[1] < pvalue(cte) & pci[2] > pvalue(cte)) > > ctag <- conover_test(ftime ~ tyretype, data = failure, + alternative = "less") > > # one-sided asymptotic p-value, page 195 > stopifnot(isequal(round(pvalue(ctag), 4), 0.0633)) > > cteg <- conover_test(ftime ~ tyretype, data = failure, + alternative = "less", distribution = "exact") > > # one-sided exact p-value, page 195 > stopifnot(isequal(round(pvalue(cteg), 4), 0.0634)) > > ctMC <- conover_test(ftime ~ tyretype, data = failure, + alternative = "less", distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ctMC), "conf.int") > > # one-sided approximated p-value, page 195 > stopifnot(pci[1] < pvalue(cteg) & pci[2] > pvalue(cteg)) > > > ### StatXact 6 manual, 413 > load("lungcancer.rda") > > ### NOTE: StatXact uses another tie handling method, see Callaert (2003) > ### and the examples below (at the end of this file) > > lta <- logrank_test(Surv(time, cens) ~ group, data = lungcancer, ties = "aver") > > # test statistic, page 415 > isequal(round(statistic(lta), 3), 2.949) # no "-", StatXact 9, p. 215 [1] TRUE > > # two-sided asymptotic p-value, page 415 > isequal(round(pvalue(lta), 4), 0.0032) [1] TRUE > > lte <- logrank_test(Surv(time, cens) ~ group, data = lungcancer, + distribution = "exact") > > # two-sided exact p-value, page 415 > stopifnot(isequal(round(pvalue(lte), 4), 0.0010)) > > ltel <- logrank_test(Surv(time, cens) ~ group, data = lungcancer, + distribution = "exact", alternative = "great") > > # one-sided exact p-value, page 415 > stopifnot(isequal(round(pvalue(ltel), 4), 0.0010)) > > # two-sided approximated p-value > ltMC <- logrank_test(Surv(time, cens) ~ group, data = lungcancer, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ltMC), "conf.int") > > stopifnot(pci[1] < pvalue(lte) & pci[2] > pvalue(lte)) > > # one-sided approximated p-value > ltMC <- logrank_test(Surv(time, cens) ~ group, data = lungcancer, + alternative = "greater", + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ltMC), "conf.int") > > stopifnot(pci[1] < pvalue(ltel) & pci[2] > pvalue(ltel)) > > > ### StatXact manual 6, page 418 > srv <- data.frame(time = c(3, 5, 7, 8, 18, 12, 19, 20, 20, 33), + event = c(0, 0, 1, 1, 1, 1, 1, 1, 0, 0), + treatment = gl(2, 5), + gender = factor(c("M", "M", "F", "F", "M", + "M", "M", "F", "F", "F"))) > > # page 419 > logrank_test(Surv(time, event) ~ treatment | gender, data = srv, + ties = "average") Asymptotic Two-Sample Logrank Test data: Surv(time, event) by treatment (1, 2) stratified by gender Z = -1.5446, p-value = 0.1224 alternative hypothesis: true theta is not equal to 1 > > pvalue(logrank_test(Surv(time, event) ~ treatment | gender, data = srv, + ties = "average", + distribution = approximate(nresample = 10000))) [1] 0.1582 99 percent confidence interval: 0.1489148 0.1678126 > > > ### StatXact 6 manual, page 448 > hypnosis <- data.frame(skin = c(23, 58, 11, 24, 34, + 23, 53, 10, 20, 40, + 23, 54, 22, 21, 22), + subject = factor(rep(c(1, 2, 3), rep(5, 3))), + treatment = factor(rep(1:5, 3), + labels = c("Fear", "Happiness", "Depression", + "Calmness", "Agitation"))) > > fta <- friedman_test(skin ~ treatment | subject, data = hypnosis) > > # test statistic, page 449 > stopifnot(isequal(round(statistic(fta), 3), 9.153)) > > # two-sided asymptotical p-value, page 449 > stopifnot(isequal(round(pvalue(fta), 4), 0.0574)) > > # two-sided approximated p-value > > ftMC <- friedman_test(skin ~ treatment | subject, data = hypnosis, + distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ftMC), "conf.int") > > stopifnot(pci[1] < 0.0268 & pci[2] > 0.0268) > > > ### StatXact 6 manual, page 458 > > analgesic_eff <- data.frame(response = factor( + c(0, 1, 1, + 0, 1, 1, + 1, 0, 1, + 0, 0, 0, + 0, 0, 1, + 0, 1, 1, + 1, 0, 1, + 0, 0, 1, + 0, 0, 0, + 0, 0, 1, + 0, 1, 0, + 0, 0, 1), labels = c("success", "failure")), + treatment = factor(rep(c("Placebo", "Aspirin", "NewDrug"), 12)), + subject = factor(rep(1:12, rep(3, 12)))) > > bta <- mh_test(response ~ treatment | subject, data = analgesic_eff) > > # asymptotic p-value, page 459 (no frame, see text!) > stopifnot(isequal(round(pvalue(bta), 3), 0.02)) > > # approximative p-value > btMC <- mh_test(response ~ treatment | subject, + data = analgesic_eff, distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(btMC), "conf.int") > stopifnot(pci[1] < 0.026 & pci[2] > 0.026) > > > ### StatXact 6 manual, page 467, Page test > cotton <- data.frame(strength = c(7.46, 7.17, 7.76, 8.14, 7.63, + 7.68, 7.57, 7.73, 8.15, 8.00, + 7.21, 7.80, 7.74, 7.87, 7.93), + potash = ordered(rep(c(144, 108, 72, 54, 36), 3), + levels = c(144, 108, 72, 54, 36)), + block = factor(rep(1:3, rep(5, 3)))) > > fta <- friedman_test(strength ~ potash | block, data = cotton) > > # test statistic, page 467 > stopifnot(isequal(round(statistic(fta), 3), 2.656)) > > # asymptotical p-value, page 467 > stopifnot(isequal(round(pvalue(fta), 4), 0.0079)) > > # approximate p-value > ftMC <- friedman_test(strength ~ potash | block, data = cotton, + distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(ftMC), "conf.int") > stopifnot(pci[1] < 0.005 & pci[2] > 0.005) > > class(cotton$potash) <- "factor" > > # approximate p-value > ftMC <- friedman_test(strength ~ potash | block, data = cotton, + distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(ftMC), "conf.int") > stopifnot(pci[1] < 0.0376 & pci[2] > 0.0376) > > > ### StatXact 6 manual, page 486 > tox <- data.frame(response = c(0, 1, 8, 10, + 0, 0, 3, 3, 8, + 5, 6, 7, 14, 14, + 1, 1, 6, 7, 7, 7, 8, 8, 10, + 7, 10, 11, 12, 13), + drug = factor(paste("Drug", c(rep(1, 4), + rep(2, 5), + rep(3, 5), + rep(4, 9), + rep(5, 5))))) > > > mta <- independence_test(response ~ drug, data = tox, + ytrafo = function(data) + trafo(data, numeric_trafo = median_trafo), + teststat = "quad") > > # StatXact reports results of scaled Pearson statistic! > a <- factor(tox$response <= 7) > mta <- chisq_test(a ~ drug, data = tox) > > # test statistic, page 487 > stopifnot(isequal(round(statistic(mta), 3), 4.317)) > > # asymptotic p-value, page 487 > stopifnot(isequal(round(pvalue(mta), 4), 0.3648)) > > mtMC <- independence_test(response ~ drug, data = tox, + ytrafo = function(data) + trafo(data, numeric_trafo = median_trafo), + teststat = "quad", + distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(mtMC), "conf.int") > pvalue(mtMC) [1] 0.4336 99 percent confidence interval: 0.4208191 0.4464444 > > stopifnot(pci[1] < 0.4289 & pci[2] > 0.4289) > > > kta <- kruskal_test(response ~ drug, data = tox) > > # test statistic, page 493 > stopifnot(isequal(round(statistic(kta), 3), 9.415)) > > # asymptotic p-value, page 493 > stopifnot(isequal(round(pvalue(kta), 4), 0.0515)) > > pta <- oneway_test(response ~ drug, data = tox, teststat = "quad") > > # test statistic, page 508 > stopifnot(isequal(round(statistic(pta), 2), 10.98)) > > # asymptotic p-value, page 508 > stopifnot(isequal(round(pvalue(pta), 4), 0.0268)) > > > ### StatXact 6 manual, page 509 > brain <- data.frame(time = c(4, 5, 9, 12, 20, 25, 30, + 1, 4, 9, 12, 15, 23, 30, + 3, 7, 14, 20, 27, 30, 32, 50, + 5, 15, 20, 31, 39, 47, 55, 67), + event = c(1, 1, 1, 1, 0, 1, 0, + 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 0, 0, + 1, 1, 1, 1, 1, 1, 0,0), + treatment = factor(paste("trt", + c(rep(1, 7), rep(2, 7), + rep(3, 8), rep(3, 8))))) > > lta <- logrank_test(Surv(time, event) ~ treatment, data = brain, + ties.method = "average-scores") > > # test statistic, page 516 > isequal(round(statistic(lta), 3), 5.012) [1] 4.298 [1] 5.012 [1] FALSE > > # asymptotic p-value, page 516 > isequal(round(pvalue(lta), 4), 0.1709) [1] 0.1166 [1] 0.1709 [1] FALSE > > > ### StatXact 6 manual, page 519 > oring <- data.frame(temp = c(66, 67, 67, 67, 68, 68, 70, 70, 72, + 73, 75, 76, 76, 78, 79, 80, 81, + 57, 58, 63, 70, 70, + 75, + 53), + incidents = ordered(c(rep("None", 17), rep("One", 5), + rep("Two", 1), rep("Three", 1)), + levels = c("None", "One", "Two", "Three"))) > > pta <- oneway_test(temp ~ incidents, data = oring) > > # test statistic, page 524 > stopifnot(isequal(round(statistic(pta), 3), -2.698)) > > stopifnot(isequal(round(pvalue(pta), 4), 0.007)) > > > brain$treatment <- ordered(brain$treatment) > > lta <- logrank_test(Surv(time, event) ~ treatment, data = brain, + ties.method = "average-scores") > > # test statistic, page 536 > isequal(round(statistic(lta), 3), 1.773) [1] 1.445 [1] 1.773 [1] FALSE > > # asymptotic p-value, page 536 > isequal(round(pvalue(lta), 4), 0.0762) [1] 0.1485 [1] 0.0762 [1] FALSE > > > ### StatXact 6 manual, ... > a <- as.table(matrix(c(5, 5, 9, 1), nrow = 2, + dimnames = list(response = c("Response", "No Response"), + drug = c("A", "B")))) > > # exact p-value, page 671 > stopifnot(isequal(round(fisher.test(a)$p.value, 4), 0.1409)) > > cta <- chisq_test(a) > > # test statistic, page 673 > stopifnot(isequal(round(statistic(cta), 2), 3.81)) > > # asymptotic p-value, page 673 > stopifnot(isequal(round(pvalue(cta), 3), 0.051)) > > ctMC <- chisq_test(a, distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ctMC), "conf.int") > > stopifnot(pci[1] < 0.1409 & pci[2] > 0.1409) > > ctMC <- cmh_test(a, distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ctMC), "conf.int") > > stopifnot(pci[1] < 0.1409 & pci[2] > 0.1409) > > > ### StatXact 6 manual, page 793 > csom <- as.table(matrix(c(17066, 48, 14464, 38, 788, 5, 126, 1, 37, 1), nrow = 2, + dimnames = list(Malformation = c("Absent", "Present"), + Alcohol = c("0", "<1", "1-2", "3-5", ">=6")))) > > lta <- lbl_test(csom, scores = list(Alcohol = 0:4)) > > # test statistic, page 796 > stopifnot(isequal(round(statistic(lta), 3), -1.352)) > > # asymptotic p-value, page 796 > stopifnot(isequal(round(pvalue(lta), 4), 0.1764)) > stopifnot(isequal(round(pvalue(lbl_test(csom)), 4), + round(prop.trend.test(csom[2,], colSums(csom))$p.value, 4))) > > ltMC <- lbl_test(csom, distribution = approximate(nresample = 100)) > > pci <- attr(pvalue(ltMC), "conf.int") > stopifnot(pci[1] < 0.179 & pci[2] > 0.179) > > lta <- lbl_test(csom, scores = list(Alcohol = c(0, 0.5, 1.5, 4, 7))) > > # test statistic, page 807 > stopifnot(isequal(round(statistic(lta), 3), -2.563)) > > # asymptotic p-value, page 807 > stopifnot(isequal(round(pvalue(lta), 4), 0.0104)) > stopifnot(isequal(round(pvalue(lbl_test(csom, + scores = list(Alcohol = c(0, 0.5, 1.5, 4, 7)))), 4), + round(prop.trend.test(csom[2,], colSums(csom), + score = c(0, 0.5, 1.5, 4, 7))$p.value, 4))) > > > ### StatXact manual 6, page 797 > load("LANCET.rda") > > lta <- lbl_test(lancet) > > # test statistic, page 799 > stopifnot(isequal(round(statistic(lta), 3), 4.335)) > > # asymptotic p-value, page 799 > stopifnot(isequal(round(pvalue(lta), 4), 0)) > > ltMC <- lbl_test(lancet, distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(ltMC), "conf.int") > stopifnot(pci[1] < 0.0000045 & pci[2] > 0.0000045) > > > ### StatXact manual, page 808 > tumor <- data.frame(number = c( 0, 0, 0, 1, + 7, 10, 6, 8, + 0, 1, 0, 1, + 11, 9, 13, 14, + 1, 1, 1, 2, + 29, 26, 28, 20), + tumor = factor(rep(c(rep("Present", 4), rep("Absent", 4)), 3), + levels = c("Present", "Absent")), + dose = ordered(rep(paste(c(0, 1, 5, 50), "units"), 6), + levels = paste(c(0, 1, 5, 50), "units")), + stratum = factor(rep(paste("Stratum", 3:5), rep(8, 3)))) > > lta <- lbl_test(tumor ~ dose | stratum, data = tumor, weights = ~ number, + scores = list(dose = c(0, 1, 5, 50))) > > # test statistic, page 812 > stopifnot(isequal(round(statistic(lta), 3), 1.739)) > > # asymptotic p-value, page 812 > stopifnot(isequal(round(pvalue(lta), 3), 0.082)) > > lta <- lbl_test(xtabs(number ~ dose + tumor + stratum, data = tumor), + scores = list(dose = c(0, 1, 5, 50))) > > # test statistic, page 812 > stopifnot(isequal(round(statistic(lta), 3), 1.739)) > > # asymptotic p-value, page 812 > stopifnot(isequal(round(pvalue(lta), 3), 0.082)) > > ltMC <- lbl_test(xtabs(number ~ dose + tumor + stratum, data = tumor), + scores = list(dose = c(0, 1, 5, 50)), + distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(ltMC), "conf.int") > stopifnot(pci[1] < 0.0769 & pci[2] > 0.0769) > > > ### StatXact 6 manual, page 832 > endo <- as.table(matrix(c(6, 9, 9, 12, 2, 4, 2, 1, 3, 2, 3, 2, 1, 1, 1, 1), + nrow = 4, dimnames = + list(Cases = c("0", "0.1-0.299", "0.3-0.625", ">0.625"), + Controls = c("0", "0.1-0.299", "0.3-0.625", ">0.625")))) > > bta <- mh_test(endo, scores = list(response = c(0, 0.2, 0.5125, 0.7))) > > # test statistic, page 837 > stopifnot(isequal(round(statistic(bta), 3), 3.735)) > > # asymptotic p-value, page 837 > stopifnot(isequal(round(pvalue(bta), 4), 0.0002)) > > btMC <- mh_test(endo, scores = list(response = c(0, 0.2, 0.5125, 0.7)), + distribution = approximate(nresample = 10000)) > > print(pvalue(btMC)) [1] <1e-04 99 percent confidence interval: 0.0000000000 0.0005296914 > pci <- attr(pvalue(btMC), "conf.int") > (pci[1] < 0.0001 & pci[2] > 0.0001) [1] TRUE > > > ### StatXact 6 manual, page 944 > oral_lesions <- as.table(matrix(c(0, 8, 0, 0, 0, 0, 0, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 0, 0, + 0, 8, 0, 0, 0, 0, 0, 1, 1), ncol = 3, + dimnames = list(site = paste("site", 1:9), + region = c("Kerala", "Gujarat", "Andhra")))) > > cta <- chisq_test(oral_lesions) > > # test statistic, page 947 > stopifnot(isequal(round(statistic(cta), 1), 22.1)) > > # asymptotic p-value, page 947 > stopifnot(isequal(round(pvalue(cta), 2), 0.14)) > > # approximate p-value, page 947 > ctMC <- cmh_test(oral_lesions, distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(ctMC), "conf.int") > stopifnot(pci[1] < 0.0269 & pci[2] > 0.0269) > > ### StatXact 6 manual, 984 > dr <- as.table(matrix(c(100, 18, 50, 50, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1), nrow = 4, + dimnames = list(dose = paste((1:4)*100, "mg"), + tox = c("mild", "moderate", "severe", "death")))) > > lta <- lbl_test(dr, scores = list(dose = (1:4)*100, tox = 1:4)) > > # teststatistic, page 993 > stopifnot(isequal(round(statistic(lta), 3), 1.807)) > > # asymptotical p-value, page 993 > stopifnot(isequal(round(pvalue(lta), 4), 0.0708)) > > ltMC <- lbl_test(dr, distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(ltMC), "conf.int") > > stopifnot(pci[1] < 0.0792 & pci[2] > 0.0792) > > > lta <- lbl_test(dr, scores = list(tox = c(1, 3, 9, 27))) > > # teststatistic, page 993 > stopifnot(isequal(round(statistic(lta), 3), 1.734)) > > # asymptotical p-value, page 993 > stopifnot(isequal(round(pvalue(lta), 4), 0.0828)) > > ltMC <- lbl_test(dr, scores = list(tox = c(1, 3, 9, 27)), + distribution = approximate(nresample = 10000)) > > pci <- attr(pvalue(ltMC), "conf.int") > stopifnot(pci[1] < 0.078 & pci[2] > 0.078) > > > ### StatXact 6 manual, page 1012 > load("army.rda") > > cta <- cmh_test(army) > > # teststatistic, page 1014 > stopifnot(isequal(round(statistic(cta), 2), 25.18)) > > # asymptotical p-value, page 1014 > stopifnot(isequal(round(pvalue(cta), 6), 0.002774)) > > > ### > data(jobsatisfaction) > > cta <- cmh_test(jobsatisfaction) > > # teststatistic, page 1017 > stopifnot(isequal(round(statistic(cta), 1), 10.2)) > > # asymptotical p-value, page 1017 > stopifnot(isequal(round(pvalue(cta), 4), 0.3345)) > > # only Job.Satisfaction ordered > cta <- cmh_test(jobsatisfaction, + scores = list(Job.Satisfaction = c(1, 3, 4, 5))) > > # teststatistic, page 1018, is _wrong_ (but StatXact itself > # computes the correct result) > # (isequal(round(statistic(cta), 3), 9.226)) > # Agresti, 2002, Table 7.12, page 297 > stopifnot(isequal(round(statistic(cta), 4), 9.0342)) > > # asymptotical p-value, page 1018, is _wrong_ (but StatXact itself > # computes the correct result) > # (isequal(round(pvalue(cta), 4), 0.02643)) > # Agresti, 2002, Table 7.12, page 297 > stopifnot(isequal(round(pvalue(cta), 4), 0.0288)) > > > lta <- lbl_test(jobsatisfaction, + scores = list(Job.Satisfaction = c(1, 3, 4, 5), + Income = c(3, 10, 20, 35))) > > # teststatistic, page 1020 > stopifnot(isequal(round(statistic(lta), 3), 2.481)) > > # asymptotical p-value, page 1020 > stopifnot(isequal(round(pvalue(lta), 5), 0.01309)) > > lta <- cmh_test(jobsatisfaction, + scores = list(Job.Satisfaction = c(1, 3, 4, 5), + Income = c(3, 10, 20, 35))) > > # teststatistic, page 1020 > stopifnot(isequal(round(statistic(lta), 3), 2.481)) > > # asymptotical p-value, page 1020 > stopifnot(isequal(round(pvalue(lta), 5), 0.01309)) > > > ### --------------------------------------------------------- ### > > # some additional checks (always add new tests at the end because of the RNG's) > > lta <- logrank_test(Surv(time, event) ~ treatment | gender, data = srv) > stopifnot(isequal(round(pvalue(lta), 4), 0.1224)) > > ### example from Callaert (2003), AmStat 57, 214-217 > exdata <- data.frame(time = c(1, 1, 5, 6, 6, 6, 6, 2, 2, 2, 3, 4, 4, 5, 5), + event = rep(TRUE, 15), + group = factor(c(rep(0, 7), rep(1, 8)))) > p <- pvalue(logrank_test(Surv(time, event) ~ group, data = exdata, + distribution = exact())) > stopifnot(isequal(round(p, 4), 0.0505)) > p <- pvalue(logrank_test(Surv(time, event) ~ group, data = exdata, + distribution = exact(), ties = "average")) > stopifnot(isequal(round(p, 4), 0.0468)) > > > ### symmetry problem, page 288 > load("AIDS.rda") > > tmp <- data.frame(y = c(AIDS$post, AIDS$pre), + x = gl(2, nrow(AIDS)), + block = factor(rep(1:nrow(AIDS), 2))) > > wsa <- wilcoxsign_test(post ~ pre, data = AIDS, distribution = "asymptotic", + zero.method = "Wilcoxon") > wsa2 <- wilcoxsign_test(y ~ x | block, data = tmp, distribution = "asymptotic", + zero.method = "Wilcoxon") > > stopifnot(all.equal(statistic(wsa), statistic(wsa2))) > > # asymptotic p-value, page 290 > stopifnot(isequal(round(statistic(wsa), 3), -2.896)) > > # asymptotic p-value, page 290 > stopifnot(isequal(round(pvalue(wsa), 4), 0.0038)) > > wsa <- wilcoxsign_test(post ~ pre, data = AIDS, distribution = "asymptotic", alternative = "less", + zero.method = "Wilcoxon") > > # asymptotic p-value, page 290 > stopifnot(isequal(round(statistic(wsa), 3), -2.896)) > > # asymptotic p-value, page 290 > stopifnot(isequal(round(pvalue(wsa), 4), 0.0019)) > > wse <- wilcoxsign_test(post ~ pre, data = AIDS, distribution = "exact", + zero.method = "Wilcoxon") > > # exact p-value, page 290 > stopifnot(isequal(round(pvalue(wse), 4), 0.0021)) > > wsa <- wilcoxsign_test(post ~ pre, data = AIDS, distribution = "exact", alternative = "less", + zero.method = "Wilcoxon") > > # exact one-sided p-value, page 290 > stopifnot(isequal(round(pvalue(wsa), 4), 0.0011)) > > ### permutation test > diff <- AIDS$pre - AIDS$post > diff <- diff[abs(diff) > 0] > y <- as.vector(t(cbind(abs(diff) * (diff > 0), abs(diff) * (diff <= 0)))) > block <- gl(length(diff), 2) > x <- factor(rep(c("pos", "neg"), length(diff))) > > sta <- symmetry_test(y ~ x | block) > > # asymptotic p-value, page 300 > stopifnot(isequal(round(statistic(sta), 3), -1.707)) > stopifnot(isequal(round(pvalue(sta), 4), 0.0878)) > > # approximate p-value, page 300 > stMC <- symmetry_test(y ~ x | block, distribution = approximate(nresample = 10000)) > pci <- attr(pvalue(stMC), "conf.int") > stopifnot(pci[1] < 0.0021 & pci[2] > 0.0021) > > > ### sign test > sa <- sign_test(post ~ pre, data = AIDS) > sa2 <- sign_test(y ~ x | block, data = tmp) > > stopifnot(all.equal(statistic(sa), statistic(sa2))) > > # asymptotic p-value, page 107 SX9 manual > stopifnot(isequal(round(statistic(sa), 3), -3)) > > # asymptotic p-value, page 107 SX9 manual > stopifnot(isequal(round(pvalue(sa), 4), 0.0027)) > > sa <- sign_test(post ~ pre, data = AIDS, alternative = "less") > > # asymptotic p-value, page 107 SX9 manual > stopifnot(isequal(round(pvalue(sa), 4), 0.0013)) > > se <- sign_test(post ~ pre, data = AIDS, distribution = "exact") > > # exact p-value, page 107 SX9 manual > stopifnot(isequal(round(pvalue(se), 4), 0.0042)) > > sa <- sign_test(post ~ pre, data = AIDS, distribution = "exact", + alternative = "less") > > # exact one-sided p-value, page 107 SX9 manual > stopifnot(isequal(round(pvalue(sa), 4), 0.0021)) > > proc.time() user system elapsed 1.96 0.14 2.09
Generated by dwww version 1.15 on Sun Jun 16 15:11:25 CEST 2024.