On Thu, 26 Apr 2012, Christopher Desjardins wrote:
Hi,
I am trying to replicate Lambert (1992)'s simulation with zero-inflated
Poisson models. The citation is here:
@article{lambert1992zero,
Author = {Lambert, D.},
Journal = {Technometrics},
Pages = {1--14},
Publisher = {JSTOR},
Title = {Zero-inflated {P}oisson regression, with an application to defects
in manufacturing},
Year = {1992}}
Specifically I am trying to recreate Table 2. But my estimates for Gamma
are not working out. Any ideas why?
Your implementation of the inverse logit link is wrong, it should be
exp(x)/(1 + exp(x)) and not exp(x)/exp(1 + x). In R I never code this by
hand but either use qlogis()/plogis() or make.link("logit").
Your setting resulting in an almost constant probability of zero inflation
and hence gamma was completely off.
Below is my cleaned up code which nicely replicates the results for n =
100. The other two would require additional work because one would need to
catch non-convergence etc.
hth,
Z
## data-generating process
dgp <- function(nobs = 100) {
gamma <- c(-1.5, 2)
beta <- c(1.5, -2)
x <- seq(0, 1, length.out = nobs)
lambda <- exp(beta[1] + beta[2] * x)
p <- plogis(gamma[1] + gamma[2] * x)
y <- ifelse(runif(nobs) <= p, 0, rpois(nobs, lambda = lambda))
data.frame(y = y, x = x)
}
## simulation of coefficients and standard errors
sim <- function(nrep = 2000, ...) {
onesim <- function(i) {
d <- dgp(...)
m <- zeroinfl(y ~ x, data = d)
c(coef(m), sqrt(diag(vcov(m))))
}
t(sapply(1:nrep, onesim))
}
## run simulation #3
library("pscl")
set.seed(1)
cfse <- sim(nobs = 100)
## mean coefficient estimates
apply(cfse[, 1:4], 2, mean)
## median coefficient estimates
apply(cfse[, 1:4], 2, median)
## sd of coefficient estimates
apply(cfse[, 1:4], 2, sd)
## mean of the standard deviation estimates from observed information
apply(cfse[, 5:8], 2, mean)
Please cc me on a reply!
Thanks,
Chris
## ZIP simulations based on Lambert (1992)'s conditions
# Empty workspace
rm(list=ls())
# Load zero-inflation package
library(pscl)
# Simulation conditions #3
NumSimulations=2000
X=seq(from=0,to=1,length.out=N)
Model.Matrix=cbind(rep(1,length(X)),X)
Gamma=c(-1.5,2)
Beta=c(1.5,-2)
P=exp(Model.Matrix%*%Gamma)/exp(1+Model.Matrix%*%Gamma)
Lambda=exp(Model.Matrix%*%Beta)
CoefSimulations=matrix(nrow=NumSimulations,ncol=2*dim(Model.Matrix)[2])
for(i in 1 : NumSimulations){
Lambda.Draw=rpois(N,Lambda)
U=runif(N)
Y=ifelse(U<=P,0,Lambda.Draw)
CoefSimulations[i,]=coef(zeroinfl(Y~X|X))
}
# What were the estimates?
colMeans(CoefSimulations) # My gamma estimates aren't even close ...
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
