zero-truncated Poisson

Jeff,

The zero truncated Poisson distribution can be described as a
probability function conditional on Y>0:

Pr(Y=y|Y>0) = Pr(Y=y) / (1-Pr(Y=0)), y=1,2,3,...

This leads to this function 'nll' for the negative log-likelihood:

nll <- function(theta, y, X) {
    lambda <- exp(drop(X %*% theta))
    -sum(dpois(y, lambda, log=TRUE) - log(1-exp(-lambda)))
}

where theta is the vector of parameters, y is the observation vector,
X is the design matrix.

A simple simulation is useful to check that the theory works. I used
'optim' with the default Nelder-Mead simplex algorithm, initial values
are based in non-truncated Poisson estimates to speed up convergence:

theta <- c(1, 2)
X <- model.matrix(~rnorm(100))
y <- rpois(100,exp(X %*% theta))
X <- X[y>0,]
y <- y[y>0]
inits <- coef(glm(y ~ X-1, family=poisson))
res <- optim(par=inits, fn=nll, y=y, X=X)
cbind(theta=theta, theta.hat=res$par)

Here I used the data you provided, note that 'hessian=TRUE' for some
subsequent SE calculations:

inits <- coef(glm(taken ~ mass + year, blja, family=poisson))
y <- blja$taken
X <- model.matrix(~ mass + year, blja)
res <- optim(par=inits, fn=nll, y=y, X=X, hessian=TRUE)

Vector of parameters that minimize 'nll' (maximize the log-likelihood)
given the data:

cf <- res$par

Standard error (note that due to minimization, I used the Hessian
instead of -H):

se <- sqrt(diag(solve(res$hessian)))

z and p values and the final output table:

zstat <- cf/se
pval <- 2 * pnorm(-abs(zstat))
out <- cbind(cf, se, zstat, pval)
colnames(out) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")

printCoefmat(out, digits = 4, signif.legend = FALSE)

For your data set, this gives:

            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.76392    0.10659  -7.167 7.66e-13 ***
mass         0.24008    0.01262  19.028  < 2e-16 ***
yeartwo     -0.29006    0.11974  -2.423   0.0154 *

Cheers,

Peter

P?ter S?lymos
Alberta Biodiversity Monitoring Institute
and Boreal Avian Modelling project
Department of Biological Sciences
CW 405, Biological Sciences Bldg
University of Alberta
Edmonton, Alberta, T6G 2E9, Canada
Phone: 780.492.8534
Fax: 780.492.7635
email <- paste("solymos", "ualberta.ca", sep = "@")
http://www.abmi.ca
http://www.borealbirds.ca
http://sites.google.com/site/psolymos

On Fri, Apr 8, 2011 at 6:40 PM, Stratford, Jeffrey

zero-truncated Poisson

Thread (6 messages)