Vector generalized linear model

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.

We suppose that the response or outcome or the dependent variable(s), $\backslash boldsymbol\{y\}=(y\_1,\backslash ldots,y\_\{Q\_1\})^T$, are assumed to be generated from a particular distribution. Most distributions are univariate, so that $Q\_1=1$, and an example of $Q\_1=2$ is the bivariate normal distribution.

Sometimes we write our data as $(\backslash boldsymbol\{x\}\_\{i\},w\_i,\backslash boldsymbol\{y\}\_\{i\})$

for $i=1,\backslash ldots,n$. Each of the n observations are considered to be

independent.

Then $\backslash boldsymbol\{y\}\_i\; =\; (y\_\{i1\},\backslash ldots,y\_\{iQ\_1\})^T$.

The $w\_\{i\}$ are known positive prior weights, and often $w\_\{i\}=1$.

The explanatory or independent variables are written $\backslash boldsymbol\{x\}=(x\_1,\backslash ldots,x\_p)^T$,

or when i is needed, as $\backslash boldsymbol\{x\}\_i\; =\; (x\_\{i1\},\backslash ldots,x\_\{ip\})^T$.

Usually there is an intercept, in which case

$x\_1\; =\; 1$ or $x\_\{i1\}\; =\; 1$.

Actually, the VGLM framework allows for S responses, each of dimension $Q\_1$.

In the above S = 1. Hence the dimension of $\backslash boldsymbol\{y\}\_\{i\}$ is more generally $Q\; =\; S\; \backslash times\; Q\_1$. 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.

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 $\backslash ell$, first derivatives $\backslash partial\; \backslash ell\; /\; \backslash partial\; \backslash theta\_j$ and expected information matrix that can be computed. The model is required to satisfy the usual MLE regularity conditions.

: 2. Linear predictors $\backslash eta\_j$ described below to model each parameter $\backslash theta\_j$, $j=1,\backslash ldots,M.$

: 3. Link functions $g\_j$ such that $\backslash theta\_j\; =\; g\_j^\{-1\}(\backslash eta\_j).$

: 4. Constraint matrices $\backslash boldsymbol\{H\}\_k$ for $k=1,\backslash ldots,p,$ each of full column-rank and known.

Each linear predictor is a quantity which incorporates

information about the independent variables into the model.

