### R code from vignette source 'vignettes/SSPA/inst/doc/SSPA.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: SSPA.Rnw:45-46 ################################################### options(width=60) ################################################### ### code chunk number 3: simulated-data-generation ################################################### library(SSPA) library(genefilter) set.seed(12345) m <- 5000 J <- 10 pi0 <- 0.8 m0 <- as.integer(m*pi0) mu <- rbitri(m - m0, a = log2(1.2), b = log2(4), m = log2(2)) data <- simdat(mu, m=m, pi0=pi0, J=J, noise=0.01) statistics <- rowttests(data, factor(rep(c(0, 1), each=J)))$statistic ################################################### ### code chunk number 4: simulated-data-plotting ################################################### pdD <- pilotData(statistics = statistics, samplesize = sqrt(1/(1/J +1/J)), distribution="norm") pdD plot(pdD) ################################################### ### code chunk number 5: simulated-data-deconvolution ################################################### ssD <- sampleSize(pdD, control=list(from=-6, to=6)) ssD plot(ssD, panel = function(x, y, ...) { panel.xyplot(x, y) panel.curve(dbitri(x), lwd=2, lty=2, n=500) }, ylim=c(0, 0.6)) ################################################### ### code chunk number 6: simulated-data-power ################################################### Jpred <- seq(10, 20, by=2) N <- sqrt(Jpred/2) pwrD <- predictpower(ssD, samplesizes=N, alpha=0.05) matplot(Jpred, pwrD, type="b", pch=16, ylim=c(0, 1), ylab="predicted power", xlab="sample size (per group)") grid() ################################################### ### code chunk number 7: simulated-data-effectsize ################################################### pdC <- pilotData(statistics = statistics, samplesize = sqrt(2*J), distribution="t", df=2*J-2) ssC <- sampleSize(pdC, method="congrad", control=list(from=-6, to=6, resolution=250)) plot(ssC, panel = function(x, y, ...) { panel.xyplot(x, y) panel.curve(2*dbitri(2*x), lwd=2, lty=2, n=500) }, xlim=c(-2,2), ylim=c(0, 1.2)) ################################################### ### code chunk number 8: tikohonov ################################################### ssT <- sampleSize(pdC, method="tikhonov", control=list(resolution=250, scale="pdfstat", lambda = 10^seq(-10, 10, length=50), verbose=FALSE, modelselection="lcurve")) plot(ssT, panel = function(x, y, ...) { panel.xyplot(x, y, type="b") panel.curve(2*dbitri(2*x), lwd=2, lty=2, n=500) }, xlim=c(-2,2), ylim=c(0, 1.2)) ################################################### ### code chunk number 9: nutrigenomics-effectsize ################################################### library(lattice) data(Nutrigenomics) str(Nutrigenomics) pd <- apply(Nutrigenomics, 2, function(x) pilotData(statistics=x[-1], samplesize=sqrt(x[1]), distribution="norm")) ss <- lapply(pd, sampleSize, control=list(pi0Method="Storey", a=0, resolution=2^10, verbose=FALSE)) ##ss <- lapply(pd, sampleSize, ## method = "congrad", ## control=list(verbose=FALSE, resolution=2^10, from=-10, to=10)) compounds <- c("Wy14,643", "fenofibrate", "trilinolenin (C18:3)", "Wy14,643", "fenofibrate") exposure <- c("5 Days", "6 Hours") effectsize <- data.frame(exposure = factor(rep(rep(exposure, c(2, 3)), each=1024)), compound = factor(rep(compounds, each=1024)), lambda = as.vector(sapply(ss, function(x)x@lambda)), theta = as.vector(sapply(ss, function(x)x@theta))) print(xyplot(lambda~theta|exposure, group=compound, data=effectsize, type=c("g", "l"), layout=c(1,2), lwd=2, xlab="effect size", ylab="", auto.key=list(columns=3, lines=TRUE, points=FALSE, cex=0.7))) ################################################### ### code chunk number 10: nutrigenomics-power ################################################### samplesize <- seq(2, 8) averagepower <- data.frame(power = as.vector(sapply(ss, function(x) as.numeric(predictpower(x, samplesize=sqrt(samplesize))))), exposure = factor(rep(rep(exposure, c(2, 3)), each=length(samplesize))), compound = factor(rep(compounds, each=length(samplesize))), samplesize = rep(2*samplesize, 5)) print(xyplot(power~samplesize|exposure, group=compound, data=averagepower, type=c("g", "b"), layout=c(1,2), lwd=2, pch=16, xlab="sample size (per group)", ylab="", auto.key=list(columns=3, lines=TRUE, points=FALSE, cex=0.7))) ################################################### ### code chunk number 11: obtain-teststatistics (eval = FALSE) ################################################### ## ##files contains the full path and file names of each sample ## targets <- data.frame(files=files, ## group=rep(c("DCLK", "WT"), 4), ## description=rep(c("transgenic (Dclk1) mouse hippocampus", ## "wild-type mouse hippocampus"), 4)) ## d <- readDGE(targets) ##reading the data ## ##filter out low read counts ## cpm.d <- cpm(d) ## d <- d[rowSums(cpm.d > 1) >= 4, ] ## ## design <- model.matrix(~group, data=d$samples) ## ##estimate dispersion ## disp <- estimateGLMCommonDisp(d, design) ## disp <- estimateGLMTagwiseDisp(disp, design) ## ##fit model ## glmfit.hoen <- glmFit(d, design, dispersion = disp$tagwise.dispersion) ## ##perform likelihood ratio test ## lrt.hoen <- glmLRT(glmfit.hoen) ## ##extract results ## tbl <- topTags(lrt.hoen, n=nrow(d))[[1]] ## statistics <- tbl$LR ################################################### ### code chunk number 12: deepsage-effectsize ################################################### library(lattice) data(deepSAGE) str(deepSAGE) pd <- pilotData(statistics=deepSAGE, samplesize=8, distribution="chisq", df=1) ss <- sampleSize(pd, method="congrad", control=list(trim=c(0, 0.98), symmetric=FALSE, from=0, to=10)) pwr <- predictpower(ss, samplesize=c(8, 16, 32)) op <- par(mfcol=c(2,1), mar=c(5,4,1,1)) plot(ss@theta, ss@lambda, xlab="effect size", ylab="", type="l") plot(c(8, 16, 32), pwr, xlab="sample size", ylab="power", type="b", ylim=c(0,1)) grid(col=1) par(op) ################################################### ### code chunk number 13: sessioninfo ################################################### toLatex(sessionInfo())