inst/doc/DirichletMultinomial.R

In DirichletMultinomial: Dirichlet-Multinomial Mixture Model Machine Learning for Microbiome Data

### R code from vignette source 'DirichletMultinomial.Rnw'

###################################################
### code chunk number 1: library
###################################################
library(DirichletMultinomial)
library(lattice)
library(xtable)
library(parallel)


###################################################
### code chunk number 2: colors
###################################################
options(width=70, digits=2)
full <- FALSE
.qualitative <- DirichletMultinomial:::.qualitative
dev.off <- function(...) invisible(grDevices::dev.off(...))


###################################################
### code chunk number 3: data-input
###################################################
fl <- system.file(package="DirichletMultinomial", "extdata",
                  "Twins.csv")
count <- t(as.matrix(read.csv(fl, row.names=1)))
count[1:5, 1:3]


###################################################
### code chunk number 4: taxon-counts
###################################################
cnts <- log10(colSums(count))
pdf("taxon-counts.pdf")
densityplot(cnts, xlim=range(cnts),
            xlab="Taxon representation (log 10 count)")
dev.off()


###################################################
### code chunk number 5: fit
###################################################
if (full) {
    fit <- mclapply(1:7, dmn, count=count, verbose=TRUE)
    save(fit, file=file.path(tempdir(), "fit.rda"))
} else data(fit)
fit[[4]]


###################################################
### code chunk number 6: min-laplace
###################################################
lplc <- sapply(fit, laplace)
pdf("min-laplace.pdf")
plot(lplc, type="b", xlab="Number of Dirichlet Components",
     ylab="Model Fit")
dev.off()
(best <- fit[[which.min(lplc)]])


###################################################
### code chunk number 7: mix-weight
###################################################
mixturewt(best)
head(mixture(best), 3)


###################################################
### code chunk number 8: fitted
###################################################
pdf("fitted.pdf")
splom(log(fitted(best)))
dev.off()


###################################################
### code chunk number 9: isoMDS
###################################################



###################################################
### code chunk number 10: isoMDS-plot
###################################################



###################################################
### code chunk number 11: posterior-mean-diff
###################################################
p0 <- fitted(fit[[1]], scale=TRUE)     # scale by theta
p4 <- fitted(best, scale=TRUE)
colnames(p4) <- paste("m", 1:4, sep="")
(meandiff <- colSums(abs(p4 - as.vector(p0))))
sum(meandiff)


###################################################
### code chunk number 12: table-1
###################################################
diff <- rowSums(abs(p4 - as.vector(p0)))
o <- order(diff, decreasing=TRUE)
cdiff <- cumsum(diff[o]) / sum(diff)
df <- head(cbind(Mean=p0[o], p4[o,], diff=diff[o], cdiff), 10)


###################################################
### code chunk number 13: xtable
###################################################
xtbl <- xtable(df,
    caption="Taxonomic contributions (10 largest) to Dirichlet components.",
    label="tab:meandiff", align="lccccccc")
print(xtbl, hline.after=0, caption.placement="top")


###################################################
### code chunk number 14: heatmap-similarity
###################################################
pdf("heatmap1.pdf")
heatmapdmn(count, fit[[1]], best, 30)
dev.off()


###################################################
### code chunk number 15: twin-pheno
###################################################
fl <- system.file(package="DirichletMultinomial", "extdata",
                  "TwinStudy.t")
pheno0 <- scan(fl)
lvls <- c("Lean", "Obese", "Overwt")
pheno <- factor(lvls[pheno0 + 1], levels=lvls)
names(pheno) <- rownames(count)
table(pheno)


###################################################
### code chunk number 16: subsets
###################################################
counts <- lapply(levels(pheno), csubset, count, pheno)
sapply(counts, dim)
keep <- c("Lean", "Obese")
count <- count[pheno %in% keep,]
pheno <- factor(pheno[pheno %in% keep], levels=keep)


###################################################
### code chunk number 17: fit-several-
###################################################
if (full) {
    bestgrp <- dmngroup(count, pheno, k=1:5, verbose=TRUE, 
                        mc.preschedule=FALSE)
    save(bestgrp, file=file.path(tempdir(), "bestgrp.rda"))
} else data(bestgrp)


###################################################
### code chunk number 18: best-several
###################################################
bestgrp
lapply(bestgrp, mixturewt)
c(sapply(bestgrp, laplace),
  `Lean+Obese`=sum(sapply(bestgrp, laplace)),
  Single=laplace(best))


###################################################
### code chunk number 19: confusion
###################################################
xtabs(~pheno + predict(bestgrp, count, assign=TRUE))


###################################################
### code chunk number 20: cross-validate
###################################################
if (full) {
    ## full leave-one-out; expensive!
    xval <- cvdmngroup(nrow(count), count, c(Lean=1, Obese=3), pheno,
                       verbose=TRUE, mc.preschedule=FALSE)
    save(xval, file=file.path(tempdir(), "xval.rda"))
} else data(xval)


###################################################
### code chunk number 21: ROC-dmngroup
###################################################
bst <- roc(pheno[rownames(count)] == "Obese",
           predict(bestgrp, count)[,"Obese"])
bst$Label <- "Single"
two <- roc(pheno[rownames(xval)] == "Obese",
           xval[,"Obese"])
two$Label <- "Two group"
both <- rbind(bst, two)
pars <- list(superpose.line=list(col=.qualitative[1:2], lwd=2))
pdf("roc.pdf")
xyplot(TruePostive ~ FalsePositive, group=Label, both,
       type="l", par.settings=pars,
       auto.key=list(lines=TRUE, points=FALSE, x=.6, y=.1),
       xlab="False Positive", ylab="True Positive")
dev.off()


###################################################
### code chunk number 22: sessionInfo
###################################################
toLatex(sessionInfo())