inst/doc/gMCP.R

In gMCP: Graph Based Multiple Comparison Procedures

## ----OptionsAndLibraries, include=FALSE, message=FALSE------------------------------------------------------------------------------------
# knitr has to be loaded for 'set_parent' and CRAN checks and also for opt_chunk during build process.
library(knitr)
if (exists("opts_chunk")) {
  opts_chunk$set(concordance=TRUE)
  opts_chunk$set(tidy.opts=list(keep.blank.line=FALSE, width.cutoff=95))
  opts_chunk$set(size="footnotesize")
  opts_chunk$set(cache=TRUE)
  opts_chunk$set(autodep=TRUE)
}
library(gMCP, quietly=TRUE)
options(width=140)
options(digits=4)
gMCPReport <- function(...) {invisible(NULL)}
graphGUI <- function(...) {invisible(NULL)}
#options(prompt="> ", continue="+ ")


## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------
graph <- BonferroniHolm(3)
cat(graph2latex(graph, scale=0.7, labelTikZ="near start,fill=blue!20", scaleText=FALSE))

## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------
graph <- BonferroniHolm(3)
cat(graph2latex(graph, scale=0.7, nodeTikZ="minimum size=1.2cm", scaleText=FALSE))
cat("\\\\$\\downarrow$ reject $H_1$\\\\")
graph <- rejectNode(graph, "H1")
cat(graph2latex(graph, scale=0.7, nodeTikZ="minimum size=1.2cm", scaleText=FALSE))
cat("\\\\$\\downarrow$ reject $H_3$\\\\\\ \\\\")
graph <- rejectNode(graph, "H3")
cat(graph2latex(graph, scale=0.7, nodeTikZ="minimum size=1.2cm", scaleText=FALSE))

## ----echo=TRUE, eval=FALSE----------------------------------------------------------------------------------------------------------------
#  library(gMCP)
#  graphGUI()

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
library(gMCP)
graph <- BonferroniHolm(3)
gMCP(graph, pvalues=c(0.01,0.07,0.02), alpha=0.05)

## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------
graph <- BretzEtAl2011()
cat(graph2latex(graph, scale=0.7, scaleText=FALSE))

## ----echo=TRUE, tidy=FALSE----------------------------------------------------------------------------------------------------------------
m <- rbind(H11=c(0,   0.5, 0,   0.5, 0,   0  ),
           H21=c(1/3, 0,   1/3, 0,   1/3, 0  ),
           H31=c(0,   0.5, 0,   0,   0,   0.5),
           H12=c(0,   1,   0,   0,   0,   0  ),
           H22=c(0.5, 0,   0.5, 0,   0,   0  ),
           H32=c(0,   1,   0,   0,   0,   0  ))

graph <- matrix2graph(m)
graph <- setWeights(graph, c(1/3, 1/3, 1/3, 0, 0, 0))

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
print(graph)

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
graph@nodeAttr$X <- c(H11=100, H21=300, H31=500, H12=100, H22=300, H32=500)
graph@nodeAttr$Y <- c(H11=100, H21=100, H31=100, H12=300, H22=300, H32=300)	

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
graph <- placeNodes(graph, nrow=2)

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
cat(graph2latex(graph))

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
edgeAttr(graph, "H11", "H21", "labelX") <- 200
edgeAttr(graph, "H11", "H21", "labelY") <- 80

## ----echo=TRUE----------------------------------------------------------------------------------------------------------------------------
graphGUI("graph")

## ----echo=TRUE,tidy=FALSE-----------------------------------------------------------------------------------------------------------------
graph <- BretzEtAl2011()

# We can reject a single node:
print(rejectNode(graph, "H11"))

# Or given a vector of pvalues let the function gMCP do all the work:  
pvalues <- c(0.1, 0.008, 0.005, 0.15, 0.04, 0.006)
result <- gMCP(graph, pvalues)
print(result)

## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------
cat(graph2latex(result@graphs[[4]], scale=0.7, scaleText=FALSE))

## ----echo=TRUE, eval=FALSE----------------------------------------------------------------------------------------------------------------
#  gMCPReport(result, "Report.tex")

## ----include=FALSE------------------------------------------------------------------------------------------------------------------------
d1 <- c(1.68005156523844, 1.95566697423700, 0.00137860945822299, 0.660052238622464, 
		1.06731835721526, 0.39479303427265, -0.312462050794408, 0.323637755662837, 
		0.490976552328251, 2.34240774442652)

d2 <- c(0.507878380203451, 1.60461475524144, 2.66959621483759, 0.0358289240280020, 
		-1.13014087491324, 0.792461583741794, 0.0701657425268248, 3.15360436883856, 
		0.217669661552567, 1.23979492014026)

