Variance function

Derivation

The generalized linear model (GLM), is a generalization of ordinary regression analysis that extends to any member of the exponential family. It is particularly useful when the response variable is categorical, binary or subject to a constraint (e.g. only positive responses make sense). A quick summary of the components of a GLM are summarized on this page, but for more details and information see the page on generalized linear models.

A GLM consists of three main ingredients:

1. Random Component: a distribution of y from the exponential family, ${\displaystyle E[y\mid X]=\mu }$

2. Linear predictor: ${\displaystyle \eta =XB=\sum _{j=1}^{p}X_{ij}^{T}B_{j}}$

3. Link function: ${\displaystyle \eta =g(\mu ),\mu =g^{-1}(\eta )}$

First it is important to derive a couple key properties of the exponential family.

Any random variable ${\displaystyle {\textit {y}}}$ in the exponential family has a probability density function of the form,

${\displaystyle f(y,\theta ,\phi )=\exp \left({\frac {y\theta -b(\theta )}{\phi }}-c(y,\phi )\right)}$

with loglikelihood,

${\displaystyle \ell (\theta ,y,\phi )=\log(f(y,\theta ,\phi ))={\frac {y\theta -b(\theta )}{\phi }}-c(y,\phi )}$

Here, ${\displaystyle \theta }$ is the canonical parameter and the parameter of interest, and ${\displaystyle \phi }$ is a nuisance parameter which plays a role in the variance. We use the Bartlett's Identities to derive a general expression for the variance function. The first and second Bartlett results ensures that under suitable conditions (see Leibniz integral rule), for a density function dependent on ${\displaystyle \theta ,f_{\theta }()}$,

${\displaystyle \operatorname {E} _{\theta }\left[{\frac {\partial }{\partial \theta }}\log(f_{\theta }(y))\right]=0}$

${\displaystyle \operatorname {Var} _{\theta }\left[{\frac {\partial }{\partial \theta }}\log(f_{\theta }(y))\right]+\operatorname {E} _{\theta }\left[{\frac {\partial ^{2}}{\partial \theta ^{2}}}\log(f_{\theta }(y))\right]=0}$

These identities lead to simple calculations of the expected value and variance of any random variable ${\displaystyle {\textit {y}}}$ in the exponential family ${\displaystyle E_{\theta }[y],Var_{\theta }[y]}$.

Expected value of Y: Taking the first derivative with respect to ${\displaystyle \theta }$ of the log of the density in the exponential family form described above, we have

