6
$\begingroup$

I would like to understand how the weight argument of glm.nb is affecting the likelihood function. I understand that glm.nb find the MLE in an alternating iteration process where for a given theta the GLM is fitted using glm() with family=negative.binomial(theta=theta.this.iteration) and for fixed means the theta parameter is estimated using score and information iterations.

Regarding the weight input, the help(glm.nb) only says that it is passed to glm() (See below).

formula, data, weights, subset, na.action, start, etastart, mustart, control, method, model, x, y, contrasts, ...
arguments for the glm() function. Note that these exclude family and offset (but offset() can be used).

So I am assuming that the weights in glm.nb is directly passed to glm() with family=negative.binomial(theta=theta.this.iteration).

Regarding the weights, help(glm) says that the dispersion parameter is rescaled as $\phi=\phi^*/w_i$ (see below):

Non-NULL weights can be used to indicate that different observations have different dispersions (with the values in weights being inversely proportional to the dispersions); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations. For a binomial GLM prior weights are used to give the number of trials when the response is the proportion of successes: they would rarely be used for a Poisson GLM.

In my understanding, the negative binomial distribution with given theta (denoted $\theta$) belongs to exponential family with a dispersion parameter = 1 (is my understanding correct?). To see this,

Suppose a random variable $y_i$ has a negative binomial distribution with mean $\mu_i$ and variance $\mu_i + \frac{1}{\theta}\mu_i^2$ then its density is:

$$ f(y_i;\mu_i,\theta) = \frac{\Gamma(\theta + y_i)}{ \Gamma(\theta) y_i ! } \left(\frac{\mu_i}{\mu_i+\theta}\right)^{y_i} \left(\frac{\theta}{\mu_i+\theta}\right)^{\theta}. $$

If $\theta$ is know, then this negative binomial distribution is from the exponential family as:

\begin{align*} f(y_i;\mu_i,\theta) &= \exp\left( y_i\log \left(\frac{\mu_i}{\mu_i+\theta}\right) + \theta \log \left(\frac{\theta}{\mu_i+\theta}\right) + \log\left[ \frac{\Gamma(\theta + y_i)}{ \Gamma(\theta) y_i ! } \right] \right)\\ &= \exp\left( \frac{m_i}{\phi}\left(y_i\kappa_i-b(\kappa_i)\right)+c(y_i;\phi) \right) \end{align*} where: \begin{align*} \kappa_i&=\log\left(\frac{\mu_i}{\mu_i+\theta} \right)\\ \phi&=1\\ m_i&=1\\ b(\kappa_i) &= - \theta\log(1-e^{\kappa_i})\\ c(y_i;\phi) &=\log\left[ \frac{\Gamma(\theta + y_i)}{ \Gamma(\theta) y_i ! } \right] \end{align*}

so the dispersion parameter $\phi=1$!! Therefore the weighted negative binomial density can be obtained by simply letting $\phi=\phi^*/w_i$:

\begin{align} f(y_i;\mu_i,\theta,w_i) &= \exp \left(w_i \frac{m_i}{\phi^*} (y_i\kappa_i - b(\kappa_i)) + c(y_i;\phi^*/w_i)\right) \\ &= \exp\left( w_iy_i\log \left(\frac{\mu_i}{\mu_i+\theta}\right) + w_i\theta \log \left(\frac{\theta }{\mu_i+\theta }\right) + \log\left[ \frac{\Gamma(\theta + y_i)}{ \Gamma(\theta) y_i ! } \right] \right)\\ \end{align} The last term does NOT depend on $w_i$ because $c(y_i;\phi)$ function does not depend on $\phi$ which leads to $c(y_i;\phi^*/w_i)=c(y_i;\phi)$.

So I originally believed that glm.nb is iteratively maximizing the above weighted likelihood $f(y_i;\mu_i,\theta,w_i)$. In the step where $\theta$ is maximized for given $\mu_i$, glm.nb use internal function named theta.ml. By closely looking at this function, I found that it uses a score function defined as:

  score <- function(n, th, mu, y, w) sum(w * (digamma(th + 
    y) - digamma(th) + log(th) + 1 - log(th + mu) - (y + 
    th)/(mu + th)))

This looks to be the score function of the likelihood:

\begin{align} f^R(y_i;\mu_i,\theta,w_i)= \exp\left( w_iy_i\log \left(\frac{\mu_i}{\mu_i+\theta}\right) + w_i\theta \log \left(\frac{\theta }{\mu_i+\theta }\right) + w_i \log\left[ \frac{\Gamma(\theta + y_i)}{ \Gamma(\theta) y_i ! } \right] \right),\\ \end{align} which is different from $f(y_i;\mu_i,\theta,w_i)$.

Is glm.nb maximizing $f(y_i;\mu_i,\theta,w_i)$ or $^Rf(y_i;\mu_i,\theta,w_i)$?

$\endgroup$
  • $\begingroup$ Have you supplied trace=TRUE? One result from binomial (and poisson and other) likelihoods is that glm(y ~ x, family=binomial, weights=w) and glm(rep.int(y, w) ~ rep.int(x, w), family=binomial) produce identical results up to the deviance when w is frequency weights. If the same held for glm.nb, I'd be inclined to say it's a likelihood of the $f^R$ form. $\endgroup$ – AdamO Mar 21 '18 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.