d3 <- c(-1.31499425534849, 1.62201370145649, 0.89391826766116, 0.845473572033649, 
		2.17912435223573, 1.07521368050267, 0.791598289847664, 1.58537210294519, 
		-0.079778759456515, 0.97295072606043)

est <- c(0.860382, 0.9161474, 0.9732953)
s <- c(0.8759528, 1.291310, 0.8570892)
pval <- c(0.01260, 0.05154, 0.02124)/2

df <- 9
# Statistics:
st <- qt(pval/2, df=df, lower.tail=FALSE)
# Estimates:
est <- st*s/sqrt(10)

## ----confint, echo=TRUE, tidy=FALSE-------------------------------------------------------------------------------------------------------
# Estimates:
est <- c("H1"=0.860382, "H2"=0.9161474, "H3"=0.9732953)
# Sample standard deviations:
ssd <- c("H1"=0.8759528, "H2"=1.291310, "H3"=0.8570892)

pval <- c(0.01260, 0.05154, 0.02124)/2

simConfint(BonferroniHolm(3), pvalues=pval, 
		confint=function(node, alpha) {
			c(est[node]-qt(1-alpha,df=9)*ssd[node]/sqrt(10), Inf)
		}, estimates=est, alpha=0.025, mu=0, alternative="greater")

# Note that the sample standard deviations in the following call
# will be calculated from the pvalues and estimates.
simConfint(BonferroniHolm(3), pvalues=pval, 
		confint="t", df=9, estimates=est, alpha=0.025, alternative="greater")

## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------

cat(graph2latex(parallelGatekeeping(), nodeTikZ="minimum size=1.2cm", tikzEnv=FALSE))
cat(graph2latex(improvedParallelGatekeeping(), nodeTikZ="minimum size=1.2cm", tikzEnv=FALSE, offset=c(300, 0), nodeR=27))


## ----echo=TRUE, tidy=FALSE----------------------------------------------------------------------------------------------------------------
m <- rbind(H1=c(0,           0,           0.5,           0.5          ),
           H2=c(0,           0,           0.5,           0.5          ),
           H3=c("\\epsilon", 0,           0,             "1-\\epsilon"),
           H4=c(0,           "\\epsilon", "1-\\epsilon", 0            ))

graph <- matrix2graph(m)
#graph <- improvedParallelGatekeeping()
graph
substituteEps(graph, eps=0.001)

gMCP(graph, pvalues=c(0.02, 0.04, 0.01, 0.02), eps=0.001)

## ----tidy=FALSE---------------------------------------------------------------------------------------------------------------------------
m <- rbind(H1=c(0, 0, 1, 0, 0),
           H2=c(0, 0, 1, 0, 0),
           H3=c(0, 0, 0, 0.9999, 1e-04),
           H4=c(0, 1, 0, 0, 0),
           H5=c(0, 0, 0, 0, 0))
weights <- c(1, 0, 0, 0, 0)
subgraph1 <- new("graphMCP", m=m, weights=weights)

m <- rbind(H1=c(0, 0, 1, 0, 0),
           H2=c(0, 0, 1, 0, 0),
           H3=c(0, 0, 0, 1e-04, 0.9999),
           H4=c(0, 0, 0, 0, 0),
           H5=c(1, 0, 0, 0, 0))
weights <- c(0, 1, 0, 0, 0)
subgraph2 <- new("graphMCP", m=m, weights=weights)

weights <- c(0.5, 0.5)
graph <- new("entangledMCP", subgraphs=list(subgraph1, subgraph2), weights=weights)

## ----eval=TRUE, tidy=FALSE----------------------------------------------------------------------------------------------------------------

cr <- rbind(H1=c(1   , 0.5 , 0.3 , 0.15),
            H2=c(0.5 , 1   , 0.15, 0.3 ),
            H3=c(0.3 , 0.15, 1   , 0.5 ),
            H4=c(0.15, 0.3 , 0.5 , 1   ))


## ----eval=FALSE---------------------------------------------------------------------------------------------------------------------------
#  
#  library(bindata)
#  
#  n1 <- 20
#  n2 <- 1e4
#  
#  pvals <- t(replicate(n2,
#                       sapply(colSums(rmvbin(n1, margprob = c(0.35, 0.4, 0.25, 0.3), bincorr = cr)),
#                              function(x, ...) {binom.test(x, ...)$p.value}, n=n1, alternative="less")
#                       ))
#  

## ----include=FALSE------------------------------------------------------------------------------------------------------------------------
load("pvals.RData")

## ----eval=TRUE----------------------------------------------------------------------------------------------------------------------------

graph <- generalSuccessive(gamma=0, delta=0)
out <- graphTest(pvalues=pvals, graph = graph)
extractPower(out)