${\displaystyle {\frac {\partial }{\partial \theta }}\log(f(y,\theta ,\phi ))={\frac {\partial }{\partial \theta }}\left[{\frac {y\theta -b(\theta )}{\phi }}-c(y,\phi )\right]={\frac {y-b'(\theta )}{\phi }}}$

Then taking the expected value and setting it equal to zero leads to,

${\displaystyle \operatorname {E} _{\theta }\left[{\frac {y-b'(\theta )}{\phi }}\right]={\frac {\operatorname {E} _{\theta }[y]-b'(\theta )}{\phi }}=0}$

${\displaystyle \operatorname {E} _{\theta }[y]=b'(\theta )}$

Variance of Y: To compute the variance we use the second Bartlett identity,

${\displaystyle \operatorname {Var} _{\theta }\left[{\frac {\partial }{\partial \theta }}\left({\frac {y\theta -b(\theta )}{\phi }}-c(y,\phi )\right)\right]+\operatorname {E} _{\theta }\left[{\frac {\partial ^{2}}{\partial \theta ^{2}}}\left({\frac {y\theta -b(\theta )}{\phi }}-c(y,\phi )\right)\right]=0}$

${\displaystyle \operatorname {Var} _{\theta }\left[{\frac {y-b'(\theta )}{\phi }}\right]+\operatorname {E} _{\theta }\left[{\frac {-b''(\theta )}{\phi }}\right]=0}$

${\displaystyle \operatorname {Var} _{\theta }\left[y\right]=b''(\theta )\phi }$

We have now a relationship between ${\displaystyle \mu }$ and ${\displaystyle \theta }$, namely

${\displaystyle \mu =b'(\theta )}$ and ${\displaystyle \theta =b'^{-1}(\mu )}$, which allows for a relationship between ${\displaystyle \mu }$ and the variance,

${\displaystyle V(\theta )=b''(\theta )={\text{the part of the variance that depends on }}\theta }$

${\displaystyle \operatorname {V} (\mu )=b''(b'^{-1}(\mu )).\,}$

Note that because ${\displaystyle \operatorname {Var} _{\theta }\left[y\right]>0,b''(\theta )>0}$, then ${\displaystyle b':\theta \rightarrow \mu }$ is invertible. We derive the variance function for a few common distributions.

Example – normal

The normal distribution is a special case where the variance function is a constant. Let ${\displaystyle y\sim N(\mu ,\sigma ^{2})}$ then we put the density function of y in the form of the exponential family described above:

${\displaystyle f(y)=\exp \left({\frac {y\mu -{\frac {\mu ^{2}}{2}}}{\sigma ^{2}}}-{\frac {y^{2}}{2\sigma ^{2}}}-{\frac {1}{2}}\ln {2\pi \sigma ^{2}}\right)}$

where

${\displaystyle \theta =\mu ,}$

${\displaystyle b(\theta )={\frac {\mu ^{2}}{2}},}$

${\displaystyle \phi =\sigma ^{2},}$

${\displaystyle c(y,\phi )=-{\frac {y^{2}}{2\sigma ^{2}}}-{\frac {1}{2}}\ln {2\pi \sigma ^{2}}}$

To calculate the variance function ${\displaystyle V(\mu )}$, we first express ${\displaystyle \theta }$ as a function of ${\displaystyle \mu }$. Then we transform ${\displaystyle V(\theta )}$ into a function of ${\displaystyle \mu }$

${\displaystyle \theta =\mu }$

${\displaystyle b'(\theta )=\theta =\operatorname {E} [y]=\mu }$

${\displaystyle V(\theta )=b''(\theta )=1}$

Therefore, the variance function is constant.

Example – Bernoulli

Let ${\displaystyle y\sim {\text{Bernoulli}}(p)}$, then we express the density of the Bernoulli distribution in exponential family form,

${\displaystyle f(y)=\exp \left(y\ln {\frac {p}{1-p}}+\ln(1-p)\right)}$

${\displaystyle \theta =\ln {\frac {p}{1-p}}=}$ logit(p), which gives us ${\displaystyle p={\frac {e^{\theta }}{1+e^{\theta }}}=}$ expit${\displaystyle (\theta )}$

${\displaystyle b(\theta )=\ln(1+e^{\theta })}$ and

${\displaystyle b'(\theta )={\frac {e^{\theta }}{1+e^{\theta }}}=}$ expit${\displaystyle (\theta )=p=\mu }$

${\displaystyle b''(\theta )={\frac {e^{\theta }}{1+e^{\theta }}}-\left({\frac {e^{\theta }}{1+e^{\theta }}}\right)^{2}}$

This give us

${\displaystyle V(\mu )=\mu (1-\mu )}$

Example – Poisson

Let ${\displaystyle y\sim {\text{Poisson}}(\lambda )}$, then we express the density of the Poisson distribution in exponential family form,

${\displaystyle f(y)=\exp(y\ln \lambda -\ln \lambda )}$

${\displaystyle \theta =\ln \lambda =}$ which gives us ${\displaystyle \lambda =e^{\theta }}$

${\displaystyle b(\theta )=e^{\theta }}$ and

${\displaystyle b'(\theta )=e^{\theta }=\lambda =\mu }$

${\displaystyle b''(\theta )=e^{\theta }=\mu }$

This give us

${\displaystyle V(\mu )=\mu }$

Here we see the central property of Poisson data, that the variance is equal to the mean.

Example – Gamma

The Gamma distribution and density function can be expressed under different parametrizations. We will use the form of the gamma with parameters ${\displaystyle (\mu ,\nu )}$

${\displaystyle f_{\mu ,\nu }(y)={\frac {1}{\Gamma (\nu )y}}\left({\frac {\nu y}{\mu }}\right)^{\nu }e^{-{\frac {\nu y}{\mu }}}}$

Then in exponential family form we have

${\displaystyle f_{\mu ,\nu }(y)=\exp \left({\frac {-{\frac {1}{\mu }}y+\ln({\frac {1}{\mu }})}{\frac {1}{\nu }}}+\ln \left({\frac {\nu ^{\nu }y^{\nu -1}}{\Gamma (\nu )}}\right)\right)}$

${\displaystyle \theta ={\frac {-1}{\mu }}\rightarrow \mu ={\frac {-1}{\theta }}}$

${\displaystyle \phi ={\frac {1}{\nu }}}$

${\displaystyle b(\theta )=-\ln(-\theta )}$

${\displaystyle b'(\theta )={\frac {-1}{\theta }}={\frac {-1}{\frac {-1}{\mu }}}=\mu }$

${\displaystyle b''(\theta )={\frac {1}{\theta ^{2}}}=\mu ^{2}}$

And we have ${\displaystyle V(\mu )=\mu ^{2}}$

Application – weighted least squares

A very important application of the variance function is its use in parameter estimation and inference when the response variable is of the required exponential family form as well as in some cases when it is not (which we will discuss in quasi-likelihood). Weighted least squares (WLS) is a special case of generalized least squares. Each term in the WLS criterion includes a weight that determines that the influence each observation has on the final parameter estimates. As in regular least squares, the goal is to estimate the unknown parameters in the regression function by finding values for parameter estimates that minimize the sum of the squared deviations between the observed responses and the functional portion of the model.

While WLS assumes independence of observations it does not assume equal variance and is therefore a solution for parameter estimation in the presence of heteroscedasticity. The Gauss–Markov theorem and Aitken demonstrate that the best linear unbiased estimator (BLUE), the unbiased estimator with minimum variance, has each weight equal to the reciprocal of the variance of the measurement.

In the GLM framework, our goal is to estimate parameters ${\displaystyle \beta }$, where ${\displaystyle Z=g(E[y\mid X])=X\beta }$. Therefore, we would like to minimize ${\displaystyle (Z-XB)^{T}W(Z-XB)}$ and if we define the weight matrix W as

${\displaystyle \underbrace {W} _{n\times n}={\begin{bmatrix}{\frac {1}{\phi V(\mu _{1})g'(\mu _{1})^{2}}}&0&\cdots &0&0\\0&{\frac {1}{\phi V(\mu _{2})g'(\mu _{2})^{2}}}&0&\cdots &0\\\vdots &\vdots &\vdots &\vdots &0\\\vdots &\vdots &\vdots &\vdots &0\\0&\cdots &\cdots &0&{\frac {1}{\phi V(\mu _{n})g'(\mu _{n})^{2}}}\end{bmatrix}},}$

where ${\displaystyle \phi ,V(\mu ),g(\mu )}$ are defined in the previous section, it allows for iteratively reweighted least squares (IRLS) estimation of the parameters. See the section on iteratively reweighted least squares for more derivation and information.

Also, important to note is that when the weight matrix is of the form described here, minimizing the expression ${\displaystyle (Z-XB)^{T}W(Z-XB)}$ also minimizes the Pearson distance. See Distance correlation for more.

The matrix W falls right out of the estimating equations for estimation of ${\displaystyle \beta }$. Maximum likelihood estimation for each parameter ${\displaystyle \beta _{r},1\leq r\leq p}$, requires

${\displaystyle \sum _{i=1}^{n}{\frac {\partial l_{i}}{\partial \beta _{r}}}=0}$, where ${\displaystyle \operatorname {l} (\theta ,y,\phi )=\log(\operatorname {f} (y,\theta ,\phi ))={\frac {y\theta -b(\theta )}{\phi }}-c(y,\phi )}$ is the log-likelihood.

Looking at a single observation we have,

${\displaystyle {\frac {\partial l}{\partial \beta _{r}}}={\frac {\partial l}{\partial \theta }}{\frac {\partial \theta }{\partial \mu }}{\frac {\partial \mu }{\partial \eta }}{\frac {\partial \eta }{\partial \beta _{r}}}}$

${\displaystyle {\frac {\partial \eta }{\partial \beta _{r}}}=x_{r}}$

${\displaystyle {\frac {\partial l}{\partial \theta }}={\frac {y-b'(\theta )}{\phi }}={\frac {y-\mu }{\phi }}}$

${\displaystyle {\frac {\partial \theta }{\partial \mu }}={\frac {\partial b'^{-1}(\mu )}{\mu }}={\frac {1}{b''(b'(\mu ))}}={\frac {1}{V(\mu )}}}$

This gives us

${\displaystyle {\frac {\partial l}{\partial \beta _{r}}}={\frac {y-\mu }{\phi V(\mu )}}{\frac {\partial \mu }{\partial \eta }}x_{r}}$, and noting that

${\displaystyle {\frac {\partial \eta }{\partial \mu }}=g'(\mu )}$ we have that

${\displaystyle {\frac {\partial l}{\partial \beta _{r}}}=(y-\mu )W{\frac {\partial \eta }{\partial \mu }}x_{r}}$

The Hessian matrix is determined in a similar manner and can be shown to be,

${\displaystyle H=X^{T}(y-\mu )\left[{\frac {\partial }{\beta _{s}}}W{\frac {\partial }{\beta _{r}}}\right]-X^{T}WX}$

Noticing that the Fisher Information (FI),

${\displaystyle {\text{FI}}=-E[H]=X^{T}WX}$, allows for asymptotic approximation of ${\displaystyle {\hat {\beta }}}$

${\displaystyle {\hat {\beta }}\sim N_{p}(\beta ,(X^{T}WX)^{-1})}$, and hence inference can be performed.

Application – quasi-likelihood

Because most features of GLMs only depend on the first two moments of the distribution, rather than the entire distribution, the quasi-likelihood can be developed by just specifying a link function and a variance function. That is, we need to specify

• the link function, ${\displaystyle E[y]=\mu =g^{-1}(\eta )}$
• the variance function, ${\displaystyle V(\mu )}$, where the ${\displaystyle \operatorname {Var} _{\theta }(y)=\sigma ^{2}V(\mu )}$

With a specified variance function and link function we can develop, as alternatives to the log-likelihood function, the score function, and the Fisher information, a quasi-likelihood, a quasi-score, and the quasi-information. This allows for full inference of ${\displaystyle \beta }$.

Quasi-likelihood (QL)

Though called a quasi-likelihood, this is in fact a quasi-log-likelihood. The QL for one observation is

${\displaystyle Q_{i}(\mu _{i},y_{i})=\int _{y_{i}}^{\mu _{i}}{\frac {y_{i}-t}{\sigma ^{2}V(t)}}\,dt}$

And therefore the QL for all n observations is

${\displaystyle Q(\mu ,y)=\sum _{i=1}^{n}Q_{i}(\mu _{i},y_{i})=\sum _{i=1}^{n}\int _{y_{i}}^{\mu _{i}}{\frac {y-t}{\sigma ^{2}V(t)}}\,dt}$

From the QL we have the quasi-score

Quasi-score (QS)

Recall the score function, U, for data with log-likelihood ${\displaystyle \operatorname {l} (\mu \mid y)}$ is

${\displaystyle U={\frac {\partial l}{d\mu }}.}$

We obtain the quasi-score in an identical manner,

${\displaystyle U={\frac {y-\mu }{\sigma ^{2}V(\mu )}}}$

Noting that, for one observation the score is

${\displaystyle {\frac {\partial Q}{\partial \mu }}={\frac {y-\mu }{\sigma ^{2}V(\mu )}}}$

The first two Bartlett equations are satisfied for the quasi-score, namely

${\displaystyle E[U]=0}$

and

${\displaystyle \operatorname {Cov} (U)+E\left[{\frac {\partial U}{\partial \mu }}\right]=0.}$

In addition, the quasi-score is linear in y.

Ultimately the goal is to find information about the parameters of interest ${\displaystyle \beta }$. Both the QS and the QL are actually functions of ${\displaystyle \beta }$. Recall, ${\displaystyle \mu =g^{-1}(\eta )}$, and ${\displaystyle \eta =X\beta }$, therefore,

${\displaystyle \mu =g^{-1}(X\beta ).}$

Quasi-information (QI)

The quasi-information, is similar to the Fisher information,

${\displaystyle i_{b}=-\operatorname {E} \left[{\frac {\partial U}{\partial \beta }}\right]}$

QL, QS, QI as functions of ${\displaystyle \beta }$

The QL, QS and QI all provide the building blocks for inference about the parameters of interest and therefore it is important to express the QL, QS and QI all as functions of ${\displaystyle \beta }$.

Recalling again that ${\displaystyle \mu =g^{-1}(X\beta )}$, we derive the expressions for QL, QS and QI parametrized under ${\displaystyle \beta }$.

Quasi-likelihood in ${\displaystyle \beta }$,

${\displaystyle Q(\beta ,y)=\int _{y}^{\mu (\beta )}{\frac {y-t}{\sigma ^{2}V(t)}}\,dt}$

The QS as a function of ${\displaystyle \beta }$ is therefore

${\displaystyle U_{j}(\beta _{j})={\frac {\partial }{\partial \beta _{j}}}Q(\beta ,y)=\sum _{i=1}^{n}{\frac {\partial \mu _{i}}{\partial \beta _{j}}}{\frac {y_{i}-\mu _{i}(\beta _{j})}{\sigma ^{2}V(\mu _{i})}}}$

${\displaystyle U(\beta )={\begin{bmatrix}U_{1}(\beta )\\U_{2}(\beta )\\\vdots \\\vdots \\U_{p}(\beta )\end{bmatrix}}=D^{T}V^{-1}{\frac {(y-\mu )}{\sigma ^{2}}}}$

Where,

${\displaystyle \underbrace {D} _{n\times p}={\begin{bmatrix}{\frac {\partial \mu _{1}}{\partial \beta _{1}}}&\cdots &\cdots &{\frac {\partial \mu _{1}}{\partial \beta _{p}}}\\{\frac {\partial \mu _{2}}{\partial \beta _{1}}}&\cdots &\cdots &{\frac {\partial \mu _{2}}{\partial \beta _{p}}}\\\vdots \\\vdots \\{\frac {\partial \mu _{m}}{\partial \beta _{1}}}&\cdots &\cdots &{\frac {\partial \mu _{m}}{\partial \beta _{p}}}\end{bmatrix}}\underbrace {V} _{n\times n}=\operatorname {diag} (V(\mu _{1}),V(\mu _{2}),\ldots ,\ldots ,V(\mu _{n}))}$

The quasi-information matrix in ${\displaystyle \beta }$ is,

${\displaystyle i_{b}=-{\frac {\partial U}{\partial \beta }}=\operatorname {Cov} (U(\beta ))={\frac {D^{T}V^{-1}D}{\sigma ^{2}}}}$

Obtaining the score function and the information of ${\displaystyle \beta }$ allows for parameter estimation and inference in a similar manner as described in Application – weighted least squares.