#
# zip-sim.q
#
# ----------------------------------------------------------------------
#
# PURPOSE:  simulate overdispersed count data (Zero-Inflated Poisson)
#
# AUTHOR:  P. Heagerty
#
# DATE:  23Jan2001
#        25Jan2007
#
# ----------------------------------------------------------------------
Nsim <- 1000
p <- 3
save.beta <- matrix( NA, Nsim, 3 )
save.se.poisson <- matrix( NA, Nsim, 3 )
save.se.scale <- matrix( NA, Nsim, 3 )
save.se.sandwich <- matrix( NA, Nsim, 3 )
save.se.empirical <- matrix( NA, Nsim, 3 )
save.alpha <- rep( NA, Nsim )
save.phi <- rep( NA, Nsim )
#
save.beta.nb <- matrix( NA, Nsim, 3 )
save.se.nb <- matrix( NA, Nsim, 3 )
save.se.emp.nb <- matrix( NA, Nsim, 3 )
save.alpha.nb <- rep( NA, Nsim )
# ----------------------------------------------------------------------
nobs <- 150
gamma <- c( -1, 2, 1 )
gamma <- c( 0, 2, 1 )
delta <- c( log(0.6), log(0.80/0.60), 0 )
# ----------------------------------------------------------------------
#
library(MASS)
#
for( j in 1:Nsim ){
#x1 <- runif(nobs)
x1 <- c(1:nobs)/nobs
x2 <- c( rep(0,nobs/2), rep(1,nobs/2) )
X <- cbind( 1, x1, x2 )
#
p <- exp( as.vector( X%*%delta ) )
lambda <- exp( as.vector( X%*%gamma ) )
#
#####
#
nu <- rbinom( n=rep(1,nobs), size=rep(1,nobs), prob=p )
z <- rpois( n=rep(1,nobs), lambda = lambda )
#
y <- nu * z
#
simdata <- data.frame(
		y = y,
		x1 = x1,
                x2 = x2  )
#
##### fit the model
#
options( contrasts=c("contr.treatment","contr.helmert") )
#
glmfit <- glm( y ~ x1 + x2, family=poisson, data=simdata )
#
nbfit <- glm.nb( y~ x1 + x2, data=simdata )
#
##### moment estimates
#
#
mu <- fitted( glmfit )
#
resids <- residuals( glmfit, type="pearson" )
#
phi <- sum( resids^2 )/( nobs - ncol(X) )
##### (real value):  phi <- mean( 1 + ((1-p)/p)*mu )
#
alpha <- mean( (resids^2 - 1)/mu )
##### (real value):  alpha <- mean( (1-p)/p )
#
#####
##### estimate model information
#####
#
X <- cbind( 1, x1, x2 )
#
XtWX <- t(X)%*%( (mu)*X )
I1inv <- solve( XtWX )
XtWVWX <- t(X)%*%( (mu+alpha*mu^2)*X )
#
est.fnx <- X*(y-mu)
UUt <- t(est.fnx)%*%est.fnx
#
se.model <- sqrt( diag( I1inv ) )
se.scale <- sqrt( phi ) * se.model
se.sandwich <- sqrt( diag( I1inv%*%XtWVWX%*%I1inv ) )
se.empirical <- sqrt( diag( I1inv%*%UUt%*%I1inv ) )
#
bhat <- glmfit$coef
#
#####
##### sandwich for NB model too
#####
mu <- fitted( nbfit )
Ainv <- summary(nbfit)$cov.scaled
theta <- nbfit$theta
est.fnx <- X*(1/(1+mu/theta))*(y-mu)
UUt <- t(est.fnx)%*%est.fnx
se.emp.nb <- sqrt( diag( Ainv%*%UUt%*%Ainv ) )
#
save.beta[ j, ] <- bhat
save.se.poisson[ j , ] <- se.model
save.se.scale[ j , ] <- se.scale
save.se.sandwich[ j , ] <- se.sandwich
save.se.empirical[ j , ] <- se.empirical
save.alpha[ j ] <- alpha
save.phi[ j ] <- phi
#
save.beta.nb[ j , ] <- nbfit$coef
save.se.nb[ j , ] <- sqrt( diag( summary(nbfit)$cov.scaled ) )
save.se.emp.nb[j ,] <- se.emp.nb
save.alpha.nb[ j ] <- 1/nbfit$theta 
}
#
#
# end-of-file...