R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), 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("maxstat-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('maxstat')
Loading required package: exactRankTests
Loading required package: mvtnorm
Loading required package: survival
Loading required package: splines
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "DLBCL"
> 
> ### * DLBCL
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: DLBCL
> ### Title: Diffuse large B-cell lymphoma
> ### Aliases: DLBCL
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> data(DLBCL)
> 
> # compute the cutpoint and plot the empirical process 
> 
> mod <- maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL, smethod="LogRank")
> 
> print(mod)

	Maximally selected LogRank statistics using none

data:  Surv(time, cens) by MGE 
M = 3.1772, p-value = NA
sample estimates:
estimated cutpoint 
         0.1860526 

> 
> ## Not run: 
> ##D   # postscript("statDLBCL.ps", horizontal=F, width=8, height=8)
> ##D   pdf("statDLBCL.pdf", width=8, height=8)
> ## End(Not run)
> par(mai=c(1.0196235, 1.0196235, 0.8196973, 0.4198450))
> plot(mod, cex.lab=1.6, cex.axis=1.6, xlab="Mean gene expression",lwd=2)
> ## Not run: 
> ##D   dev.off()
> ## End(Not run)
> 
> # significance of the cutpoint
> # limiting distribution
> 
> maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL,
+              smethod="LogRank", pmethod="Lau92", iscores=TRUE)

	Maximally selected LogRank statistics using Lau92

data:  Surv(time, cens) by MGE 
M = 3.171, p-value = 0.03611
sample estimates:
estimated cutpoint 
         0.1860526 

> 
> # improved Bonferroni inequality, plot with significance bound
> 
> maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL,
+              smethod="LogRank", pmethod="Lau94", iscores=TRUE)

	Maximally selected LogRank statistics using Lau94

data:  Surv(time, cens) by MGE 
M = 3.171, p-value = 0.02435
sample estimates:
estimated cutpoint 
         0.1860526 

> 
> mod <- maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL, smethod="LogRank",
+                     pmethod="Lau94", alpha=0.05)
> plot(mod, xlab="Mean gene expression")
> 
> ## Not run: 
> ##D #  postscript(file="RNewsStat.ps",horizontal=F, width=8, height=8)
> ##D    pdf("RNewsStat.pdf", width=8, height=8)
> ##D 
> ## End(Not run)
> par(mai=c(1.0196235, 1.0196235, 0.8196973, 0.4198450))
> plot(mod, xlab="Mean gene expression", cex.lab=1.6, cex.axis=1.6)
> ## Not run: 
> ##D   dev.off()
> ## End(Not run)
> 
> # small sample solution Hothorn & Lausen
> 
> maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL,
+              smethod="LogRank", pmethod="HL")

	Maximally selected LogRank statistics using HL

data:  Surv(time, cens) by MGE 
M = 3.171, p-value = 0.02218
sample estimates:
estimated cutpoint 
         0.1860526 

> 
> # normal approximation
> 
> maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL,
+              smethod="LogRank", pmethod="exactGauss", iscores=TRUE,
+              abseps=0.01)

	Maximally selected LogRank statistics using exactGauss

data:  Surv(time, cens) by MGE 
M = 3.171, p-value = 0.01030
sample estimates:
estimated cutpoint 
         0.1860526 

> 
> # conditional Monte-Carlo
> 
> maxstat.test(Surv(time, cens) ~ MGE, data=DLBCL,
+              smethod="LogRank", pmethod="condMC", B = 9999) 

	Maximally selected LogRank statistics using condMC

data:  Surv(time, cens) by MGE 
M = 3.1772, p-value = 0.0112
sample estimates:
estimated cutpoint 
         0.1860526 

