ExamineConvest.R

library(limma)

# ============= 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)
 #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=50)

pi0 <- c() for(i in seq(1,100, by=1)) {

 pi0[i] <-   convest(pvalues, niter=i)

}

1. Interestingly the convergence rate varies if m is small

plot(pi0[1:100], type="l")

plot(diff(pi0), type="l")

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

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

1. Weighted mean

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

convest(pvalues, niter=50) pi0[100]