GLM Fundamentals

This page covers the statistical theory behind the GLM tab. See that page for usage instructions.

Model Formulation

GLM (the generalized linear model) extends the normal linear model to the exponential family of distributions defined in the next section, in order to handle response variables that a normal distribution does not fit—binary outcomes, counts, and positive continuous values—within a regression framework. It was introduced by Nelder & Wedderburn (1972). A GLM is defined by three components:

Distribution family: The response variable $Y$ follows a distribution in the exponential family
Linear predictor: $\eta = X\beta$ (a linear combination of predictors)
Link function: A monotonic function $g$ such that $\eta = g(\mu)$ , connecting the linear predictor to the mean $\mu = E[Y]$

Exponential Family

The distribution of the response changes with the kind of data: Binomial for binary outcomes, Poisson for counts. GLM avoids rebuilding the estimation theory for each distribution by restricting the choice to the exponential family. Within this family, the differences between distributions reduce to a single function expressing the mean-variance relationship, and replacing that function lets the same estimation algorithm (IRLS) cover every family.

A family of distributions is called an exponential family if its density (or mass) function can be written as:

f(y \mid \theta, \phi) = \exp\!\left\{\frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi)\right\}

Here, $\theta$ is called the natural (canonical) parameter, $\phi$ the dispersion parameter, and $b(\theta)$ the log-partition function; $a(\phi)$ is a known function determined by the family (see the table below). From this definition, the mean and variance follow as derivatives of the log-partition function $b(\theta)$ :

$E[Y] = b'(\theta) = \mu$
$\operatorname{Var}(Y) = b''(\theta) \cdot a(\phi)$

Rewriting $b''(\theta)$ as a function of $\mu$ gives the variance function $V(\mu)$ , which expresses the mean-variance relationship $\operatorname{Var}(Y) = V(\mu) \cdot a(\phi)$ through a single function. The per-family forms of $V(\mu)$ are tabulated in Variance Functions and Overdispersion.

Exponential family parameters for each distribution family:

Family	$\theta$ (natural parameter)	$a(\phi)$	$b(\theta)$	$c(y, \phi)$
Gaussian	$\mu$	$\phi$	$\theta^2/2$	$-\dfrac{y^2}{2\phi} - \dfrac{\log(2\pi\phi)}{2}$
Binomial	$\log\!\bigl(\mu/(1-\mu)\bigr)$	$1/n_i$	$\log(1+e^\theta)$	$\log\binom{n_i}{k_i}$
Poisson	$\log\mu$	$1$	$e^\theta$	$-\log(y!)$
Gamma	$-1/\mu$	$\phi$	$-\log(-\theta)$	$(1/\phi - 1)\log y + (1/\phi)\log(1/\phi) - \log\Gamma(1/\phi)$
Negative Binomial	$\log\!\bigl(\mu/(\mu+r)\bigr)$	$1$	$-r\log(1-e^\theta)$	$\log\Gamma(y+r) - \log\Gamma(r) - \log(y!)$

In the Binomial row, $y_i$ is the proportion of successes ( $y_i = k_i/n_i$ , $0 \le y_i \le 1$ ), $k_i$ is the number of successes at observation $i$ , $n_i$ is the number of trials at observation $i$ (not the sample size $n$ but the per-observation trial count), and $\mu$ is the success probability. When $n_i=1$ , it reduces to the Bernoulli distribution
The Negative Binomial $r$ is displayed as $\theta$ in the MIDAS UI; on this page we use $r$ to avoid confusion with the exponential family natural parameter. The Negative Binomial belongs to the exponential family only when $r$ is known. In MIDAS's automatic estimation mode, $r$ is estimated by maximizing the profile likelihood $L_p(r) = \max_\beta L(\beta, r)$ (with $\beta$ profiled out) in an outer loop (see GLM usage). The standard errors for $\hat\beta$ reported in this mode are computed from the information matrix with $r = \hat r$ treated as known, so uncertainty in $r$ is not reflected

Choosing the Link Function

The link function is a monotonic function $\eta = g(\mu)$ connecting the linear predictor $\eta$ to the expected value $\mu$ of the response. A link function satisfying $g(\mu) = \theta$ (the natural parameter) is called the canonical link. While the distribution family is mostly determined by the type of the response, the link function leaves room for choice. For how the family is chosen, see Choosing a distribution family.