## ----eval=FALSE---------------------------------------------------------------------------------------------------------------------------
#  # Other example:
#  
#  library(copula)
#  cop <- mvdc(normalCopula(0.75, dim=4), c("norm", "norm", "exp", "exp"), list(list(mean=0.5), list(mean=1), list(rate = 2), list(rate = 3)))
#  x <- rMvdc(250, cop)
#  # ...
#  
#  # Boot example?
#  # Permutation test example?
#  #

## ----echo=TRUE, tidy=FALSE----------------------------------------------------------------------------------------------------------------

# Bonferroni adjustment
G <- diag(2)
weights <- c(0.5,0.5)
corMat <- diag(2)+matrix(1,2,2)
theta <- c(1,2)
calcPower(weights, alpha=0.025, G, theta, corMat)
calcPower(weights, alpha=0.025, G, 2*theta, 2*corMat)


## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------

cat(graph2latex(generalSuccessive(), scale=0.7, scaleText=FALSE))


## ----echo=TRUE, tidy=FALSE----------------------------------------------------------------------------------------------------------------
graph <- generalSuccessive()
graph

## ----echo=FALSE,results='asis'------------------------------------------------------------------------------------------------------------

cat(graph2latex(result@graphs[[3]], pvalues=pvalues, scale=0.7, scaleText=FALSE))


## ----echo=TRUE, message=FALSE, warning=FALSE----------------------------------------------------------------------------------------------
data(hydroquinone)

pvalues <- c() 
x <- hydroquinone$micronuclei[hydroquinone$group=="C-"]
for (dose in c("30 mg/kg", "50 mg/kg", "75 mg/kg", "100 mg/kg", "C+")) {
	y <- hydroquinone$micronuclei[hydroquinone$group==dose]
	result <- wilcox.test(x, y, alternative="less", correct=TRUE)
	pvalues <- c(result$p.value,  pvalues)
}
pvalues

library(coin, quietly=TRUE)
pvalues <- c() 
for (dose in c("30 mg/kg", "50 mg/kg", "75 mg/kg", "100 mg/kg", "C+")) {
	subdata <- droplevels(hydroquinone[hydroquinone$group %in% c("C-", dose),])
	result <- wilcox_test(micronuclei ~ group, data=subdata, distribution="exact")
	pvalues <- c(pvalue(result), pvalues)
}
pvalues

## ----echo=FALSE,fig.width=8,fig.height=6--------------------------------------------------------------------------------------------------
data(hydroquinone)
boxplot(micronuclei~group, data=hydroquinone)

## ----echo=FALSE, results='asis'-----------------------------------------------------------------------------------------------------------

graphs <- list(BonferroniHolm(4),
               parallelGatekeeping(),
               improvedParallelGatekeeping(),
               BretzEtAl2011(),
               HungEtWang2010(),
               HuqueAloshEtBhore2011(),
               HommelEtAl2007(),
               HommelEtAl2007Simple(),
               MaurerEtAl1995(),
               improvedFallbackI(weights=rep(1/3, 3)),
               improvedFallbackII(weights=rep(1/3, 3)),
               #cycleGraph(nodes=paste("H",1:4,sep=""), weights=rep(1/4, 4)),
               fixedSequence(5),
               fallback(weights=rep(1/4, 4)),
               generalSuccessive(weights = c(1/2, 1/2)),
               simpleSuccessiveI(),
               simpleSuccessiveII(),
               truncatedHolm(),
               BauerEtAl2001(),
               BretzEtAl2009a(),
               BretzEtAl2009b(),
               BretzEtAl2009c(),
               FerberTimeDose2011(times=5, doses=3, w=1/2),
               Ferber2011(),
               WangTing2014()
               #Entangled1Maurer2012(),
               #Entangled2Maurer2012()
)

make.url <- function (x) {
    return(gsub("(https?://[^ ]+)", "\\\\url{\\1}", x))
}

texify <- function(x, linebreak="\n\n") {
  x <- make.url(x)
  x <- gsub(linebreak, "\\\\\\\\\n", x)
  x <- gsub("&", "\\\\&", x)
  x <- gsub("\\\\epsilon", "$\\\\epsilon$", x)
  x <- gsub("\\\\tau", "$\\\\tau$", x)
  x <- gsub("\\\\nu", "$\\\\nu$", x)
  return(x)
}

# first.line("a\nb\nc")
first.line <- function(x) {
  return(strsplit(x, split="\n")[[1]][1])
}

count <-0
for (g in graphs) {
  descr <- attr(g, "descr")
  cat(graph2latex(g, fig=TRUE, scaleText=TRUE, scale=0.7, fig.caption=texify(descr, linebreak = "\n"), fig.caption.short=texify(first.line(descr))))
  count <- count + 1
  if(count%%6==0) cat("\\cleardoublepage\n")
}