https://github.com/cran/coin
Tip revision: f2c846c830000634333e7f8ba60468d503e80753 authored by Torsten Hothorn on 30 November 2005, 00:00:00 UTC
version 0.4-1
version 0.4-1
Tip revision: f2c846c
coin-Ex.Rout.save
R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.2.0 (2005-10-06 r35749)
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)
+ delayedAssign("T", stop("T used instead of TRUE"),
+ assign.env = .CheckExEnv)
+ delayedAssign("F", stop("F used instead of FALSE"),
+ assign.env = .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-Ex.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)
> library('coin')
Loading required package: survival
Loading required package: splines
Loading required package: mvtnorm
Loading required package: modeltools
>
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "ContingencyTests"
>
> ### * ContingencyTests
>
> flush(stderr()); flush(stdout())
>
> ### Name: ContingencyTests
> ### Title: Independence in Three-Way 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 = approximate(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)
Asymptotic 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)
Asymptotic 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)))
Asymptotic 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 = approximate(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
>
> ### Smoking and HDL cholesterin status
> ### (from Jeong, Jhun and Kim, 2005, CSDA 48, 623-631, Table 2)
> smokingHDL <- as.table(
+ matrix(c(15, 8, 11, 5,
+ 3, 4, 6, 1,
+ 6, 7, 15, 11,
+ 1, 2, 3, 5), ncol = 4,
+ dimnames = list(smoking = c("none", "< 5", "< 10", ">=10"),
+ HDL = c("normal", "low", "borderline", "abnormal"))
+ ))
> ### use interval mid-points as scores for smoking
> lbl_test(smokingHDL, scores = list(smoking = c(0, 2.5, 7.5, 15)))
Asymptotic Linear-by-Linear Association Test
data: HDL (ordered) by groups none < < 5 < < 10 < >=10
T = 9.8557, df = 1, p-value = 0.001693
>
>
>
>
> 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")
12-26 Weeks 30
>
> ### its expectation
> expectation(wt)
12-26 Weeks
40
>
> ### and variance
> covariance(wt)
12-26 Weeks
12-26 Weeks 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)
Asymptotic 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 = approximate(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 (1999), 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 = approximate(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 = approximate(B = 90000))
+
+ ### global p-value
+ print(pvalue(NDWD))
+
+ ### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
+ print(pvalue(NDWD, method = "single-step"))
+ }
Loading required package: multcomp
[1] 0.0004888889
99 percent confidence interval:
0.0003199269 0.0007126921
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 <- "MarginalHomogeneityTest"
>
> ### * MarginalHomogeneityTest
>
> flush(stderr()); flush(stdout())
>
> ### Name: MarginalHomogeneityTest
> ### Title: Marginal Homogeneity 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)
Asymptotic Marginal-Homogenity Test
data: response by
groups ExtramaritalSex, PremaritalSex
stratified by block
T = 271.9222, df = 3, p-value < 2.2e-16
>
> ### and as ordinal
> mh_test(PreExSex, scores = list(response = 1:length(opinions)))
Asymptotic Marginal-Homogenity Test for Ordered Data
data: response (ordered) by
groups ExtramaritalSex, PremaritalSex
stratified by block
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)
Asymptotic 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)
>
> ### plot density
> plot(support(at), dens, type = "s")
>
> ### 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)
Asymptotic 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 (1999), 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 (1999), 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 = approximate(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")
>
> ### points to a difference in location
> pvalue(ltmax, method = "single-step")
Location Scale
Control 0.0067991 0.6189816
>
> ### Funny: We could have used a simple Bonferroni procedure
> ### since the correlation between the Wilcoxon and Ansari-Bradley
> ### test statistics is zero
> covariance(ltmax)
Control:Location Control:Scale
Control:Location 85 0
Control:Scale 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)
Asymptotic 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
>
>
> ### asymptotic tests for carcinoma data
> data(ocarcinoma, package = "coin")
> surv_test(Surv(time, event) ~ stadium, data = ocarcinoma)
Asymptotic Logrank Test
data: Surv(time, event) by groups II, IIA
T = -2.3373, p-value = 0.01942
alternative hypothesis: two.sided
> survdiff(Surv(time, event) ~ stadium, data = ocarcinoma)
Call:
survdiff(formula = Surv(time, event) ~ stadium, data = ocarcinoma)
N Observed Expected (O-E)^2/E (O-E)^2/V
stadium=II 15 6 11.3 2.51 5.57
stadium=IIA 20 16 10.7 2.67 5.57
Chisq= 5.6 on 1 degrees of freedom, p= 0.0183
>
> ### example data given in Callaert (2003)
> 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 = 0.0523
> survdiff(Surv(time, event) ~ group, data = exdata)
Call:
survdiff(formula = Surv(time, event) ~ group, data = exdata)
N Observed Expected (O-E)^2/E (O-E)^2/V
group=0 7 7 9.84 0.82 3.76
group=1 8 8 5.16 1.56 3.76
Chisq= 3.8 on 1 degrees of freedom, p= 0.0523
> ### p = 0.0505
> surv_test(Surv(time, event) ~ group, data = exdata,
+ distribution = exact())
Exact Logrank Test
data: Surv(time, event) by groups 0, 1
T = -1.9201, p-value = 0.0505
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)
Asymptotic 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 (1999), 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) %*% t(contrMat(table(x), "Tukey"))
+ ),
+ 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, method = "single-step"))
+ }
Loading required package: multcomp
[1] 0.003455914
99 percent confidence interval:
0.003205723 0.003706106
Round Out-Narrow Angle 0.623924061
Wide Angle-Narrow Angle 0.053791966
Wide Angle-Round Out 0.003398077
>
> ### 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 for ordered alternatives
> ft <- friedman_test(strength ~ potash | block, data = sc)
> ft
Asymptotic 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 = approximate(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
> ### 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)
x y.1 y.2 y.3 y.4 y.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))
x <= 0.39 x <= -0.621 x <= 1.125 x <= -0.045 x <= -0.016 x <= 0.944
1 0 0 0 0 0 0
2 1 0 1 0 0 1
3 1 1 1 1 1 1
4 1 1 1 1 1 1
5 0 0 1 0 0 0
6 1 0 1 1 1 1
7 1 0 1 0 1 1
8 0 0 1 0 0 1
9 0 0 1 0 0 1
10 0 0 1 0 0 1
x <= 0.821 x <= 0.594
1 0 0
2 1 1
3 1 1
4 1 1
5 0 0
6 1 1
7 1 1
8 0 0
9 1 0
10 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
[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 <- "alzheimer"
>
> ### * alzheimer
>
> flush(stderr()); flush(stdout())
>
> ### Name: alzheimer
> ### Title: Smoking and Alzheimer's Disease
> ### Aliases: alzheimer
> ### Keywords: datasets
>
> ### ** Examples
>
> data("alzheimer")
>
> ### Cochran-Mantel-Haenszel test
> cmh_test(alzheimer)
Asymptotic Generalised Cochran-Mantel-Haenszel Test
data: disease by
groups None, <10, 10-20, >20
stratified by gender
T = 23.3163, df = 6, p-value = 0.0006972
>
> ### Linear-by-Linear Association test
> cmh_test(alzheimer, scores = list(smoking = c(0, 5, 15, 25)))
Asymptotic Linear-by-Linear Association Test
data: disease by
groups None < <10 < 10-20 < >20
stratified by gender
T = 10.3366, df = 2, p-value = 0.005694
> statistic(cmh_test(alzheimer, scores = list(smoking = c(0, 5, 15, 25))),
+ "standardized")
Alzheimer's Other dementias Other diagnoses
0*None + 5*<10 + 15*10-2 + 25*>20 -2.334202 3.087989 -0.6459272
>
>
>
>
> 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 Variance / Covariance of Linear
> ### Statistics
> ### Aliases: expectation expectation-methods
> ### expectation,IndependenceTest-method
> ### expectation,IndependenceTestStatistic-method
> ### expectation,ScalarIndependenceTestStatistic-method
> ### expectation,MaxTypeIndependenceTestStatistic-method
> ### expectation,QuadTypeIndependenceTestStatistic-method covariance
> ### covariance-methods covariance,CovarianceMatrix-method
> ### covariance,IndependenceTest-method
> ### covariance,IndependenceTestStatistic-method
> ### covariance,ScalarIndependenceTestStatistic-method
> ### covariance,MaxTypeIndependenceTestStatistic-method
> ### covariance,QuadTypeIndependenceTestStatistic-method variance
> ### variance-methods variance,Variance-method
> ### variance,CovarianceMatrix-method variance,IndependenceTest-method
> ### variance,IndependenceTestStatistic-method
> ### variance,ScalarIndependenceTestStatistic-method
> ### variance,MaxTypeIndependenceTestStatistic-method
> ### variance,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:1 2:1 3:1 1:2 2:2 3:2 1:3 2:3
0.6666667 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667
3:3
0.6666667
>
> ### covariance
> Vl <- covariance(ct)
> Vl
1:1 2:1 3:1 1:2 2:2 3:2
1:1 0.35555556 -0.17777778 -0.17777778 -0.17777778 0.08888889 0.08888889
2:1 -0.17777778 0.35555556 -0.17777778 0.08888889 -0.17777778 0.08888889
3:1 -0.17777778 -0.17777778 0.35555556 0.08888889 0.08888889 -0.17777778
1:2 -0.17777778 0.08888889 0.08888889 0.35555556 -0.17777778 -0.17777778
2:2 0.08888889 -0.17777778 0.08888889 -0.17777778 0.35555556 -0.17777778
3:2 0.08888889 0.08888889 -0.17777778 -0.17777778 -0.17777778 0.35555556
1:3 -0.17777778 0.08888889 0.08888889 -0.17777778 0.08888889 0.08888889
2:3 0.08888889 -0.17777778 0.08888889 0.08888889 -0.17777778 0.08888889
3:3 0.08888889 0.08888889 -0.17777778 0.08888889 0.08888889 -0.17777778
1:3 2:3 3:3
1:1 -0.17777778 0.08888889 0.08888889
2:1 0.08888889 -0.17777778 0.08888889
3:1 0.08888889 0.08888889 -0.17777778
1:2 -0.17777778 0.08888889 0.08888889
2:2 0.08888889 -0.17777778 0.08888889
3:2 0.08888889 0.08888889 -0.17777778
1:3 0.35555556 -0.17777778 -0.17777778
2:3 -0.17777778 0.35555556 -0.17777778
3:3 -0.17777778 -0.17777778 0.35555556
>
> ### the standardized contingency table (hard way)
> (l - El) / sqrt(variance(ct))
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: Surv(time, event) 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: Surv(time, event) by groups Control, RIT
T = 3.2215, p-value = 0.0001588
alternative hypothesis: two.sided
>
> ### stratified logrank test
> surv_test(Surv(time, event) ~ group | histology, data = glioma,
+ distribution = approximate(B = 10000))
Approximative Logrank Test
data: Surv(time, event) by
groups Control, RIT
stratified by histology
T = 3.6704, 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)
Asymptotic 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 <- "mercuryfish"
>
> ### * mercuryfish
>
> flush(stderr()); flush(stdout())
>
> ### Name: mercuryfish
> ### Title: Chromosomal Effects of Mercury Contaminated Fish Consumption
> ### Aliases: mercuryfish
> ### Keywords: datasets
>
> ### ** Examples
>
> data("mercuryfish")
>
> coherence <- function(data) {
+ x <- as.matrix(data)
+ matrix(apply(x, 1, function(y)
+ sum(colSums(t(x) < y) == ncol(x)) -
+ sum(colSums(t(x) > y) == ncol(x))), ncol = 1)
+ }
>
> ### POSET-test
> poset <- independence_test(mercury + abnormal + ccells ~ group, data =
+ mercuryfish, ytrafo = coherence)
>
> ### linear statistic (T in Rosenbaum's, 1994, notation)
> statistic(poset, "linear")
control -237
>
> ### expectation
> expectation(poset)
control
0
>
> ### variance (Rosenbaum, 1994, uses the unconditional approach)
> covariance(poset)
control
control 3097.954
>
> ### the standardized statistic
> statistic(poset)
control
-4.258051
>
> ### and asymptotic p-value
> pvalue(poset)
[1] 2.062169e-05
>
> ### exact p-value
> independence_test(mercury + abnormal + ccells ~ group, data =
+ mercuryfish, ytrafo = coherence, distribution = "exact")
Exact General Independence Test
data: mercury, abnormal, ccells by groups control, exposed
T = -4.2581, p-value = 4.486e-06
alternative hypothesis: two.sided
>
> ### multivariate analysis
> mvtest <- independence_test(mercury + abnormal + ccells ~ group,
+ data = mercuryfish)
>
> ### global p-value
> pvalue(mvtest)
[1] 0.007050273
99 percent confidence interval:
0.006381220 0.007719326
>
> ### adjusted univariate p-value
> pvalue(mvtest, method = "single-step")
mercury abnormal ccells
control 0.007191719 0.01705445 0.03836982
>
>
>
>
> 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 = approximate(B = 10000),
+ alternative = "less", ytrafo = function(data) trafo(data,
+ numeric_trafo = consal_trafo))
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)
II
-2.337284
>
> ### p-value
> pvalue(lrt)
[1] 0.01819756
>
>
>
>
> cleanEx(); ..nameEx <- "photocar"
>
> ### * photocar
>
> flush(stderr()); flush(stdout())
>
> ### Name: photocar
> ### Title: Multiple Dosing Photococarcinogenicity Experiment
> ### Aliases: photocar
> ### Keywords: datasets
>
> ### ** Examples
>
>
> data("photocar")
>
> layout(matrix(1:3, ncol = 3))
> plot(survfit(Surv(time, event) ~ group, data = photocar), xmax = 50, lty =
+ 1:3, main = "Survival Time")
> legend("bottomleft", lty = 1:3, levels(photocar$group), bty = "n")
> plot(survfit(Surv(dmin, tumor) ~ group, data = photocar), xmax = 50, lty =
+ 1:3, main = "Time to First Tumor")
> legend("bottomleft", lty = 1:3, levels(photocar$group), bty = "n")
> boxplot(ntumor ~ group, data = photocar, main = "Number of Tumors")
>
> ### global test (all three responses)
> fm <- Surv(time, event) + Surv(dmin, tumor) + ntumor ~ group
> it <- independence_test(fm, data = photocar,
+ distribution = approximate(B = 10000))
> pvalue(it)
[1] 0
99 percent confidence interval:
0.0000000000 0.0005296914
>
> ### why was the global null hypothesis rejected?
> statistic(it, "standardized")
Surv(time, event) Surv(dmin, tumor) ntumor
A -2.327338 -2.178704 0.2642120
B -4.750336 -4.106039 0.1509783
C 7.077674 6.284743 -0.4151904
> pvalue(it, "single-step")
Surv(time, event) Surv(dmin, tumor) ntumor
A 0.129 0.1850 0.9999
B 0.000 0.0003 1.0000
C 0.000 0.0000 0.9988
>
>
>
>
> 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
>
> ### bivariate 2-sample problem
> df <- data.frame(y1 = rnorm(20) + c(rep(0, 10), rep(1, 10)),
+ y2 = rnorm(20),
+ x = gl(2, 10))
>
> it <- independence_test(y1 + y2 ~ x, data = df,
+ distribution = approximate(B = 9999))
> pvalue(it, method = "single-step")
y1 y2
1 0.0110011 0.9998
> pvalue(it, method = "step-down")
y1 y2
1 0.0110011 0.9789979
> pvalue(it, method = "discrete")
y1 y2
1 0.01097084 0.999559
>
>
>
>
> 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)
Asymptotic 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)
Asymptotic Maxstat Test
data: Surv(RFS, event) by SPF
T = 2.3966, p-value = 0.1582
sample estimates:
cutpoint in SPF
107
>
> ### reproduce the test statistic reported in Hothorn & Lausen (2003)
> maxstat_test(Surv(RFS, event) ~ SPF, data = sphase,
+ ytrafo = function(data) trafo(data, surv_trafo = function(x)
+ logrank_trafo(x, ties.method = "HL")))
Asymptotic Maxstat Test
data: Surv(RFS, event) by SPF
T = 2.4033, p-value = 0.1557
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
> ### 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)
Asymptotic 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: 46.51 1.08 49.81 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')
