SpinevsConvest.R

Code snipits and programs written in R, S or S-PLUS library(limma) library(siggenes)

simNorm2 <- function(m=1000, pi0=0.8, mu0=0, s1=1, s2=3) {

 m1 <- round(m * (1-pi0)) 

1 m0 <- m - m1

 means <- rep(0,m)
 sigma <- c(rep(s2,m1), rep(s1,m0))
 DEflag <- as.logical(sigma - min(sigma))
 X <- rnorm(m, means, sigma)
 pvalues <- 2*(1-pnorm(abs(X),0,s1))
 return(pvalues)

}

pvalues <- simNorm2(m=1000, pi0=0.25)

convest(pvalues, niter=150, doreport=TRUE, file="conv.txt") conv <- read.table("conv.txt", header=TRUE, as.is=TRUE)

1. Interestingly the convergence rate varies if m is small

plot(conv$pi0, type="l")

lambdaseq <- seq(0, 0.95, by=0.05) pi0.est(pvalues, lambdaseq, ncs.value="paper", ncs.weights=1-lambdaseq)$p0

splitPlots <- function(pvalues, m=1000) {

 lambdaseq <- seq(0,(m-1)/m, by=0.01)
 pi0lambda <- tapply(lambdaseq, lambdaseq, function(x){sum(pvalues>x)/(m*(1-x))})
 splinefit <- smooth.spline(lambdaseq, pi0lambda,  w=1-lambdaseq, df=3)$y
 print(paste("Estimate1=", splinefit[100]))
 par(mfrow=c(1,2))
 hist(pvalues, breaks=20, prob=T)
 abline(h=1, lty=1)
 plot(lambdaseq[1:99], pi0lambda[1:99])
 lines(lambdaseq, splinefit, lty=2)

}

splitPlots(pvalues)