#
# sample-size.q
#
#----------------------------------------------------------------------
#
mu0 <- 5.2
mu1 <- mu0 * 0.75
#
alpha <- 1/0.4391
alpha <- 2.5
#
scale <- 13.5
#
n <- 800 + c(1:20)*50
#
Zalpha <- qnorm(0.025)
Zbeta <- qnorm(0.80)
#
#####
##### Method #1 -- get the variance of beta1 using delta method
#####
#
##### negative binomial (NB2) variance
#
vvv.1 <- (2/n)*( mu0 + alpha*mu0^2)/(mu0^2) +
       (2/n)*( mu1 + alpha*mu1^2)/(mu1^2) 
#
##### scale overdispersion
#
vvv.2 <- (2/n)*( 1/mu0 + 1/mu1 )*scale   
#
ThePower.1 <- pnorm( Zalpha - log(0.75)/sqrt(vvv.1) )
ThePower.2 <- pnorm( Zalpha - log(0.75)/sqrt(vvv.2) )
#
print( cbind( n, ThePower.1, ThePower.2 ) )
method1.results <- cbind( n, ThePower.1, ThePower.2 )
#
#
#####
##### Method #2 -- use EE methods to compute s.e.
#####
#
##### consider precision for 1 TX and 1 CONTROL subject
#
X <- cbind( 1, c(0,1) )
#
muVec <- c( mu0, mu1 )
#
Amat <- t(X)%*%diag( muVec )%*%X
#
##### negative binomial (NB2) variance
#
Bmat.1 <- t(X)%*%diag( muVec + alpha*muVec^2 )%*%X
#
##### scale overdispersion
#
Bmat.2 <- t(X)%*%diag( scale * muVec )%*%X
#
sandwich.1 <- solve(Amat)%*%Bmat.1%*%solve(Amat)
sandwich.2 <- solve(Amat)%*%Bmat.2%*%solve(Amat)
#
vvv.1 <- (2/n)*sandwich.1[2,2]
vvv.2 <- (2/n)*sandwich.2[2,2]
#
ThePower.1 <- pnorm( Zalpha - log(0.75)/sqrt(vvv.1) )
ThePower.2 <- pnorm( Zalpha - log(0.75)/sqrt(vvv.2) )
#
print( cbind( n, ThePower.1, ThePower.2 ) )
method2.results <- cbind( n, ThePower.1, ThePower.2 )
#
##### compare!
#
print( round( cbind( method1.results, method2.results ), 3 ) )
#
# end-of-file...