In statistics, the variance function is a smooth function that depicts the variance of a random quantity as a function of its mean. The variance function is a measure of heteroscedasticity and plays a large role in many settings of statistical modelling. It is a main ingredient in the generalized linear model framework and a tool used in non-parametric regression,[1] semiparametric regression[1] and functional data analysis.[2] In parametric modeling, variance functions take on a parametric form and explicitly describe the relationship between the variance and the mean of a random quantity. In a non-parametric setting, the variance function is assumed to be a smooth function.
In a regression model setting, the goal is to establish whether or not a relationship exists between a response variable and a set of predictor variables. Further, if a relationship does exist, the goal is then to be able to describe this relationship as best as possible. A main assumption in linear regression is constant variance or (homoscedasticity), meaning that different response variables have the same variance in their errors, at every predictor level. This assumption works well when the response variable and the predictor variable are jointly normal. As we will see later, the variance function in the Normal setting is constant; however, we must find a way to quantify heteroscedasticity (non-constant variance) in the absence of joint Normality.
When it is likely that the response follows a distribution that is a member of the exponential family, a generalized linear model may be more appropriate to use, and moreover, when we wish not to force a parametric model onto our data, a non-parametric regression approach can be useful. The importance of being able to model the variance as a function of the mean lies in improved inference (in a parametric setting), and estimation of the regression function in general, for any setting.
Variance functions play a very important role in parameter estimation and inference. In general, maximum likelihood estimation requires that a likelihood function be defined. This requirement then implies that one must first specify the distribution of the response variables observed. However, to define a quasi-likelihood, one need only specify a relationship between the mean and the variance of the observations to then be able to use the quasi-likelihood function for estimation.[3] Quasi-likelihood estimation is particularly useful when there is overdispersion. Overdispersion occurs when there is more variability in the data than there should otherwise be expected according to the assumed distribution of the data.
In summary, to ensure efficient inference of the regression parameters and the regression function, the heteroscedasticity must be accounted for. Variance functions quantify the relationship between the variance and the mean of the observed data and hence play a significant role in regression estimation and inference.
The variance function and its applications come up in many areas of statistical analysis. A very important use of this function is in the framework of generalized linear models and non-parametric regression.
When a member of the exponential family has been specified, the variance function can easily be derived.[4] The general form of the variance function is presented under the exponential family context, as well as specific forms for Normal, Bernoulli, Poisson, and Gamma. In addition, we describe the applications and use of variance functions in maximum likelihood estimation and quasi-likelihood estimation.
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,
E[y\midX]=\mu
2. Linear predictor:
η=XB=
p | |
\sum | |
j=1 |
T | |
X | |
ij |
Bj
3. Link function:
η=g(\mu),\mu=g-1(η)
First it is important to derive a couple key properties of the exponential family.
Any random variable
it{y}
f(y,\theta,\phi)=\exp\left(
y\theta-b(\theta) | |
\phi |
-c(y,\phi)\right)
\ell(\theta,y,\phi)=log(f(y,\theta,\phi))=
y\theta-b(\theta) | |
\phi |
-c(y,\phi)
\theta
\phi
\theta,f\theta
\operatorname{E} | ||||
|
log(f\theta(y))\right]=0
\operatorname{Var} | ||||
|
log(f\theta(y))\right]+\operatorname{E}
|
log(f\theta(y))\right]=0
These identities lead to simple calculations of the expected value and variance of any random variable
it{y}
E\theta[y],Var\theta[y]
Expected value of Y:Taking the first derivative with respect to
\theta
\partial | |
\partial\theta |
log(f(y,\theta,\phi))=
\partial | \left[ | |
\partial\theta |
y\theta-b(\theta) | |
\phi |
-c(y,\phi)\right]=
y-b'(\theta) | |
\phi |
\operatorname{E} | ||||
|
\right]=
\operatorname{E | |
\theta[y]-b'(\theta)}{\phi}=0 |
\operatorname{E}\theta[y]=b'(\theta)
Variance of Y:To compute the variance we use the second Bartlett identity,
\operatorname{Var} | \left( | ||||
|
y\theta-b(\theta) | |
\phi |
-
c(y,\phi)\right)\right]+\operatorname{E} | \left( | ||||
|
y\theta-b(\theta) | |
\phi |
-c(y,\phi)\right)\right]=0
\operatorname{Var} | ||||
|
\right]+\operatorname{E}\theta\left[
-b''(\theta) | |
\phi |
\right]=0
\operatorname{Var}\theta\left[y\right]=b''(\theta)\phi
We have now a relationship between
\mu
\theta
\mu=b'(\theta)
\theta=b'-1(\mu)
\mu
V(\theta)=b''(\theta)=thepartofthevariancethatdependson\theta
\operatorname{V}(\mu)=b''(b'-1(\mu)).
\operatorname{Var}\theta\left[y\right]>0,b''(\theta)>0
b':\theta → \mu
The normal distribution is a special case where the variance function is a constant. Let
y\simN(\mu,\sigma2)
f(y)=\exp\left(
| ||||||
\sigma2 |
-
y2 | |
2\sigma2 |
-
1 | |
2 |
ln{2\pi\sigma2}\right)
where
\theta=\mu,
b(\theta)=
\mu2 | |
2 |
,
\phi=\sigma2,
c(y,\phi)=-
y2 | |
2\sigma2 |
-
1 | |
2 |
ln{2\pi\sigma2}
To calculate the variance function
V(\mu)
\theta
\mu
V(\theta)
\mu
\theta=\mu
b'(\theta)=\theta=\operatorname{E}[y]=\mu
V(\theta)=b''(\theta)=1
Therefore, the variance function is constant.
Let
y\simBernoulli(p)
f(y)=\exp\left(yln
p | |
1-p |
+ln(1-p)\right)
\theta=ln
p | |
1-p |
=
p=
e\theta | |
1+e\theta |
=
(\theta)
b(\theta)=ln(1+e\theta)
b'(\theta)=
e\theta | |
1+e\theta |
=
(\theta)=p=\mu
b''(\theta)=
e\theta | |
1+e\theta |
-\left(
e\theta | |
1+e\theta |
\right)2
This give us
V(\mu)=\mu(1-\mu)
Let
y\simPoisson(λ)
f(y)=\exp(ylnλ-lnλ)
\theta=lnλ=
λ=e\theta
b(\theta)=e\theta
b'(\theta)=e\theta=λ=\mu
b''(\theta)=e\theta=\mu
This give us
V(\mu)=\mu
Here we see the central property of Poisson data, that the variance is equal to the mean.
The Gamma distribution and density function can be expressed under different parametrizations. We will use the form of the gamma with parameters
(\mu,\nu)
f\mu,\nu(y)=
1 | \left( | |
\Gamma(\nu)y |
\nuy | |
\mu |
\right)\nu
| ||||
e |
f\mu,\nu(y)=\exp\left(
| ||||||||
|
+ln\left(
\nu\nuy\nu-1 | |
\Gamma(\nu) |
\right)\right)
\theta=
-1 | |
\mu |
→ \mu=
-1 | |
\theta |
\phi=
1 | |
\nu |
b(\theta)=-ln(-\theta)
b'(\theta)=
-1 | |
\theta |
=
-1 | |||
|
=\mu
b''(\theta)=
1 | |
\theta2 |
=\mu2
And we have
V(\mu)=\mu2
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
\beta
Z=g(E[y\midX])=X\beta
(Z-XB)TW(Z-XB)
\underbrace{W}n=\begin{bmatrix}
1 | ||||||||
|
&0& … &0&0\\ 0&
1 | ||||||||
|
&0& … &0\\ \vdots&\vdots&\vdots&\vdots&0\\ \vdots&\vdots&\vdots&\vdots&0\\ 0& … & … &0&
1 | ||||||||
|
\end{bmatrix},
\phi,V(\mu),g(\mu)
Also, important to note is that when the weight matrix is of the form described here, minimizing the expression
(Z-XB)TW(Z-XB)
The matrix W falls right out of the estimating equations for estimation of
\beta
\betar,1\leqr\leqp
n | |
\sum | |
i=1 |
\partialli | |
\partial\betar |
=0
\operatorname{l}(\theta,y,\phi)=log(\operatorname{f}(y,\theta,\phi))=
y\theta-b(\theta) | |
\phi |
-c(y,\phi)
\partiall | |
\partial\betar |
=
\partiall | |
\partial\theta |
\partial\theta | |
\partial\mu |
\partial\mu | |
\partialη |
\partialη | |
\partial\betar |
\partialη | |
\partial\betar |
=xr
\partiall | |
\partial\theta |
=
y-b'(\theta) | = | |
\phi |
y-\mu | |
\phi |
\partial\theta | |
\partial\mu |
=
\partialb'-1(\mu) | |
\mu |
=
1 | |
b''(b'(\mu)) |
=
1 | |
V(\mu) |
This gives us
\partiall | = | |
\partial\betar |
y-\mu | |
\phiV(\mu) |
\partial\mu | |
\partialη |
xr
\partialη | |
\partial\mu |
=g'(\mu)
\partiall | |
\partial\betar |
=(y-\mu)W
\partialη | |
\partial\mu |
xr
The Hessian matrix is determined in a similar manner and can be shown to be,
H=
| |||||
X | W |
\partial | |
\betar |
\right]-XTWX
Noticing that the Fisher Information (FI),
FI=-E[H]=XTWX
\hat{\beta}
\hat{\beta}\sim
TWX) | |
N | |
p(\beta,(X |
-1)
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
E[y]=\mu=g-1(η)
V(\mu)
\operatorname{Var}\theta(y)=\sigma2V(\mu)
\beta
Quasi-likelihood (QL)
Though called a quasi-likelihood, this is in fact a quasi-log-likelihood. The QL for one observation is
Qi(\mui,yi)=
\mui | |
\int | |
yi |
yi-t | |
\sigma2V(t) |
dt
Q(\mu,y)=
n | |
\sum | |
i=1 |
Qi(\mui,yi)=
n | |
\sum | |
i=1 |
\mui | |
\int | |
yi |
y-t | |
\sigma2V(t) |
dt
From the QL we have the quasi-score
Quasi-score (QS)
Recall the score function, U, for data with log-likelihood
\operatorname{l}(\mu\midy)
U=
\partiall | |
d\mu |
.
We obtain the quasi-score in an identical manner,
U=
y-\mu | |
\sigma2V(\mu) |
Noting that, for one observation the score is
\partialQ | |
\partial\mu |
=
y-\mu | |
\sigma2V(\mu) |
The first two Bartlett equations are satisfied for the quasi-score, namely
E[U]=0
\operatorname{Cov}(U)+E\left[
\partialU | |
\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
\beta
\beta
\mu=g-1(η)
η=X\beta
\mu=g-1(X\beta).
Quasi-information (QI)
The quasi-information, is similar to the Fisher information,
ib=-\operatorname{E}\left[
\partialU | |
\partial\beta |
\right]
QL, QS, QI as functions of
\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
\beta
Recalling again that
\mu=g-1(X\beta)
\beta
Quasi-likelihood in
\beta
Q(\beta,y)=
\mu(\beta) | |
\int | |
y |
y-t | |
\sigma2V(t) |
dt
The QS as a function of
\beta
Uj(\betaj)=
\partial | |
\partial\betaj |
Q(\beta,y)=
n | |
\sum | |
i=1 |
\partial\mui | |
\partial\betaj |
yi-\mui(\betaj) | ||||||
|
U(\beta)=\begin{bmatrix}U1(\beta)\\ U2(\beta)\\ \vdots\\ \vdots\\ Up(\beta) \end{bmatrix}=DTV-1
(y-\mu) | |
\sigma2 |
Where,
\underbrace{D}n=\begin{bmatrix}
\partial\mu1 | & … & … & | |
\partial\beta1 |
\partial\mu1 | \\ | |
\partial\betap |
\partial\mu2 | & … & … & | |
\partial\beta1 |
\partial\mu2 | \\ \vdots\\ \vdots\\ | |
\partial\betap |
\partial\mum | & … & … & | |
\partial\beta1 |
\partial\mum | |
\partial\betap |
\end{bmatrix} \underbrace{V}n x =\operatorname{diag}(V(\mu1),V(\mu2),\ldots,\ldots,V(\mun))
The quasi-information matrix in
\beta
ib=-
\partialU | |
\partial\beta |
=\operatorname{Cov}(U(\beta))=
DTV-1D | |
\sigma2 |
Obtaining the score function and the information of
\beta
Non-parametric estimation of the variance function and its importance, has been discussed widely in the literature[5] [6] [7] In non-parametric regression analysis, the goal is to express the expected value of your response variable(y) as a function of your predictors (X). That is we are looking to estimate a mean function,
g(x)=\operatorname{E}[y\midX=x]
g(x)
gv(x)=\operatorname{Var}(Y\midX=x)
gv(x)=\operatorname{Var}(Y\midX=x)=\operatorname{E}[y2\midX=x]-\left[\operatorname{E}[y\midX=x]\right]2
An example is detailed in the pictures to the right. The goal of the project was to determine (among other things) whether or not the predictor, number of years in the major leagues (baseball), had an effect on the response, salary, a player made. An initial scatter plot of the data indicates that there is heteroscedasticity in the data as the variance is not constant at each level of the predictor. Because we can visually detect the non-constant variance, it useful now to plot
gv(x)=\operatorname{Var}(Y\midX=x)=\operatorname{E}[y2\midX=x]-\left[\operatorname{E}[y\midX=x]\right]2
\operatorname{E}[y2\midX=x]
\left[\operatorname{E}[y\midX=x]\right]2
. Peter McCullagh. Nelder, John . John Nelder . Generalized Linear Models . London: Chapman and Hall . 1989 . second. 0-412-31760-5 .