> 
> # survival analysis and plotting like in Alizadeh et al. (2000)
> 
> if(require(survival, quietly = TRUE)) {
+ 
+   splitGEG <- rep(1, nrow(DLBCL))
+   DLBCL <- cbind(DLBCL, splitGEG)
+   DLBCL$splitGEG[DLBCL$GEG == "Activated B-like"] <- 0
+ 
+   plot(survfit(Surv(time, cens) ~ splitGEG, data=DLBCL),
+        xlab="Survival time in month", ylab="Probability")
+ 
+   text(90, 0.7, "GC B-like")
+   text(60, 0.3, "Activated B-like")
+ 
+   splitIPI <- rep(1, nrow(DLBCL))
+   DLBCL <- cbind(DLBCL, splitIPI)
+   DLBCL$splitIPI[DLBCL$IPI <= 2] <- 0
+ 
+   plot(survfit(Surv(time, cens) ~ splitIPI, data=DLBCL),
+        xlab="Survival time in month", ylab="Probability")
+ 
+   text(90, 0.7, "Low clinical risk")
+   text(60, 0.25, "High clinical risk")
+ 
+   # survival analysis using the cutpoint 
+ 
+   splitMGE <- rep(1, nrow(DLBCL))
+   DLBCL <- cbind(DLBCL, splitMGE)
+   DLBCL$splitMGE[DLBCL$MGE <= mod$estimate] <- 0
+ 
+   ## Not run: 
+ ##D    # postscript("survDLBCL.ps",horizontal=F, width=8, height=8)
+ ##D     pdf("survDLBCL.pdf", width=8, height=8)
+ ##D 
+ ##D   
+ ## End(Not run)
+   par(mai=c(1.0196235, 1.0196235, 0.8196973, 0.4198450))
+ 
+   plot(survfit(Surv(time, cens) ~ splitMGE, data=DLBCL),
+        xlab = "Survival time in month",
+        ylab="Probability", cex.lab=1.6, cex.axis=1.6, lwd=2)
+ 
+   text(90, 0.9, expression("Mean gene expression" > 0.186), cex=1.6)   
+   text(90, 0.45, expression("Mean gene expression" <= 0.186 ), cex=1.6)   
+ 
+   ## Not run: 
+ ##D     dev.off()
+ ##D   
+ ## End(Not run)
+ }
> 
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "corrmsrs"
> 
> ### * corrmsrs
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: corrmsrs
> ### Title: Correlation Matrix
> ### Aliases: corrmsrs
> ### Keywords: misc
> 
> ### ** Examples
> 
> 
> # matrix of hypothetical prognostic factors
> 
> X <- matrix(rnorm(30), ncol=3) 
> 
> # this function
> 
> print(system.time(a <- corrmsrs(X, minprop=0, maxprop=0.999)))
[1] 0 0 0 0 0
> 
> # coded by just typing the definition of the correlation
> 
> testcorr <- function(X) {
+   wh <- function(cut, x)
+     which(x <= cut)
+   index <- function(x) {
+     ux <- unique(x)
+     ux <- ux[ux < max(ux)]
+     lapply(ux, wh, x = x)
+   }
+   a <- unlist(test <- apply(X, 2, index), recursive=FALSE)
+   cnull <- rep(0, nrow(X))
+   mycorr <- diag(length(a))
+   for (i in 1:(length(a)-1)) {
+     for (j in (i+1):length(a)) {
+       cone <- cnull
+       cone[a[[i]]] <- 1
+       ctwo <- cnull
+       ctwo[a[[j]]] <- 1
+       sone <- sqrt(sum((cone - mean(cone))^2))
+       stwo <- sqrt(sum((ctwo - mean(ctwo))^2))
+       tcorr <- sum((cone - mean(cone))*(ctwo - mean(ctwo)))
+       tcorr <- tcorr/(sone * stwo)
+       mycorr[i,j] <- tcorr
+     }
+   }
+   mycorr
+ }
> 
> print(system.time(tc <- testcorr(X)))
[1] 0.07 0.00 0.07 0.00 0.00
> tc <- tc + t(tc)
> diag(tc) <- 1
> stopifnot(all.equal(tc, a))
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "hohnloser"
> 
> ### * hohnloser
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hohnloser
> ### Title: Left ventricular ejection fraction of patients with malignant
> ###   ventricular tachyarrhythmias.
> ### Aliases: hohnloser
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> data(hohnloser)
> 
> # limiting distribution
> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser, 
+ smethod="LogRank", pmethod="Lau92")

	Maximally selected LogRank statistics using Lau92

data:  Surv(month, cens) by EF 
M = 3.5691, p-value = 0.01065
sample estimates:
estimated cutpoint 
                39 

> 
> # with integer valued scores for comparison
> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser, 
+ smethod="LogRank", pmethod="Lau92", iscores=TRUE)

	Maximally selected LogRank statistics using Lau92

data:  Surv(month, cens) by EF 
M = 3.5639, p-value = 0.01083
sample estimates:
estimated cutpoint 
                39 

> 
> # improved Bonferroni inequality
> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser,
+ smethod="LogRank", pmethod="Lau94")

	Maximally selected LogRank statistics using Lau94

data:  Surv(month, cens) by EF 
M = 3.5691, p-value = 0.005453
sample estimates:
estimated cutpoint 
                39 

> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser,
+ smethod="LogRank", pmethod="Lau94", iscores=TRUE)

	Maximally selected LogRank statistics using Lau94

