View source for Random coefficient models


{| class="wikitable"
!Copyright notice <!-- don't remove! --> 
|-
| This article ''Linear mixed models (=Random Coefficient Models)'' was adapted from an original article by Nicholas Tibor Longford, which appeared in ''StatProb: The Encyclopedia Sponsored by Statistics and Probability Societies''. The original article ([<nowiki>http://statprob.com/encyclopedia/RandomCoefficientModels.html</nowiki> StatProb Source], Local Files: [[Media:random_coefficient_models.pdf|pdf]] | [[Media:random_coefficient_models.tex|tex]]) is copyrighted by the author(s), the article has been donated to ''Encyclopedia of Mathematics'', and its further issues are under ''Creative Commons Attribution Share-Alike License'. All pages from StatProb are contained in the [[:Category:Statprob|Category StatProb]].
|-
|}

{{MSC|62J10}}


<!-- \documentclass[10pt]{article} 
     \usepackage{epsf} -->

$ \newcommand{\refb}{\ref} $

$ \newcommand{\bbeta}{ { {\mathbf \beta}}} $
$ \newcommand{\bdelta}{ { {\mathbf \delta}}} $
$ \newcommand{\bomega}{ { {\mathbf \omega}}} $
$ \newcommand{\bOmega}{ { {\mathbf \Omega}}} $
$ \newcommand{\bSigma}{ { {\mathbf \Sigma}}} $
$ \newcommand{\tr}{\mbox{tr}} $
$ \newcommand{\tra}{^{\top}} $
$ \newcommand{\pri}{^{\prime}} $
$ \newcommand{\mE}{\mbox {E}} $
$ \newcommand{\bG}{ {\bf G}} $
$ \newcommand{\bH}{ {\bf H}} $
$ \newcommand{\bI}{ {\bf I}} $
$ \newcommand{\bUU}{ {\bf U}} $
$ \newcommand{\bV}{ {\bf V}} $
$ \newcommand{\bW}{ {\bf W}} $
$ \newcommand{\bX}{ {\bf X}} $
$ \newcommand{\bZ}{ {\bf Z}} $
$ \newcommand{\bY}{ {\bf Y}} $
$ \newcommand{\be}{ {\bf e}} $
$ \newcommand{\bi}{ {\bf i}} $
$ \newcommand{\bq}{ {\bf q}} $
$ \newcommand{\bs}{ {\bf s}} $
$ \newcommand{\buu}{ {\bf u}} $
$ \newcommand{\bx}{ {\bf x}} $
$ \newcommand{\by}{ {\bf y}} $
$ \newcommand{\bz}{ {\bf z}} $
$ \newcommand{\bzer}{ {\bf 0}} $
$ \newcommand{\var}{\mbox {var}} $
$ \newcommand{\cor}{\mbox {cor}} $

<!-- \begin{document} -->

<center>'''Random Coefficient Models\footnote {N. T. Longford, SNTL and Departament d'Economia i Empresa, Universitat Pompeu Fabra, Ramon Trias Fargas 25--27, 08005 Barcelona, Spain; email:  NTL@sntl.co.uk'''</center>
}

<center> Nicholas T. Longford 
\\ SNTL and UPF, Barcelona, Spain </center>

<!-- \maketitle -->

=== Summary ===
Random coefficient models are intended for settings with two or more sources
of random variation. The widest range of applications is found for them
when observational units form natural clusters, such that the units within
a cluster are more similar than units in general.
Models for independent observations have to be extended to allow
for within- and between-cluster variation.
\vspace*{1cm}

<!-- \noindent -->

Keywords: ''Analysis of variance; borrowing strength; clusters;
correlation structure; empirical Bayes analysis;
longitudinal analysis; maximum likelihood; ordinary regression.''

\clearpage

Independence of the observations is a key assumption of many standard
statistical methods, such as analysis of variance (ANOVA) and ordinary
regression, and some of its extensions.
Common examples of data structures that do not fit into such a framework
arise in longitudinal analysis, in which observations are made on subjects
at subject-specific sequences of time points, and in studies that involve
subjects (units) occurring naturally in clusters, such as individuals within
families, school\-children within classrooms, employees within companies,
and the like. The assumption of independence of the observational units
is not tenable when observations within a cluster tend to be more similar
than observations in general. Such similarity
can be conveniently represented by a positive correlation (dependence).

This article describes an adaptation of the ordinary regression
for clustered observations. Such observations require two indices,
one for elements within clusters, $i=1, \ldots, n_{j\,}$,
and another for clusters, $j=1, \ldots, m$.
Thus, we have $n = n_1 + \cdots + n_m$ elementary units and $m$ clusters.
The ordinary regression model
$$ \label{ML1}
y_{ij} \,=\, \bx_{ij\,}\bbeta + \varepsilon_{ij} \,,
\tag{1}$$

with the usual assumptions of normality, independence and equal variance
(homoscedasticity) of the deviations $\varepsilon_{ij\,}$,
$\varepsilon_{ij} \sim {\cal{N}}(0, \sigma^2)$, i.i.d.,
implies that the regressions within the clusters $j$
have a common vector of coefficients $\bbeta$.
This restriction can be relaxed by allowing the regressions to differ
in their intercepts. A practical way of defining such a model is
by the equation
$$ \label{ML2}
y_{ij} \,=\, \bx_{ij\,} \bbeta + \delta_j + \varepsilon_{ij} \,,
\tag{2}$$

where $\delta_{j\,}$, $j=1, \ldots, m$,
are a random sample from a centred normal distribution,
$\delta_j \sim {\cal{N}}(0, \sigma^2_{\rm B})$, i.i.d., independent
from the $\varepsilon$'s.
This differs from the model for analysis of covariance (ANCOVA) only
by the status of the deviations $\delta_{j\,}$. In ANCOVA, they are fixed
(constant across hypothetical replications), whereas in \ref{ML2}
they are random.

In the model in \ref{ML2}, the within-cluster regressions are parallel
--- their intercepts are $\beta_0 + \delta_{j\,}$,
but the coefficients on all the other variables in $\bx$ are common
to the clusters.
A more appealing interpretation of the model is that observations
in a cluster are correlated,
$$
\cor \left ( y_{i_1, j\,}, y_{i_2, j} \right ) \,=\,
\frac{\sigma^2_{\rm B}} {\sigma^2 + \sigma^2_{\rm B}} \,,
$$
because they share the same deviation $\delta_{j\,}$.
Further relaxation of how the within-cluster regressions differ
is attained by allowing some (or all) the regression slopes to be
specific to the clusters.
We select a set of variables in $\bx$, denoted by $\bz$, and assume
that the regressions with respect to these variables differ across
the clusters, but are constant with respect to the remaining variables;
$$ \label{ML3}
y_{ij} \,=\, \bx_{ij\,} \bbeta + \bz_{ij\,}\bdelta_j + \varepsilon_{ij} \,,
\tag{3}$$

where $\bdelta_{j\,}$, $j=1, \ldots, m$, are a random sample from a multivariate normal distribution
${\cal{N}}(\bzer, \bSigma_{\rm B})$, independent from the $\varepsilon$'s.
We say that the variables in $\bz$ are associated with (cluster-level)
variation. The variance of an observation $y_{ij\,}$, without conditioning
on the cluster $j$, is
$$
\var \left ( y_{ij} \right ) \,=\, \sigma^2 +
\bx_{ij} \bSigma_{\rm B\,} \bx_{ij} \tra \,.
$$
We refer to $\sigma^2$ and $\bz_{ij} \bSigma_{\rm B\,} \bz_{ij}\tra$
as the ''variance components'' (at the elementary and cluster
levels, respectively). The principle of invariance
with respect to linear transformations of $\bz$ implies that
the intercept should always be included in $\bz$,
unless $\bz$ is empty, as in the model in \ref{ML1}.
The function $V(\bz) = \bz \bSigma_{\rm B} \bz\tra$,
over the feasible values of $\bz$, defines the ''pattern of variation'',
and it can be described by its behaviour
(local minima, points of inflection, and the like).
By way of an example, suppose $\bz$ contains the intercept and
a single variable $z$. Denote the variances in $\bSigma_{\rm B}$
by $\sigma^2_0$ and $\sigma^2_{\rm z\,}$, and the covariance by
$\sigma_{0\rm z\,}$. Then
$$ \label{ML4}
V(\bz) \,=\, \sigma^2_0 + 2z \sigma_{0\rm z} + z^2 \sigma^2_{\rm z} \,,
\tag{4}$$

and this quadratic function has a unique minimum at
$z^{\ast} = -\sigma_{0\rm z}/\sigma_{\rm z\,}^2$, unless $\sigma^2_{\rm z} =0$,
in which case we revert to the model in \ref{ML2} in which $V(\bz)$
is constant. Figure 1 illustrates four patterns of variation
on examples with a single covariate.


<span id="Fig1">
[[File:FigMiLo.png|center|frame|Figure 1:
Patterns of variation in random coefficient models with a single
covariate. Thick dashes mark the average regression ($\bx\bbeta$)
and thin solid lines the within-cluster regressions
($\bx\bbeta + \bz\bdelta_j$).

In the megaphone pattern, the matrix $\bSigma_{\rm B}$ is singular,
with $\sigma^2_{\rm z} > 0$.([[Media:FigMiLo.ps|ps]]) ]]
</span>


The model in \ref{ML3} is fitted by maximum likelihood (ML)
which maximizes the log-likelihood function1
$$ \label{Fisc1}
l \left (\bbeta, \sigma^2, \bSigma_{\rm B} \right ) \,=\,
-\frac{1}{2} \sum_{j=1}^m
\left [ \log \left \{ \det \left ( \bV_j \right ) \right \} \,+\,
\left ( \by_j - \bX_j \bbeta \right ) \tra \bV_j^{-1}
\left ( \by_j - \bX_j \bbeta \right ) \right ] \,,
\tag{5}$$

where $\bV_j$ is the variance matrix of the observations in cluster $j$,
$\by_j$ the vector of the outcomes for the observations in cluster $j$,
and $\bX_j$ the corresponding regression design matrix
formed by vertical stacking of the rows $\bx_{ij\,}$, $i=1, \ldots, n_{j\,}$.
The variation design matrices $\bZ_{j\,}$,
$j = 1, \ldots, m$, are defined similarly; with them,
$\bV_j = \sigma^2 \bI_{n_j} + \bZ_j \bSigma_{\rm B} \bZ_j \tra$,
where $\bI_{n_j}$ is the $n_j \times n_j$ identity matrix.
The Fisher scoring algorithm for maximising the log-likelihood
in \ref{Fisc1} is described in the Appendix;
for details and applications, see see Longford (1993),
and for an alternative method Goldstein (2000).
These and other algorithms are implemented
in most standard statistical packages.
A key to their effective implementation are closed-form expressions
for the inverse and determinant of patterned matrices (Harville, 1997).

Model selection entails two tasks,
selecting a set of variables to form $\bx$
and selecting its subset to form $\bz$.
The variables in $\bx$ can be defined for elements or clusters;
the latter can be defined as being constant within clusters.
Inclusion of cluster-level variables in $\bz$ does not have
an interpretation in terms of varying regression coefficients,
so associating them with variation is in most contexts not meaningful.
However, the identity in \ref{ML4} and its generalisations for
$\bSigma_{\rm B}$ with more than two rows and columns indicate
that $\bz$ can be used for modelling variance heterogeneity.
The likelihood ratio test statistic and various information criteria
can be used for selecting among alternative models,
so long as one is a submodel of the other; that is,
the variables in both $\bx$ and $\bz$ of one model are subsets
of (or coincide with) their counterparts in the other model.

Random coefficients can be applied to a range of models much wider
than ordinary regression.
In principle, we can conceive any ''basis model'',
characterized by a vector of parameters, which applies to every cluster.
A subset of these parameters is constant across the clusters
and the remainder varies according to a model for cluster-level variation.
The latter model need not be a multivariate normal distribution,
although suitable alternatives to it are difficult to identify.
The basis model itself can be complex, such as a random coefficient model
itself.
This gives rise to three- or, generally, ''multilevel models'',
in which elements are clustered within two-level units,
these units in three-level units, and so on.

Generalized linear mixed models (GLMM) have generalized linear models
(McCullagh and Nelder, 1989) as their basis; see Pinheiro and Bates (2000).
For cluster $j$ we posit the model
$$
g \left \{ \mE(\by_j\,|\,\bdelta_j) \right \} \,=\,
\bX_j \bbeta + \bZ_j \bdelta_j \,,
$$
with a (monotone) link function $g$,
which transforms the expected outcomes to a linear scale;
the outcomes (elements of $\by_j$) are conditionally independent
given $\bdelta_j$ and have a specified distribution,
such as binary or gamma.
Some advantage, in both interpretation and computing, is gained
by using ''canonical'' link functions, for which the conditional
likelihood has a compact set of sufficient statistics.
For example, the logistic link, $g(p) = \log(p)-\log(1-p)$,
is the canonical link for binary outcomes, and the logarithm
is the canonical link for Poisson outcomes. Models with the
(standard) normality assumptions correspond to the link $g(y)=y$.

Without conditioning, the likelihood for non-normally distributed
$\by_j$ has an intractable form.
An established approach to fitting GLMM maximises a (tractable)
normal-like approximation to the log-likelihood. This algorithm
can be described as an iteratively reweighted version of the
Fisher scoring (or another) algorithm for fitting
the (normal) linear mixed model.
With the advent of modern computing,
using numerical quadrature and other methods for numerical integration
has because feasible,
especially for fitting models with univariate deviations $\delta_{j\,}$.
An alternative framework for GLMM and a computational algorithm
for fitting them are constructed by Lee and Nelder (2001).

Random coefficient models are well suited for analysing surveys
in which clusters arise naturally as a consequence of the organisation
(design) of the survey and the way the studied population is structured.
They can be applied also in settings in which multiple observations
are made on subjects, as in longitudinal studies (Molenberghs and
Verbeke, 2000). In some settings it is contentious as to whether
the clusters should be regarded as fixed or random.
For example, small-area estimation (Rao, 2003) is concerned
with inferences about districts or another partition of a country
when some (or all) districts are represented in the analysed national survey
by small subsamples. In one perspective, district-level quantities,
such as their means of a variable, should be regarded as fixed
because they are the inferential targets,
fixed across hypothetical replications.
When they are assumed to be random the (random coefficient) models
are often more parsimonious than their fixed-effects (ANCOVA) counterparts,
because the number of parameters involved does not depend
on the number of clusters.

''Borrowing strength'' (Robbins, 1955,
Efron and Morris, 1972) is a general principle for efficient inference
about each cluster (district) by exploiting the similarity of the clusters.
It is the foundation of the empirical Bayes analysis, in which
the between-cluster variance matrix plays a role similar to the Bayes
prior for the within-cluster regression coefficients. The qualifier
`empirical' refers to using a data-based estimator $\widehat{\bSigma}_{\rm B}$
in place of the unknown $\bSigma_{\rm B\,}$.


====References====
{|
|-
|valign="top"|{{Ref|1}}||valign="top"|  Efron, B., and Morris, C. N. (1972). Limiting the risk of Bayes and empirical Bayes estimators --- Part II: empirical Bayes case. ''Journal of the American Statistical Association'' '''67''', 1286--1289.
|-
|valign="top"|{{Ref|2}}||valign="top"| Goldstein, H. (2000). ''Multilevel Statistial Models.'' 2nd Edition. London: Edward Arnold.
|-
|valign="top"|{{Ref|3}}||valign="top"| Harville, D. (1997). ''Matrix Algebra from a Statistician's Perspective.'' New York: Springer-Verlag.
|-
|valign="top"|{{Ref|4}}||valign="top"| Lee, Y., and Nelder, J. A. (2001). Hierarchical generalised linear models: a synthesis of generalised linear models, random effect models and structured dispersions. ''Biometrika'' '''88''', 987--1004.
|-
|valign="top"|{{Ref|5}}||valign="top"| Longford, N. T. (1993). ''Random Coefficient Models.'' Oxford: Oxford University Press.
|-
|valign="top"|{{Ref|6}}||valign="top"| Magnus, J. R., and Neudecker, H. (1988). ''Matrix Differential Calculus with Applications in Statistics and Econometrics.'' New York: Wiley.
|-
|valign="top"|{{Ref|7}}||valign="top"| McCullagh, P., and Nelder, J. A. (1989). ''Generalized Linear Models.'' 2nd Edition. London: Chapman and Hall.
|-
|valign="top"|{{Ref|8}}||valign="top"| Pinheiro, J. C., and Bates, D. M. (2000). ''Mixed-Effects Models in  S and  Splus.'' New York: Springer-Verlag.
|-
|valign="top"|{{Ref|9}}||valign="top"| Rao, J. N. K. (2003). ''Small Area Estimation.'' New York: Wiley and Sons.
|-
|valign="top"|{{Ref|10}}||valign="top"| Robbins, H. (1955). An empirical Bayes approach to statistics. ''Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability,'' '''1''', 157--164. Berkeley, CA: University of California Press.
|-
|valign="top"|{{Ref|11}}||valign="top"| Verbeke, G., and Molenberghs, G. (2000). ''Linear Mixed Models for Longitudinal Data.'' New York: Springer-Verlag.
|-
|}
*{1cm}

== Appendix. Fisher scoring algorithm ==
This Appendix describes a method for fitting a random coefficient
model by maximum likelihood.
We prefer to use the scaled variance matrices
$\bW_j = \sigma^{-2} \bV_j$ and $\bOmega = \sigma^{-2} \bSigma_{\rm B\,}$,
so that $\bW_j = \bI_{n_d} + \bZ_j \bOmega \bZ_j \tra$ does not depend
on $\sigma^2$.
The log-likelihood in \ref{Fisc1} is equal to
$$ \label{Fisc2}
l \left (\bbeta, \sigma^2, \bOmega \right ) \,=\,
-\frac{1}{2} \sum_{j=1}^m
\left [ n \log \left ( \sigma^2 \right ) +
\log \left \{ \det \left ( \bW_j \right ) \right \} \,+\,
\frac{1}{\sigma^2} \, \be_j \tra \bW_j^{-1} \be_j \right ] \,.
\tag{6}$$

where $\be_j = \by_j - \bX_j \bbeta$ is the vector of residuals
for cluster $j$.
We have the following closed-form expressions
for the inverse and determinant of $\bW_{j\,}$:
\begin{eqnarray} \label{FiscI}
\bW_j^{-1} &=& \bI_{n_d} - \bZ_j \bOmega \bG_j^{-1} \bZ_j \tra
\nonumber 
\\
\det \left ( \bW_j \right ) &=& \sigma^{2n_d} \det \left ( \bG_j \right ) \,,
\end{eqnarray}
where $\bG_j = \bI_{r} + \bZ_j \tra \bZ_j \bOmega$.
All the matrices $\bG_j$ have the same dimension, $r \times r$, where
$r$ is the number of variables in $\bZ$.

Assuming that the log-likelihood $l$ has a single maximum,
it can be found as the root of the score vector.
By matrix differentiation we obtain
$$
\frac{\partial l}{\partial \bbeta} \,=\,
\frac{1}{\sigma^2} \sum_{j=1}^m \bX_j \tra \bW_j^{-1} \be_j \,,
$$
and so the maximum likelihood estimator of $\bbeta$ is
the generalised least-squares estimator
$$ \label{FiscB}
\hat{\bbeta} \,=\,
\left ( \sum_{j=1}^m \bX_j \tra \hat{\bW}_j^{-1} \bX_j \right )^{-1}
\sum_{j=1}^m \bX_j \tra \hat{\bW}_j^{-1} \by_j \,.
\tag{7}$$

The matrices $\hat{\bW}_j$ are the estimated versions of $\bW_{j\,}$,
with estimator $\hat{\bOmega}$ substituted for $\bOmega$.
Estimation of $\bOmega$ is described below.

The elementary-level (residual) variance $\sigma^2$ is estimated
by the root of its score,
$$
\frac{\partial l}{\partial \sigma^2} \,=\,
-\frac{1}{2} \left ( \frac{n}{\sigma^2} - \frac{1}{\sigma^4}
\sum_{j=1}^m \be_j \tra \bW_j^{-1} \be_j \right ) \,,
$$
which is
$$
\hat{\sigma}^2 \,=\, \frac{1}{n} \sum_{j=1}^m
\be_j \tra \hat{\bW}_j^{-1} \be_j \,.
$$
The elements of $\bOmega$ are estimated by the Fisher scoring algorithm.
Let these elements, without any redundancy, comprise the vector $\bomega$.
In most applications, $\bomega$ comprises the $\frac{1}{2} r \times (r+1)$
unique elements, $r$ variances and $\frac{1}{2} r \times (r-1)$ covariances.
The Fisher scoring algorithm proceeds by iterations
that update the estimate of $\omega$,
based on the score vector $\bs$ and the expected information matrix $\bH$
evaluated at the current solution.
The vector $\bs$ and matrix $\bH$ are derived by matrix differentiation.
See Magnus and Neudecker (1988) and Harville (1997) for background.

Let $\bi_h$ be the ($r \times 1$) indicator vector for element $h$.
For example, when $r = 5$, $\bi_2 = (0,1,0,0,0)\tra$.
The element of $\bomega$ that corresponds to the (scaled) covariance
$(h, h\pri)$ in $\bOmega$ can be expressed
as
$$
\omega \,=\, \bi_h \tra \bOmega \bi_{h\pri} \,,
$$
and $\partial \bOmega /\partial \omega = \bi_{h\,}\bi_{h\pri}\tra +
\bi_{h\pri\,}\bi_{h}\tra$.
If we replace each (scaled) variance in $\bomega$ by its half,
then this expression holds also for these half-variance parameters.
With this parametrisation, noting that
$\partial \bW_j /\partial \omega =
\bZ_{j\,} \partial \bOmega/ \partial \omega_{\,} \bZ_j\tra$, we have
\begin{eqnarray*}
\frac{\partial l}{\partial \omega} &=&
- \frac{1}{2} \sum_{j=1}^m \left \{
\tr \left ( \bW_j^{-1} \frac{\partial{\bW_j}}{\partial \omega} \right )
- \frac{1}{\sigma^2} \, \be_j \tra \bW_j^{-1}
\frac{\partial{\bW_j}}{\partial \omega} \bW_j^{-1} \be_j \right \}

\\ &=&
\sum_{j=1}^m \left ( - \bi_h \tra bUU_{j\,} \bi_{h\pri}
+ \frac{1}{\sigma^2}\, \buu_j \tra \bi_h \, \buu_j \tra \bi_{h\pri}
\right )

\\ &=&
\sum_{j=1}^m \left ( - U_{j,hh\pri} + \frac{1}{\sigma^2}\, u_{j,h\,}u_{j,h\pri}
\right ) \,,
\end{eqnarray*}
where $bUU_j = \bZ_j \tra \bW_j^{-1} \bZ_{j\,}$,
$\buu_j = \bZ_j \tra \bW_j^{-1} \be_{j\,}$,
$U_{j,hh\pri}$ is the $(h,h\pri)$-element of $bUU_j$
and $u_{j,h}$ the $h$-element of $\buu_{j\,}$.
Multiplying by a vector $\bi_h$ amounts to extracting an element.
Thus, evaluation of the score $\partial l/\partial \omega$
requires calculation of quadratic forms $\bq \tra \bW_j\tra\bZ_{j\,}$.
From \ref{FiscI} we have, for an arbitrary $n_d \times 1$ vector $\bq$,
the identity
$$
\bZ_j \tra \bW_j^{-1} \bq \,=\, \bG_j ^{-1} \bZ_j \tra \bq \,,
$$
so we do not have to form the matrices $\bW_j$ and do not have to
invert any matrices of large size.
Further differentiation yields the expression
{<!-- \arraycolsep 5pt  -->
\begin{eqnarray*}
\frac{\partial ^2 l}{\partial \omega_{1\,} \partial \omega_2} &=&
\sum_{j=1}^m \left [ \bi_{h_1} \tra bUU_{j\,}
\frac{\partial \bOmega}{\partial \omega_2} \,
bUU_j \bi_{h\pri_1} \,-\, \frac{2}{\sigma^2} \left \{
\buu_j \tra \bi_{h_1\,} \buu_j \tra \frac{\partial \bOmega}{\partial \omega_2} \,
bUU_j \bi_{h\pri_1} +
\buu_j \tra \bi_{h\pri_1\,} \buu_j \tra
\frac{\partial \bOmega}{\partial \omega_2} \,
bUU_j \bi_{h_1} \right \} \right ]

\\ &=&
\sum_{j=1}^m \left \{ U_{j,h_1h_2\,}U_{j,h\pri_1h\pri_2\,} +
U_{j,h\pri_1h_2\,}U_{j,h_1h\pri_2\,} \right.

\\ && - \left. \! \frac{2}{\sigma^2}\,
\left ( u_{j,h_1\,}u_{j,h_2\,} U_{j,h\pri_1h\pri_2} +
u_{j,h_1\,}u_{j,h\pri_2\,} U_{j,h\pri_1h_2} +
u_{j,h\pri_1\,}u_{j,h_2\,} U_{j,h_1h\pri_2} +
u_{j,h\pri_1\,}u_{j,h\pri_2\,} U_{j,h_1h_2} \right ) \right \}
\end{eqnarray*}}
\vspace*{-6mm}

<!-- \noindent -->

for $\omega_1$ and $\omega_2$ associated with the respective
elements $(h_{1\,}, h\pri_1)$ and $(h_{2\,}, h\pri_2)$ of $\bOmega$.
The Hessian matrix, its negative expectation, comprises elements
\begin{eqnarray*}
H \left ( \omega_{1\,},\omega_2 \right ) &=& \sum_{j=1}^m \left (
\bi_{h_1} \tra bUU_j \bi_{h_2} \,
\bi_{h\pri_1} \tra bUU_j \bi_{h\pri_2} \,+\,
\bi_{h\pri_1} \tra bUU_j \bi_{h_2} \,
\bi_{h_1} \tra bUU_j \bi_{h\pri_2} \right )

\\ &=&
\sum_{j=1}^m \left ( U_{j,h_1h_2\,}U_{j,h\pri_1h\pri_2} +
U_{j,h\pri_1h_2\,}U_{j,h_1h\pri_2} \right ) \,,
\end{eqnarray*}
which are cluster-level totals of the cross-products of various
elements of $bUU_{j\,}$.
In each iteration $t$, the estimate of $\bomega$ is updated as
$$
\hat{\bomega}^{(t)} \,=\, \hat{\bomega}^{(t-1)} +
r \hat{\bH}_t^{-1} \hat{\bs}_t \,,
$$
where $r=1$, unless the updated matrix $\bOmega^{(t)}$ is not
positive definite.
One way to avoid this is to keep halving $r$ until the updated
matrix $\bOmega^{(t)}$ is positive definite. Having to do so
in many (or all) iterations is usually a sign of having included
too many variables in $\bZ$, and the model should be revised
accordingly.
An alternative approach estimates the Cholesky (or another)
decomposition of $\bOmega$. The iterations are terminated when
the norm of the updating vector $\bH^{-1}\bs$ is sufficiently small.

\vspace*{10mm}

<!-- \noindent -->

Based on an article from Lovric, Miodrag (2011),
''International Encyclopedia of Statistical Science.''
Heidelberg: Springer Science + Business Media, LLC.
<!-- \end{document} -->

<references />

[[Category:Statprob]]