ExamineConvest.R

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

1. Load local copy of convest

source("convest.R")

1. ------------- Generating a mixture of normally distributed data ------------- #

simNorm <- function(m=1000, pi0=0.8, mu0=0, mu1=3, sigma=1) {

 m1 <- round(m * (1-pi0)) 
 m0 <- m - m1
 means <- c(rep(mu1,m1), rep(mu0,m0))
 DEflag <- as.logical(means)
 X <- rnorm(m, means, sigma)
 pvalues <- 2*(1-pnorm(abs(X),0,sigma))
 return(pvalues)

}

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

 m1 <- round(m * (1-pi0)) 
 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)

}

1. ---------------------------------------------------------------------------- #

1. convest arguement doreport=TRUE provides iteration information

pvalues <- simNorm(m=500, pi0=0.25)

1. 1) Illustrating that modified convest function is identical to limma version

limma::convest(pvalues, niter=50, doreport=TRUE) convest(pvalues, niter=50, doreport=TRUE)

1. 2) Writing report information to file

convest(pvalues, niter=100, 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") plot(diff(conv$pi0), type="l")

1. Weighted average of convest iterations after removing burnin period
2. Raw

mean(conv$pi0[-(1:10)])

1. Weighted mean

weighted.mean(conv$pi0[-(1:10)], w=seq(0,1, length=length(conv$pi0[-(1:10)])))

1. 3) Tolerance adjustment

convest(pvalues, niter=300, doreport=TRUE, file="conv.txt", tol=1e-6) conv <- read.table("conv.txt", header=TRUE, as.is=TRUE)

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