data:  Surv(month, cens) by EF 
M = 3.5639, p-value = 0.005556
sample estimates:
estimated cutpoint 
                39 

> 
> # small sample solution by Hothorn & Lausen
> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser,
+ smethod="LogRank", pmethod="HL")

	Maximally selected LogRank statistics using HL

data:  Surv(month, cens) by EF 
M = 3.5639, p-value = 0.00667
sample estimates:
estimated cutpoint 
                39 

> 
> # normal approximation
> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser,
+ smethod="LogRank", pmethod="exactGauss")
Warning in pexactgauss(STATISTIC, x, minprop, maxprop, ...) : 
	 Completion with error > abseps

	Maximally selected LogRank statistics using exactGauss

data:  Surv(month, cens) by EF 
M = 3.5691, p-value = 0.003808
sample estimates:
estimated cutpoint 
                39 

> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser,
+ smethod="LogRank", pmethod="exactGauss", iscores=TRUE)

	Maximally selected LogRank statistics using exactGauss

data:  Surv(month, cens) by EF 
M = 3.5639, p-value = 0.002379
sample estimates:
estimated cutpoint 
                39 

> 
> # conditional Monte-Carlo
> 
> maxstat.test(Surv(month, cens) ~ EF, data=hohnloser,
+ smethod="LogRank", pmethod="condMC", B = 9999)

	Maximally selected LogRank statistics using condMC

data:  Surv(month, cens) by EF 
M = 3.5691, p-value = 0.0042
sample estimates:
estimated cutpoint 
                39 

> 
> 
> 
> 
> cleanEx(); ..nameEx <- "maxstat.test"
> 
> ### * maxstat.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: maxstat.test
> ### Title: Maximally Selected Rank and Statistics
> ### Aliases: maxstat.test maxstat.test.data.frame maxstat.test.default
> ###   maxstat
> ### Keywords: htest
> 
> ### ** Examples
> 
> 
> x <- sort(runif(20))
> y <- c(rnorm(10), rnorm(10, 2))
> mydata <- data.frame(cbind(x,y))
> 
> mod <- maxstat.test(y ~ x, data=mydata, smethod="Wilcoxon", pmethod="HL",
+                     minprop=0.25, maxprop=0.75, alpha=0.05)
> print(mod)

	Maximally selected Wilcoxon statistics using HL

data:  y by x 
M = 2.9481, p-value = 0.01603
sample estimates:
estimated cutpoint 
         0.5728534 

> plot(mod)
> 
> # adjusted for more than one prognostic factor.
> 
> data(DLBCL)
> 
> mstat <- maxstat.test(Surv(time, cens) ~ IPI + MGE, data=DLBCL, 
+                       smethod="LogRank", pmethod="exactGauss", 
+                       abseps=0.01)
Warning in pexactgauss(STATISTIC, x, minprop, maxprop, ...) : 
	 Completion with error > abseps
