### R code from vignette source 'vignettes/edd/inst/doc/HOWTO-edd.Rnw' ################################################### ### code chunk number 1: getGolub ################################################### library(edd) library(golubEsets) library(xtable) data(Golub_Merge) ################################################### ### code chunk number 2: filterGolub ################################################### madvec <- apply(exprs(Golub_Merge),1,mad) minvec <- apply(exprs(Golub_Merge),1,min) keep <- (madvec > median(madvec)) & (minvec > 300) gmfilt <- Golub_Merge[keep==TRUE,] ################################################### ### code chunk number 3: splitGolub ################################################### ALL <- gmfilt$ALL.AML=="ALL" gall <- gmfilt[,ALL==TRUE] gaml <- gmfilt[,ALL==FALSE] show(gall) ################################################### ### code chunk number 4: eddByType ################################################### set.seed(12345) alldists <- edd(gall, meth="nnet", size=10, decay=.2) amldists <- edd(gaml, meth="nnet", size=10, decay=.2) ################################################### ### code chunk number 5: HOWTO-edd.Rnw:105-107 ################################################### greo <- match(sapply(eddDistList,tag),names(table(alldists))) if (any(ii <- is.na(greo))) greo = greo[-which(ii)] ################################################### ### code chunk number 6: HOWTO-edd.Rnw:112-116 ################################################### gn <- row.names(exprs(gall)) alld <- as.character(alldists)[1:5] names(alld) <- gn[1:5] print(alld) ################################################### ### code chunk number 7: lookKNN ################################################### set.seed(123) alldistsKNN <- edd(gall, meth="knn", k=1, l=0) alldistsTEST <- edd(gall, meth="test", thresh=.3) ################################################### ### code chunk number 8: compareKNNvNN ################################################### cap <- "Comparison of distribution shape classification by nnet (rows) and by knn (columns) methods in edd." print(try(xtable( latEDtable(table(alldists,alldistsKNN),reorder=greo),digits=rep(0,length(table(alldists))+1),caption=cap, label="conc1"))) ################################################### ### code chunk number 9: testTable ################################################### print(table(alldistsTEST)) ################################################### ### code chunk number 10: ALLtable ################################################### cap <- "Frequencies of distributional shapes in filtered ALL data." print(xtable(latEDtable( table(alldists),reorder=greo),digits=rep(0,length(table(alldists))+1),caption=cap,label="marg1")) ################################################### ### code chunk number 11: compare ################################################### mymat <- matrix(c(1,2),nr=1) layout(mymat)#,heights=c(lcm(6),lcm(6)),widths=c(lcm(6),lcm(6)),respect=TRUE) barplot(table(alldists),las=2) barplot(table(amldists),las=2) ################################################### ### code chunk number 12: discordTable ################################################### cap <- "Rows are gene-specific distribution shapes for ALL, columns for AML, and cell entries are counts of genes." print(xtable(latEDtable(table(alldists,amldists),reord=greo),cap=cap, label="disco1")) ################################################### ### code chunk number 13: findBiMod ################################################### print((1:540)[alldists==".75N(0,1)+.25N(4,1)" & amldists=="N(0,1)"][1:5]) ################################################### ### code chunk number 14: lookD87593 ################################################### par(mfrow=c(2,2)) plotED(MIXN1,data=exprs(gall)[65,]) hist(exprs(gall)[65,],xlim=c(0,3500)) plotED(N01,data=exprs(gaml)[65,]) hist(exprs(gaml)[65,],xlim=c(0,3500)) ################################################### ### code chunk number 15: namesEDL ################################################### names(eddDistList) ################################################### ### code chunk number 16: plotRefs ################################################### par(mfrow=c(4,2)) for (i in 1:8) plotED(eddDistList[[i]]) ################################################### ### code chunk number 17: addMIXN3 ################################################### #> median(rmixnorm(10000,.75,0,1,6,1)) #[1] 0.4337298 #> mad(rmixnorm(10000,.75,0,1,6,1)) #[1] 1.550768 -- these are resistant stats! MIXN3 <- new("eddDist", stub="mixnorm", parms=c(p1=.75,m1=0,s1=1,m2=6,s2=1), median=.43, mad=1.55, tag=".75N(0,1)+.25N(6,1)", plotlim=c(-3,11), latexTag="$\\frac{3}{4}\\Phi+\\frac{1}{4}\\Phi_{6,1}$") eddDistList[["MIXN3"]] <- MIXN3 set.seed(12345) alldists2 <- edd(gall, meth="nnet", size=10, decay=.2) print(alldists2[65]) ################################################### ### code chunk number 18: checkFit3 ################################################### plotED(MIXN3,data=exprs(gall)[65,]) ################################################### ### code chunk number 19: checkFit1 ################################################### plotED(MIXN1,data=exprs(gall)[65,])