LegoCondInf.R
In coin: Conditional Inference Procedures in a Permutation Test Framework

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

###################################################
### code chunk number 1: setup
###################################################
options(width = 65, prompt = "R> ", continue = "   ")
require("coin")
set.seed(290875)
### get rid of the NAMESPACE
#load(file.path(.find.package("coin"), "R", "all.rda"))
anonymous <- FALSE


###################################################
### code chunk number 2: authors
###################################################
if(!anonymous) {
    cat("\\author{Torsten Hothorn$^1$, Kurt Hornik$^2$, \\\\
            Mark A. van de Wiel$^3$ and Achim Zeileis$^2$}\n")
} else  {
    cat("\\author{TAS MS05-239, Revision}\n")
}


###################################################
### code chunk number 3: affil
###################################################
if(!anonymous)
    cat("\\noindent$^1$ Institut f\\\"ur Medizininformatik, Biometrie und Epidemiologie\\\\
           Friedrich-Alexander-Universit\\\"at Erlangen-N\\\"urnberg\\\\
           Waldstra{\\ss}e 6, D-91054 Erlangen, Germany \\\\
           \\texttt{Torsten.Hothorn@R-project.org}
         \\newline

         \\noindent$^2$ Department f\\\"ur Statistik und Mathematik,
            Wirtschaftsuniversit\\\"at Wien \\\\
            Augasse 2-6, A-1090 Wien, Austria \\\\
            \\texttt{Kurt.Hornik@R-project.org} \\\\
            \\texttt{Achim.Zeileis@R-project.org}
         \\newline

         \\noindent$^3$ Department of Mathematics, Vrije Universiteit \\\\
                        De Boelelaan 1081a, 1081 HV Amsterdam, The Netherlands \\\\
            \\texttt{mark.vdwiel@vumc.nl}
         \\newline\n")


###################################################
### code chunk number 4: coincite
###################################################
if (anonymous) {
    cat(" \\citep{PKG:coina} ")
} else {
    cat(" \\citep{PKG:coin} ")
}


###################################################
### code chunk number 5: alpha-data-figure
###################################################
n <- table(alpha$alength)
par(cex.lab = 1.3, cex.axis = 1.3)
boxplot(elevel ~ alength, data = alpha, ylab = "Expression Level",
        xlab = "NACP-REP1 Allele Length", varwidth = TRUE)
axis(3, at = 1:3, labels = paste("n = ", n))
rankif <- function(data) trafo(data, numeric_trafo = rank_trafo)


###################################################
### code chunk number 6: alpha-kruskal
###################################################
kruskal.test(elevel ~ alength, data = alpha)


###################################################
### code chunk number 7: alpha-kruskal
###################################################
independence_test(elevel ~ alength, data = alpha, ytrafo = rank_trafo, teststat = "quadratic")


###################################################
### code chunk number 8: mpoints
###################################################
mpoints <- function(x) c(2, 7, 11)[unlist(x)]


###################################################
### code chunk number 9: alpha-kruskal-ordered
###################################################
independence_test(elevel ~ alength, data = alpha, ytrafo = rank_trafo, xtrafo = mpoints)


###################################################
### code chunk number 10: alzheimer-demographics
###################################################
total <- nrow(alzheimer)
stopifnot(total == 538)
male <- sum(alzheimer$gender == "Male")
stopifnot(male == 200)
female <- sum(alzheimer$gender == "Female")
stopifnot(female == 338)
disease <- table(alzheimer$disease)
smoked <- sum(alzheimer$smoking != "None")
atab <- xtabs(~ smoking + + disease + gender, data = alzheimer)
### there is a discrepancy between Table 1 (32% smokers of 117 women
### suffering from other diagnoses) and Table 4 (63% non-smokers).
### We used the data as given in Table 4.


###################################################
### code chunk number 11: alzheimer-tab
###################################################
x <- t(atab[,,"Female"])
lines <- paste(paste(dimnames(x)$disease, " & "),
               paste(apply(x, 1, function(l) paste(l, collapse = " & ")), "\\\\"))
for (i in 1:length(lines)) cat(lines[i], "\n")


###################################################
### code chunk number 12: alzheimer-tab
###################################################
x <- t(atab[,,"Male"])
lines <- paste(paste(dimnames(x)$disease, " & "),
               paste(apply(x, 1, function(l) paste(l, collapse = " & ")), "\\\\"))
for (i in 1:length(lines)) cat(lines[i], "\n")


###################################################
### code chunk number 13: alzheimer-plot
###################################################
layout(matrix(1:2, ncol = 2))
spineplot(disease ~ smoking, data = alzheimer, subset = gender == "Male",
main = "Male", xlab = "Smoking", ylab = "Disease", tol = 1)
spineplot(disease ~ smoking, data = alzheimer, subset = gender == "Female",
main = "Female", xlab = "Smoking", ylab = "Disease", tol = 1)


###################################################
### code chunk number 14: alzheimer-mantelhaen
###################################################
it_alz <- independence_test(disease ~ smoking | gender, data = alzheimer,
                            teststat = "quadratic")
it_alz


###################################################
### code chunk number 15: alzheimer-statistic
###################################################
statistic(it_alz, type = "linear")


###################################################
### code chunk number 16: alzheimer-men
###################################################
females <- alzheimer$gender == "Female"
males <- alzheimer$gender == "Male"
pvalue(independence_test(disease ~ smoking, data = alzheimer,
       subset = females, teststat = "quadratic"))
pvalue(independence_test(disease ~ smoking, data = alzheimer,
       subset = males, teststat = "quadratic"))


###################################################
### code chunk number 17: alzheimer-max
###################################################
it_alzmax <- independence_test(disease ~ smoking, data = alzheimer,
       subset = males, teststat = "maximum")
it_alzmax


###################################################
### code chunk number 18: alzheimer-maxstat
###################################################
statistic(it_alzmax, type = "standardized")


###################################################
### code chunk number 19: alzheimer-qperm
###################################################
qperm(it_alzmax, 0.95)


###################################################
### code chunk number 20: alzheimer-MTP
###################################################
pvalue(it_alzmax, method = "single-step")


###################################################
### code chunk number 21: photocar-plot
###################################################
par(cex.lab = 1.3, cex.axis = 1.3)
layout(matrix(1:3, ncol = 3))
plot(survfit(Surv(time, event) ~ group, data = photocar), xmax = 50,
     xlab = "Survival Time (in weeks)", ylab = "Probability",
     lty = 1:3)
     legend("bottomleft", lty = 1:3, levels(photocar$group), bty = "n")
plot(survfit(Surv(dmin, tumor) ~ group, data = photocar), xmax = 50,
     xlab = "Time to First Tumor (in weeks)", ylab = "Probability",
     lty = 1:3)
     legend("bottomleft", lty = 1:3, levels(photocar$group), bty = "n")
boxplot(ntumor ~ group, data = photocar,
        ylab = "Number of Tumors", xlab = "Treatment Group",
        varwidth = TRUE)


###################################################
### code chunk number 22: photocar-global
###################################################
it_ph <- independence_test(Surv(time, event) + Surv(dmin, tumor) + ntumor ~ group,
                           data = photocar)
it_ph


###################################################
### code chunk number 23: photocar-linear
###################################################
statistic(it_ph, type = "linear")


###################################################
### code chunk number 24: photocar-stand
###################################################
statistic(it_ph, type = "standardized")


###################################################
### code chunk number 25: photocar-stand (eval = FALSE)
###################################################
## pvalue(it_ph, method = "single-step")


###################################################
### code chunk number 26: photocar-stand
###################################################
round(pvalue(it_ph, method = "single-step"), 5)


###################################################
### code chunk number 27: mercuryfish-plot
###################################################
par(cex.lab = 1.3, cex.axis = 1.3)
layout(matrix(1:3, ncol = 3))
boxplot(I(log(mercury)) ~ group, data = mercuryfish,
        ylab = "Mercury Blood Level (in logs)", varwidth = TRUE)
boxplot(abnormal ~ group, data = mercuryfish,
        ylab = "Abnormal Cells (in %)", varwidth = TRUE)
boxplot(ccells ~ group, data = mercuryfish,
        ylab = "Chromosome Aberrations (in %)", varwidth = TRUE)


###################################################
### code chunk number 28: mercurysfish-score
###################################################
coherence <- function(data) {
    x <- t(as.matrix(data))
    apply(x, 2, function(y)
        sum(colSums(x < y) == nrow(x)) - sum(colSums(x > y) == nrow(x)))
}


###################################################
### code chunk number 29: mercuryfish-poset
###################################################
poset <- independence_test(mercury + abnormal + ccells ~ group,
    data = mercuryfish, ytrafo = coherence, distribution = exact())


###################################################
### code chunk number 30: mercuryfish-pvalue
###################################################
pvalue(poset)


###################################################
### code chunk number 31: mercuryfish-ppermplot
###################################################
par(cex.lab = 1.3, cex.axis = 1.1)
ite <- poset
ita <- independence_test(mercury + abnormal + ccells ~ group, data =
                           mercuryfish, ytrafo = coherence)
site <- support(ite)
layout(matrix(1:2, ncol = 2))
site <- site[site <= qperm(ite, 0.1) & site > -3]
pite <- sapply(site, function(x) pperm(ite, x))
pita <- sapply(site, function(x) pperm(ita, x))

plot(site, pite, type = "S", ylab = "Probability", xlab = "Standardized Statistic")
lines(site, pita, lty = 3)
legend("topleft", lty = c(1,3), legend = c("Conditional Distribution",
"Approximation"), bty = "n")

site <- support(ite)
site <- site[site >= qperm(ite, 0.9) & site < 3]
pite <- sapply(site, function(x) pperm(ite, x))
pita <- sapply(site, function(x) pperm(ita, x))

plot(site, pite, type = "S", ylab = "Probability", xlab = "Standardized Statistic")
lines(site, pita, lty = 3)