inst/doc/MultipleComparisons.R

In PairViz: Visualization using Graph Traversal

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(digits=4)

## -----------------------------------------------------------------------------
suppressPackageStartupMessages(library(PairViz))
data(cancer)  # Need this step to load the data
str(cancer)   # Summary of structure of the data
# We can separate the survival times by which organ is affected
organs <- with(cancer, split(Survival, Organ))
# And record their names for use later!
organNames <- names(organs)
# the structure of the organs data
str(organs)

## ---- fig.align="center", fig.width=6, fig.height=5---------------------------

library(colorspace)
cols <- rainbow_hcl(5, c = 50)  # choose chromaticity of 50 to dull colours
boxplot(organs, col=cols, 
        ylab="Survival time", 
        main="Cancer treated by vitamin C")

## ---- fig.align="center", fig.width=6, fig.height=5---------------------------
# Split the data
sqrtOrgans <- with(cancer, split(sqrt(Survival), Organ))
boxplot(sqrtOrgans, col=cols, 
        ylab=expression(sqrt("Survival time")), 
        main="Cancer treated by vitamin C")

## -----------------------------------------------------------------------------
ord <- eulerian(5)
ord

## ---- fig.align="center", fig.width=7.5, fig.height=5-------------------------
boxplot(sqrtOrgans[ord], col=cols[ord], 
        ylab=expression(sqrt("Survival time")), 
        main="Cancer treated by vitamin C", cex.axis=.6)

## -----------------------------------------------------------------------------
ordHam <-  hpaths(5, matrix = FALSE)
ordHam

## -----------------------------------------------------------------------------
# Get the test results
test <- with(cancer,
             pairwise.t.test(sqrt(Survival), Organ))
pvals <- test$p.value
pvals

## -----------------------------------------------------------------------------
# First construct a vector, removing NAs.
weights <- pvals[!is.na(pvals)]
weights <-edge2dist(weights)

## -----------------------------------------------------------------------------
weights <- as.matrix(weights)
rownames(weights) <- organNames
colnames(weights)<- rownames(weights)
weights

## -----------------------------------------------------------------------------
g <- mk_complete_graph(weights)

## ----fig.width=5, fig.height=5, fig.align='center'----------------------------
requireNamespace("igraph")
igplot <- function(g,weights=FALSE,layout=igraph::layout_in_circle, 
                   vertex.size=60, vertex.color="lightblue",...){
    g <- igraph::graph_from_graphnel(as(g, "graphNEL"))
    op <- par(mar=c(1,1,1,1))
    if (weights){
      ew <- round(igraph::get.edge.attribute(g,"weight"),3) 
      igraph::plot.igraph(g,layout=layout,edge.label=ew,vertex.size=vertex.size,vertex.color=vertex.color,...)
    }
    else
    igraph::plot.igraph(g,layout=layout,vertex.size=vertex.size,vertex.color=vertex.color,...)
    par(op)
}
igplot(g,weights=TRUE,edge.label.color="black")

## ----fig.width=5, fig.height=5,fig.align='center', eval=FALSE-----------------
#  library(Rgraphviz)
#  ew <- round(unlist(edgeWeights(g)),3)
#  ew <- ew[setdiff(seq(along=ew), removedEdges(g))]
#  names(ew) <- edgeNames(g)
#  plot(g,  "circo",edgeAttrs=list(label=ew))

## -----------------------------------------------------------------------------
low2highEulord <- eulerian(weights); colnames(weights)[low2highEulord]
## or equivalently
eulerian(g)

## -----------------------------------------------------------------------------
bestHam <- order_best(weights)
colnames(weights)[bestHam]

## ---- fig.align="center", fig.width=7.5, fig.height=5-------------------------
boxplot(sqrtOrgans[low2highEulord], col=cols[ord], 
        ylab=expression(sqrt("Survival time")), 
        main="Cancer treated by vitamin C", cex.axis=.6)

## ----fig.width=5, fig.height=5, align='center'--------------------------------
aovOrgans <-  aov(sqrt(Survival) ~ Organ,data=cancer)
TukeyHSD(aovOrgans,conf.level = 0.95)

## -----------------------------------------------------------------------------
tuk <-TukeyHSD(aovOrgans,conf.level = 0.95)
ptuk <- tuk$Organ[,"p adj"]
dtuk <- as.matrix(edge2dist(ptuk))
rownames(dtuk)<- colnames(dtuk)<- organNames
g <- mk_complete_graph(weights)
eulerian(dtuk)

