[R] Negative Binomial Predictions

[Achim Zeileis]

On Wed, 1 Oct 2008, Donald Catanzaro, PhD wrote:

> Good Day All,
>
> I have a negative binomial model which I have developed using the MASS 
> library.  I now would like to develop some predictions from it. 
> Running the predict.glm (stats library) using type="response" gives me a 
> non-integer value which was rather puzzling.  I would like to confirm that 
> this is actually the mean predicted value of the probability mass function as 
> opposed to the most likely value.  I am afraid that reading the help file for 
> predict.glm either does not state this or I don't understand what it states 
> (which of course could always be the case)

It gives the conditional mean, i.e., the inverse link applied to the 
linear predictor x'b.

> Given that this is a negative binomial model, the mean is often times to the 
> right of the most likely value, so I'd like to ask how one would go about 
> predicting the most likely value.

You can compute the expected probabilities for each count using the 
predprob() method in package "pscl", or by hand using the dnbinom() 
function:

## load data and fit model
library("pscl")
data("bioChemists")
fm <- glm.nb(art ~ ., data = bioChemists)

## compute expected probabilites for y = 0
## via predprob()
pp <- predprob(fm)
pp[,1]
## by hand
dnbinom(0, mu = fitted(fm), size = fm$theta)

## compute the count with the highest probability for each observation
apply(pp, 1, which.max) - 1
## or the median
apply(pp, 1, function(x) max(which(cumsum(x) < 0.5)))

hth,
Z

> Thanks in advance !
>
> -- 
>
> -Don 
> Don Catanzaro, PhD                  Landscape Ecologist
> dgcatanzaro at gmail.com               16144 Sigmond Lane
> 479-751-3616                        Lowell, AR 72745
>
> ______________________________________________
> R-help at 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.
>
>