The first criterion is the effect scale: whether the coefficients are to be read as differences in the response or as ratios. With the Identity link, coefficients directly represent differences in the response. With the Log link, addition on the linear predictor becomes multiplication on the mean, so $\exp(\beta)$ reads as a ratio—a rate ratio (IRR) for Poisson and Negative Binomial. With Logit, $\exp(\beta)$ is an odds ratio (OR). For Logit and Log, the coefficient table shows the $\exp(\hat\beta)$ column corresponding to this interpretation.

Link Function	Formula	Coefficient Interpretation	Canonical Link For
Identity	$\eta = \mu$	$\beta$ is a difference in the response	Gaussian
Logit	$\eta = \log\!\bigl(\mu / (1 - \mu)\bigr)$	$\exp(\beta)$ is an odds ratio	Binomial
Log	$\eta = \log(\mu)$	$\exp(\beta)$ is a ratio (multiplicative effect)	Poisson
Inverse	$\eta = 1/\mu$	$\beta$ is a change on the $1/\mu$ scale	Gamma
Probit	$\eta = \Phi^{-1}(\mu)$	$\beta$ is a change in $\Phi^{-1}(\mu)$ (the z-value of the probability)	—

What most strongly supports a particular effect scale is knowledge of how the data arise. If each factor scales the rate by a constant multiple, effects add up on the log scale. If the response can be thought of as an unobserved continuous variable crossing a threshold, coefficients read as effects on the scale of that variable, and assuming a normal distribution for it gives Probit (see the threshold model in GLMM Fundamentals).

The canonical link is the default choice when the effect scale does not settle the question. Since $\eta = \theta$ , $X'y$ becomes a sufficient statistic for $\beta$ , and the log-likelihood is concave in $\beta$ . When the design matrix $X$ has full rank and the MLE exists, the solution is unique and IRLS converges stably¹. These are properties of estimation convenience; being canonical does not make the model correct.

Non-canonical links forfeit these properties but may be chosen when the effect scale takes priority. For example, the canonical link for Gamma is Inverse ( $\eta = 1/\mu$ ), which puts coefficients on a $1/\mu$ scale that is hard to interpret. The Log link ( $\exp(\beta)$ as a multiplicative effect) is more commonly used in practice.

A link chosen for the effect scale is not guaranteed to be compatible with the family's domain. For example, Poisson + Identity, chosen to read rate differences, produces a fitted mean outside the domain (positive values) when the linear predictor goes negative for some observation. Whether this happens is determined by the data, not by the combination itself; when it does, MIDAS does not force the mean into range and aborts the fit with an error (see the notes in GLM usage).

Parameter Estimation (IRLS)

GLM parameters $\beta$ are estimated by maximum likelihood. Under regularity conditions (differentiability of the log-likelihood, the true parameter being an interior point of the parameter space, etc.), the estimator is consistent, asymptotically normal, and asymptotically efficient. In general no closed-form solution exists, so IRLS (Iteratively Reweighted Least Squares) is used².

At each iteration, working weights $W$ and an adjusted dependent variable $z$ are computed, then the weighted least squares problem:

