inst/doc/Implementation.R

In coin: Conditional Inference Procedures in a Permutation Test Framework

### R code from vignette source 'Implementation.Rnw'

###################################################
### code chunk number 1: coin-setup
###################################################
options(prompt = "R> ", continue = "+  ")
library("coin")
library("e1071")
set.seed(290875)

### extract slots of a class
c2t <- function(x) {

    classdef <- getClassDef(x)

    extends <- names(classdef@contains)[1]
    if (!is.null(extends)) {
        eslots <- names(getClassDef(extends)@slots)
        slots <- classdef@slots[!names(classdef@slots) %in% eslots]
    } else {
        slots <- classdef@slots
    }

    RET <- cbind(names(slots), slots)
    attr(RET, "contains") <- extends
    attr(RET, "name") <- x
    class(RET) <- "c2t"
    RET
}

### pretty printing
toLatex.c2t <- function(object, center = TRUE, ...) {

    RET <- c()

    if (center) RET <- c(RET, "\\begin{center}")

    ### class name
    RET <- c(RET, "\\begin{tabular}{ll}",
                  paste("\\multicolumn{2}{l}{Class \\Rclass{",
                        attr(object, "name"), "}} \\\\", sep = ""))

    ### extends?
    if (!is.null(attr(object, "contains")))
        RET <- c(RET, paste("\\multicolumn{2}{l}{Contains \\Rclass{",
                            attr(object, "contains"), "}} \\\\", sep = ""))

    ### slots
    RET <- c(RET, " & \\\\", "Slot & Class \\\\ \\hline ",
             apply(object, 1, function(x) {
                 x <- cbind(paste("\\code{", x[1], "}", sep = ""),
                            paste("\\Rclass{", x[2], "}", sep = ""))
                 paste(paste(x, collapse = " & "), "\\\\")
             }),
             "\\hline")
    RET <- c(RET, "\\end{tabular}")

    if (center) RET <- c(RET, "\\end{center}")

    class(RET) <- "Latex"
    return(RET)
}


###################################################
### code chunk number 2: Ex
###################################################
library("coin")
data("rotarod", package = "coin")
independence_test(time ~ group, data = rotarod,
  ytrafo = rank_trafo, distribution = exact())


###################################################
### code chunk number 3: IndependenceProblem
###################################################
toLatex(c2t("IndependenceProblem"))


###################################################
### code chunk number 4: Ex-IndependenceProblem
###################################################
ip <- new("IndependenceProblem",
  y = rotarod["time"], x = rotarod["group"])


###################################################
### code chunk number 5: IndependenceTestProblem
###################################################
toLatex(c2t("IndependenceTestProblem"))


###################################################
### code chunk number 6: Ex-IndependenceTestProblem
###################################################
itp <- new("IndependenceTestProblem", ip, ytrafo = rank_trafo)


###################################################
### code chunk number 7: IndependenceLinearStatistic
###################################################
toLatex(c2t("IndependenceLinearStatistic"))


###################################################
### code chunk number 8: Ex-IndependenceLinearStatistic
###################################################
ils <- new("IndependenceLinearStatistic", itp)


###################################################
### code chunk number 9: Ex-IndependenceLinearStatistic-statistic
###################################################
statistic(ils, type = "linear")


###################################################
### code chunk number 10: Ex-IndependenceLinearStatistic-statistic
###################################################
expectation(ils)
variance(ils)


###################################################
### code chunk number 11: IndependenceTestStatistic
###################################################
toLatex(c2t("IndependenceTestStatistic"))


###################################################
### code chunk number 12: ScalarIndependenceTestStatistic
###################################################
toLatex(c2t("ScalarIndependenceTestStatistic"))


###################################################
### code chunk number 13: Ex-ScalarIndependenceTestStatistic
###################################################
sits <- new("ScalarIndependenceTestStatistic", ils,
  alternative = "two.sided")
statistic(sits, type = "standardized")


###################################################
### code chunk number 14: MaxTypeIndependenceTestStatistic
###################################################
toLatex(c2t("MaxTypeIndependenceTestStatistic"))


###################################################
### code chunk number 15: QuadTypeIndependenceTestStatistic
###################################################
toLatex(c2t("QuadTypeIndependenceTestStatistic"))


