The HPGENSELECT Procedure

Log-Likelihood Functions

The HPGENSELECT procedure forms the log-likelihood functions of the various models as

\[  L(\bmu ;\mb{y}) = \sum _{i=1}^{n} f_ i \,  l(\mu _ i;y_ i,w_ i)  \]

where $l(\mu _ i;y_ i,w_ i)$ is the log-likelihood contribution of the ith observation that has weight $w_ i$, and $f_ i$ is the value of the frequency variable. For the determination of $w_ i$ and $f_ i$, see the WEIGHT and FREQ statements. The individual log-likelihood contributions for the various distributions are as follows.

In the following, the mean parameter $\mu _ i$ for each observation i is related to the regression parameters $\bbeta _ i$ through the linear predictor $\eta _ i=\mb{x}_ i^\prime \bbeta $ by

\[  \mu _ i = g^{-1}(\eta _ i) \]

where g is the link function.

There are two link functions and linear predictors that are associated with zero-inflated Poisson and zero-inflated negative binomial distributions: one for the zero-inflation probability $\omega $, and another for the parameter $\lambda $, which is the Poisson or negative binomial mean if there is no zero-inflation. Each of these parameters is related to regression parameters through an individual link function,

\begin{align*}  \eta _ i & = \mb{x}_ i^\prime \bbeta \\ \kappa _ i & = \mb{z}_ i^\prime \bgamma \\ \lambda _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ \omega _ i(\bgamma ) & = h^{-1}(\kappa _ i) \\ \end{align*}

where h is one of the following link functions that are associated with binary data: complementary log-log, log-log, logit, or probit. These link functions are also shown in Table 7.9.

Binary Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith binary observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i) & = y_ i \log \{ \mu _ i\}  + (1-y_ i)\log \{ 1-\mu _ i\}  \end{align*}

Here, $\mu _ i$ is the probability of an event, and the variable $y_ i$ takes on the value 1 for an event and the value 0 for a non-event. The inverse link function $g^{-1}(\cdot )$ maps from the scale of the linear predictor $\eta _ i$ to the scale of the mean. For example, for the logit link (the default),

\[  \mu _ i(\bbeta ) = \frac{\exp \{ \eta _ i\} }{1+\exp \{ \eta _ i\} }  \]

You can control which binary outcome in your data is modeled as the event by specifying the response-options in the MODEL statement, and you can choose the link function by specifying the LINK= option in the MODEL statement.

If a WEIGHT statement is specified and $w_ i$ denotes the weight for the current observation, the log-likelihood function is computed as

\[  l(\mu _ i(\bbeta );y_ i,w_ i) = w_ i l(\mu _ i(\bbeta );y_ i)  \]

Binomial Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith binomial observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i,w_ i) & = w_ i \left( y_ i \log \{ \mu _ i\}  + (n_ i - y_ i) \log \{ 1-\mu _ i\}  \right) \\ & + w_ i \left( \log \{ \Gamma (n_ i+1)\}  - \log \{ \Gamma (y_ i+1)\}  - \log \{ \Gamma (n_ i-y_ i+1)\} \right) \end{align*}

where $y_ i$ and $n_ i$ are the values of the events and trials of the ith observation, respectively. $\mu _ i$ measures the probability of events (successes) in the underlying Bernoulli distribution whose aggregate follows the binomial distribution.

Gamma Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i,w_ i) & = \frac{w_ i}{\phi } \log \left( \frac{w_ i y_ i}{\phi \mu _ i} \right) - \frac{w_ i y_ i}{\phi \mu _ i} - \log (y_ i) - \log \left( \Gamma \left( \frac{w_ i}{\phi } \right) \right) \end{align*}

For the gamma distribution, $\nu =\frac{1}{\phi }$ is the estimated dispersion parameter that is displayed in the output.

Inverse Gaussian Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i,w_ i) & = -\frac{1}{2} \left[ \frac{w_ i(y_ i-\mu _ i)^2}{y_ i \mu ^2 \phi } + \log \left( \frac{\phi y_ i^3}{w_ i} \right) + \log (2 \pi ) \right] \end{align*}

where $\phi $ is the dispersion parameter.

Multinomial Distribution

The multinomial distribution that is modeled by the HPGENSELECT procedure is a generalization of the binary distribution; it is the distribution of a single draw from a discrete distribution with J possible values. The log-likelihood function for the ith observation is

\[  l(\bmu _ i;\mb{y}_ i,w_ i) = w_ i \sum _{j=1}^{J} y_{ij}\log \{ \mu _{ij}\}   \]

In this expression, J denotes the number of response categories (the number of possible outcomes) and $\mu _{ij}$ is the probability that the ith observation takes on the response value that is associated with category j. The category probabilities must satisfy

\[  \sum _{j=1}^{J} \mu _ j = 1  \]

and the constraint is satisfied by modeling $J-1$ categories. In models that have ordered response categories, the probabilities are expressed in cumulative form, so that the last category is redundant. In generalized logit models (multinomial models that have unordered categories), one category is chosen as the reference category and the linear predictor in the reference category is set to 0.

