Nothing
## ----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)
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.