In statistics, the class of vector generalized linear models (VGLMs) was proposed to enlarge the scope of models catered for by generalized linear models (GLMs). In particular, VGLMs allow for response variables outside the classical exponential family and for more than one parameter. Each parameter (not necessarily a mean) can be transformed by a link function. The VGLM framework is also large enough to naturally accommodate multiple responses; these are several independent responses each coming from a particular statistical distribution with possibly different parameter values.
Vector generalized linear models are described in detail in Yee (2015). The central algorithm adopted is the iteratively reweighted least squares method, for maximum likelihood estimation of usually all the model parameters. In particular, Fisher scoring is implemented by such, which, for most models, uses the first and expected second derivatives of the log-likelihood function.
Motivation
GLMs essentially cover one-parameter models from the classical exponential family, and include 3 of the most important statistical regression models: the linear model, Poisson regression for counts, and logistic regression for binary responses. However, the exponential family is far too limiting for regular data analysis. For example, for counts, zero-inflation, zero-truncation and overdispersion are regularly encountered, and the makeshift adaptations made to the binomial and Poisson models in the form of quasi-binomial and quasi-Poisson can be argued as being ad hoc and unsatisfactory. But the VGLM framework readily handles models such as zero-inflated Poisson regression, zero-altered Poisson (hurdle) regression, positive-Poisson regression, and negative binomial regression. As another example, for the linear model, the variance of a normal distribution is relegated as a scale parameter and it is treated often as a nuisance parameter (if it is considered as a parameter at all). But the VGLM framework allows the variance to be modelled using covariates.
As a whole, one can loosely think of VGLMs as GLMs that handle many models outside the classical exponential family and are not restricted to estimating a single mean. During estimation, rather than using weighted least squares during IRLS, one uses generalized least squares to handle the correlation between the M linear predictors.
Data and notation
We suppose that the response or outcome or the dependent variable(s), , are assumed to be generated from a particular distribution. Most distributions are univariate, so that , and an example of is the bivariate normal distribution.
Sometimes we write our data as for . Each of the n observations are considered to be independent. Then . The are known positive prior weights, and often .
The explanatory or independent variables are written , or when i is needed, as . Usually there is an intercept, in which case or .
Actually, the VGLM framework allows for S responses, each of dimension . In the above S = 1. Hence the dimension of is more generally . One handles S responses by code such as vglm(cbind(y1, y2, y3) ~ x2 + x3, ..., data = mydata) for S = 3. To simplify things, most of this article has S = 1.
Model components
The VGLM usually consists of four elements: 1. A probability density function or probability mass function from some statistical distribution which has a log-likelihood , first derivatives and expected information matrix that can be computed. The model is required to satisfy the usual MLE regularity conditions. 2. Linear predictors described below to model each parameter , 3. Link functions such that 4. Constraint matrices for each of full column-rank and known.
Linear predictors
Each linear predictor is a quantity which incorporates information about the independent variables into the model. The symbol (Greek "eta") denotes a linear predictor and a subscript j is used to denote the jth one. It relates the jth parameter to the explanatory variables, and is expressed as linear combinations (thus, "linear") of unknown parameters i.e., of regression coefficients .
The jth parameter, , of the distribution depends on the independent variables, through
Let be the vector of all the linear predictors. (For convenience we always let be of dimension M). Thus all the covariates comprising potentially affect all the parameters through the linear predictors . Later, we will allow the linear predictors to be generalized to additive predictors, which is the sum of smooth functions of each and each function is estimated from the data.
Link functions
Each link function provides the relationship between a linear predictor and a parameter of the distribution. There are many commonly used link functions, and their choice can be somewhat arbitrary. It makes sense to try to match the domain of the link function to the range of the distribution's parameter value. Notice above that the allows a different link function for each parameter. They have similar properties as with generalized linear models, for example, common link functions include the logit link for parameters in , and the log link for positive parameters. The VGAM package has function identitylink() for parameters that can assume both positive and negative values.
Constraint matrices
More generally, the VGLM framework allows for any linear constraints between the regression coefficients of each linear predictors. For example, we may want to set some to be equal to 0, or constraint some of them to be equal. We have
where the are the constraint matrices. Each constraint matrix is known and prespecified, and has M rows, and between 1 and M columns. The elements of constraint matrices are finite-valued, and often they are just 0 or 1. For example, the value 0 effectively omits that element while a 1 includes it. It is common for some models to have a parallelism assumption, which means that for , and for some models, for too. The special case when for all is known as trivial constraints; all the regression coefficients are estimated and are unrelated. And is known as an intercept-only parameter if the jth row of all the are equal to for , i.e., equals an intercept only. Intercept-only parameters are thus modelled as simply as possible, as a scalar.
The unknown parameters, , are typically estimated by the method of maximum likelihood. All the regression coefficients may be put into a matrix as follows:
The xij facility
With even more generally, one can allow the value of a variable to have a different value for each . For example, if each linear predictor is for a different time point then one might have a time-varying covariate. For example, in discrete choice models, one has conditional logit models, nested logit models, generalized logit models, and the like, to distinguish between certain variants and fit a multinomial logit model to, e.g., transport choices. A variable such as cost differs depending on the choice, for example, taxi is more expensive than bus, which is more expensive than walking. The xij facility of VGAM allows one to generalize to .
The most general formula is
Here the is an optional offset; which translates to be a matrix in practice. The VGAM package has an xij argument that allows the successive elements of the diagonal matrix to be inputted.
Software
Yee (2015) describes an R package implementation in the called VGAM. Currently this software fits approximately 150 models/distributions. The central modelling functions are vglm() and vgam(). The family argument is assigned a VGAM family function, e.g., family = negbinomial for negative binomial regression, family = poissonff for Poisson regression, family = propodds for the proportional odd model or cumulative logit model for ordinal categorical regression.
Fitting
Maximum likelihood
We are maximizing a log-likelihood
where the are positive and known prior weights. The maximum likelihood estimates can be found using an iteratively reweighted least squares algorithm using Fisher's scoring method, with updates of the form:
where is the Fisher information matrix at iteration a. It is also called the expected information matrix, or EIM.
VLM For the computation, the (small) model matrix constructed from the RHS of the formula in vglm() and the constraint matrices are combined to form a big model matrix. The IRLS is applied to this big X. This matrix is known as the VLM matrix, since the vector linear model is the underlying least squares problem being solved. A VLM is a weighted multivariate regression where the variance-covariance matrix for each row of the response matrix is not necessarily the same, and is known. (In classical multivariate regression, all the errors have the same variance-covariance matrix, and it is unknown). In particular, the VLM minimizes the weighted sum of squares
This quantity is minimized at each IRLS iteration. The working responses (also known as pseudo-response and adjusted dependent vectors) are
where the are known as working weights or working weight matrices. They are symmetric and positive-definite. Using the EIM helps ensure that they are all positive-definite (and not just the sum of them) over much of the parameter space. In contrast, using Newton–Raphson would mean the observed information matrices would be used, and these tend to be positive-definite in a smaller subset of the parameter space.
Computationally, the Cholesky decomposition is used to invert the working weight matrices and to convert the overall generalized least squares problem into an ordinary least squares problem.
Examples
Generalized linear models
Of course, all generalized linear models are a special cases of VGLMs. But we often estimate all parameters by full maximum likelihood estimation rather than using the method of moments for the scale parameter.
Ordered categorical response If the response variable is an ordinal measurement with M + 1 levels, then one may fit a model function of the form:
where
for Different links g lead to proportional odds models or ordered probit models, e.g., the VGAM family function cumulative(link = probit) assigns a probit link to the cumulative probabilities, therefore this model is also called the cumulative probit model. In general they are called cumulative link models.
For categorical and multinomial distributions, the fitted values are an (M + 1)-vector of probabilities, with the property that all probabilities add up to 1. Each probability indicates the likelihood of occurrence of one of the M + 1 possible values.
Unordered categorical response
If the response variable is a nominal measurement, or the data do not satisfy the assumptions of an ordered model, then one may fit a model of the following form:
for The above link is sometimes called the multilogit link, and the model is called the multinomial logit model. It is common to choose the first or the last level of the response as the reference or baseline group; the above uses the last level. The VGAM family function multinomial() fits the above model, and it has an argument called refLevel that can be assigned the level used for as the reference group.
Count data
Classical GLM theory performs Poisson regression for count data. The link is typically the logarithm, which is known as the canonical link. The variance function is proportional to the mean:
where the dispersion parameter is typically fixed at exactly one. When it is not, the resulting quasi-likelihood model is often described as Poisson with overdispersion, or quasi-Poisson; then is commonly estimated by the method-of-moments and as such, confidence intervals for are difficult to obtain.