Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,359 characters · 15 sections · 76 citation commands
Weighted-Average Least Squares for Negative Binomial Regression
\thispagestyle{empty}
\startlist[main]{lof} \startlist[main]{lot}
\pagestyle{main}
In many empirical applications, model uncertainty emerges for a variety of reasons. For example, competing theories exist that can describe the data, or different assumptions are imposed on the data-generating process (DGP). The two most common approaches for dealing with model uncertainty are model selection and model averaging. In model selection, the user selects the best performing model according to an estimation criterion and then carries out inference based on the chosen model. This approach is problematic because the uncertainty in the initial model selection step is often ignored, which could lead to overly confident decisions and predictions steel2020ma. In contrast, model averaging accounts for model uncertainty by averaging over a set of candidate models, typically aiming at improving predictive accuracy ando2014ma.
As datasets become larger, researchers commonly find themselves in high-dimensional settings with many potential covariates to model their response variable. Choosing appropriate regressors is particularly difficult in these situations because the number of candidate models grows exponentially with the number of regressors, i.e.\ for $k$ regressors, $2^k$ different subsets exist that may be considered as candidates. For the same reason, managing the model space and computational resources is key to applying model averaging procedures in the presence of covariate uncertainty. Bayesian model averaging (BMA) provides two general approaches: 1.\ Markov chain Monte Carlo methods (MCMC) and 2.\ non-MCMC approximation methods, see e.g.\ hoeting1999bma for an early overview. A common solution adopted in frequentist model averaging (FMA), e.g.\ in zhang2016optmaglm, is to prescreen for a viable set of models. In contrast, weighted-average least squares (WALS), first proposed by magnus2010growth for the linear regression model and then extended by deluca2018glm to generalized linear models (GLMs), omits a preselection of models by combining Bayesian and frequentist aspects and, especially, leveraging a semiorthogonal transformation of the regressors allowing for fast computation times. Earlier work by heumann2010logit generalizes WALS to logistic regression using a similar transformation as in deluca2018glm.
Most of the literature, particularly in economics, has focused on model averaging for linear regression models. However, many interesting applications require nonlinear models, e.g.\ classification, count data modeling and survival analysis. The negative binomial (NB) distribution, especially of type 2 (NB2), is a popular distribution featuring overdispersion for count data regression, see e.g.\ cameron1986doctor and cameron1988health for applications in health economics. Notably, the NB2 regression model is not a GLM when its dispersion parameter is estimated from the data. deb2002rand extend it to hurdle and finite mixture models and greene2008nbp develops a more general form, called NBP, which encompasses the NB of type 1 and 2.
Despite its wide application, very limited literature exists on model averaging methods for the NB regression model that jointly estimate the regression coefficients and the dispersion parameter. One of the few open-source packages for model averaging is BMA by raftery2020bma, which currently supports BMA for GLMs and survival models. Hence, it is only able to fit an NB2 with pre-specified dispersion parameter, which is a GLM.
In this paper, I extend WALS GLM by deluca2018glm to the NB2 regression model (WALS NB) to account for covariate uncertainty in the specification of the linear predictor. WALS is particularly well suited as it elegantly circumvents a preselection of models by transforming the regressors, allowing me to focus on the averaging procedure. Analogous to deluca2018glm, I first derive the one-step maximum likelihood estimator based on a Taylor expansion of the NB2 log-likelihood function and then employ a transformation akin to the semiorthogonal transformation used in WALS GLM.
At the time of writing, the asymptotic distribution of the WALS estimator for GLMs is still an open research topic and its variance estimator has been a subject of debate. Recent work by deluca2022sampling proposes a new estimator for the variance of WALS in the linear regression model instead of the Bayesian posterior variance that has traditionally been used. deluca2023interval further analyze the confidence and prediction intervals of WALS in the linear model and propose a new simulation-based method that corrects for bias in the WALS estimator. In contrast, this work focuses on the predictive power of model averaging and leaves the challenging issue of inference (after model averaging) for future research. Model averaging estimators typically improve the predictive accuracy compared to using a single model. For example, in an early application of BMA, madigan1994occam find that BMA achieves better logarithmic predictive score than any single model. Moreover, min1993bma show that the expected squared error loss of predictive mean forecasts is always minimized by BMA, if the data-generating model is included in the model space considered for averaging. In this paper, I compare the proposed WALS NB method to traditional maximum likelihood (ML) estimation of the NB2 regression model in a simulation experiment using the classical precision measure, root mean squared error (RMSE), and scoring rules gneiting2007scores as measures for the distributional fit. Finally, the method is also compared to the lasso estimator wang2016pencount in an empirical application on modeling doctor visits. Both the simulation experiment and the empirical application show that WALS NB improves on the ML estimator in sparse situations with few observations and many covariates. In the latter, its fit is competitive with lasso while being computationally more efficient.
The setup and derivation of WALS NB mostly follow the steps in deluca2018glm for WALS GLM. Assume that data $y_i, i = 1, 2, \dotsc, n$, are conditionally independent given $k$-dimensional regressors $x_{i}$ and follow an NB2 distribution with mean $\mu_i$ and dispersion parameter $\rho$, i.e.\ $y_{i} | x_{i} \sim \mathrm{NB2}(\mu_i, \rho)$. As in the standard GLM setup, I model the mean using an inverse link function $h$ on $\mu_i := \mu(\eta(\beta, x_i)) = h(\eta(\beta, x_i))$ with linear predictor $\eta_i := \eta(\beta, x_i) = x_{i}^{\top} \beta$ and regression coefficients $\beta$. The NB2 distribution has the probability mass function
where $\Gamma$ is the gamma function, and its conditional variance is given by
A distribution from the exponential family has the following density
where $b$ and $l$ are known functions. Typical formulations as in e.g.\ fahrmeir2013regression include a dispersion parameter which, without loss of generality, I set equal to one. Moreover, the following two identities hold for the mean and variance
For WALS estimation, I rewrite the NB2 probability mass function into a similar form as the exponential family with a $\log$-link on $\rho$ by using
Thus, the probability mass function becomes
where I dropped the dependence of $\theta$ on $\mu$ and $\rho$, and of $\rho$ on $\alpha$ for notational brevity. From the last line, we can identify the following building blocks of the exponential family:
Thus, for fixed $\rho$, the NB2 is a member of the exponential family and leads to a GLM. However, $\rho$ is estimated from the data in the WALS procedure and, hence, the underlying model is not a GLM anymore. Furthermore, I separate $l(y_i, \rho)$ into two terms using
so the NB2 probability mass function can be rewritten as
which will simplify the derivation of the WALS estimator later.
I allow for uncertainty in the specification of the linear predictor while assuming that the (conditional) probability mass function of $y_i$ and the inverse link $h$ are correctly specified. First, collect over all observations $n$ the response $y_i$ to an $n$-vector $y$ and the regressors $x_i$ to an $n \times k$ matrix $X$ that contains $x_{i}^{\top}$ as $i$th row. Then, partition the regressors into focus and auxiliary regressors $X = (X_{1}, X_{2})$, where $X_{p}$ is an $n \times k_{p}$ matrix with $i$th row equal to $x_{ip}^{\top}, p = 1, 2$, and $k_1 + k_2 = k$. Further, let $\beta = (\beta_{1}^{\top}, \beta_{2}^{\top})^{\top}$ so the linear predictor can be expressed as $\eta_i = x_{i1}^{\top} \beta_1 + x_{i2}^{\top} \beta_2$. Stacking the linear predictors over all $n$ observations then gives the vector $\eta(\beta) = X_1 \beta_1 + X_2 \beta_2$.
Consider averaging over models containing all focus regressors $X_1$ but arbitrary subsets of the $k_2$ auxiliary regressors in $X_2$, which leads to a total of $2^{k_2}$ possible models. The $j$th model is represented by the restriction $R_{j}^{\top} \beta_{2} = 0$, where $R_{j}$ denotes a $k_2 \times r_j$ matrix of rank $0 \leq r_j \leq k_2$, such that $R_{j}^{\top} = (I_{r_j}, 0)$ or column-permutations thereof. Thus, the matrix $R_j$ specifies which auxiliary regressors are excluded from the $j$th model and its rank $r_j$ denotes the number of excluded auxiliary regressors. Note that $0$ represents a scalar, vector or matrix filled with zeroes of matching dimension unless otherwise stated. For example, $0$ in $R_{j}^{\top} = (I_{r_j}, 0)$ is an $r_{j} \times (k_2 - r_{j})$ matrix.
I start with the classical maximum likelihood estimator of the NB2 regression model. Under conditional independence, the (conditional) $\log$-likelihood is
where $b_i := b(\theta_i, \rho)$ and $a_i := a(y_i, \rho)$. In the following, I will generally omit the dependence of $\theta$, $\mu$, $\eta$, $\rho$, $b$ and $a$ on their parameters to reduce clutter. Moreover, only the $\log$-link is considered for the mean (and dispersion) parameter, i.e.\ $h(\eta_{i}) = \exp(\eta_{i})$, but the general notation using $h$ is retained in many places below to facilitate comparisons with WALS GLM by deluca2018glm and to allow easier extension of the method to other link functions in the future.
The score functions follow:
with
where $\dig(x) := \partial \log \Gamma(x) / \partial x$ is the digamma function. Furthermore, let $H(\beta, \alpha)$ be the negative Hessian of the log-likelihood, which is composed of several submatrices that are listed below. The first components are
where
The next submatrices are defined as
They further simplify thanks to
so $H_{p\alpha}(\beta, \alpha)$ may be rewritten as
using $c_{i} := c(\eta_{i}, \rho) = \partial^2 \theta_i / \partial \eta_{i} \partial \rho$. Finally, the last part is
with
where $\trig(x) := \partial^2 \log \Gamma(x) / \partial x^2$ is the trigamma function, and
The ML estimator for the $j$th model solves the following constrained optimization problem
As a first step towards the solution, I construct the Lagrangian
where $\nu_{j}$ denotes the $r_{j}$-vector of Lagrange multipliers. Setting the first derivatives equal to zero yields the system of nonlinear equations
Following deluca2018glm, I consider a one-step ML estimator that approximates the solution of the system. In contrast to iterative procedures such as Newton-Raphson, which are typically used for solving nonlinear equation systems, the one-step ML estimator admits closed-form expressions.
In the remainder of the paper, I assume that all necessary conditions for the algebraic manipulations, e.g.\ rank conditions on the regressor matrix $X$, are satisfied. Detailed proofs are found in (ref).
I expand the estimating equations of (ref) (except for $R_{j}^{\top} \beta_2 = 0$) around starting values $\bar{\beta} = (\bar{\beta}_{1}^{\top}, \bar{\beta}_{2}^{\top})^{\top}$ and $\bar{\alpha}$. Further, $\bar{\rho} = \rho(\bar{\alpha})$, since the mapping from $\alpha$ to $\rho$ is strictly monotonic (log-link). Using a first-order Taylor expansion and ignoring the remainder term yields
where $\bs_{p} := s_{p}(\bbeta, \balpha), \bH_{pq} := H_{pq}(\bbeta, \balpha)$ and $\bH_{p\alpha} = H_{p \alpha}(\bbeta, \balpha)$, $p = 1, 2$. In the following, all quantities evaluated at $(\bbeta, \balpha)$ are denoted by a bar and the approximations in (ref) are treated as equalities for a simpler notation.
First, consider the unrestricted model with $R_u = 0$. Define the data transformations
which involve the $n \times n$ matrices
and the $n$-vectors $\bmu := \mu(\bar{\eta}) := \left(h(\bar{\eta}_{1}), h(\bar{\eta}_{2}), \dotsc, h(\bar{\eta}_{n})\right)^{\top}$ and $\bar{\eta} := \left(\bar{\eta}_{1}, \bar{\eta}_{2}, \dotsc, \bar{\eta}_{n}\right)^{\top} = X_1 \bbeta_1 + X_2 \bbeta_2$ with $\bar{\eta}_{i} := \eta(\bar{\beta}, x_i)$. Using the log-link further guarantees $\rank(\bPsi) = n$ because $ \psi(\bar{\eta}_{i}, \bar{\rho}, y_i) = \bar{\mu}_{i} \bar{\rho} (y_i + \bar{\rho}) / (\mu_{i} + \bar{\rho})^{2} > 0 $ since $\bar{\mu}_{i} > 0$, $\bar{\rho} > 0$ and $y_{i} \geq 0$ for all $i$. Moreover, define
with $n$-vectors
where $\bk_{i} := k(\bar{\eta}_{i}, \bar{\rho}, y_i)$, $\bkappa_{i} := \kappa(\bar{\eta}_{i}, \bar{\rho}, y_i)$ and $\ones := (1, \dotsc, 1)^{\top}$ is an $n$-vector filled with ones. Notice the slight abuse in notation, where $\mu$, $k$ and $\kappa$ are vector-valued functions here, whereas they were scalar-valued in the sections before. Furthermore, let
Then, the solution to the linearized system of likelihood equations (ref) can be expressed in closed form as
where
is a symmetric matrix. In contrast to deluca2018glm, $\bM_1$ is not idempotent anymore due to the rank-1 perturbation in $I_n + \bar{g} \bepsilon \bPsi^{-1/2} \bq \bq^{\top} \bPsi^{-1/2}$, which is a consequence of the additional dispersion parameter $\rho$ in the NB2 model compared to GLMs.
Likewise, consider the general one-step ML estimator for the $j$th model. Define the symmetric and idempotent $k_2 \times k_2$ matrix
the $k_1 \times k_2$ matrix
and the following transformation of the unrestricted one-step ML estimator $\tbeta_{2u}$
Then, analogous to Proposition 1 of deluca2018glm, I obtain the one-step ML estimator for the $j$th model in the following proposition.
The WALS NB estimator relies on a preliminary transformation of the auxiliary regressors to reduce the computational burden, akin to WALS for the linear regression model magnus2010growth and GLMs deluca2018glm.
First, scale the focus regressors by defining
with the $k_1 \times k_1$ diagonal matrix $\bDelta_1 := \diag\left(\xmat{\bX_1}{\bX_1} / n \right)^{-1/2}$ such that $\diag\left( \xmat{\bZ_1}{\bZ_1} / n \right) = (1, \dotsc, 1)$. The only purpose of the transformation is to improve the numerical accuracy by normalizing all regressors to be the same scale in $\bZ_1$ deluca2018glm. It further implies
so scaling by $\bDelta_1$ has no effect on $\bM_1$. Next, transform the auxiliary regressors by
where
and I assumed $\bX_{2}^{\top} \bM_1 \bX_2$ to be positive definite so $\bXi^{1/2}$ exists. Furthermore, the $k_2 \times k_2$ diagonal matrix $\bDelta_{2} := \diag\left(\bX_{2}^{\top} \bM_{1} \bX_2 / n \right)^{-1/2}$ is chosen such that $\diag(\bXi) = (1, \dotsc, 1)$. Unlike the matrix $\bDelta_1$, the transformation by $\bDelta_2$ serves the dual purpose of improving numerical accuracy and making the WALS NB estimator equivariant to scale transformations of the auxiliary regressors. Otherwise it would be only scale equivariant for the focus regressors deluca2011stata.
Notice that combining (ref) and (ref) leads to
In contrast to deluca2018glm, $\bM_1 \bZ_2 / \sqrt{n}$ is not semiorthogonal\footnote{Semiorthogonality is defined as $A A^{\top} = I$ or $A^{\top} A = I$ for a general (non-square) matrix $A$ zhang2017matrix.} anymore, since $\bM_1$ is not idempotent. The transformation further implies
Using (ref) and (ref) I can show for the unrestricted model that
so the linear predictor stays the same for the unrestricted model. Therefore, all the quantities that only depend on $\bar{\alpha}$ and indirectly on $\bar{\beta}$ via
where $\bgamma_{1} = \bDelta_{1}^{-1} \bar{\beta}_1$ and $\bgamma_{2} = \bXi^{1/2} \bDelta_{2}^{-1} \bar{\beta}_2$, remain the same (e.g.\ $\bar{\mu}$, $\bPsi$, $\bQ$, $\bar{g}$, $\bar{t}$, \ldots), as they do not depend on $\bar{\beta}$ directly. Note that deluca2018glm suggest using the fully iterated unrestricted ML estimates as starting values $\bar{\beta}$ and $\balpha$. In this case, the starting value of the dispersion parameter $\balpha_{Z}$ for the transformed regressors $Z$ is identical to $\balpha$ for the original regressors $X$ since the estimated conditional means are equal, i.e.\ $h(\eta(\bar{\beta}, x_i)) = h(\eta(\bgamma, z_i))$ for all $i$, where $z_{i}^{\top}$ is the $i$th row vector of $Z = (Z_1, Z_2)$.
It follows from (ref) using (ref) that the one-step ML estimators for the $j$th transformed model are given by
where the fully restricted and unrestricted estimators are
with
Exploiting $R_{j}^{\top}R_{j} = I_{r_j}$, the following terms simplify
Analogous to $\tvartheta$, using (ref) yields
As a direct consequence of (ref), both $P_{j}$ and $W_{j}$ become nonrandom projection matrices that are different from $\bP_j$ and $\bW_j$ used for the estimation with the untransformed regressors. Furthermore, $W_{j}$ reduces to a diagonal matrix with $k_2 - r_j$ ones and $r_j$ zeros on its main diagonal. The $h$th diagonal element of $W_{j}$ is zero, when the $h$th component of $\gamma_{2}$ is constrained to be zero in the $j$th model. Otherwise, the $h$th component is one. Combining this observation with $\tgamma_{2j}$ from (ref), it follows that all models that include the $h$th column of $Z_2$ as regressor will have the same estimator for the $h$th component, namely the $h$th component of $\tgamma_{2u}$.
Note that the $j$th model for the transformed regressors is generally not equivalent to the $j$th model of the untransformed regressors because the restriction in (ref) differs. The exceptions are the unrestricted model $u$ and the fully restricted model $r$, where the restriction is irrelevant:
This implies that $\tgamma_{2j} \neq \bXi^{1/2} \bDelta_{2}^{-1} \tbeta_{2j}$ for $j \notin \{u, r\}$ and $k_2 \geq 2$ auxiliary regressors. For $k_2 = 1$, there exist only two models: 1.\ the unrestricted and 2.\ the fully restricted model, so $j \in \{u, r\}$. The results are summarized in (ref).
Consider the model averaging estimators of $\gamma_1$, $\gamma_{2}$ and $\alpha$
where $\lambda_j$ are data-dependent model weights satisfying the restrictions
Note that the regularity condition $\lambda_j = \lambda_j(\sqrt{n} \tgamma_{2u})$ is equivalent to the condition on the model weights used by hjort2003ma.
From (ref) I get
where $W = \sum_{j=1}^{2^{k_2}} \lambda_j W_{j}$ is a diagonal matrix with entries $w_{h} \in [0, 1]$, because $W_{j}$ is a diagonal matrix with entries $w_{j,h} \in \{0, 1\}$, $h = 1, 2, \dotsc, k_2$ (notice the slight abuse of notation: $h$ is used as an index here and does not refer to the inverse link). Next, I can transform $\hat{\alpha}$ to an estimate for $\rho$ by applying the inverse of the log-link, i.e.\ $\hat{\rho} = \exp(\hat{\alpha})$. Furthermore, using $\hat{\gamma}_1$ and $\hat{\gamma}_2$, the WALS estimators of the original parameters $\beta_1$ and $\beta_{2}$ are given by
The final step in completing the WALS NB model averaging estimator is to estimate the model weights $\lambda_j$. However, notice that both $\hat{\gamma}_1$ and $\hat{\alpha}$ can be expressed as functions of $\hat{\gamma}_2$. Therefore, it is sufficient to find an expression for $\hat{\gamma}_2$ instead of directly estimating the weights $\lambda_j$. Similar to deluca2018glm, I construct $\hat{\gamma}_{2}$ as a Bayesian shrinkage estimator by exploiting the approximate normality and independence of $\tgamma_{2u}$ under the local misspecification framework hjort2003ma. First, let the auxiliary parameters be $\beta_{2} = \delta / \sqrt{n}$, where $\delta$ is an unknown constant vector that represents the departure of the DGP from the unrestricted model. Then, if the fully iterated ML estimator of the unrestricted model is used as starting values $\bar{\beta}_1$, $\bar{\beta}_2$ and $\bar{\alpha}$ and mild regularity conditions are assumed, I can show that
in large samples, where $\gamma_{2n} = d / \sqrt{n}$, $d$ = $\Xi^{1/2} \Delta_{2}^{-1} \delta$, $\Xi = \plim \bXi$ and $\Delta_{2} = \plim \bDelta_{2}$ (see the supplementary materials\ for more details). Further, consider $\hat{\gamma}_2$ from (ref) and assume analogously to deluca2018glm that each diagonal element $w_h, h = 1, 2, \dotsc, k_2$, of $W$ only depends on the $h$th component $\sqrt{n} \tgamma_{2u,h}$ of $\sqrt{n} \tilde{\gamma}_{2u}$. Then, (ref) implies that the components of $\hat{\gamma}_{2}$ are also approximately independent. This assumption further simplifies the estimation problem by reducing the $k_2$-dimensional problem of estimating $\hat{\gamma}_{2}$ to $k_2$ times a one-dimensional problem of estimating each element of $\hat{\gamma}_{2}$. Moreover, $\hat{\gamma}_{2,h}$ is a shrunken version of $\tgamma_{2u,h}$ because $0 \leq w_h \leq 1$, therefore, $\hat{\gamma}_{2}$ is a shrinkage estimator of $\gamma_{2n}$.
The previous two observations suggest that the Bayesian posterior mean is a suitable shrinkage estimator for $\sqrt{n} \gamma_{2n,h}$. Thus, the $h$th component of the WALS NB estimator $\hat{\gamma}_{2}$ follows as
where $\sqrt{n} \tgamma_{2u,h} \approx \normal(d_h, 1)$ with prior mean $d_h$, which is the $h$th element of $d$ and is assumed to have a symmetric and unimodal prior $f$ (see section 9 of magnus2016wals for more details on the prior and the estimation). Notice again that $\hat{\gamma}_{2,h}$ lies between 0 and the `observed data' $\tgamma_{2u,h}$.
magnus2016wals require the desirable properties of robustness\footnote{A prior $\pi(\gamma)$ is robust if the posterior mean $m(x)$ based on $\pi$ satisfies $x - m(x) \rightarrow 0$ as $x \rightarrow \infty$.}, neutrality\footnote{A prior $\pi(\gamma)$ is neutral if the prior median of $\gamma$ is zero and the prior median of $|\gamma|$ is one.} and minimax regret\footnote{Regret is defined as difference between risk and the infimum of risk, where risk is defined as expected squared loss.} for the prior $f$, which further motivates the use of the Bayesian posterior mean as the shrinkage estimator in $\hat{\gamma}_{2,h}$. The reflected Weibull, under suitable parameter values, is a prior that fulfills all the properties mentioned above. In contrast, the Laplace prior is neutral but not robust magnus2016wals. However, it admits a closed-form expression for the posterior mean in (ref) magnus2010growth and therefore calculating the posterior mean under the Laplace prior is computationally less complex than under the reflected Weibull, which requires numerical integration.
In order to compare the performance of WALS NB with other methods, I first need to define performance metrics. The classical performance measure for regression is the RMSE, which is given by
where $\hat{\mu}_i$ is the predicted mean for observation $i$. However, I would like to evaluate the fit of an entire distribution and not only the expectation. Traditional measures used in machine learning such as (R)MSE only focus on point predictions, i.e.\ the conditional expectation of the fitted distribution, in relation to the observed values and do not make judgment on other aspects of the fitted distribution. czado2009predictive recommend scoring rules for evaluation of count data models, which have also been used in kolassa2016count. WALS NB and all other methods considered in this paper fit an entire (conditional) distribution for each individual that allows probabilistic predictions/forecasts, which is exactly the scenario for which scoring rules provide quality assessment gneiting2007scores. For count data, a probabilistic forecast is a predictive probability distribution $\hat{P}$ on the set of nonnegative integers $\mathbb{N}_{0}$ czado2009predictive.
Following czado2009predictive, I take scoring rules to be penalties I wish to minimize. Specifically, the penalty $s(\hat{P}, y)$ is incurred when the forecaster quotes predictive distribution $\hat{P}$ and count $y$ is realized. Moreover, let $s(\hat{P}, Q)$ denote the expected value of $s(\hat{P}, \cdot)$ under distribution $Q$
In practice, the average over suitable pairs $(\hat{P}, y)$ is used:
where $\hat{P}_i$ refers to the $i$th predictive distribution and $y_i$ the $i$th observed count. In the simulation experiment and empirical application of (ref), respectively, scores will always refer to a suitable average.
Suppose the forecaster has predictive distribution $Q$ available. Then the forecaster has no incentive to predict any $\hat{P} \neq Q$ and is encouraged to quote her true belief, $\hat{P} = Q$, if the scoring rule is strictly proper. Strict propriety is defined by
with equality if and only if $\hat{P} = Q$, and encourages honest quotes (czado2009predictive; gneiting2007scores). If $s(Q, Q) \leq s(\hat{P}, Q)$ for all $\hat{P}$ and $Q$, then the scoring rule is only proper. Since only strict propriety ensures that both calibration (consistency with actual realizations) and sharpness (concentration of the predictive distribution) of the predictive distribution are addressed winkler1996scores, I exclusively use strictly proper scoring rules.
czado2009predictive propose a number of strictly proper scoring rules for count data. It is a priori unclear which scoring rule to use unless there is a unique and clearly defined underlying decision problem. Since probabilistic forecasts often have many uses, it is appropriate to use a variety of scores to take advantage of their differing emphases czado2009predictive. In this paper, I use the logarithmic (log), Brier and spherical score, which I briefly summarize here: Let $\hat{p}_y := \hat{P}(Y = y)$ denote the probability mass at count $y$ (for continuous distributions it is the density at $y$), then the log score is defined as
The sum of log scores corresponds to the negative log-likelihood. Further define
where the infinite sum may be truncated if no closed-form expression exists. The quadratic score, also called Brier score, is then
The spherical score uses the same components differently:
The aim is to compare the performance of WALS NB with the traditional ML estimator of the NB2 regression model in a controlled environment. The DGP is inspired by the local misspecification framework so I can assess the influence of varying numbers of focus and auxiliary regressors.
The dependent count variable is sampled from an NB2 using a log-link, i.e.\
for $i = 1, 2, \dotsc, n$, where $x_i = (x_{i1}^{\top}, x_{i2}^{\top})^{\top}$ is a random vector of dimension $k = k_1 + k_2$ composed of $k_1$ focus regressors $x_{i1}$ and $k_2$ auxiliary regressors $x_{i2}$. Analogously, the coefficient vector is separated into two parts: $\beta = (\beta_{1}^{\top}, \beta_{2}^{\top})^{\top}$. Inspired by the simulation experiments in zhang2019inference and deluca2023interval, who compare confidence and prediction intervals of model averaging methods for the linear regression model, I choose the regressors to be multivariate normal because it allows me to analyze the effect of the correlation between the regressors on the performance of the methods. For simplicity, I specify each element of $x_i$ to have variance $1$ and pairwise correlation $b$, i.e.\ \[ \Sigma_{k} =
. \] The same offset $\alpha = \log(3)$ is used in all experiments such that the DGP produces reasonable counts, see (ref) for a visualization of a training set from a specific run.
Moreover, the regression coefficients are generated as follows: Define the vectors $\bar{\beta}_{1} := (\bar{\beta}_{1,1}, \bar{\beta}_{1,2}, \dotsc, \bar{\beta}_{1,10})^{\top}$ and $\bar{\beta}_{2} := (\bar{\beta}_{2,1}, \bar{\beta}_{2,2}, \dotsc, \bar{\beta}_{2,100})^{\top}$. Then, the maximum number of regression coefficients $k_1 = 10$ and $k_2 = 100$ are randomly sampled once according to the following rules:
The simulations then only take the first $k_1$ and $k_2$ values from these vectors as regression coefficients $\beta_1$ and $\beta_2$. For example, in the setting $k_1 = 5, k_2 = 10$, $\beta_1 = (\bar{\beta}_{1,1}, \bar{\beta}_{1,2}, \dotsc, \bar{\beta}_{1,5})^{\top}$ and $\beta_{2} = (\bar{\beta}_{2,1}, \bar{\beta}_{2,2}, \dotsc, \bar{\beta}_{2,10})^{\top}$. Hence, the magnitude of the elements in $\beta_2$ is much smaller than in $\beta_1$, therefore the regressors $x_{i2}$ are considered auxiliary regressors and the main variation is driven by $x_{i1}$. (ref) and (ref) in (ref) show the entries of $\bar{\beta}_1$ and $\bar{\beta}_2$, respectively.
All values of the parameters used in the experiment are summarized in (ref). A total of 1080 scenarios consisting of all combinations of the parameters are simulated for $R = 300$ runs each.
I compare six different procedures that are named according to the pattern `method-specification'. The two methods are called `walsNB', which estimates the NB2 regression model using WALS NB, and `ML', which uses maximum likelihood. For WALS NB procedures, the Weibull prior is used as it theoretically provides the best tradeoff between robustness and regret, for more details see magnus2016wals. The results for other priors are expected to be quite similar as WALS for the linear regression model has empirically shown to be relatively insensitive to the choice of the prior deluca2022sampling.\footnote{I also conducted the simulation experiment using the Laplace prior and the results are similar to the ones using the Weibull prior.}
The procedures considered are
The second WALS NB specification, walsNB-aux, is included to analyze the extent to which prior information about the focus regressors in walsNB-dgp affects performance. Ideally, including $x_{i1}$ as focus regressors in the procedure should improve performance as they are the covariates that dominate and should therefore be included in all submodels of WALS NB. However, in walsNB-aux their coefficients are also subjected to the regularization of the Bayesian estimation step, which may improve performance. Thus, a priori it is unclear which model will dominate.
All ML specifications are estimated using a log-link for the mean parameter, while the dispersion parameter is estimated directly without a link (default setting). Moreover, I increase the maximum number of iterations for both the alternation process between IRLS and ML estimation of $\rho$ and the IRLS algorithm itself from the default setting of 25 to 2500 to increase the odds for convergence. The remaining settings, e.g.\ convergence criteria, are left at their default values.
Moreover, the WALS NB specifications use a $\log$-link for the mean and dispersion parameter following the DGP (ref) and are initialized using the ML estimates of the unrestricted model, which are given by the ML-U procedure (using the increased maximum number of iterations as described above). This initialization is recommended by deluca2018glm as it produces lower RMSE for the WALS GLM estimator in their Monte Carlo simulations deluca2018glm compared to using the estimates of the fully restricted model as starting values. It further ensures that (ref) approximately holds for the one-step estimators of the auxiliary regression coefficients, which I exploit in the Bayesian estimation step to reduce the $k_2$-dimensional posterior mean estimation to $k_2$ one-dimensional problems (see (ref)).
Finally, the Weibull prior for all WALS specifications uses the parameters recommended in magnus2016wals, which are minimax regret solutions for the normal location problem.
In order to compare the performance of the procedures, I follow the benchmark experiments framework of hothorn2005design and more specifically the `Simulation Problem'. The simulation is structured to emulate the typical use of the methods: For each scenario and run, a training sample of size $n$ is generated, where all procedures are applied and performance criteria are computed on an independently generated validation set that is fixed in size to $n_e = 4000$ to avoid any variation due to its size. I do not employ hypothesis tests to check if the performance differences are significant because the simulation experiment itself is already computationally intensive due to the large number of parameter settings.
I consider the RMSE as a classical precision measure for regression and additionally log, Brier and spherical scores to assess the distributional fit as described in (ref). For the scoring rules, the average is taken over the validation sample as in (ref). Further, I truncate the infinite sum in $|| \hat{p} ||$ from (ref) used in the Brier and spherical score at the count $r = 150$ because the response typically does not exceed 150 and it would be meaningless to extrapolate beyond the observed data. See (ref) for a visualization of the training data of a single simulation run of the setting $k_1 = 10, k_2 = 100, \rho = 1, b = 0$ (this setting should maximize the range of $y$, since $k_1 = 10$ and $k_2 = 100$ allow for the largest possible means $\mu_i$ and $\rho = 1$ maximizes variance), where the response only ranges between 0 and 71.
In the following, only the scenarios highlighting the differences between the procedures are discussed. For more results, see the supplementary materials.
First, I analyze how the procedures behave when the number of regressors is varied. In the following plots, the points represent the mean validation metric over all successful runs of the experiment, i.e.\ $R = 300$ if the method never fails to converge. The total number of failed runs is given below the corresponding points in the plots and the shaded area displays the interquartile range (the box of a boxplot) of the validation metric.
(ref) shows that walsNB-aux performs similarly to walsNB-dgp in terms of mean validation RMSE when we vary $k_2$ with fixed $k_1 = 1$, $\rho = 1$ and $b = 0$, because most of the regressors are auxiliary and the former includes all regressors as auxiliary. Moreover, both WALS NB specifications outperform ML-U on average when $k_2 \geq 20$ and $n$ is small ($n \leq 2000$). In fact, WALS NB specifications show lower mean validation RMSE than all ML specifications in these scenarios, except for ML-focus that contains only the focus regressors. For $k_2 = 100$ and $n < 2000$ the `typical' performance of walsNB-dgp and walsNB-aux is also better than ML-U and ML-AC as their interquartile ranges do not overlap. On the other hand, when $k_2$ is small and/or $n$ is large, their interquartile ranges are similar. The largest difference in mean RMSE is observed at $n = 500$, $k_1 = 1$, $k_2 = 100$ where walsNB-aux exhibits around 8.9% lower mean RMSE than ML-U. ML-focus performs the best in all scenarios, especially when $k_2$ is large and $n$ small.
Therefore, if we know the focus regressors, then ML-focus yields the best fit in very sparse situations with few observations. Otherwise, walsNB-dgp and walsNB-aux are better than using all regressors in the large regression model ML-U. The outperformance in walsNB-dgp and walsNB-aux compared to ML-U is likely due to the reduced variance thanks to the Bayesian regularization step, which typically reduces variance and leads to lower RMSE via the bias-variance trade-off. In reality it is unlikely that we can exactly identify which regressors are the focus regressors, so walsNB-aux offers a great alternative that does not require variable selection. In all scenarios, it performs at least as well as ML-U but better when the data is sparse and few observations are available. For large $n$ or small $k_2$, all procedures fit the data equally well as their mean validation RMSE converges to the RMSE of the oracle.
When I vary $k_1$ with fixed $k_2 = 1$, $\rho = 1$ and $b = 0$ in (ref), the picture changes. In all scenarios, ML-AC returns the highest RMSE and the remaining specifications perform similarly as their interquartile ranges overlap. Increasing $k_1$ shifts the `RMSE-curve' up for all procedures, including the oracle, while retaining their relative order. This behavior is explained by the form of the variance of the NB2 distribution in (ref). The more focus regressors with large coefficients are included, the more likely it is that the conditional expectation $\mu_i$ is large, which increases the variation of the response $y_i$ since the conditional variance is monotonically increasing in $\mu_i$. Thus, even if we could exactly estimate the true $\beta_1$ and $\beta_2$, the RMSE would still increase due to the increased conditional variance, which is demonstrated by the behavior of the oracle.
The same patterns hold for the validation log score. Firstly, (ref) shows that WALS NB specifications generally perform better than ML specifications in terms of log score, when the number of auxiliary regressors is high compared to the number of focus regressors and few observations are available. The exception is again ML-focus, which performs the best across all scenarios. The largest difference in mean log score between walsNB-aux and ML-U is realized at $k_1 = 1$, $k_2 = 100$ and $n = 500$ where the mean log score of walsNB-aux is around 3.9% lower. Moreover, the typical performance of walsNB-aux is also better than ML-U in this scenario as their interquartile ranges do not overlap. For large $n$ the distributional fit of all models is similar because the mean log scores converge to the log score of the oracle.
Secondly, similar to the results for RMSE in (ref), I find a small upwards shift of the mean validation log scores in (ref) when increasing $k_1$ given $k_2 = 1$. As expected, the distributional fit of ML-AC, which only includes the auxiliary regressors, is the worst among the procedures when $k_1 \geq 5$. Finally, the interquartile range of all models except ML-AC overlap, so their performance in terms of log score typically does not differ. The relative ranking of the procedures for the Brier and spherical score is the same as for the log score, so their results are only shown in the supplementary materials.
In summary, the WALS NB specifications generally outperform ML-U in terms of RMSE and log score when the number of auxiliary regressors is very large relative to the number of focus regressors and when the number of observations is small. This is in line with the results from abadie2019regularized for the pretest estimator, which is the predecessor of the WALS estimator: The authors consider a `Spike and Normal' process for noisy estimates $\hat{\mathcal{X}}_1, \hat{\mathcal{X}}_2, \dotsc, \hat{\mathcal{X}}_k$ of e.g.\ the coefficients from a linear regression model abadie2019regularized: The estimates $\hat{\mathcal{X}}_j$ are assumed to follow $\hat{\mathcal{X}}_j \sim \normal(m_j, s_j^2)$ for $j = 1, 2, \dotsc, k$, e.g.\ $\hat{\mathcal{X}}_j$ are elements of the ordinary least squares (OLS) estimator in a linear regression model with homoskedastic normal error terms. In this setup, the mean $m_j$ can be regarded as the true value of the regression coefficient that is estimated as $\hat{\mathcal{X}}_j$. The idea is that regularized estimators such as lasso, ridge, and pretest modify the OLS estimator $\hat{\mathcal{X}}_j$. Furthermore, the mean $m_j$ is set to zero (spike) with a fixed probability $p$, and with probability $1-p$ the coefficient follows $m_j \sim \normal(m_0, s_{0}^2)$ for all $j$. Under this setting, the authors show that the pretest estimator exhibits smaller integrated risk (integrated expected squared error over the space of distributions of the data distribution, see abadie2019regularized for details) than lasso and ridge, when the process is very sparse, i.e.\ $p$ is high and $m_0$ is large, so many coefficients are set to zero and the non-zero coefficients are far away from zero. The results further agree with the Monte Carlo simulations of deluca2023interval for WALS in the linear regression model: The authors find that the ratio of the MSE of OLS relative to the MSE of WALS increases when the number of auxiliary regressors $k_2$ becomes larger. Moreover, for all $k_2$, the ratio decreases when the sample size $n$ increases. Both observations are in line with (ref), where walsNB-dgp dominates in terms of (R)MSE compared to the unrestricted estimator ML-U when $n$ is small and $k_2$ is large.
I fixed $k_1 = k_2 = 5$ so that varying the correlation $b$ affects the correlation within focus and auxiliary regressors, as well as the correlation between focus and auxiliary regressors. In contrast, if I had set $k_1 = k_2 = 1$, only the correlation between focus and auxiliary regressors would be modified.
(ref) shows that for fixed $b = 0$ and $k_1 = k_2 = 5$, all procedures yield similar mean validation RMSE across all $\rho$, except for ML-AC, which exhibits much higher values compared to the other methods. Note that the mean validation RMSE generally decreases for all procedures, even the oracle, when $\rho$ increases. This is due to the fact that higher $\rho$ leads to less overdispersion, i.e.\ lower conditional variance, resulting in lower RMSE for all procedures.
Increasing the correlation $b$ between all regressors for fixed $\rho = 1$ and $k_1 = k_2 = 5$ in (ref), the mean validation RMSE shifts down for all procedures, especially for ML-AC as it only includes the auxiliary regressors and a constant. The larger the correlation, the better it can compensate the lack of focus regressors. Generally, the choice of regressors matters less for prediction when the regressors are highly correlated as each of them will contain similar information for the prediction task. The remaining procedures perform very similarly when increasing $b$ and converge to the mean validation RMSE of the oracle for large $n$. Except for ML-AC in the cases with $b < 0.9$, the typical RMSE of the procedures are comparable as the interquartile ranges overlap and have similar widths in all scenarios.
The results for the mean validation log score when varying $\rho$ and $b$ with fixed $k_1 = k_2 = 5$ in (ref) and (ref) are qualitatively the same as for the mean validation RMSE. Interestingly, I also observe a downward shift in the mean validation log score for all procedures and $n$ when I increase $\rho$. The argument used to explain the downward shift for the mean validation RMSE, namely that the variance around the conditional mean is lower the higher $\rho$, does not hold anymore since less overdispersion does not necessarily lead to lower log scores. Intuitively, less overdispersion leads to less variation around the conditional mean that could allow for a more precise estimation of the conditional mean, resulting in an improved distributional fit and, hence, a lower log score.
Finally, (ref) shows the mean validation log scores varying $n$ and the correlation $b$. Similar to the RMSE, the mean validation log scores generally decrease across all $n$, when $b$ is increased. The reduction is especially large for ML-AC due to the same reasons as for the RMSE in (ref). The remaining procedures perform very similarly: Their mean validation log scores are similar and converge to the oracle when $n$ is large and their interquartile ranges overlap.
The relative ranking of the procedures for the Brier and spherical score is similar to that for the log score, so their plots are only shown in the supplementary materials.
The aim of the empirical illustration is to compare the predictive performance of WALS NB with ML and lasso estimation of the NB2 regression model on real data, and to check whether the observations from the simulation experiment translate to a real-world application.
I use the cross-sectional data set called `DoctorVisits', which derives from the 1977-1978 Australian Health Survey and was analyzed in cameron1986doctor and mullahy1997hetero. The dataset contains $n = 5190$ observations from individuals over 18 years of age on twelve variables, including the response visits, which describes the number of doctor visits in a two-week period before the interview. It further provides explanatory variables such as income and age, as well as health-related variables like recent illnesses and health insurance coverage. The data are available via the \proglang{R} package \pkg{AER} kleiber2008aer as DoctorVisits based on the original from the Journal of Applied Econometrics Data Archive.\footnote{\url{https://www.journaldata.zbw.eu/dataset/heterogeneity-excess-zeros-and-the-structure-of-count-data-models}}
(ref) and (ref) in (ref) provide a description and summary statistics of the variables in the DoctorVisits dataset. Further, (ref) shows a visualization of the response visits, which clearly exhibits overdispersion and will be modeled using regression models for count data. For the computation of the Brier and spherical score, I truncate the infinite sum in $||\hat{p}||$ from (ref) at the largest observed count of the dataset, which is 9.
Inspired by the applications of rupp2012mlatom and faber2020qml in quantum chemistry, I apply $K$-fold cross-validation (CV) to produce `learning curves' that allow me to compare the performance of the procedures for different sizes of the training set. (ref) illustrates the process for generating a $K$-fold cross-validated learning curve for any evaluation metric and procedure.
Note that only the training sets $\mathcal{T}_{l,k}$ vary in size but the validation sets $\mathcal{V}_{k}$ remain the same. Similar to meek2002learncurve, the training sets $\mathcal{T}_{l, k}$ are nested, i.e.\ $\mathcal{T}_{l, k} \subset \mathcal{T}_{l + 1, k}$ for $l = 1,2 \dotsc, L-1$. For all experiments below, I use $K = 10$ folds.
I compare procedures that differ in the estimator and specification of the mean, where the choices for the latter are inspired by the applications in cameron1986doctor. The different combinations of estimator and specification are named following the pattern: `estimator-specification'. Again, `walsNB' and `ML' represent the WALS NB and ML estimator, respectively, while `lasso' estimates the NB2 regression model using the lasso estimator of wang2016pencount (see (ref) for details). I consider a total of six estimator-specification combinations:
All procedures are fitted using a log-link for the mean parameter. WALS NB procedures use the Laplace prior because the Weibull led to numerical instabilities in some small subsamples resulting from the numerical integration procedure required for computing the posterior mean of the auxiliary regression coefficients in (ref). The parameters of the Laplace prior are taken from magnus2016wals, which are minimax regret solutions for the normal location problem. The remaining settings for WALS NB and ML specifications are retained from the simulation experiment of (ref). Notably, all WALS NB specifications use the unrestricted ML estimator for NB2 as starting values for the regression coefficients and the dispersion parameter. By unrestricted, I refer to the unrestricted model given the covariates that are included in the specification, i.e.\ ML-main for walsNB-main and ML-int for, both, walsNB-main-focus and walsNB-int.
The lasso specification `lasso-int' performs tuning (maximizing 10-fold CV log-likelihood) in the training set $\mathcal{T}_{l,k}$ of each fold $k$ (and each training set size $t_l$), as recommended by hothorn2005design who include tuning and final model fit in the estimation process. This is sensible, as tuning of the regularization parameter is key to the performance of lasso. Different values of the regularization parameter correspond to different levels of regularization and the regressors included in the model may also differ. Moreover, the method also standardizes the regressors in the training set of each fold to have zero mean and unit variance before estimation (i.e., it uses the estimated means and variances of the regressors in the subsample and not over the entire dataset).
In (ref) and (ref) we observe that all WALS NB specifications except for walsNB-main-focus outperform the ML specifications in terms of RMSE for all numbers of training observations. The differences are particularly large for small training sets, e.g.\ for $t_{l} = 500$ the CV RMSE of walsNB-int is almost 19% smaller than ML-int. Except for $t_{l} < 1500$, walsNB-int and walsNB-main also outperform the lasso specification lasso-int. Further, note that walsNB-int outperforms walsNB-main-focus although the only difference between the two procedures is that the latter specifies some of the covariates as focus regressors. This observation seems to contradict the results of the simulation experiment, where walsNB-dgp and walsNB-aux perform very similarly even though the latter considers all covariates as auxiliary regressors and the former considers part of them as focus regressors. However, walsNB-dgp chooses the same focus regressors as the DGP of the simulation, which is unlikely in empirical applications.
The CV log scores in (ref) and (ref) show that all procedures perform similarly in terms of the distributional fit. Moreover, the curves decrease as I increase the number of training observations and flatten at about 2000 observations. This shows that the methods are able to `learn' more (i.e.\ improve the fit in terms of log score), when more training observations are available but stop `learning' at some point, i.e.\ when the curves flatten.
The other metrics for distributional fit, Brier and spherical score, show qualitatively similar results but the curves are flatter, hence the results are only shown in the supplementary materials. This further underlines that the distributional fit of the methods does not improve drastically when the dataset becomes larger.
Note that WALS NB specifications are computationally less demanding than lasso, while performing similarly in terms of CV RMSE and log score. They do not require any tuning unlike lasso, which performs an `internal' 10-fold CV to choose the optimal regularization parameter. Consequently, the fitting times of WALS NB are typically shorter than those of lasso and competitive with the ML specifications. Of course, one may change the parameters of the fitting algorithm of lasso to improve the computing time. However, it should not result in better performance metrics as the current setup already favors lasso: It allows many iterations in the fitting algorithm and, thus, a high chance for convergence.
In conclusion, all WALS NB specifications, except for walsNB-main-focus, perform better than the ML specifications in terms of RMSE, while metrics for the distributional fit such as log, Brier and spherical score are similar or minimally worse than the ML specifications. Moreover, the RMSE is similar or slightly lower than for lasso at large $t_{l}$, while demanding less computational resources as WALS NB does not require tuning by CV.
This paper extends the WALS approach to NB2 regression models (WALS NB) for count data based on WALS GLM of deluca2018glm and compares its predictive performance to the traditional ML and lasso estimator in simulated and real count datasets using the classical measure RMSE and strictly proper scoring rules.
In the simulation experiment, WALS NB outperforms the ML estimator in very sparse situations, i.e.\ where the number of auxiliary regressors is large and the number of training observations is small. When increasing the number of training observations, WALS NB and the unrestricted ML estimator converge in all performance metrics. Interestingly, whether WALS NB includes all regressors as auxiliary regressors or parts of them as focus does not change the results substantially. This shows that specifying all regressors as auxiliary is a reasonable choice if no prior information is available on the importance of the individual regressors. Moreover, it highlights that the regularized Bayesian estimation of the coefficients of the auxiliary regressors is key for the performance of WALS NB.
The empirical illustration emphasizes the results found in the simulation experiment: For small training sets, WALS NB using all covariates as auxiliary regressors outperforms all ML specifications in terms of RMSE while yielding a comparable distributional fit measured by strictly proper scores. Only the lasso estimator yields lower RMSE for small training sets. However, it is more computationally demanding than WALS NB due to the additional 10-fold CV that is run for determining the optimal regularization parameter. This highlights an important advantage of WALS compared to other model averaging techniques: low computational costs. Moreover, WALS NB using all covariates as auxiliary regressors outperformed all other specifications of WALS NB. Thus, if only the predictive power is of concern, WALS NB is a viable alternative to established estimation methods for the NB2 regression model that is easy to specify (choose all regressors as auxiliary), regularized and computationally efficient.
For future research, it would be interesting to generalize WALS to hurdle or zero-inflation models to handle count data with excess zeros. Thus far, WALS has been limited to univariate response variables, therefore extending the methodology to multivariate outcomes would allow a larger variety of applications, such as joint modeling of related count processes. Lastly, an investigation of the large sample properties of WALS could improve our understanding of statistical inference after model averaging.
\addcontentsline{toc}{section}{Acknowledgements}
The scientific computing center sciCORE (\url{https://scicore.unibas.ch/}) at the University of Basel provided me with valuable computing resources. I would also like to thank Christian Kleiber for our helpful discussions.
\stoplist[main]{lof} \stoplist[main]{lot}
\startlist[appendix]{lof} \startlist[appendix]{lot} \addappheadtotoc