Solved – Negative Binomial dispersion parameter in SPSS

generalized linear modelnegative-binomial-distributionoverdispersionspss

I am working in accident prediction modeling and I'm using SPSS Generalized Linear Model procedure with Negative Binomial distribution and Log Link Function.

Does anyone know the form in wich SPSS returns the dispersion parameter?
Is it the form Var(x)=mu+dispersion*mu^2?
All i know is that the dispersion parameter is equal to 0 for a Poisson distribution, where Var(x)=mu.

Many thanks in advance

Best Answer

Yes you have the form right, specifically the variance for the $j$th observation is:

$$\text{Var}(x_j) = \mu_j(1 + au_j)$$

Stata and SPSS have the same defaults, so I will refer to Stata's documentation (see page 5) for a nice walkthrough. This is the NB2 model from Cameron and Trivedi (which I have not read, but is cited as such in the Long and Freese Regression Models for Categorical Dependent Variables Using Stata).

It can be tricky when reading texts. Sometimes people refer to the dispersion as the entire thing in the parentheses (which changes per the mean estimate of the variable), but other times they only refer to $a$ (which is what the software reports). I've written a blog post translating between a few of the regularly used notations for negative binomial parameterizations.

Related Solutions

Negative Binomial Distribution – Dispersion Parameter of Negbin Distribution

For a general exponential family, we have the variance in the following form:

$$Var(Y_{i})=\phi h(E[Y_{i}])$$

for some function $h(.)$. Using the wikipedia definition of negative binomial we have a pdf of:

$$p(Y_{i}=y|r,p)={r+y-1 \choose y}p^{y}(1-p)^{r}\;\;\;\;\;\;y=0,1,2,\dots$$

And this has expectation $E[Y_{i}]=\frac{pr}{1-p}$ and a variance equal to $Var[Y_{i}]=\frac{pr}{(1-p)^{2}}$. Note that this cannot be written in the usual form for a generalised linear model, but has the form:

$$Var[Y_{i}]=E[Y_{i}]+\frac{1}{r}E[Y_{i}]^{2}$$

And as such it can be seen as taking the function $h(x)=x+\frac{1}{r}x^{2}$. hence dispersion equal to "1". although neg binomial technically not a member of exponential family (as it is a mixture of exponential family, similar to student distribution).

Solved – Negative binomial regression in R allowing for correlation between dispersion & regression coefficients

I haven't found another R package which does this, but I have written code which, based on the maximum likelihood estimates of a model fitted with glm.nb, calculates the full variance covariance matrix using the observed information matrix.

Comparing to values from SAS this appears to match, but if anyone spots an error or finds that it does not match the variance covariance matrix from SAS or Stata, please add a comment to this answer.

glm.nb.cov <- function(mod) {
  #given a model fitted by glm.nb in MASS, this function returns a variance covariance matrix for the
  #regression coefficients and dispersion parameter, without assuming independence between these
  #note that the model must have been fitted with x=TRUE argument so that design matrix is available

  #formulae based on p23-p24 of http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf
  #and http://www.math.mcgill.ca/~dstephens/523/Papers/Lawless-1987-CJS.pdf

  k <- mod$theta
  #p is number of regression coefficients
  p <- dim(vcov(mod))[1]

  #construct observed information matrix
  obsInfo <- array(0, dim=c(p+1, p+1))

  #first calculate top left part for regression coefficients
  for (i in 1:p) {
    for (j in 1:p) {
      obsInfo[i,j] <- sum( (1+mod$y/mod$theta)*mod$fitted.values*mod$x[,i]*mod$x[,j] / (1+mod$fitted.values/mod$theta)^2  )
    }
  }

  #information for dispersion parameter
  obsInfo[(p+1),(p+1)] <- -sum(trigamma(mod$theta+mod$y) - trigamma(mod$theta) -
                                 1/(mod$fitted.values+mod$theta) + (mod$theta+mod$y)/(mod$theta+mod$fitted.values)^2 - 
                                 1/(mod$fitted.values+mod$theta) + 1/mod$theta)

  #covariance between regression coefficients and dispersion
  for (i in 1:p) {
    obsInfo[(p+1),i] <- -sum(((mod$y-mod$fitted.values) * mod$fitted.values / ( (mod$theta+mod$fitted.values)^2 )) * mod$x[,i] )
    obsInfo[i,(p+1)] <- obsInfo[(p+1),i]
  }

  #return variance covariance matrix
  solve(obsInfo)
}

Best Answer

Related Solutions

Negative Binomial Distribution – Dispersion Parameter of Negbin Distribution

Solved – Negative binomial regression in R allowing for correlation between dispersion & regression coefficients

Related Question