ZINB power analysis

Power analysis for a zero-inflated negative binomial model.

library("knitr")
library("emdbook") ## for rzinbinom
library("pscl")    ## for zeroinfl
library("splines") ## interpSpline, backSpline
## simulate data
##' @param beta0 baseline log-counts
##' @param beta1 difference in log-counts
##' @param k NB overdispersion
##' @param pz zero-inflation
##' @param n1 sample size of group 1
##' @param n2 sample size of group 2
simfun <- function(beta1=0.2,beta0=1,k=2.5,pz=0.3,n1=50,n2=n1) {
    mu <- exp(rep(c(beta0,beta0+beta1),c(n1,n2)))
    data.frame(f=factor(rep(1:2,c(n1,n2))),
               y=rzinbinom(n1+n2,mu=mu,zprob=pz,size=k))
}
## calc. p-value of difference
sumfun <- function(d) {
    m1 <- zeroinfl(y~f | 1, data=d)
    m2 <- zeroinfl(y~1 | 1, data=d)
    ## LRT: too bad zeroinfl doesn't have an anova() method ...
    unname(pchisq(2*(logLik(m1)-logLik(m2)),df=1,lower.tail=FALSE))
}
## calculate power
powfun <- function(beta1=0.2,n1=50,...,nsim=200,alpha=0.05) {
    mean(replicate(nsim,sumfun(simfun(beta1=beta1,n1=n1,...)))<alpha)
}

Basic tests:

Plot sim results

plot(log(1+y)~f,simfun(beta1=1))

p-value for a single sim

sumfun(simfun(beta1=0.5))
## 'log Lik.' 0.0003119597 (df=3)

Power

powfun(beta1=0.7,n1=200)
## [1] 1

Power as a function of effect size (log-ratio difference):

svec <- seq(0.1,2,by=0.1)
set.seed(101)
powvec1 <- sapply(svec,powfun)
maxpow <- which(powvec1==1)[1]
ss <- interpSpline(svec[1:maxpow],powvec1[1:maxpow])
bs <-backSpline(ss)
par(las=1,bty="l")
plot(svec,powvec1,xlab="effect size (log-ratio)",ylab="power")
svec2 <- seq(0.1,svec[maxpow],length=101)
lines(predict(ss,svec2))
abline(h=0.8,lty=2)
crit.effsize <- predict(bs,0.8)$y
abline(v=crit.effsize,lty=2)

Effect size for 80% power = 0.719

nvec <- round(10^seq(1,3.5,length=15))
set.seed(102)
powvec2 <- sapply(nvec,function(x) powfun(n1=x))
##diff(powvec2)
pow.OK <- -(1:3)
ss2 <- interpSpline(nvec[pow.OK],powvec2[pow.OK])
bs2 <-backSpline(ss2)
par(las=1,bty="l")
plot(nvec,powvec2,xlab="sample size (n1)",ylab="power",log="x")
nvec2 <- seq(10,3000,length=101)
lines(predict(ss2,nvec2))
abline(h=0.8,lty=2)
crit.n <- predict(bs2,0.8)$y
abline(v=crit.n,lty=2)
abline(h=0.8,lty=2)

For the default effect size (0.2), the required (half!) sample size for 80% power is 839