Negative Binomial Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i,w_ i) & = y_ i\log \left(\frac{k\mu }{w_ i} \right) - (y_ i+w_ i/k)\log \left(1+\frac{k\mu }{w_ i} \right) + \log \left(\frac{\Gamma (y_ i+w_ i/k)}{\Gamma (y_ i+1)\Gamma (w_ i/k)}\right) \end{align*}

where k is the negative binomial dispersion parameter that is displayed in the output.

Normal Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i,w_ i) & = -\frac{1}{2} \left[ \frac{w_ i(y_ i-\mu _ i)^2}{\phi } + \log \left( \frac{\phi }{w_ i} \right) + \log (2 \pi ) \right] \end{align*}

where $\phi $ is the dispersion parameter.

Poisson Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\mu _ i(\bbeta );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \mu _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ l(\mu _ i(\bbeta );y_ i,w_ i) & = w_ i[y_ i \log (\mu _ i) - \mu _ i - \log (y_ i!) ] \end{align*}

Tweedie Distribution

The Tweedie distribution does not in general have a closed form log-likelihood function in terms of the mean, dispersion, and power parameters. The form of the log likelihood is

\[  L(\bmu ;\mb{y}) = \sum _{i=1}^{n} f_ i \,  l(\mu _ i;y_ i,w_ i)  \]

where

\[  l(\mu _ i,y_ i,w_ i) = \log (f(y_ i;\mu _ i,p,\frac{\phi }{w_ i}))  \]

and $f(y,\mu ,p,\phi )$ is the Tweedie probability distribution, which is described in the section Tweedie Distribution. Evaluation of the Tweedie log-likelihood for model fitting is performed numerically as described in Dunn and Smyth (2005, 2008).

Quasi-likelihood

The extended quasi-likelihood (EQL) is constructed according to the definition of McCullagh and Nelder (1989, Chapter 9) as

\[  Q_ p(\mb{y}, \bmu , \phi , p) = \sum _ i q(y_ i, \mu _ i, \phi , p)  \]

where the contribution from an observation is

\[  q(y_ i, \mu _ i, \phi , p) = -0.5 \log (2 \pi \frac{\phi }{w_ i} y_ i^ p) - w_ i \left( \frac{y_ i^{2-p}-(2-p)y_ i \mu _ i^{1-p} + (1-p) \mu _ i^{2-p}}{(1-p)(1-p)} \right) / \phi  \]

where $1<p<2$. This EQL is used in computing initial values for the iterative maximization of the Tweedie log likelihood, as specified using the OPTMETHOD= Tweedie option in Table 7.5. If you specify the OPTMETHOD=EQL Tweedie-optimization-option in Table 7.5, then the parameter estimates are computed by using the EQL instead of the log likelihood.

Zero-Inflated Negative Binomial Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\lambda _ i(\bbeta ),\omega _ i(\bgamma );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \kappa _ i & = \mb{z}_ i^\prime \bgamma \\ \lambda _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ \omega _ i(\bgamma ) & = h^{-1}(\kappa _ i)\\ l(\mu _ i(\bbeta ),\omega _ i(\bgamma );y_ i,w_ i) & = \left\{  \begin{array}{ll} \log [\omega _ i + (1-\omega _ i)(1+\frac{k}{w_ i}\lambda )^{-\frac{1}{k}}] &  y_ i=0 \\ \log (1-\omega _ i)+ y_ i\log \left(\frac{k\lambda }{w_ i} \right) & \\ -(y_ i+\frac{w_ i}{k})\log \left(1+\frac{k\lambda }{w_ i} \right) & \\ +\log \left(\frac{\Gamma (y_ i+\frac{w_ i}{k})}{\Gamma (y_ i+1)\Gamma (\frac{w_ i}{k})}\right) &  y_ i>0 \\ \end{array} \right. \end{align*}

where k is the zero-inflated negative binomial dispersion parameter that is displayed in the output.

Zero-Inflated Poisson Distribution

The HPGENSELECT procedure computes the log-likelihood function $l(\lambda _ i(\bbeta ),\omega _ i(\bgamma );y_ i)$ for the ith observation as

\begin{align*}  \eta _ i & = \mb{x}_ i’\bbeta \\ \kappa _ i & = \mb{z}_ i^\prime \bgamma \\ \lambda _ i(\bbeta ) & = g^{-1}(\eta _ i) \\ \omega _ i(\bgamma ) & = h^{-1}(\kappa _ i)\\ l(\mu _ i(\bbeta ),\omega _ i(\bgamma );y_ i,w_ i) & = \left\{  \begin{array}{ll} w_ i\log [\omega _ i + (1-\omega _ i)\exp (-\lambda _ i)] &  y_ i=0 \\ w_ i[\log (1-\omega _ i)+y_ i \log (\lambda _ i) - \lambda _ i-\log (y_ i!)] &  y_ i>0 \\ \end{array} \right. \end{align*}