The symbol $\backslash eta\_j$ (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

$\backslash eta\_j$ is expressed as linear combinations (thus, "linear")

of unknown parameters $\backslash boldsymbol\{\backslash beta\}\_j,$

i.e., of regression coefficients $\backslash beta\_\{(j)k\}$.

The jth parameter, $\backslash theta\_j$, of the distribution depends on the

independent variables, $\backslash boldsymbol\{x\},$ through

: $$

g_j( \theta_j) = \eta_j = \boldsymbol{\beta}_j^T \boldsymbol{x}.

Let $\backslash boldsymbol\{\backslash eta\}\; =\; (\backslash eta\_1,\; \backslash ldots,\; \backslash eta\_M)^T$ be the vector of

all the linear predictors. (For convenience we always let $\backslash boldsymbol\{\backslash eta\}$

be of dimension M).

Thus all the covariates comprising $\backslash boldsymbol\{x\}$ potentially affect all the parameters through the linear predictors $\backslash eta\_j$. Later, we will allow the linear predictors to be generalized to additive predictors, which is the sum of smooth functions of each $x\_k$ and each function is estimated from the data.

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 $g\_j$ 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 $(0,\; 1)$,

and the log link for positive parameters. The VGAM package has function identitylink() for parameters that can assume both positive and negative values.

More generally, the VGLM framework allows for any linear constraints between the regression coefficients $\backslash beta\_\{(j)k\}$ 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

: $\backslash boldsymbol\{\backslash eta\}\; =$

\sum_{k=1}^p \, \boldsymbol{\beta}^{}_{(k)} \, x_{k} =

\sum_{k=1}^p \, \boldsymbol{H}_k \; \boldsymbol{\beta}^{*}_{(k)} \, x_{k}

where the $\backslash boldsymbol\{H\}\_k$ 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

$\backslash boldsymbol\{H\}\_k\; =\; \backslash boldsymbol\{1\}\_M$ for $k=2,\backslash ldots,p$, and

for some models, for $k=1$ too.

The special case when $\backslash boldsymbol\{H\}\_k\; =\; \backslash boldsymbol\{I\}\_M$ for

all $k=1,\backslash ldots,p$ is known as trivial constraints; all the

regression coefficients are estimated and are unrelated.

And $\backslash theta\_j$ is known as an intercept-only parameter

if the jth row of all the $\backslash boldsymbol\{H\}\_k\; =$ are equal to $\backslash boldsymbol\{0\}^T$ for $k=2,\backslash ldots,p$, i.e., $\backslash eta\_j\; =\; \backslash beta\_\{(j)1\}^\{*\}$ equals an intercept only. Intercept-only parameters are thus modelled as simply as possible, as a scalar.

The unknown parameters, $\backslash boldsymbol\{\backslash beta\}^\{*\}\; =\; (\backslash boldsymbol\{\backslash beta\}\_\{(1)\}^\{*T\},\backslash ldots,\backslash boldsymbol\{\backslash beta\}\_\{(p)\}^\{*T\})^T$,

are typically estimated by the method of maximum likelihood.

All the regression coefficients may be put into a matrix as follows:

: $$

\boldsymbol{\eta}_i =

\boldsymbol{B}^T \boldsymbol{x}_i =

\begin{pmatrix}

\boldsymbol{\beta}^{T}_{1} \, \boldsymbol{x}_i \\

\vdots \\

\boldsymbol{\beta}^{T}_{M} \, \boldsymbol{x}_i \\

\end{pmatrix} =

\left( \boldsymbol{\beta}^{}_{(1)}, \ldots, \boldsymbol{\beta}^{}_{(p)} \right) \; \boldsymbol{x}_i.

With even more generally, one can allow the value of a variable $x\_k$

to have a different value for each $\backslash eta\_j$.

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 $\backslash eta\_j(\backslash boldsymbol\{x\}\_i)$

to $\backslash eta\_j(\backslash boldsymbol\{x\}\_\{ij\})$.

The most general formula is

: $\backslash boldsymbol\{\backslash eta\}\_i\; =$

\boldsymbol{o}_i +

\sum_{k=1}^{p} \, diag(x_{ik1},\ldots,x_{ikM}) \,

\mathbf{H}_k \, \boldsymbol{\beta}_{(k)}^{*}.

Here the $\backslash boldsymbol\{o\}\_i$ is an optional offset; which translates

to be a $n\; \backslash times\; M$ matrix in practice.

The VGAM package has an xij argument that allows

the successive elements of the diagonal matrix to be inputted.

Yee (2015) describes an R package

implementation in the

called VGAM.{{cite web |url=http://cran.R-project.org/package=VGAM|

title=Vector Generalized Linear Models |

date=2016-01-18}}

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.

We are maximizing a log-likelihood

:$$

\ell =

\sum_{i=1}^{n} \, w_i \, \ell_i,

where the $w\_i$ 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:

:$$

\boldsymbol\beta^{(a+1)} =

\boldsymbol\beta^{(a)} +

\boldsymbol{\mathcal{I}}^{-1}(\boldsymbol\beta^{(a)}) \, \,

\mathbf{u}(\boldsymbol\beta^{(a)}),

where $\backslash boldsymbol\{\backslash mathcal\{I\}\}(\backslash boldsymbol\backslash beta^\{(a)\})$ is

the Fisher information matrix at iteration a.

It is also called the expected information matrix, or EIM.

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

:$$

\mathrm{ResSS} =

\sum_{i=1}^{n} \;

w_i

\left\{ \mathbf{z}^{(a-1)}_{i} - \boldsymbol{\eta}^{(a-1)}_{i} \right\}^T

\mathbf{W}_{i}^{(a-1)}

\left\{ \mathbf{z}^{(a-1)}_{i} - \boldsymbol{\eta}^{(a-1)}_{i} \right\}

This quantity is minimized at each IRLS iteration.

The working responses (also known as pseudo-response and adjusted dependent vectors) are

:$$

\mathbf{z}_i = \boldsymbol{\eta}_i + \mathbf{W}_i^{-1} \mathbf{u}_i,

where the $\backslash mathbf\{W\}\_i$ 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.

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.

If the response variable is an ordinal measurement with M + 1 levels, then one may fit a model function of the form:

:$g(\backslash theta\_j)\; =\; \backslash eta\_j$   where $\backslash theta\_j\; =\; \backslash mathrm\{Pr\}(Y\; \backslash leq\; j),$

for $j=1,\backslash ldots,M.$

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.

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:

:$$

\log \left[ \frac{Pr(Y = j)}{\mathrm{Pr}(Y = M+1)} \right] = \eta_j,

for $j=1,\backslash ldots,M.$ 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.

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:

:$\backslash operatorname\{Var\}(Y\_\{i\})\; =\; \backslash tau\backslash mu\_\{i\},\backslash ,$

where the dispersion parameter $\backslash tau$ 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 $\backslash tau$ is commonly estimated by the method-of-moments and as such,

confidence intervals for $\backslash tau$ are difficult to obtain.

In contrast, VGLMs offer a much richer set of models to handle overdispersion with respect to the Poisson, e.g., the negative binomial distribution and several variants thereof. Another count regression model is the generalized Poisson distribution. Other possible models are the zeta distribution and the Zipf distribution.

RR-VGLMs are VGLMs where a subset of

the B matrix is of a lower rank.

Without loss of generality, suppose that $\backslash boldsymbol\{x\}=(\backslash boldsymbol\{x\}\_1^T,\; \backslash boldsymbol\{x\}\_2^T)^T$ is a partition of the covariate vector. Then the part of the B matrix corresponding to $\backslash boldsymbol\{x\}\_2$ is of the form $\backslash boldsymbol\{A\}\; \backslash boldsymbol\{C\}^T$ where $\backslash boldsymbol\{A\}$ and

$\backslash boldsymbol\{C\}$ are thin matrices (i.e., with R columns), e.g., vectors if the rank R = 1. RR-VGLMs potentially offer several advantages when applied to certain

models and data sets. Firstly, if M and p are large then the number of regression coefficients

that are estimated by VGLMs is large ($M\; \backslash times\; p$). Then RR-VGLMs can reduce the number of estimated regression coefficients enormously if R is low, e.g., R = 1

or R = 2. An example of a model where this is particularly useful is the RR-multinomial logit model, also known as the stereotype model.

Secondly, $\backslash boldsymbol\{\backslash nu\}\; =\; \backslash boldsymbol\{C\}^T\; \backslash boldsymbol\{x\}\_2\; =\; (\backslash nu\_1,\backslash ldots,\backslash nu\_R)^T$

is an R-vector of latent variables, and often these can be usefully interpreted.

If R = 1 then we can write $\backslash nu\; =\; \backslash boldsymbol\{c\}^T\; \backslash boldsymbol\{x\}\_2$

so that the latent variable comprises loadings on the explanatory variables.

It may be seen that RR-VGLMs take optimal linear combinations of the $\backslash boldsymbol\{x\}\_2$

and then a VGLM is fitted to the explanatory variables $(\backslash boldsymbol\{x\}\_1,\; \backslash boldsymbol\{\backslash nu\})$. Thirdly, a biplot can be produced if R = 2 , and this allows the model to be visualized.

It can be shown that RR-VGLMs are simply VGLMs where the constraint matrices for

the variables in $\backslash boldsymbol\{x\}\_2$ are unknown and to be estimated.

It then transpires that $\backslash boldsymbol\{H\}\_k\; =\; \backslash boldsymbol\{A\}$ for

such variables.

RR-VGLMs can be estimated by an alternating algorithm which fixes $\backslash boldsymbol\{A\}$

and estimates $\backslash boldsymbol\{C\},$ and then fixes $\backslash boldsymbol\{C\}$ and estimates $\backslash boldsymbol\{A\}$, etc.

In practice, some uniqueness constraints are needed for $\backslash boldsymbol\{A\}$

and/or $\backslash boldsymbol\{C\}$. In VGAM, the rrvglm() function uses corner constraints by default, which means that the top R rows of $\backslash boldsymbol\{A\}$ is set to $\backslash boldsymbol\{I\}\_R$. RR-VGLMs were proposed in 2003.{{cite journal |last1= Yee|first1=T. W. |last2= Hastie|first2=T. J.|

year=2003 | title=Reduced-rank vector generalized linear models |

journal=Statistical Modelling |volume=3 |issue=1 |pages=15–41 |doi=10.1191/1471082x03st045oa|citeseerx=10.1.1.36.3700 |s2cid=122810408 }}

A special case of RR-VGLMs is when R = 1 and M = 2. This is dimension reduction from 2 parameters to 1 parameter. Then it can be shown that

: $$

\theta_2 = g_2^{-1}

\left( t_1 + a_{21} \cdot g_1(\theta_1) \right),

where elements $t\_1$ and $a\_\{21\}$ are estimated. Equivalently,

: $$

\eta_2 = t_1 + a_{21} \cdot \eta_1.

This formula provides a coupling of $\backslash eta\_1$ and $\backslash eta\_2$. It induces a relationship between two parameters of a model that can be useful, e.g., for modelling a mean-variance relationship. Sometimes there is some choice of link functions, therefore it offers a little flexibility when coupling the two parameters, e.g., a logit, probit, cauchit or cloglog link for parameters in the unit interval. The above formula is particularly useful for the negative binomial distribution, so that the RR-NB has variance function

: $$

\operatorname{Var}(Y\mid\boldsymbol{x}) = \mu(\boldsymbol{x}) + \delta_1 \, \mu(\boldsymbol{x})^{\delta_2}.

This has been called the NB-P variant by some authors. The $\backslash delta\_1$ and $\backslash delta\_2$ are estimated, and it is also possible to obtain approximate confidence intervals for them too.

Incidentally, several other useful NB variants can also be fitted, with the help of selecting the right combination of constraint matrices. For example, NB − 1, NB − 2 (negbinomial() default), NB − H; see Yee (2014){{cite journal |last= Yee|first=T. W. |year=1996 |

title=Reduced-rank vector generalized linear models with two linear predictors |

journal=Computational Statistics & Data Analysis |

volume=71 | pages=889–902|doi=10.1016/j.csda.2013.01.012 }} and Table 11.3 of Yee (2015).

The subclass of row-column interaction models

(RCIMs) has also been proposed; these are a special type of RR-VGLM.

RCIMs apply only to a matrix Y response and there are

no explicit explanatory variables $\backslash boldsymbol\{x\}$.

Instead, indicator variables for each row and column are explicitly set up, and an order-R

interaction of the form $\backslash boldsymbol\{A\}\; \backslash boldsymbol\{C\}^T$ is allowed.

Special cases of this type of model include the Goodman RC association model

and the quasi-variances methodology as implemented by the qvcalc R package.

RCIMs can be defined as a RR-VGLM applied to Y with

:$$

g_1(\theta_1) \equiv

\eta_{1ij} = \beta_0 + \alpha_i + \gamma_j + \sum_{r=1}^R c_{ir} \, a_{jr} .

For the Goodman RC association model, we have $\backslash eta\_\{1ij\}=\backslash log\; \backslash mu\_\{ij\},$ so that

if R = 0 then it is a Poisson regression fitted to a matrix of counts with row effects and column effects; this has a similar idea to a no-interaction two-way ANOVA model.

Another example of a RCIM is if $g\_1$ is the identity link and the parameter is the median and the model corresponds to an asymmetric Laplace distribution; then a no-interaction RCIM is similar to a technique called median polish.

In VGAM, rcim() and grc() functions fit the above models.

And also Yee and Hadi (2014)

{{cite journal |last1= Yee|first1=T. W. |last2= Hadi|first2=A. F. |year=2014 |title=Row-column interaction models, with an R implementation |journal=Computational Statistics |volume=29 |issue=6 |pages=1427–1445 |doi=10.1007/s00180-014-0499-9|s2cid=253724333 }}

show that RCIMs can be used to fit unconstrained quadratic ordination

models to species data; this is an example of indirect gradient analysis in

ordination (a topic in statistical ecology).

Vector generalized additive models (VGAMs) are a major

extension to VGLMs in which the linear predictor $\backslash eta\_j$ is not restricted to be

linear in the covariates $x\_k$ but is the

sum of smoothing functions applied to the $x\_k$:

:$\backslash boldsymbol\{\backslash eta\}(\backslash boldsymbol\{x\})\; =$

\boldsymbol{H}_1 \, \boldsymbol{\beta}_{(1)}^{*} +

\boldsymbol{H}_2 \, \boldsymbol{f}_{(2)}^{*}(x_2) +

\boldsymbol{H}_3 \, \boldsymbol{f}_{(3)}^{*}(x_3) + \cdots \,\!

where $\backslash boldsymbol\{f\}\_\{(k)\}^\{*\}(x\_k)\; =\; (f\_\{(1)k\}^\{*\}(x\_k),f\_\{(2)k\}^\{*\}(x\_k),\backslash ldots)^T.$

These are M additive predictors.

Each smooth function $f\_\{(j)k\}^\{*\}$ is estimated from the data.

Thus VGLMs are model-driven while VGAMs are data-driven.

Currently, only smoothing splines are implemented in the VGAM package.

For M > 1 they are actually vector splines, which estimate the component functions

in $f\_\{(j)k\}^\{*\}(x\_k)$ simultaneously.

Of course, one could use regression splines with VGLMs.

The motivation behind VGAMs is similar to

that of

Hastie and Tibshirani (1990)

{{cite book|last1 = Hastie |first1= T. J. |last2=Tibshirani | first2=R. J.|

title = Generalized Additive Models|

publisher = London:Chapman and Hall|

year = 1990}}

and

Wood (2017).

{{cite book|last = Wood |first= S. N. |

title = Generalized Additive Models: An Introduction with R|

publisher = London: Chapman and Hall|

year = 2017|

edition = second|

isbn=9781498728331}}

VGAMs were proposed in 1996

.

{{cite journal |last1= Yee|first1=T. W. |last2= Wild|first2=C. J. |year=1996 |title=Vector generalized additive models |journal=Journal of the Royal Statistical Society, Series B |volume=58 |issue=3 |pages=481–493|doi=10.1111/j.2517-6161.1996.tb02095.x }}

Currently, work is being done to estimate VGAMs using P-splines

of Eilers and Marx (1996)

.

{{cite journal |last1= Eilers|first1=P. H. C. |last2= Marx|first2=B. D. |year=1996 |

title=Flexible smoothing with B-splines and penalties |

journal=Statistical Science |volume=11 |issue=2 |pages=89–121 |doi=10.1214/ss/1038425655|citeseerx=10.1.1.47.4521}}

This allows for several advantages over using smoothing splines and vector backfitting, such as the

ability to perform automatic smoothing parameter selection easier.

These add on a quadratic in the latent variable to the RR-VGLM class.

The result is a bell-shaped curve can be fitted to each response, as

a function of the latent variable.

For R = 2, one has bell-shaped surfaces as a function of the 2

latent variables---somewhat similar to a

bivariate normal distribution.

Particular applications of QRR-VGLMs can be found in ecology,

in a field of multivariate analysis called ordination.

As a specific rank-1 example of a QRR-VGLM,

consider Poisson data with S species.

The model for Species s is the Poisson regression

:$$

\log\,\mu_s(\nu) =

\eta_s(\nu) =

\beta_{(s)1} + \beta_{(s)2} \, \nu + \beta_{(s)3} \, \nu^2 =

\alpha_s - \frac{1}{2} \left( \frac{\nu -u_s}{t_s} \right)^2,

for $s=1,\backslash ldots,S$. The right-most parameterization which uses the symbols $\backslash alpha\_s,$ $u\_s,$ $t\_s,$ has particular ecological meaning, because they relate to the species abundance, optimum and tolerance respectively. For example, the tolerance is a measure of niche width, and a large value means that that species can live in a wide range of environments. In the above equation, one would need in order

to obtain a bell-shaped curve.

QRR-VGLMs fit Gaussian ordination models by maximum likelihood estimation, and

they are an example of direct gradient analysis.

The cqo() function in the VGAM package currently

calls optim() to search for the optimal

$\backslash boldsymbol\{C\}$, and given that, it is easy to calculate

the site scores and fit a suitable generalized linear model to that.

The function is named after the acronym CQO, which stands for

constrained quadratic ordination: the constrained is for direct

gradient analysis (there are environmental variables, and a linear combination

of these is taken as the latent variable) and the quadratic is for the

quadratic form in the latent variables $\backslash boldsymbol\{\backslash nu\}$

on the $\backslash boldsymbol\{\backslash eta\}$ scale.

Unfortunately QRR-VGLMs are sensitive to outliers in both the response

and explanatory variables, as well as being computationally expensive, and

may give a local solution rather than a global solution.

QRR-VGLMs were proposed in 2004.

{{cite journal |last= Yee|first=T. W. |year=2004 |

title=A new technique for maximum-likelihood canonical Gaussian ordination |

journal=Ecological Monographs |volume=74 |issue=4 |pages=685–701 |doi=10.1890/03-0078}}

• generalized linear models
• R (software)
• Regression analysis
• Statistical model
• Natural exponential family

{{reflist}}

• {{cite book | author=Hilbe, Joseph | title = Negative Binomial Regression | publisher = Cambridge: Cambridge University Press | year = 2011 |edition=2nd |isbn=978-0-521-19815-8}}

{{statistics|correlation}}

Category:Actuarial science

*

Category:Regression models