###################################################
### code chunk number 16: PValue
###################################################
toLatex(c2t("PValue"))


###################################################
### code chunk number 17: NullDistribution
###################################################
toLatex(c2t("NullDistribution"))


###################################################
### code chunk number 18: Ex-NullDistribution-pvalue
###################################################
end <- ExactNullDistribution(sits)
pvalue(end, statistic(sits))
qperm(end, 0.95)


###################################################
### code chunk number 19: IndependenceTest
###################################################
toLatex(c2t("IndependenceTest"))


###################################################
### code chunk number 20: IndependenceTest
###################################################
new("IndependenceTest", statistic = sits, distribution = end)


###################################################
### code chunk number 21: Ex-distribution
###################################################
set.seed(2908)
correxample <- data.frame(x = rnorm(7), y = rnorm(7))
sexact <- function(object) {
  x <- object@xtrans
  y <- object@ytrans
  perms <- permutations(nrow(x))
  pstats <- apply(perms, 1, function(p) sum(x[p,] * y))
  pstats <- (pstats - as.vector(expectation(object))) /
      sqrt(as.vector(variance(object)))
  p <- function(q) 1 - mean(pstats > q)
  new("PValue", name = "Exact Distribution", p = p, pvalue = p)
}


###################################################
### code chunk number 22: Ex-distribution
###################################################
independence_test(y ~ x, data = correxample, alternative = "less",
  distribution = sexact)


###################################################
### code chunk number 23: Implementation.Rnw:867-868
###################################################
mood_score <- function(y) (rank_trafo(y) - (sum(!is.na(y)) + 1) / 2)^2


###################################################
### code chunk number 24: Implementation.Rnw:872-881
###################################################
ip <- new("IndependenceProblem",
  y = rotarod["time"], x = rotarod["group"])
itp <- new("IndependenceTestProblem", ip,
  ytrafo = mood_score)
ils <- new("IndependenceLinearStatistic", itp)
sits <- new("ScalarIndependenceTestStatistic", ils,
  alternative = "two.sided")
new("ScalarIndependenceTest", statistic = sits,
  distribution = ExactNullDistribution(sits, algorithm = "split-up"))


###################################################
### code chunk number 25: Implementation.Rnw:885-887
###################################################
independence_test(time ~ group, data = rotarod, ytrafo = mood_score,
  distribution = exact(algorithm = "split-up"))


###################################################
### code chunk number 26: js
###################################################
data("jobsatisfaction", package = "coin")
js <- jobsatisfaction
dimnames(js)[[2]] <- c("VeryDiss", "LitSat", "ModSat", "VerySat")
ftable(Job.Satisfaction ~ Gender + Income, data = js)


###################################################
### code chunk number 27: js-plot
###################################################
library("vcd")
cotabplot(js,
  split_vertical = TRUE, spacing = spacing_highlighting,
  gp = gpar(fill = rev(gray.colors(4))),
  labeling_args = list(rot_labels = 0, varnames = FALSE,
    just_labels = c("center", "right")),
  panel_args = list(margins = c(3, 1, 2, 3.5)))


###################################################
### code chunk number 28: jobsatisfaction-it
###################################################
it <- independence_test(js, teststat = "quadratic",
  distribution = asymptotic())
it


###################################################
### code chunk number 29: jobsatisfaction-T
###################################################
statistic(it, type = "linear")


###################################################
### code chunk number 30: jobsatisfaction-margin
###################################################
margin.table(js, 1:2)


###################################################
### code chunk number 31: jobsatisfaction-stat
###################################################
statistic(it, type = "standardized")


###################################################
### code chunk number 32: jobsatisfaction-ordinal
###################################################
it <- independence_test(js, distribution = approximate(nresample = 10000),
  scores = list(Job.Satisfaction = 1:4, Income = 1:4))
pvalue(it)


###################################################
### code chunk number 33: jobsatisfaction-max
###################################################
independence_test(js, teststat = "maximum")


###################################################
### code chunk number 34: jobsatisfaction-minp
###################################################
pvalue(independence_test(js, teststat = "maximum"),
       method = "single-step")


###################################################
### code chunk number 35: coin-doxygen (eval = FALSE)
###################################################
## browseURL(system.file("documentation", "html", "index.html",
##   package = "coin"))