## ----fig.width=5, fig.height=5, align='center'--------------------------------
par(mar=c(3,8,3,3))
plot(TukeyHSD(aovOrgans,conf.level = 0.95),las=1,tcl = -.3)

## ----fig.align="center", fig.width=7.5, fig.height=5--------------------------

mc_plot(sqrtOrgans,aovOrgans,main="Pairwise comparisons of cancer types", 
        ylab="Sqrt Survival",col=cols,cex.axis=.6)

## ----fig.align="center", fig.width=7.5, fig.height=5--------------------------

suppressPackageStartupMessages(library(multcomp))
fitVitC <- glht(aovOrgans, linfct = mcp(Organ= "Tukey"))

# this gives confidence intervals without a family correction
confint(fitVitC,level=.99,calpha = univariate_calpha())

# for short, define
cifunction<- function(f, lev) 
  confint(f,level=lev,calpha = univariate_calpha())$confint

# calculate the confidence intervals
conf <- cbind(cifunction(fitVitC,.9),cifunction(fitVitC,.95)[,-1],cifunction(fitVitC,.99)[,-1])


mc_plot(sqrtOrgans,conf,path=low2highEulord,
        main="Pairwise comparisons of cancer types", 
        ylab="Sqrt Survival",col=cols,cex.axis=.6)




## -----------------------------------------------------------------------------
if (!requireNamespace("Sleuth3", quietly = TRUE)){
    install.packages("Sleuth3")
}
library(Sleuth3)
mice <- case0501
str(mice)
levels(mice$Diet)
# get rid of "/"
levels(mice$Diet) <-  c("NN85", "NR40", "NR50", "NP" ,   "RR50" ,"lopro")

## ---- fig.align="center", fig.width=6, fig.height=5---------------------------

life <- with(mice, split(Lifetime ,Diet))
cols <- rainbow_hcl(6, c = 50) 
boxplot(life, col=cols, 
        ylab="Lifetime", 
        main="Diet Restriction and Longevity")

## -----------------------------------------------------------------------------
aovMice   <- aov(Lifetime ~ Diet-1, data=mice)
fitMice <- glht(aovMice,
          linfct=c("DietNR50 - DietNN85 = 0", 
          "DietRR50  - DietNR50 = 0",
          "DietNR40  - DietNR50 = 0",
          "Dietlopro - DietNR50 = 0",
          "DietNN85  - DietNP   = 0")) 
  summary(fitMice,test=adjusted("none")) # No multiple comparison adjust.
  confint(fitMice, calpha = univariate_calpha()) # No adjustment

## -----------------------------------------------------------------------------
g <- new("graphNEL", nodes=names(life))

## ----fig.align='center', fig.width=5, fig.height=5----------------------------
fitMiceSum <- summary(fitMice,test=adjusted("none"))
pvalues <- fitMiceSum$test$pvalues
pvalues
# Extract  labels from the p-values for the edges
edgeLabs <- unlist(strsplit(names(pvalues), " - "))
edgeLabs <- matrix(substring(edgeLabs,5), nrow=2)
g <- addEdge(edgeLabs[1,], edgeLabs[2,], g,pvalues)

pos <- rbind(c(-1,0), c(0,-1), c(0,0), c(-2,0),c(1,0), c(0,1))
igplot(g, weights=TRUE, layout=pos,vertex.size=32)

## -----------------------------------------------------------------------------
eulerian(g)

## ----fig.align='center', fig.width=5, fig.height=5----------------------------
g1 <- addEdge("NR40","NP",g,1)
g1 <- addEdge("lopro","RR50",g1,1)
igplot(g1, weights=TRUE, layout=pos,vertex.size=32)
eulerian(g1)

## ---- fig.align="center", fig.width=6, fig.height=5---------------------------

eul <- eulerian(g1)
# make eul numeric
eul <- match(eul, names(life))

fitMice1 <- glht(aovMice, linfct = mcp(Diet= "Tukey"))

# need to construct the confidence intervals for all pairs
conf <- cbind(cifunction(fitMice1,.9),cifunction(fitMice1,.95)[,-1],cifunction(fitMice1,.99)[,-1])
# these comparisons are not relevant
conf[c(1,4,5,7,8,9,13,14,15),]<- NA

mc_plot(life,conf,path=eul,
        main="Diet Restriction and Longevity", 
        ylab="Lifetime",col=cols,cex.axis=.6)