Warning in maxstat(y = c(77.4, 3.4, 71.3, 69.6, 51.2, 3.2, 8.3, 102.4, 89.8,  : 
	 pvmnorm:  Completion with error > abseps
> plot(mstat)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "pLausen92"
> 
> ### * pLausen92
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pLausen92
> ### Title: Approximating Maximally Selected Statistics
> ### Aliases: pLausen92 qLausen92
> ### Keywords: distribution
> 
> ### ** Examples
> 
> 
> # Compute quantiles. Should be equal to Table 2 in Lausen and Schumacher
> 
> data(LausenTab2)
> 
> a <- rev(c(0.01, 0.025, 0.05, 0.1))
> prop <- rbind(c(0.25, 0.75), c(0.4, 0.6), c(0.1, 0.9), c(0.4, 0.9))
> Quant <- matrix(rep(0, length(a)*nrow(prop)), nrow=length(a)) 
> 
> for (i in 1:length(a)) {                                            
+   for (j in 1:nrow(prop)) {                            
+     Quant[i,j] <- qLausen92(a[i], minprop=prop[j,1], maxprop=prop[j,2]) 
+   }
+ }
> 
> Quant <- round(Quant, 3)
> rownames(Quant) <- a
> colnames(Quant) <- c("A2575", "A46", "A19", "A49")
> Quant <- as.data.frame(Quant)
> rownames(LausenTab2) <- a
> 
> Quant
      A2575   A46   A19   A49
0.1   2.539 2.263 2.784 2.597
0.05  2.829 2.560 3.054 2.883
0.025 3.088 2.828 3.297 3.138
0.01  3.395 3.149 3.588 3.442
> 
> LausenTab2
      A2575   A46   A19   A49
0.1   2.539 2.263 2.784 2.597
0.05  2.829 2.560 3.054 2.883
0.025 3.088 2.828 3.297 3.138
0.01  3.395 3.149 3.588 3.442
> 
> if(!all.equal(LausenTab2, Quant)) stop("error checking pLausen92")
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "pLausen94"
> 
> ### * pLausen94
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pLausen94
> ### Title: Approximating Maximally Selected Statistics
> ### Aliases: pLausen94 qLausen94
> ### Keywords: distribution
> 
> ### ** Examples
> 
> 
> p <- pLausen94(2.5, 20, 0.25, 0.75)
> 
> # Lausen 94, page 489
> 
> if (round(p, 3) != 0.073) stop("error checking pLausen94")
> 
> # the same
> 
> p2 <- pLausen94(2.5, 200, 0.25, 0.75, m=seq(from=50, to=150, by=10))
> 
> stopifnot(all.equal(round(p,3), round(p2,3)))
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "pexactgauss"
> 
> ### * pexactgauss
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pexactgauss
> ### Title: Computing Maximally Selected Gauss Statistics
> ### Aliases: pexactgauss qexactgauss
> ### Keywords: distribution
> 
> ### ** Examples
> 
> 
> x <- rnorm(20)
> 
> pexact <- pexactgauss(2.5, x, abseps=0.01)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.maxtest"
> 
> ### * plot.maxtest
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.maxtest
> ### Title: Print and Plot Standardized Statistics
> ### Aliases: plot.maxtest print.maxtest plot.mmaxtest print.mmaxtest
> ### Keywords: htest
> 
> ### ** Examples
> 
> 
> x <- sort(runif(20))
> y <- rbinom(20, 1, 0.5)
> mydata <- data.frame(c(x,y))
> 
> mod <- maxstat.test(y ~ x, data=mydata, smethod="Median", 
+                     pmethod="HL", alpha=0.05)
> print(mod)

	Maximally selected Median statistics using HL

data:  y by x 
M = 1.5533, p-value = 1
sample estimates:
estimated cutpoint 
         0.7176185 

> plot(mod)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "pmaxstat"
> 
> ### * pmaxstat
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pmaxstat
> ### Title: Approximating Maximally Selected Statistics
> ### Aliases: pmaxstat qmaxstat
> ### Keywords: distribution
> 
> ### ** Examples
> 
> 
> pmaxstat(2.5, 1:20, 5:15)
[1] 0.09403538
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "sphase"
> 
> ### * sphase
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sphase
> ### Title: S-phase fraction of tumor cells
> ### Aliases: sphase
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(sphase)
> maxstat.test(Surv(RFS, cens) ~ SPF, data=sphase, smethod="LogRank",
+ pmethod="Lau94")

	Maximally selected LogRank statistics using Lau94

data:  Surv(RFS, cens) by SPF 
M = 2.4033, p-value = 0.2727
sample estimates:
estimated cutpoint 
               107 

> maxstat.test(Surv(RFS, cens) ~ SPF, data=sphase, smethod="LogRank",
+ pmethod="Lau94", iscores=TRUE)

	Maximally selected LogRank statistics using Lau94

data:  Surv(RFS, cens) by SPF 
M = 2.4017, p-value = 0.2738
sample estimates:
estimated cutpoint 
               107 

> maxstat.test(Surv(RFS, cens) ~ SPF, data=sphase, smethod="LogRank",
+ pmethod="HL")

	Maximally selected LogRank statistics using HL

data:  Surv(RFS, cens) by SPF 
M = 2.4017, p-value = 0.5494
sample estimates:
estimated cutpoint 
               107 

> maxstat.test(Surv(RFS, cens) ~ SPF, data=sphase, smethod="LogRank",
+ pmethod="condMC", B = 9999)

	Maximally selected LogRank statistics using condMC

data:  Surv(RFS, cens) by SPF 
M = 2.4033, p-value = 0.1534
sample estimates:
estimated cutpoint 
               107 

> plot(maxstat.test(Surv(RFS, cens) ~ SPF, data=sphase, smethod="LogRank"))
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "treepipit"
> 
> ### * treepipit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: treepipit
> ### Title: Tree Pipit Data
> ### Aliases: treepipit
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(treepipit)
> mod <- maxstat.test(counts ~ coverstorey, data = treepipit, 
+                     smethod = "Data", pmethod = "HL", minprop = 0.2,
+                     maxprop = 0.8)
> print(mod)

	Maximally selected Data statistics using HL

data:  counts by coverstorey 
M = 4.3139, p-value = 0.0001141
sample estimates:
estimated cutpoint 
                40 

> plot(mod)
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  10.58 0.72 11.49 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')