ExamineConvest.R

Code snipits and programs written in R, S or S-PLUS Sys.putenv("DISPLAY"=":0")

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)) 

1 m0 <- m - m1

 means <- rep(0,m)
 #c(rep(mu1,m1), rep(mu0,m0))
 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. ---------------------------------------------------------------------------- #
2. Multtest does not allow a max of 1 (NA) going into smooth.spline

1. convest arguement doreport=TRUE provides iteration information

pvalues <- simNorm2(m=500)

convest(pvalues, niter=50, 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

1. 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. Checking versions

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