\hat\beta^{(t+1)} = (X'W^{(t)}X)^{-1}X'W^{(t)}z^{(t)}

is solved to update $\beta$ . $W$ is a diagonal matrix, and its $i$ -th diagonal entry $W_{ii}$ and the $i$ -th component $z_i$ of $z$ are computed from the current $\hat\mu^{(t)}$ and the link function as:

W_{ii} = \frac{1}{V(\mu_i)\,(d\eta/d\mu)_i^2}, \qquad z_i = \eta_i + (y_i - \mu_i)\,\Bigl(\frac{d\eta}{d\mu}\Bigr)_i

where $V(\mu)$ is the variance function and $d\eta/d\mu$ is the derivative of the link function. For Grouped Binomial, the variance of observation $i$ is $V(\mu_i)/n_i$ , as given by $a(\phi) = 1/n_i$ in the exponential family table, so $W_{ii}$ is additionally multiplied by the trial count $n_i$ . See Nelder & Wedderburn (1972) for the original formulation of IRLS for GLMs. Iteration stops when the maximum absolute change in coefficients falls below the convergence threshold.

With the canonical link, the concavity of the log-likelihood ensures stable convergence. Non-canonical links may lead to slower convergence or convergence failure.

Coefficient Standard Errors and Confidence Intervals

The confidence intervals for the coefficients reported by the GLM tab are Wald intervals $\hat\beta_j \pm c \cdot \operatorname{SE}(\hat\beta_j)$ . The standard error $\operatorname{SE}(\hat\beta_j)$ is the square root of the $j$ -th diagonal entry of the estimated covariance matrix $\hat\phi\,(X'\hat{W}X)^{-1}$ , where $\hat W$ is the IRLS weight matrix $W$ evaluated at the converged $\hat\mu$ .

The estimated dispersion parameter $\hat\phi$ is determined per family. Poisson, Binomial, and Negative Binomial with estimated $r$ use $\hat\phi = 1$ . Gaussian uses the residual deviance divided by $n-p$ ; since the Gaussian deviance is $\sum_i (y_i - \hat\mu_i)^2$ , this value equals $\text{Pearson }\chi^2/(n-p)$ . Gamma and Negative Binomial with fixed $r$ use $\text{Pearson }\chi^2/(n-p)$ . The deviance-based estimator is not used for Gamma because it is not consistent there (McCullagh & Nelder, 1989). What estimating $\hat\phi$ with fixed $r$ means is discussed in Variance Functions and Overdispersion.

The quantile $c$ is taken from the $t$ distribution with $n-p$ degrees of freedom for families whose $\hat\phi$ is estimated from the data, and from the standard normal distribution for families with $\hat\phi = 1$ . For the families that use the $t$ distribution, $c$ reflects the uncertainty in estimating $\hat\phi$ .

For Gaussian + Identity, the coverage probability of this interval (the probability that the interval contains the true value) matches the nominal level in finite samples. For all other family/link combinations, the interval rests on the asymptotic normality of the maximum likelihood estimator, so its coverage probability matches the nominal level only in large samples.

Variance Functions and Overdispersion

As described in the Exponential Family section, the variance function $V(\mu) = b''(\theta)$ is the second derivative of the log-partition function rewritten in terms of $\mu$ . Through the relationship $\operatorname{Var}(Y) = V(\mu) \cdot a(\phi)$ , it defines the mean-variance relationship for each family.

Family	$V(\mu)$	$a(\phi)$	$\operatorname{Var}(Y)$
Gaussian	$1$	$\phi$	$\phi$ (= $\sigma^2$ )
Binomial	$\mu(1 - \mu)$	$1/n_i$	$\mu(1-\mu)/n_i$
Poisson	$\mu$	$1$	$\mu$
Gamma	$\mu^2$	$\phi$	$\mu^2 \phi$
Negative Binomial	$\mu + \mu^2/r$	$1$	$\mu + \mu^2/r$

Poisson and Binomial assume a dispersion parameter $\phi = 1$ . When the actual data variance exceeds this assumption, the condition is called overdispersion. Overdispersion leads to underestimated standard errors and confidence intervals that are too narrow. To diagnose overdispersion, use the Deviance/df ratio shown in the Deviance Goodness-of-Fit section of the GLM Diagnostics tab, opened via View Diagnostics. Under the assumption, this ratio should be close to 1, so a value far from 1 suggests overdispersion. Note that the ratio fluctuates more around 1 when $n - p$ is small.

When overdispersion is detected with Poisson data, switching to Negative Binomial adds a $\mu^2/r$ term to the variance, explicitly modeling the extra dispersion. When $r$ is estimated, overdispersion is absorbed into $r$ , so $\phi = 1$ . When $r$ is fixed, residual overdispersion beyond the fixed $r$ is estimated as $\hat\phi = \text{Pearson }\chi^2/(n-p)$ and reflected in standard errors and confidence intervals.

However, for binary data with $n_i = 1$ (logistic regression), each observation follows Bernoulli $(\mu_i)$ , and once the mean $\mu_i$ is fixed the marginal variance $\mu_i(1-\mu_i)$ is determined as well. There is no degree of freedom in the per-observation variance, so there is simply nothing to compare against to say "the data variance exceeds the theoretical variance." This is why Pearson $\chi^2$ and deviance cannot detect overdispersion at the individual level. This does not mean overdispersion is absent — only that it cannot be detected from the same data; extra dispersion arising from clusters or repeated measurements (e.g., treating multiple patients from the same hospital as independent) can still exist and must be handled separately (see the glossary). Classical overdispersion diagnostics and remedies are meaningful only for grouped Binomial data with $n_i > 1$ .

For grouped Binomial overdispersion, MIDAS does not currently support quasi-binomial or Beta-Binomial alternatives. If the extra dispersion arises from cluster structure, introducing random effects via GLMM is an option. When overdispersion is suspected, check the estimated dispersion parameter and consider that standard errors and confidence intervals may be underestimated.

Confidence Intervals for the Mean Response and Prediction Intervals

Mathematical background for the confidence intervals for the mean response and the prediction intervals for new observations computed by the GLM prediction feature.

In the formulas below, $\hat\phi$ is the estimated dispersion parameter defined in Coefficient Standard Errors and Confidence Intervals, and the quantile $c$ follows the rule in that section. $h_i = x_\text{new}' (X'\hat W X)^{-1} x_\text{new}$ is the leverage of the prediction point, measuring how far its predictor values are from the center of the training data.

The confidence interval for the mean response is the Wald interval on the link scale, $\hat\eta \pm c\sqrt{\hat\phi\, h_i}$ , transformed to the response scale via $g^{-1}$ . Because $g$ is monotonic, the transformation does not change the coverage probability of the interval. With a non-linear link, the transformed interval is not symmetric around $\hat\mu$ . The accuracy of this interval is the same as for the coefficient confidence intervals.

The prediction interval covers the variability of a new observation $Y_\text{new}$ . This variability cannot be captured by transforming a link-scale interval, so the computation depends on the family. All families except Gaussian + Identity use a plug-in method. A plug-in method treats the estimated parameters as if they were the true values and computes the interval from those values. Unlike a confidence interval, it does not account for parameter estimation uncertainty.

Gaussian + Identity: The analytical formula $\hat\mu \pm t_{n-p} \sqrt{\hat\phi(1 + h_i)}$ accounts for both the variance of a new observation ( $\hat\phi$ ) and the estimation uncertainty of the mean ( $\hat\phi \cdot h_i$ )
Gaussian + non-Identity link: A plug-in method is used: $\hat\mu \pm t_{n-p} \sqrt{\hat\phi}$ . The non-linear link transformation prevents exact incorporation of estimation uncertainty on the $\mu$ scale in closed form. As a result, prediction intervals for this combination are a simplified form that does not reflect estimation uncertainty, meaning predictions at the center of the data and at extrapolation points receive the same interval width. Building the interval on the link scale and back-transforming, as well as a first-order delta-method approximation, both weaken coverage guarantees under non-linear links or in small samples, so MIDAS uses neither
Poisson, Binomial, Gamma, Negative Binomial: A plug-in quantile method computes the quantiles of the fitted distribution directly:
- Poisson: quantiles of Poisson $(\hat\mu)$
- Binomial: quantiles of Binomial $(n_\text{new}, \hat\mu)$ with success probability $\hat\mu$ and trial count $n_\text{new}$ (for grouped Binomial, each row's value of the training Trials column in the prediction dataset; $1$ for binary data), divided by $n_\text{new}$ so the interval is reported on the success-proportion scale
- Gamma: quantiles of the Gamma distribution with mean $\hat\mu$ , shape $\alpha = 1/\hat\phi$ , and scale $\hat\phi \cdot \hat\mu$
- Negative Binomial: quantiles of the Negative Binomial distribution with mean $\hat\mu$ and $r$ ( $\hat r$ in automatic estimation mode, the specified value in fixed mode)

For discrete distributions (Poisson, Binomial, Negative Binomial), quantiles are rounded conservatively to the smallest integer $k$ satisfying $P(X \le k) \ge \alpha$ . Under the plug-in distribution, the coverage probability therefore meets or exceeds the nominal level. For individual Binomial data ( $n_i = 1$ ), the only quantile candidates are $\{0, 1\}$ , limiting the informativeness of the interval.

The plug-in methods do not account for parameter estimation uncertainty, so the actual coverage probability may fall below the stated confidence level, particularly in small samples or for predictions far from the observed data range.

References

Nelder, J. A., & Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society: Series A, 135(3), 370-384. https://www.jstor.org/stable/2344614
McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall.

The exception where the MLE does not exist is complete separation—when a linear combination of predictors perfectly separates the response—where the log-likelihood remains concave and $X$ full-rank but has no finite maximum. Binary logistic regression is the canonical example, and similar cases arise in other discrete-response models such as multinomial logit. MIDAS detects separation and reports a warning (see the convergence issues section in GLM usage). ↩
The exception is Gaussian + Identity: $V(\mu)=1$ and $d\eta/d\mu=1$ make $W=I$ and $z=y$ in the formulas in the main text, so the weights do not depend on the data and IRLS reaches the OLS solution $\hat\beta = (X'X)^{-1}X'y$ in a single iteration from any starting point. ↩

Also available as a Markdown file.