Nothing
### 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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.