R : Copyright 2004, The R Foundation for Statistical Computing
Version 2.0.1 (2004-11-15), ISBN 3-900051-07-0
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 a HTML browser interface to help.
Type 'q()' to quit R.
> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
>
> assign("cleanEx",
+ function(env = .GlobalEnv) {
+ rm(list = ls(envir = env, all.names = TRUE), envir = env)
+ RNGkind("default", "default")
+ set.seed(1)
+ options(warn = 1)
+ assign("T", delay(stop("T used instead of TRUE")),
+ pos = .CheckExEnv)
+ assign("F", delay(stop("F used instead of FALSE")),
+ pos = .CheckExEnv)
+ sch <- search()
+ newitems <- sch[! sch %in% .oldSearch]
+ for(item in rev(newitems))
+ eval(substitute(detach(item), list(item=item)))
+ missitems <- .oldSearch[! .oldSearch %in% sch]
+ if(length(missitems))
+ warning("items ", paste(missitems, collapse=", "),
+ " have been removed from the search path")
+ },
+ env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("coin-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> library('coin')
Loading required package: survival
Loading required package: splines
Loading required package: mvtnorm
>
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "ContingencyTests"
>
> ### * ContingencyTests
>
> flush(stderr()); flush(stdout())
>
> ### Name: ContingencyTests
> ### Title: Independence in General I x K x J Contingency Tables
> ### Aliases: chisq_test chisq_test.formula chisq_test.table
> ### chisq_test.IndependenceProblem cmh_test.formula cmh_test.table
> ### cmh_test.IndependenceProblem cmh_test lbl_test.formula lbl_test.table
> ### lbl_test.IndependenceProblem lbl_test
> ### Keywords: htest
>
> ### ** Examples
>
>
> ## Don't show:
> set.seed(290875)
> ## End Don't show
>
> data(jobsatisfaction, package = "coin")
>
> ### for females only
> chisq_test(as.table(jobsatisfaction[,,"Female"]),
+ distribution = "approx", B = 9999)
Approximative Pearson's Chi-Squared Test
data: Job.Satisfaction by
groups <5000, 5000-15000, 15000-25000, >25000
T = 6.8203, p-value = 0.6891
>
> ### both Income and Job.Satisfaction unordered
> cmh_test(jobsatisfaction)
Asymptotical Generalised Cochran-Mantel-Haenszel Test
data: Job.Satisfaction by
groups <5000, 5000-15000, 15000-25000, >25000
stratified by Gender
T = 10.2001, df = 9, p-value = 0.3345
>
> ### both Income and Job.Satisfaction ordered, default scores
> lbl_test(jobsatisfaction)
Asymptotical Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
groups <5000 < 5000-15000 < 15000-25000 < >25000
stratified by Gender
T = 6.6235, df = 1, p-value = 0.01006
>
> ### both Income and Job.Satisfaction ordered, alternative scores
> lbl_test(jobsatisfaction, scores = list(Job.Satisfaction = c(1, 3, 4, 5),
+ Income = c(3, 10, 20, 35)))
Asymptotical Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
groups <5000 < 5000-15000 < 15000-25000 < >25000
stratified by Gender
T = 6.1563, df = 1, p-value = 0.01309
>
> ### the same, null distribution approximated
> cmh_test(jobsatisfaction, scores = list(Job.Satisfaction = c(1, 3, 4, 5),
+ Income = c(3, 10, 20, 35)),
+ distribution = "approx", B = 10000)
Approximative Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
groups <5000 < 5000-15000 < 15000-25000 < >25000
stratified by Gender
T = 6.1563, p-value = 0.0123
>
>
>
>
> cleanEx(); ..nameEx <- "IndependenceTest"
>
> ### * IndependenceTest
>
> flush(stderr()); flush(stdout())
>
> ### Name: IndependenceTest
> ### Title: General Independence Tests
> ### Aliases: independence_test independence_test.formula
> ### independence_test.IndependenceProblem independence_test.table
> ### Keywords: htest
>
> ### ** Examples
>
>
> data(asat, package = "coin")
>
> ### independence of asat and group via normal scores test
> independence_test(asat ~ group, data = asat,
+
+ ### exact null distribution
+ distribution = "exact",
+
+ ### one-sided test
+ alternative = "greater",
+
+ ### apply normal scores to asat$asat
+ ytrafo = function(data) trafo(data, numeric_trafo = normal_trafo),
+
+ ### indicator matrix of 1st level of group
+ xtrafo = function(data) trafo(data, factor_trafo = function(x)
+ matrix(x == levels(x)[1], ncol = 1))
+ )
Exact General Independence Test
data: asat by groups Compound, Control
T = 1.4269, p-value = 0.07809
alternative hypothesis: greater
>
> ### same as
> normal_test(asat ~ group, data = asat, distribution = "exact",
+ alternative = "greater")
Exact Normal Quantile (van der Waerden) Test
data: asat by groups Compound, Control
T = 1.4269, p-value = 0.07809
alternative hypothesis: true mu is greater than 0
>
>
>
>
> cleanEx(); ..nameEx <- "LocationTests"
>
> ### * LocationTests
>
> flush(stderr()); flush(stdout())
>
> ### Name: LocationTests
> ### Title: Independent Two- and K-Sample Location Tests
> ### Aliases: oneway_test oneway_test.formula
> ### oneway_test.IndependenceProblem wilcox_test.formula
> ### wilcox_test.IndependenceProblem wilcox_test normal_test.formula
> ### normal_test.IndependenceProblem normal_test median_test.formula
> ### median_test.IndependenceProblem median_test kruskal_test.formula
> ### kruskal_test.IndependenceProblem kruskal_test
> ### Keywords: htest
>
> ### ** Examples
>
>
> ### Tritiated Water Diffusion Across Human Chorioamnion
> ### Hollander & Wolfe (1999), Table 4.1, page 110
> water_transfer <- data.frame(
+ pd = c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46,
+ 1.15, 0.88, 0.90, 0.74, 1.21),
+ age = factor(c(rep("At term", 10), rep("12-26 Weeks", 5))))
>
> ### Wilcoxon-Mann-Whitney test, cf. Hollander & Wolfe (1999), page 111
> ### exact p-value and confidence interval for the difference in location
> ### (At term - 12-26 Weeks)
> wt <- wilcox_test(pd ~ age, data = water_transfer,
+ distribution = "exact", conf.int = TRUE)
> print(wt)
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: pd by groups 12-26 Weeks, At term
T = -1.2247, p-value = 0.2544
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
-0.76 0.15
sample estimates:
difference in location
-0.305
>
> ### extract observed Wilcoxon statistic, i.e, the sum of the
> ### ranks for age = "12-26 Weeks"
> statistic(wt, "linear")
[,1]
12-26 Weeks 30
>
> ### its expectation
> expectation(wt)
[1] 40
>
> ### and variance
> covariance(wt)
[,1]
[1,] 66.66667
>
> ### and, finally, the exact two-sided p-value
> pvalue(wt)
[1] 0.2544123
>
> ### Confidence interval for difference (12-26 Weeks - At term)
> wilcox_test(pd ~ age, data = water_transfer,
+ xtrafo = function(data)
+ trafo(data, factor_trafo = function(x) as.numeric(x == levels(x)[2])),
+ distribution = "exact", conf.int = TRUE)
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: pd by groups 12-26 Weeks, At term
T = 1.2247, p-value = 0.2544
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
-0.15 0.76
sample estimates:
difference in location
0.305
>
> ### Permutation test, asymptotic p-value
> oneway_test(pd ~ age, data = water_transfer)
Asymptotical 2-Sample Permutation Test
data: pd by groups 12-26 Weeks, At term
T = -1.5225, p-value = 0.1279
alternative hypothesis: true mu is not equal to 0
>
> ### approximate p-value (with 99% confidence interval)
> pvalue(oneway_test(pd ~ age, data = water_transfer,
+ distribution = "approx", B = 9999))
[1] 0.1311131
99 percent confidence interval:
0.1225462 0.1400337
> ### exact p-value
> pt <- oneway_test(pd ~ age, data = water_transfer, distribution = "exact")
> pvalue(pt)
[1] 0.1318681
>
> ### plot density and distribution of the standardised
> ### test statistic
> layout(matrix(1:2, nrow = 2))
> s <- support(pt)
> d <- sapply(s, function(x) dperm(pt, x))
> p <- sapply(s, function(x) pperm(pt, x))
> plot(s, d, type = "S", xlab = "Teststatistic", ylab = "Density")
> plot(s, p, type = "S", xlab = "Teststatistic", ylab = "Cumm. Probability")
>
> ### Length of YOY Gizzard Shad from Kokosing Lake, Ohio,
> ### sampled in Summer 1984, Hollander & Wolfe, Table 6.3, page 200
> 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))))
>
> ### Kruskal-Wallis test, approximate exact p-value
> kw <- kruskal_test(length ~ site, data = YOY,
+ distribution = "approx", B = 9999)
> kw
Approximative Kruskal-Wallis Test
data: length by groups I, II, III, IV
T = 22.8524, p-value < 2.2e-16
> pvalue(kw)
[1] 0
99 percent confidence interval:
0.0000000000 0.0005297444
>
> ### Nemenyi-Damico-Wolfe-Dunn test (joint ranking)
> ### Hollander & Wolfe, 1999, page 244
> ### (where Steel-Dwass results are given)
> if (require(multcomp)) {
+
+ NDWD <- oneway_test(length ~ site, data = YOY,
+ ytrafo = function(data) trafo(data, numeric_trafo = rank),
+ xtrafo = function(data) trafo(data, factor_trafo = function(x)
+ model.matrix(~x - 1) %*% t(contrMat(table(x), "Tukey"))),
+ teststat = "maxtype", distribution = "approx", B = 90000)
+
+ ### global p-value
+ print(pvalue(NDWD))
+
+ ### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
+ print(pvalue(NDWD, adjusted = TRUE))
+ }
Loading required package: multcomp
[1] 0.0004888889
99 percent confidence interval:
0.0003199269 0.0007126921
[,1]
II-I 0.9469666667
III-I 0.0089333333
IV-I 0.0070666667
III-II 0.0006777778
IV-II 0.0004777778
IV-III 0.9998777778
>
>
>
>
> cleanEx(); ..nameEx <- "MarginalHomogenityTest"
>
> ### * MarginalHomogenityTest
>
> flush(stderr()); flush(stdout())
>
> ### Name: MarginalHomogenityTest
> ### Title: Marginal Homogenity Test
> ### Aliases: mh_test mh_test.table mh_test.formula mh_test.SymmetryProblem
> ### Keywords: htest
>
> ### ** Examples
>
>
> ### Opinions on Pre- and Extramarital Sex, Agresti (2002), page 421
> opinions <- c("always wrong", "almost always wrong",
+ "wrong only sometimes", "not wrong at all")
>
> PreExSex <- as.table(matrix(c(144, 33, 84, 126,
+ 2, 4, 14, 29,
+ 0, 2, 6, 25,
+ 0, 0, 1, 5), nrow = 4,
+ dimnames = list(PremaritalSex = opinions,
+ ExtramaritalSex = opinions)))
>
> ### treating response as nominal
> mh_test(PreExSex)
Asymptotical Marginal-Homogenity Test
data: response by
groups ExtramaritalSex, PremaritalSex
stratified by
T = 271.9222, df = 3, p-value < 2.2e-16
>
> ### and as ordinal
> mh_test(PreExSex, scores = list(response = 1:length(opinions)))
Asymptotical Marginal-Homogenity Test for Ordered Data
data: response (ordered) by
groups ExtramaritalSex, PremaritalSex
stratified by
T = 270.7402, df = 1, p-value < 2.2e-16
>
>
>
>
> cleanEx(); ..nameEx <- "MaxstatTest"
>
> ### * MaxstatTest
>
> flush(stderr()); flush(stdout())
>
> ### Name: MaxstatTest
> ### Title: Maximally Selected Statistics
> ### Aliases: maxstat_test maxstat_test.formula
> ### maxstat_test.IndependenceProblem
> ### Keywords: htest
>
> ### ** Examples
>
>
> data(treepipit, package = "coin")
>
> maxstat_test(counts ~ coverstorey, data = treepipit)
Asymptotical Maxstat Test
data: counts by coverstorey
T = 4.3139, p-value = 0.0001951
sample estimates:
cutpoint in coverstorey
40
>
>
>
>
> cleanEx(); ..nameEx <- "PermutationDistribution"
>
> ### * PermutationDistribution
>
> flush(stderr()); flush(stdout())
>
> ### Name: PermutationDistribution
> ### Title: Permutation Distribution of Conditional Independence Tests
> ### Aliases: pperm pperm-methods pperm,NullDistribution-method
> ### pperm,IndependenceTest-method pperm,ScalarIndependenceTest-method
> ### pperm,MaxTypeIndependenceTest-method
> ### pperm,QuadTypeIndependenceTest-method qperm qperm-methods
> ### qperm,NullDistribution-method qperm,IndependenceTest-method
> ### qperm,ScalarIndependenceTest-method
> ### qperm,MaxTypeIndependenceTest-method
> ### qperm,QuadTypeIndependenceTest-method dperm dperm-methods
> ### dperm,NullDistribution-method dperm,IndependenceTest-method
> ### dperm,ScalarIndependenceTest-method
> ### dperm,MaxTypeIndependenceTest-method
> ### dperm,QuadTypeIndependenceTest-method support support-methods
> ### support,NullDistribution-method support,IndependenceTest-method
> ### support,ScalarIndependenceTest-method
> ### support,MaxTypeIndependenceTest-method
> ### support,QuadTypeIndependenceTest-method
> ### Keywords: methods htest
>
> ### ** Examples
>
>
> ### artificial 2-sample problem
> df <- data.frame(y = rnorm(20), x = gl(2, 10))
>
> ### Ansari-Bradley test
> at <- ansari_test(y ~ x, data = df, distribution = "exact")
>
> ### density of the exact distribution of the Ansari-Bradley statistic
> dens <- sapply(support(at), dperm, object = at)
>
> ### 95% quantile
> qperm(at, 0.95)
[1] 1.669331
>
> ### one-sided p-value
> pperm(at, statistic(at))
[1] 0.698635
>
>
>
>
> cleanEx(); ..nameEx <- "ScaleTests"
>
> ### * ScaleTests
>
> flush(stderr()); flush(stdout())
>
> ### Name: ScaleTests
> ### Title: Independent Two- and K-Sample Scale Tests
> ### Aliases: ansari_test ansari_test.formula
> ### ansari_test.IndependenceProblem fligner_test.formula
> ### fligner_test.IndependenceProblem fligner_test
> ### Keywords: htest
>
> ### ** Examples
>
>
> ### Serum Iron Determination Using Hyland Control Sera
> ### Hollander & Wolfe (1999), page 147
> sid <- data.frame(
+ serum = c(111, 107, 100, 99, 102, 106, 109, 108, 104, 99,
+ 101, 96, 97, 102, 107, 113, 116, 113, 110, 98,
+ 107, 108, 106, 98, 105, 103, 110, 105, 104,
+ 100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99),
+ method = factor(gl(2, 20), labels = c("Ramsay", "Jung-Parekh")))
>
> ### Ansari-Bradley test, asymptotical p-value
> ansari_test(serum ~ method, data = sid)
Asymptotical Ansari-Bradley Test
data: serum by groups Ramsay, Jung-Parekh
T = -1.3363, p-value = 0.1815
alternative hypothesis: true mu is not equal to 1
>
> ### exact p-value
> ansari_test(serum ~ method, data = sid, distribution = "exact")
Exact Ansari-Bradley Test
data: serum by groups Ramsay, Jung-Parekh
T = -1.3363, p-value = 0.1881
alternative hypothesis: true mu is not equal to 1
>
> ### Platelet Counts of Newborn Infants
> ### Hollander & Wolfe, Table 5.4, page 171
> platalet_counts <- data.frame(
+ counts = c(120, 124, 215, 90, 67, 95, 190, 180, 135, 399,
+ 12, 20, 112, 32, 60, 40),
+ treatment = factor(c(rep("Prednisone", 10), rep("Control", 6))))
>
> ### Lepage test, Hollander & Wolfe, page 172
> lt <- independence_test(counts ~ treatment, data = platalet_counts,
+ ytrafo = function(data) trafo(data, numeric_trafo = function(x)
+ cbind(rank(x), ansari_trafo(x))),
+ teststat = "quadtype", distribution = "approx", B = 9999)
>
> lt
Approximative General Independence Test
data: counts by groups Control, Prednisone
T = 9.3384, p-value = 0.0042
>
> ### where did the rejection come from? Use maximum statistic
> ### instead of a quadratic form
> ltmax <- independence_test(counts ~ treatment, data = platalet_counts,
+ ytrafo = function(data) trafo(data, numeric_trafo = function(x)
+ matrix(c(rank(x), ansari_trafo(x)), ncol = 2,
+ dimnames = list(1:length(x), c("Location", "Scale")))),
+ teststat = "maxtype", distribution = "approx", B = 9999)
>
> ### points to a difference in location
> pvalue(ltmax, adjusted = TRUE)
Location Scale
Control 0.00250025 0.5948595
>
> ### Funny: We could have used a simple Bonferroni procedure
> ### since the correlation between the Wilcoxon and Ansari-Bradley
> ### test statistics is zero
> covariance(ltmax)
[,1] [,2]
[1,] 85 0
[2,] 0 21
>
>
>
>
> cleanEx(); ..nameEx <- "SpearmanTest"
>
> ### * SpearmanTest
>
> flush(stderr()); flush(stdout())
>
> ### Name: SpearmanTest
> ### Title: Spearman's Test on Independence
> ### Aliases: spearman_test spearman_test.formula
> ### spearman_test.IndependenceProblem
> ### Keywords: htest
>
> ### ** Examples
>
>
> spearman_test(CONT ~ INTG, data = USJudgeRatings)
Asymptotical Spearman Correlation Test
data: CONT by INTG
T = -1.1437, p-value = 0.2527
alternative hypothesis: true mu is not equal to 0
>
>
>
>
> cleanEx(); ..nameEx <- "SurvTest"
>
> ### * SurvTest
>
> flush(stderr()); flush(stdout())
>
> ### Name: SurvTest
> ### Title: Independent Two- and K-Sample Tests for Censored Data
> ### Aliases: surv_test surv_test.formula surv_test.IndependenceProblem
> ### Keywords: htest
>
> ### ** Examples
>
>
> data(ocarcinoma, package = "coin")
>
> surv_test(Surv(time, event) ~ stadium, data = ocarcinoma,
+ distribution = "exact")
Exact Logrank Test
data: y by groups II, IIA
T = -2.3341, p-value = 0.01837
alternative hypothesis: two.sided
>
>
>
>
> cleanEx(); ..nameEx <- "SymmetryTests"
>
> ### * SymmetryTests
>
> flush(stderr()); flush(stdout())
>
> ### Name: SymmetryTests
> ### Title: Symmetry Tests
> ### Aliases: friedman_test friedman_test.formula
> ### friedman_test.SymmetryProblem wilcoxsign_test.formula
> ### wilcoxsign_test.IndependenceProblem wilcoxsign_test
> ### Keywords: htest
>
> ### ** Examples
>
>
> ### Hollander & Wolfe (1999), Table 7.1, page 274
> ### Comparison of three methods ("round out", "narrow angle", and
> ### "wide angle") for rounding first base.
> RoundingTimes <- data.frame(
+ times = c(5.40, 5.50, 5.55,
+ 5.85, 5.70, 5.75,
+ 5.20, 5.60, 5.50,
+ 5.55, 5.50, 5.40,
+ 5.90, 5.85, 5.70,
+ 5.45, 5.55, 5.60,
+ 5.40, 5.40, 5.35,
+ 5.45, 5.50, 5.35,
+ 5.25, 5.15, 5.00,
+ 5.85, 5.80, 5.70,
+ 5.25, 5.20, 5.10,
+ 5.65, 5.55, 5.45,
+ 5.60, 5.35, 5.45,
+ 5.05, 5.00, 4.95,
+ 5.50, 5.50, 5.40,
+ 5.45, 5.55, 5.50,
+ 5.55, 5.55, 5.35,
+ 5.45, 5.50, 5.55,
+ 5.50, 5.45, 5.25,
+ 5.65, 5.60, 5.40,
+ 5.70, 5.65, 5.55,
+ 6.30, 6.30, 6.25),
+ methods = factor(rep(c("Round Out", "Narrow Angle", "Wide Angle"), 22)),
+ block = factor(rep(1:22, rep(3, 22))))
>
> ### classical global test
> friedman_test(times ~ methods | block, data = RoundingTimes)
Asymptotical Friedman Test
data: times by
groups Narrow Angle, Round Out, Wide Angle
stratified by block
T = 11.1429, df = 2, p-value = 0.003805
>
> ### parallel coordinates plot
> matplot(t(matrix(RoundingTimes$times, ncol = 3, byrow = TRUE)),
+ type = "l", col = 1, lty = 1, axes = FALSE, ylab = "Time",
+ xlim = c(0.5, 3.5))
> axis(1, at = 1:3, labels = levels(RoundingTimes$methods))
> axis(2)
>
> ### where do the differences come from?
> ### Wilcoxon-Nemenyi-McDonald-Thompson test
> ### Hollander & Wolfe, page 295
> if (require(multcomp)) {
+
+ ### all pairwise comparisons
+ rtt <- symmetry_test(times ~ methods | block, data = RoundingTimes,
+ teststat = "maxtype",
+ xtrafo = function(data)
+ trafo(data, factor_trafo = function(x)
+ model.matrix(~ x - 1)
+ ),
+ ytrafo = function(data)
+ trafo(data, numeric_trafo = rank, block = RoundingTimes$block)
+ )
+
+ ### a global test, again
+ print(pvalue(rtt))
+
+ ### simultaneous P-values for all pair comparisons
+ ### Wide Angle vs. Round Out differ (Hollander and Wolfe, 1999, page 296)
+ print(pvalue(rtt, adjusted = TRUE))
+ }
Loading required package: multcomp
[1] 0.00387132
percent confidence interval:
0.003609045 0.004133595
[,1]
xNarrow Angle 0.701908847
xRound Out 0.042678917
xWide Angle 0.003811824
>
> ### Strength Index of Cotton, Hollander & Wolfe (1999), Table 7.5, page 286
> sc <- data.frame(block = factor(c(rep(1, 5), rep(2, 5), rep(3, 5))),
+ potash = ordered(rep(c(144, 108, 72, 54, 36), 3)),
+ 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))
>
> ### Page test
> ft <- friedman_test(strength ~ potash | block, data = sc)
> ft
Asymptotical Page Test
data: strength by
groups 36 < 54 < 72 < 108 < 144
stratified by block
T = 7.0533, df = 1, p-value = 0.007912
>
> ### one-sided p-value
> 1 - pnorm(sqrt(statistic(ft)))
[1] 0.003955894
>
> ### approximate null distribution via Monte-Carlo
> pvalue(friedman_test(strength ~ potash | block, data = sc,
+ distribution = "approx", B = 9999))
[1] 0.00470047
99 percent confidence interval:
0.003124471 0.006765311
>
>
>
>
> cleanEx(); ..nameEx <- "Transformations"
>
> ### * Transformations
>
> flush(stderr()); flush(stdout())
>
> ### Name: Transformations
> ### Title: Functions for Data Transformations and Score Computations
> ### Aliases: trafo id_trafo ansari_trafo fligner_trafo normal_trafo
> ### median_trafo consal_trafo maxstat_trafo logrank_trafo f_trafo
> ### Keywords: manip
>
> ### ** Examples
>
>
> ### dummy matrices, 2-sample problem (only one column)
> f_trafo(y <- gl(2, 5))
1
1 1
2 1
3 1
4 1
5 1
6 0
7 0
8 0
9 0
10 0
>
> ### K-sample problem (K columns)
> f_trafo(y <- gl(5, 2))
1 2 3 4 5
1 1 0 0 0 0
2 1 0 0 0 0
3 0 1 0 0 0
4 0 1 0 0 0
5 0 0 1 0 0
6 0 0 1 0 0
7 0 0 0 1 0
8 0 0 0 1 0
9 0 0 0 0 1
10 0 0 0 0 1
>
> ### normal scores
> normal_trafo(x <- rnorm(10))
[1] -0.6045853 -0.1141853 -1.3351777 1.3351777 0.1141853 -0.9084579
[7] 0.3487557 0.9084579 0.6045853 -0.3487557
>
> ### and now together
> trafo(data.frame(x = x, y = y), numeric_trafo = normal_trafo)
1 2 3 4 5
1 -0.6045853 1 0 0 0 0
2 -0.1141853 1 0 0 0 0
3 -1.3351777 0 1 0 0 0
4 1.3351777 0 1 0 0 0
5 0.1141853 0 0 1 0 0
6 -0.9084579 0 0 1 0 0
7 0.3487557 0 0 0 1 0
8 0.9084579 0 0 0 1 0
9 0.6045853 0 0 0 0 1
10 -0.3487557 0 0 0 0 1
attr(,"assign")
[1] 1 2 2 2 2 2
>
> ### maximally selected statistics
> maxstat_trafo(rnorm(10))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0 0 0 0 0 0 0 0
[2,] 1 0 1 0 0 1 1 1
[3,] 1 1 1 1 1 1 1 1
[4,] 1 1 1 1 1 1 1 1
[5,] 0 0 1 0 0 0 0 0
[6,] 1 0 1 1 1 1 1 1
[7,] 1 0 1 0 1 1 1 1
[8,] 0 0 1 0 0 1 0 0
[9,] 0 0 1 0 0 1 1 0
[10,] 0 0 1 0 0 1 1 1
>
> ### apply transformation blockwise (e.g. for Friedman test)
> trafo(data.frame(y = 1:20), numeric_trafo = rank, block = gl(4, 5))
[,1]
[1,] 1
[2,] 2
[3,] 3
[4,] 4
[5,] 5
[6,] 1
[7,] 2
[8,] 3
[9,] 4
[10,] 5
[11,] 1
[12,] 2
[13,] 3
[14,] 4
[15,] 5
[16,] 1
[17,] 2
[18,] 3
[19,] 4
[20,] 5
attr(,"assign")
[1] 1
>
>
>
>
> cleanEx(); ..nameEx <- "asat"
>
> ### * asat
>
> flush(stderr()); flush(stdout())
>
> ### Name: asat
> ### Title: Toxicological Study on Female Wistar Rats
> ### Aliases: asat
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(asat, package = "coin")
>
> ### proof-of-safety based on ratio of medians
> pos <- wilcox_test(I(log(asat)) ~ group, data = asat, alternative = "less",
+ conf.int = TRUE, distribution = "exact")
>
> ### one-sided confidence set. Safety cannot be concluded since the effect of
> ### the compound exceeds 20% of the control median
> exp(confint(pos)$conf.int)
[1] 0.000000 1.337778
attr(,"conf.level")
[1] 0.95
>
>
>
>
> cleanEx(); ..nameEx <- "expectation-methods"
>
> ### * expectation-methods
>
> flush(stderr()); flush(stdout())
>
> ### Name: expectation-methods
> ### Title: Extract the Expectation and Covariance of Linear Statistics
> ### Aliases: expectation expectation-methods
> ### expectation,IndependenceTest-method
> ### expectation,ScalarIndependenceTestStatistic-method
> ### expectation,MaxTypeIndependenceTestStatistic-method
> ### expectation,QuadTypeIndependenceTestStatistic-method covariance
> ### covariance-methods covariance,IndependenceTest-method
> ### covariance,ScalarIndependenceTestStatistic-method
> ### covariance,MaxTypeIndependenceTestStatistic-method
> ### covariance,QuadTypeIndependenceTestStatistic-method
> ### Keywords: methods
>
> ### ** Examples
>
>
> df <- data.frame(y = gl(3, 2), x = gl(3, 2)[sample(1:6)])
>
> ### Cochran-Mantel-Haenzel Test
> ct <- cmh_test(y ~ x, data = df)
>
> ### the linear statistic, i.e, the contingency table
> l <- statistic(ct, type = "linear")
> l
1 2 3
1 1 0 1
2 0 2 0
3 1 0 1
>
> ### expectation
> El <- expectation(ct)
> El
1 2 3
1 0.6666667 0.6666667 0.6666667
2 0.6666667 0.6666667 0.6666667
3 0.6666667 0.6666667 0.6666667
>
> ### covariance
> Vl <- covariance(ct)
> Vl
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.35555556 -0.17777778 -0.17777778 -0.17777778 0.08888889 0.08888889
[2,] -0.17777778 0.35555556 -0.17777778 0.08888889 -0.17777778 0.08888889
[3,] -0.17777778 -0.17777778 0.35555556 0.08888889 0.08888889 -0.17777778
[4,] -0.17777778 0.08888889 0.08888889 0.35555556 -0.17777778 -0.17777778
[5,] 0.08888889 -0.17777778 0.08888889 -0.17777778 0.35555556 -0.17777778
[6,] 0.08888889 0.08888889 -0.17777778 -0.17777778 -0.17777778 0.35555556
[7,] -0.17777778 0.08888889 0.08888889 -0.17777778 0.08888889 0.08888889
[8,] 0.08888889 -0.17777778 0.08888889 0.08888889 -0.17777778 0.08888889
[9,] 0.08888889 0.08888889 -0.17777778 0.08888889 0.08888889 -0.17777778
[,7] [,8] [,9]
[1,] -0.17777778 0.08888889 0.08888889
[2,] 0.08888889 -0.17777778 0.08888889
[3,] 0.08888889 0.08888889 -0.17777778
[4,] -0.17777778 0.08888889 0.08888889
[5,] 0.08888889 -0.17777778 0.08888889
[6,] 0.08888889 0.08888889 -0.17777778
[7,] 0.35555556 -0.17777778 -0.17777778
[8,] -0.17777778 0.35555556 -0.17777778
[9,] -0.17777778 -0.17777778 0.35555556
>
> ### the standardized contingency table (hard way)
> (l - El) / matrix(sqrt(diag(Vl)), ncol = nrow(El))
1 2 3
1 0.559017 -1.118034 0.559017
2 -1.118034 2.236068 -1.118034
3 0.559017 -1.118034 0.559017
>
> ### easy way
> statistic(ct, type = "standardized")
1 2 3
1 0.559017 -1.118034 0.559017
2 -1.118034 2.236068 -1.118034
3 0.559017 -1.118034 0.559017
>
>
>
>
> cleanEx(); ..nameEx <- "glioma"
>
> ### * glioma
>
> flush(stderr()); flush(stdout())
>
> ### Name: glioma
> ### Title: Malignant Glioma Pilot Study
> ### Aliases: glioma
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(glioma, package = "coin")
>
> par(mfrow=c(1,2))
>
> ### Grade III glioma
> g3 <- subset(glioma, histology == "Grade3")
>
> ### Plot Kaplan-Meier curves
> plot(survfit(Surv(time, event) ~ group, data=g3),
+ main="Grade III Glioma", lty=c(2,1),
+ legend.text=c("Control", "Treated"),
+ legend.bty=1, ylab="Probability",
+ xlab="Survival Time in Month")
>
> ### logrank test
> surv_test(Surv(time, event) ~ group, data = g3,
+ distribution = "exact")
Exact Logrank Test
data: y by groups Control, RIT
T = 2.1711, p-value = 0.02877
alternative hypothesis: two.sided
>
> ### Grade IV glioma
> gbm <- subset(glioma, histology == "GBM")
>
> ### Plot Kaplan-Meier curves
> plot(survfit(Surv(time, event) ~ group, data=gbm),
+ main="Grade IV Glioma", lty=c(2,1),
+ legend.text=c("Control", "Treated"),
+ legend.bty=1, legend.pos=1, ylab="Probability",
+ xlab="Survival Time in Month")
>
> ### logrank test
> surv_test(Surv(time, event) ~ group, data = gbm,
+ distribution = "exact")
Exact Logrank Test
data: y by groups Control, RIT
T = 3.2444, p-value = 0.0001429
alternative hypothesis: two.sided
>
> ### stratified logrank test
> surv_test(Surv(time, event) ~ group | histology, data = glioma,
+ distribution = "approx", B = 10000)
Approximative Logrank Test
data: y by
groups Control, RIT
stratified by histology
T = 3.6818, p-value < 2.2e-16
alternative hypothesis: two.sided
>
>
>
>
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "jobsatisfaction"
>
> ### * jobsatisfaction
>
> flush(stderr()); flush(stdout())
>
> ### Name: jobsatisfaction
> ### Title: Income and Job Satisfaction
> ### Aliases: jobsatisfaction
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(jobsatisfaction, package = "coin")
>
> # Generalized Cochran-Mantel-Haenzel test
> cmh_test(jobsatisfaction)
Asymptotical Generalised Cochran-Mantel-Haenszel Test
data: Job.Satisfaction by
groups <5000, 5000-15000, 15000-25000, >25000
stratified by Gender
T = 10.2001, df = 9, p-value = 0.3345
>
>
>
>
> cleanEx(); ..nameEx <- "neuropathy"
>
> ### * neuropathy
>
> flush(stderr()); flush(stdout())
>
> ### Name: neuropathy
> ### Title: Acute Painful Diabetic Neuropathy
> ### Aliases: neuropathy
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(neuropathy, package = "coin")
>
> ### compare with Table 2 of Conover & Salsburg (1988)
> oneway_test(pain ~ group, data = neuropathy, alternative = "less",
+ distribution = "exact")
Exact 2-Sample Permutation Test
data: pain by groups control, treat
T = -1.3191, p-value = 0.09589
alternative hypothesis: true mu is less than 0
>
> wilcox_test(pain ~ group, data = neuropathy, alternative = "less",
+ distribution = "exact")
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: pain by groups control, treat
T = -0.9817, p-value = 0.1654
alternative hypothesis: true mu is less than 0
>
> oneway_test(pain ~ group, data = neuropathy, distribution = "approx",
+ alternative = "less", ytrafo = function(data) trafo(data,
+ numeric_trafo = consal_trafo), B = 10000)
Approximative 2-Sample Permutation Test
data: pain by groups control, treat
T = -1.8683, p-value = 0.0284
alternative hypothesis: true mu is less than 0
>
>
>
>
> cleanEx(); ..nameEx <- "ocarcinoma"
>
> ### * ocarcinoma
>
> flush(stderr()); flush(stdout())
>
> ### Name: ocarcinoma
> ### Title: Ovarian Carcinoma
> ### Aliases: ocarcinoma
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(ocarcinoma, package = "coin")
>
> ### logrank test with exact two-sided p-value
> lrt <- surv_test(Surv(time, event) ~ stadium, data = ocarcinoma,
+ distribution = "exact")
>
> ### the test statistic
> statistic(lrt)
[1] -2.334131
>
> ### p-value
> pvalue(lrt)
[1] 0.01836702
>
>
>
>
> cleanEx(); ..nameEx <- "pvalue-methods"
>
> ### * pvalue-methods
>
> flush(stderr()); flush(stdout())
>
> ### Name: pvalue-methods
> ### Title: Extract P-Values
> ### Aliases: pvalue pvalue-methods pvalue,NullDistribution-method
> ### pvalue,IndependenceTest-method pvalue,ScalarIndependenceTest-method
> ### pvalue,MaxTypeIndependenceTest-method
> ### pvalue,QuadTypeIndependenceTest-method
> ### Keywords: methods htest
>
> ### ** Examples
>
>
> ### artificial 2-sample problem
> df <- data.frame(y = rnorm(20), x = gl(2, 10))
>
> ### Ansari-Bradley test
> at <- ansari_test(y ~ x, data = df, distribution = "exact")
>
> at
Exact Ansari-Bradley Test
data: y by groups 1, 2
T = 0.4553, p-value = 0.7102
alternative hypothesis: true mu is not equal to 1
>
> pvalue(at)
[1] 0.710234
>
>
>
>
> cleanEx(); ..nameEx <- "rotarod"
>
> ### * rotarod
>
> flush(stderr()); flush(stdout())
>
> ### Name: rotarod
> ### Title: Rotating Rats Data
> ### Aliases: rotarod
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(rotarod, package = "coin")
>
> ### Wilcoxon-Mann-Whitney Rank Sum Test
>
> ### one-sided exact (0.0186)
> wilcox_test(time ~ group, data = rotarod,
+ alternative = "greater", distribution = "exact")
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: time by groups control, treatment
T = 2.4389, p-value = 0.01863
alternative hypothesis: true mu is greater than 0
> ### two-sided exact (0.0373)
> wilcox_test(time ~ group, data = rotarod, distribution = "exact")
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: time by groups control, treatment
T = 2.4389, p-value = 0.03727
alternative hypothesis: true mu is not equal to 0
> ### two-sided asymptotical (0.0147)
> wilcox_test(time ~ group, data = rotarod)
Asymptotical Wilcoxon Mann-Whitney Rank Sum Test
data: time by groups control, treatment
T = 2.4389, p-value = 0.01473
alternative hypothesis: true mu is not equal to 0
>
>
>
>
> cleanEx(); ..nameEx <- "sphase"
>
> ### * sphase
>
> flush(stderr()); flush(stdout())
>
> ### Name: sphase
> ### Title: S-phase Fraction of Tumor Cells
> ### Aliases: sphase
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(sphase, package = "coin")
>
> maxstat_test(Surv(RFS, event) ~ SPF, data = sphase)
Asymptotical Maxstat Test
data: y by SPF
T = 2.4033, p-value = 0.1559
sample estimates:
cutpoint in SPF
107
>
>
>
>
> cleanEx(); ..nameEx <- "statistic-methods"
>
> ### * statistic-methods
>
> flush(stderr()); flush(stdout())
>
> ### Name: statistic-methods
> ### Title: Extract Test Statistics, Linear Statistics and Standardized
> ### Statistics
> ### Aliases: statistic statistic-methods statistic,IndependenceTest-method
> ### statistic,ScalarIndependenceTestStatistic-method
> ### statistic,MaxTypeIndependenceTestStatistic-method
> ### statistic,QuadTypeIndependenceTestStatistic-method
> ### Keywords: methods
>
> ### ** Examples
>
>
> df <- data.frame(y = gl(4, 5), x = gl(5, 4))
>
> ### Cochran-Mantel-Haenzel Test
> ct <- cmh_test(y ~ x, data = df)
>
> ### chisq-type statistics
> statistic(ct)
[1] 38
>
> ### the linear statistic, i.e, the contingency table
> statistic(ct, type = "linear")
1 2 3 4
1 4 0 0 0
2 1 3 0 0
3 0 2 2 0
4 0 0 3 1
5 0 0 0 4
>
> ### the same
> table(df$x, df$y)
1 2 3 4
1 4 0 0 0
2 1 3 0 0
3 0 2 2 0
4 0 0 3 1
5 0 0 0 4
>
> ### and the standardized contingency table for illustrating
> ### departures from the null-hypothesis of independence of x and y
> statistic(ct, type = "standardized")
1 2 3 4
1 3.774917 -1.258306 -1.258306 -1.258306
2 0.000000 2.516611 -1.258306 -1.258306
3 -1.258306 1.258306 1.258306 -1.258306
4 -1.258306 -1.258306 2.516611 0.000000
5 -1.258306 -1.258306 -1.258306 3.774917
>
>
>
> cleanEx(); ..nameEx <- "treepipit"
>
> ### * treepipit
>
> flush(stderr()); flush(stdout())
>
> ### Name: treepipit
> ### Title: Tree Pipit (Anthus trivialis) Forest Data
> ### Aliases: treepipit
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data(treepipit, package = "coin")
>
> maxstat_test(counts ~ age + coverstorey + coverregen + meanregen +
+ coniferous + deadtree + cbpiles + ivytree,
+ data = treepipit)
Asymptotical Maxstat Test
data: counts by
age, coverstorey, coverregen, meanregen, coniferous, deadtree, cbpiles, ivytree
T = 4.3139, p-value = 0.0008564
sample estimates:
cutpoint in coverstorey
40
>
>
>
>
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed: 20.67 0.15 20.99 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')
