\documentclass{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amscd}
\usepackage{indentfirst}
\usepackage{natbib}
\usepackage{url}
\usepackage[utf8]{inputenc}
\hyphenation{Wagenius}
\DeclareMathOperator{\var}{var}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\cl}{cl}
\newcommand{\real}{\mathbb{R}}
\newcommand{\set}[1]{\{\, #1 \,\}}
\newcommand{\inner}[1]{\langle #1 \rangle}
\newcommand{\abs}[1]{\lvert #1 \rvert}
\newcommand{\norm}[1]{\lVert #1 \rVert}
\newcommand{\fatdot}{\,\cdot\,}
\newcommand{\opand}{\mathbin{\rm and}}
\newcommand{\opor}{\mathbin{\rm or}}
% \setlength{\textheight}{\paperheight}
% \addtolength{\textheight}{ 2 in}
% \setlength{\topmargin}{0.25 pt}
% \setlength{\headheight}{0 pt}
% \setlength{\headsep}{0 pt}
% \addtolength{\textheight}{ \topmargin}
% \addtolength{\textheight}{ \topmargin}
% \addtolength{\textheight}{ \footskip}
% \setlength{\textwidth}{\paperwidth}
% \addtolength{\textwidth}{ 2 in}
% \setlength{\oddsidemargin}{0.25 in}
% \setlength{\evensidemargin}{0.25 in}
% \addtolength{\textwidth}{ \oddsidemargin}
% \addtolength{\textwidth}{ \evensidemargin}
\addtolength{\textheight}{\headsep}
\addtolength{\textheight}{\headheight}
\setlength{\headheight}{0 pt}
\setlength{\headsep}{0 pt}
\begin{document}
\vspace*{0.9375in}
\begin{center}
{\bfseries Aster Models with Random Effects
via Penalized Likelihood} \\
By \\
Charles J. Geyer, Caroline E. Ridley, Robert G. Latta, Julie R. Etterson,
and Ruth G. Shaw \\
Technical Report No.~692 \\
School of Statistics \\
University of Minnesota \\
July 30, 2012 \\
% \today
\end{center}
\thispagestyle{empty}
\cleardoublepage
\setcounter{page}{1}
\thispagestyle{empty}
\begin{abstract}
This technical report works out details of approximate maximum likelihood
estimation for aster models with random effects.
Fixed and random effects are estimated by penalized log likelihood.
Variance components are estimated by integrating out the random effects
in the Laplace approximation of the complete data likelihood following
\citet{b+c}, which can be done analytically,
and maximizing the resulting approximate missing data likelihood.
A further approximation treats the second derivative matrix of the cumulant
function of the exponential family where it appears in the approximate
missing data log likelihood as a constant (not a function of parameters).
Then first and second derivatives of the approximate missing data log
likelihood can be done analytically. Minus the second derivative matrix
of the approximate missing data log likelihood is treated as approximate
Fisher information and used to estimate standard errors.
\end{abstract}
\section{Theory} \label{sec:theo}
Aster models \citep*{gws,aster2} have attracted much recent attention.
Several researchers have raised the issue of incorporating random effects
in aster models, and we do so here.
\subsection{Complete Data Log Likelihood} \label{sec:complete}
Although we are particularly interested in aster models \citep{gws}, our
theory works for any exponential family model. The log likelihood can
be written
$$
l(\varphi) = y^T \varphi  c(\varphi),
$$
where $y$ is the canonical statistic vector,
$\varphi$ is the canonical parameter vector,
and the cumulant function $c$ satisfies
\begin{align}
\mu(\varphi) & = E_{\varphi}(y) = c'(\varphi)
\label{eq:moo}
\\
W(\varphi) & = \var_{\varphi}(y) = c''(\varphi)
\label{eq:w}
\end{align}
where $c'(\varphi)$ denotes the vector of first partial derivatives
and $c''(\varphi)$ denotes the matrix of second partial derivatives.
We assume a canonical affine submodel with random effects determined by
\begin{equation} \label{eq:canaffsub}
\varphi = a + M \alpha + Z b,
\end{equation}
where $a$ is a known vector, $M$ and $Z$ are known matrices, $b$
is a normal random vector with mean vector zero and
variance matrix $D$.
The vector $a$ is called the \emph{offset vector} and the matrices
$M$ and $Z$ are called the \emph{model matrices} for fixed and random
effects, respectively, in the terminology of the R function \texttt{glm}.
(The vector $a$ is called the \emph{origin} in the terminology of the
R function \texttt{aster}. \emph{Design matrix} is alternative
terminology for model matrix.)
The matrix $D$ is assumed to be diagonal,
so the random effects are independent random variables.
The diagonal components of $D$ are called \emph{variance components}
in the classical terminology of random effects models \citep{searle}.
Typically the components of $b$ are divided into blocks having the same
variance \citep[Section~6.1]{searle}, so there are only
a few variance components but many random effects, but nothing in this
document uses this fact.
The unknown parameter vectors are $\alpha$ and $\nu$, where $\nu$ is
the vector of variance components. Thus $D$ is a function of $\nu$,
although this is not indicated by the notation.
The ``complete data log likelihood'' (i.~e., what the log likelihood would
be if the random effect vector $b$ were observed) is
\begin{equation} \label{eq:foo}
l_c(\alpha, b, \nu)
=
l(a + M \alpha + Z b)
 \tfrac{1}{2} b^T D^{ 1} b
 \tfrac{1}{2} \log \det(D)
\end{equation}
in case none of the variance components are zero.
We deal with the case of zero variance components in Sections~\ref{sec:zero},
\ref{sec:raw}, and~\ref{sec:roots} below.
\subsection{Missing Data Likelihood}
Ideally, inference about the parameters should be based on
the \emph{missing data likelihood}, which is the complete data likelihood with
random effects $b$ integrated out
\begin{equation} \label{eq:bar}
L_m(\alpha, \nu) = \int e^{l_c(\alpha, b, \nu)} \, d b
\end{equation}
Maximum likelihood estimates (MLE) of $\alpha$ and $\nu$ are the values that
maximize \eqref{eq:bar}. However MLE are hard to find. The integral in
\eqref{eq:bar} cannot be done analytically, nor can it be done by numerical
integration except in very simple cases.
There does exist a large literature on doing such integrals
by ordinary or Markov chain Monte Carlo
\citep*{penttinen,thompson,g+t,mcmcml,promislow,shawgeyershaw,boothhobert,sung,hunteretal,okabayashigeyer,hummeletal},
but these methods take a great
deal of computing time and are difficult for ordinary users to apply.
We wish to avoid that route if at all possible.
\subsection{A Digression on Minimization} \label{sec:minimize}
The theory of constrained optimization (Section~\ref{sec:raw} below)
has a bias in favor of minimization rather than maximization.
The explication below will be simpler if we switch now from maximization
to minimization (minimizing minus the log likelihood) rather than switch
later.
\subsection{Laplace Approximation} \label{sec:laplace}
\citet{b+c} proposed to replace the integrand in \eqref{eq:bar} by
its Laplace approximation, which is proportional to
a normal probability density function
so the random effects can be integrated out analytically.
Let $b^*$ denote the result of maximizing \eqref{eq:foo} considered
as a function of $b$ for fixed $\alpha$ and $\nu$.
Then $ \log L_m(\alpha, \nu)$ is approximated by
$$
q(\alpha, \nu)
=
\tfrac{1}{2} \log \det[\kappa''(b^*)] + \kappa(b^*)
$$
where
\begin{align*}
\kappa(b)
& =
 l_c(a + M \alpha + Z b)
\\
\kappa'(b)
& =
 Z^T [ y + \mu(a + M \alpha + Z b) ] + D^{1} b
\\
\kappa''(b)
& =
Z^T W(a + M \alpha + Z b) Z + D^{1}
\end{align*}
Hence
\begin{equation} \label{eq:logpickle}
\begin{split}
q(\alpha, \nu)
& =
 l_c(\alpha, b^*, \nu)
+ \tfrac{1}{2} \log
\det\bigl[ \kappa''(b^*) \bigr]
\\
& =
 l(a + M \alpha + Z b^*)
+ \tfrac{1}{2} (b^*)^T D^{1} b^*
+ \tfrac{1}{2} \log \det(D)
\\
& \qquad
+ \tfrac{1}{2} \log
\det\bigl[ Z^T W(a + M \alpha + Z b^*) Z + D^{1} \bigr]
\\
& =
 l(a + M \alpha + Z b^*)
+ \tfrac{1}{2} (b^*)^T D^{1} b^*
\\
& \qquad
+ \tfrac{1}{2} \log
\det\bigl[ Z^T W(a + M \alpha + Z b^*) Z D + I \bigr]
\end{split}
\end{equation}
where $I$ denotes the identity matrix of the appropriate dimension
(which must be the same as the dimension of $D$ for the expression it
appears in to make sense),
where $b^*$ is a function of $\alpha$ and $\nu$
and $D$ is a function of $\nu$, although this is not indicated by the
notation, and where
the last equality uses the rule sum of logs is log of product and
the rule product of determinants is determinant of matrix product
\citep[Theorem~13.3.4]{harville}.
Since minus the log likelihood of an exponential family is a convex function
\citep[Theorem~9.1]{barndorff} and the middle term on the righthand side
of \eqref{eq:foo} is a strictly convex function, it follows that
\eqref{eq:foo} considered as a function of $b$ for fixed $\alpha$ and $\nu$
is a strictly convex function.
Moreover, this function has bounded level sets, because the first term
on the righthand side of \eqref{eq:foo} is bounded
\citep[Theorems~4 and~6]{gdor} and the second term has bounded level sets.
It follows that there is unique global minimizer
\citep[Theorems~1.9 and~2.6]{raw}.
Thus $b^*(\alpha, \nu)$ is well defined for all values
of $\alpha$ and $\nu$.
The key idea is to use \eqref{eq:logpickle} as if it were the log likelihood
for the unknown parameters ($\alpha$ and $\nu$), although it is only an
approximation. However, this is also problematic. In doing likelihood
inference using \eqref{eq:logpickle} we need first and second derivatives
of it (to calculate Fisher information), but $W$ is already
the second derivative matrix of the cumulant function, so
first derivatives of \eqref{eq:logpickle} would involve third derivatives
of the cumulant function and
second derivatives of \eqref{eq:logpickle} would involve fourth derivatives
of the cumulant function. There are no published
formulas for derivatives higher than second of the aster model cumulant
function nor does
software (the R package \texttt{aster},
\citealp{packageaster}) provide such 
the derivatives do, of course, exist because every cumulant function of
a regular exponential family is infinitely differentiable at every
point of the canonical parameter space \citep[Theorem~8.1]{barndorff} 
they are just not readily available.
\citet{b+c} noted the same problem in the context of GLMM, and proceeded as
if $W$ were a constant function of its argument, so all derivatives of $W$
were zero. This is not a bad approximation because ``in asymptopia'' the
aster model log likelihood is exactly quadratic and $W$ is a constant function,
this being a general property of likelihoods \citep{simple}. Hence we adopt
this idea too, more because we are forced to by the difficulty of
differentiating $W$ than by our belief that we are ``in asymptopia.''
This leads to the following idea. Rather than basing inference on
\eqref{eq:logpickle}, we actually use
\begin{equation} \label{eq:logpickletoo}
q(\alpha, \nu)
=
 l(a + M \alpha + Z b^*)
+ \tfrac{1}{2} (b^*)^T D^{1} b^*
+ \tfrac{1}{2} \log
\det\bigl[ Z^T \widehat{W} Z D + I \bigr]
\end{equation}
where $\widehat{W}$ is a constant matrix (not a function of $\alpha$ and
$\nu$). This makes sense for any
choice of $\widehat{W}$ that is symmetric and positive semidefinite,
but we will choose $\widehat{W}$ that are close to
$W(a + M \hat{\alpha} + Z \hat{b})$, where
$\hat{\alpha}$ and $\hat{\nu}$ are the joint minimizers
of \eqref{eq:logpickle} and $\hat{b} = b^*(\hat{\alpha}, \hat{\nu})$.
Note that \eqref{eq:logpickletoo} is a redefinition of $q(\alpha, \nu)$.
Hereafter we will no longer use the definition \eqref{eq:logpickle}.
\subsection{A Key Concept}
Introduce
\begin{equation} \label{eq:key}
p(\alpha, b, \nu)
=
 l(a + M \alpha + Z b)
+ \tfrac{1}{2} b^T D^{1} b
+ \tfrac{1}{2} \log
\det\bigl[ Z^T \widehat{W} Z D + I \bigr]
\end{equation}
where, as the lefthand side says, $\alpha$, $b$, and $\nu$ are all
free variables and, as usual, $D$ is a function of $\nu$, although
the notation does not indicate this.
Since the terms that contain $b$ are the same in both \eqref{eq:foo}
and \eqref{eq:key}, $b^*$ can also be defined as the result of
minimizing \eqref{eq:key} considered
as a function of $b$ for fixed $\alpha$ and $\nu$.
Thus \eqref{eq:logpickletoo} is a profile of \eqref{eq:key}
and $(\hat{\alpha}, \hat{b}, \hat{\nu})$ is the joint minimizer
of \eqref{eq:key}.
Since $p(\alpha, b, \nu)$ is a much simpler function than $q(\alpha, \nu)$,
the latter having no closed form expression and requiring an optimization
as part of each evaluation, it is much simpler to find
$(\hat{\alpha}, \hat{b}, \hat{\nu})$
by minimizing the former rather than the latter.
\subsection{A Digression on Partial Derivatives}
Let $f(\alpha, b, \nu)$ be a scalarvalued function of three vector
variables. We write partial derivative vectors using subscripts:
$f_\alpha(\alpha, b, \nu)$ denotes the vector of partial derivatives
with respect to components of $\alpha$. Our convention is that we take
this to be a column vector. Similarly for $f_b(\alpha, b, \nu)$.
We also use this convention for partial derivatives with respect to
single variables: $f_{\nu_k}(\alpha, b, \nu)$, which are, of course,
scalars. We use this convention for any scalarvalued function
of any number of vector variables.
We continue this convention for second partial derivatives:
$f_{\alpha b}(\alpha, b, \nu)$ denotes the matrix of partial derivatives
having $i, j$ component that is the (mixed) second partial derivative of
$f$ with respect to $\alpha_i$ and $b_j$. Thus the row dimension of
$f_{\alpha b}(\alpha, b, \nu)$ is the dimension of $\alpha$,
the column dimension is the dimension of $b$, and
$f_{b \alpha}(\alpha, b, \nu)$ is the transpose of
$f_{\alpha b}(\alpha, b, \nu)$.
This convention allows easy indication of points at which partial derivatives
are evaluated. For example, $f_{\alpha b}(\alpha, b^*, \nu)$ indicates
that $b^*$ is plugged in for $b$ in the expression for
$f_{\alpha b}(\alpha, b, \nu)$.
We also use this convention of subscripts denoting partial derivatives
with vectorvalued functions.
If $f(\alpha, b, \nu)$ is a columnvectorvalued function of vector
variables, then $f_\alpha(\alpha, b, \nu)$ denotes the matrix of
partial derivatives having $i, j$ component that is the partial derivative
of the $i$th component of $f_\alpha(\alpha, b, \nu)$ with respect to
$\alpha_j$.
Thus the row dimension of
$f_{\alpha}(\alpha, b, \nu)$ is the dimension of
$f(\alpha, b, \nu)$ and
the column dimension is the dimension of $\alpha$.
\subsection{First Derivatives}
Start with \eqref{eq:key}. Its derivatives are
\begin{equation} \label{eq:palpha}
p_\alpha(\alpha, b, \nu)
=
 M^T \bigl[ y  \mu(a + M \alpha + Z b) \bigr]
\end{equation}
\begin{equation} \label{eq:pc}
p_b(\alpha, b, \nu)
=
 Z^T \bigl[ y  \mu(a + M \alpha + Z b) \bigr] + D^{ 1} b
\end{equation}
and
\begin{equation} \label{eq:ptheta}
p_{\nu_k}(\alpha, b, \nu)
=
 \tfrac{1}{2} b^T D^{ 1} E_k D^{ 1} b
+ \tfrac{1}{2} \tr \Bigl(
\bigl[ Z^T \widehat{W} Z D + I \bigr]^{ 1}
Z^T \widehat{W} Z E_k
\Bigr)
\end{equation}
where
\begin{equation} \label{eq:eek}
E_k = D_{\nu_k}(\nu)
\end{equation}
is the diagonal matrix whose components are equal to one
if the corresponding components of $D$ are equal to $\nu_k$ by
definition (rather than by accident when some other component of $\nu$
also has the same value) and whose components are otherwise zero.
The formula for the derivative of a matrix inverse
comes from \citet[Chapter~15, Equation~8.15]{harville}.
The formula for the derivative of the log of a determinant
comes from \citet[Chapter~15, Equation~8.6]{harville}.
The estimating equation for $b^*$ can be written
\begin{equation} \label{eq:estimatingcstar}
p_b\bigl(\alpha, b^*, \nu\bigr) = 0
\end{equation}
and by the multivariate chain rule \citep[Theorem~8.15]{browder}
we have
\begin{equation} \label{eq:qalpha}
\begin{split}
q_\alpha(\alpha, \nu)
& =
p_\alpha(\alpha, b^*, \nu)
+
b^*_\alpha(\alpha, \nu)^T p_b(\alpha, b^*, \nu)
\\
& =
p_\alpha(\alpha, b^*, \nu)
\end{split}
\end{equation}
by \eqref{eq:estimatingcstar}, and
\begin{equation} \label{eq:qtheta}
\begin{split}
q_{\nu_k}(\alpha, \nu)
& =
b^*_{\nu_k}(\alpha, \nu)^T
p_b(\alpha, b^*, \nu)
+
p_{\nu_k}(\alpha, b^*, \nu)
\\
& =
p_{\nu_k}(\alpha, b^*, \nu)
\end{split}
\end{equation}
again by \eqref{eq:estimatingcstar}.
\subsection{Second Derivatives} \label{sec:second}
By the multivariate chain rule \citep[Theorem~8.15]{browder}
\begin{align*}
q_{\alpha \alpha}(\alpha, \nu)
& =
p_{\alpha \alpha}(\alpha, b^*, \nu)
+
p_{\alpha b}(\alpha, b^*, \nu) b^*_\alpha(\alpha, \nu)
\\
q_{\alpha \nu}(\alpha, \nu)
& =
p_{\alpha \nu}(\alpha, b^*, \nu)
+
p_{\alpha b}(\alpha, b^*, \nu) b^*_{\nu}(\alpha, \nu)
\\
q_{\nu \nu}(\alpha, \nu)
& =
p_{\nu \nu}(\alpha, b^*, \nu)
+
p_{\nu b}(\alpha, b^*, \nu) b^*_{\nu}(\alpha, \nu)
\end{align*}
The estimating equation \eqref{eq:estimatingcstar}
defines $b^*$ implicitly. Thus derivatives of $b^*$ are computed using
the implicit function theorem \citep[Theorem~8.29]{browder}
\begin{align}
b^*_\alpha(\alpha, \nu)
& =

p_{b b}(\alpha, b^*, \nu)^{1} p_{b \alpha}(\alpha, b^*, \nu)
\label{eq:cstaralpharewrite}
\\
b^*_{\nu}(\alpha, \nu)
& =

p_{b b}(\alpha, b^*, \nu)^{1} p_{b \nu}(\alpha, b^*, \nu)
\label{eq:cstarthetarewrite}
\end{align}
This theorem requires that $p_{b b}(\alpha, b^*, \nu)$ be invertible,
and we shall see below that it is.
Then the second derivatives above can be rewritten
\begin{align*}
q_{\alpha \alpha}(\alpha, \nu)
& =
p_{\alpha \alpha}(\alpha, b^*, \nu)

p_{\alpha b}(\alpha, b^*, \nu)
p_{b b}(\alpha, b^*, \nu)^{1} p_{b \alpha}(\alpha, b^*, \nu)
\\
q_{\alpha \nu}(\alpha, \nu)
& =
p_{\alpha \nu}(\alpha, b^*, \nu)

p_{\alpha b}(\alpha, b^*, \nu)
p_{b b}(\alpha, b^*, \nu)^{1} p_{b \nu}(\alpha, b^*, \nu)
\\
q_{\nu \nu}(\alpha, \nu)
& =
p_{\nu \nu}(\alpha, b^*, \nu)

p_{\nu b}(\alpha, b^*, \nu)
p_{b b}(\alpha, b^*, \nu)^{1} p_{b \nu}(\alpha, b^*, \nu)
\end{align*}
a particularly simple and symmetric form. If we combine all the parameters
in one vector $\psi = (\alpha, \nu)$ and write $p(\psi, b)$ instead
of $p(\alpha, b, \nu)$ we have
\begin{equation} \label{eq:psipsi}
q_{\psi \psi}(\psi)
=
p_{\psi \psi}(\psi, b^*)

p_{\psi b}\bigl(\psi, b^*\bigr)
p_{b b}\bigl(\psi, b^*\bigr)^{ 1}
p_{b \psi}\bigl(\psi, b^*\bigr)
\end{equation}
This form is familiar from the conditional variance formula
for normal distributions if
\begin{equation} \label{eq:fat}
\begin{pmatrix} \Sigma_{1 1} & \Sigma_{1 2}
\\ \Sigma_{2 1} & \Sigma_{2 2} \end{pmatrix}
\end{equation}
is the partitioned variance matrix of a partitioned normal random vector
with components $X_1$ and $X_2$, then the variance matrix of the conditional
distribution of $X_1$ given $X_2$ is
\begin{equation} \label{eq:thin}
\Sigma_{1 1}  \Sigma_{1 2} \Sigma_{2 2}^{ 1} \Sigma_{2 1}
\end{equation}
assuming that $X_2$ is nondegenerate \citep[Theorem~2.5.1]{anderson}.
Moreover, if the conditional distribution is degenerate, that is, if there
exists a nonrandom vector $v$ such that $\var(v^T X_1 \mid X_2) = 0$, then
$$
v^T X_1 = v^T \Sigma_{1 2} \Sigma_{2 2}^{ 1} X_2
$$
with probability one, assuming $X_1$ and $X_2$ have mean zero
\citep[also by][Theorem~2.5.1]{anderson}, and the joint distribution
of $X_1$ and $X_2$ is also degenerate. Thus we conclude that if
the (joint) Hessian matrix of $p$ is nonsingular, then so is the (joint)
Hessian matrix of $q$ given by \eqref{eq:psipsi}.
The remaining work for this section is deriving second derivatives of $p$
\begin{align*}
p_{\alpha \alpha}(\alpha, b, \nu)
& =
M^T W(a + M \alpha + Z b) M
\\
p_{\alpha b}(\alpha, b, \nu)
& =
M^T W(a + M \alpha + Z b) Z
\\
p_{b b}(\alpha, b, \nu)
& =
Z^T W(a + M \alpha + Z b) Z + D^{ 1}
\\
p_{\alpha \nu_k}(\alpha, b, \nu)
& =
0
\\
p_{b \nu_k}(\alpha, b, \nu)
& =
 D^{ 1} E_k D^{ 1} b
\\
p_{\nu_j \nu_k}(\alpha, b, \nu)
& =
b^T D^{ 1} E_j D^{ 1} E_k D^{ 1} b
\\
& \qquad

\tfrac{1}{2} \tr \Bigl(
\bigl[ Z^T \widehat{W} Z D + I \bigr]^{ 1}
Z^T \widehat{W} Z E_j
\\
& \qquad \qquad
\bigl[ Z^T \widehat{W} Z D + I \bigr]^{ 1}
Z^T \widehat{W} Z E_k
\Bigr)
\end{align*}
This finishes the derivation of all the derivatives we need.
Recall that in our use of the implicit function theorem we needed
$p_{b b}(\alpha, b^*, \nu)$ to be invertible. From the explicit form
given above we see that it is actually negative definite, because
$W(a + M \alpha + Z b)$ is positive semidefinite by \eqref{eq:w}.
\subsection{Zero Variance Components} \label{sec:zero}
When some variance components are zero, the corresponding diagonal components
of $D$ are zero, and the corresponding components of $b$ are zero almost surely.
However we deal with this situation, it must have the same effect as omitting
those variance components and the corresponding random effects from the model.
\citet[Section 2.3]{b+c} suggest using the MoorePenrose pseudoinverse
\citep[Chapter~20]{harville}. Let $D^+$ denote the diagonal matrix whose
diagonal components are the reciprocals of the diagonal components
provided those are nonzero and whose diagonal components are zero otherwise.
Then
\begin{equation} \label{eq:logpickletootoo}
q(\alpha, \nu)
=
 l(a + M \alpha + Z b^*)
+ \tfrac{1}{2} (b^*)^T D^+ b^*
+ \tfrac{1}{2} \log
\det\bigl[ Z^T \widehat{W} Z D + I \bigr]
\end{equation}
is an approximate log likelihood, but in the calculation of
$b^*$ constrained penalized maximum likelihood must be used:
elements of $b$ corresponding to zero variance components must be
constrained to be zero, because \eqref{eq:logpickletootoo} does
not force them to be zero.
Although this proposal \citep[Section 2.3]{b+c} does deal with the
situation where the zero variance components are somehow known, it
does not adequately deal with estimating which variance components are
zero. That is the subject of the following two sections.
\subsection{The Theory of Constrained Optimization} \label{sec:raw}
\subsubsection{Incorporating Constraints in the Objective Function}
When zero variance components arise, optimization of \eqref{eq:key}
puts us in the
realm of constrained optimization. The theory of constrained optimization
\citep{raw} has a notational bias towards minimization \citep[p.~5]{raw}.
Thus, as explained above (Section~\ref{sec:minimize}) we have switched
from maximization to minimization.
Readers who are familiar with KarushKuhnTucker\
(\citealp[Section~9.1]{fletcher}; \citealp[Section~12.2]{nocedalwright})
theory should be warned
that that theory is not adequate for the problem at hand, because the
constraint set is not a closed set and so cannot be defined in terms of
smooth constraint functions.
Thus the need for the more general theory \citep{raw}.
The theory of constrained optimization incorporates constraints in the
objective function by the simple device of defining the objective function
(for a minimization problem) to have the value $+ \infty$ off the constraint
set \citep[Section~1A]{raw}.
Since no point where the objective function has the value $+ \infty$
can minimize it, unless the the objective function has the value $+ \infty$
everywhere, which is not the case in any application, the unconstrained
minimizer of this sort of objective function is the same as the constrained
minimizer.
Thus we need to impose constraints on our key function
\eqref{eq:key}, requiring that
each component of $\nu$ be nonnegative and when any component of $\nu$ is
zero the corresponding components of $b$ are also zero. However,
the formula \eqref{eq:key} does not make sense when components of $\nu$
are zero, so we proceed differently.
\subsubsection{Lower Semicontinuous Regularization}
Since all but the middle term on the righthand side of
\eqref{eq:key} are actually defined on some neighborhood of each point
of the constraint set and differentiable at each point
of the constraint set, we only need to deal with the middle term.
It is the sum of terms of the form $b_i^2 / \nu_k$, where $\nu_k$ is
the variance of $b_i$. Thus we investigate functions of this form
\begin{equation} \label{eq:h}
h(b, \nu)
=
b^2 / \nu
\end{equation}
where, temporarily, $b$ and $\nu$ are scalars rather than vectors
(representing single components of the vectors). In case $\nu > 0$
we have derivatives
\begin{align*}
h_b(b, \nu) & = 2 b / \nu
\\
h_\nu(b, \nu) & =  b^2 / \nu^2
\\
h_{b b}(b, \nu) & = 2 / \nu
\\
h_{b \nu}(b, \nu) & =  2 b / \nu^2
\\
h_{\nu \nu}(b, \nu) & = 2 b^2 / \nu^3
\end{align*}
The Hessian matrix
$$
h''(b, \nu)
=
\begin{pmatrix}
2 / \nu &  2 b / \nu^2
\\
 2 b / \nu^2 & 2 b^2 / \nu^3
\end{pmatrix}
$$
has nonnegative determinants of its principal submatrices, since the
diagonal components are positive and $\det\bigl(h''(b, \nu)\bigr)$ is zero.
Thus the Hessian matrix is nonnegative definite
\citep[Theorem~14.9.11]{harville}, which implies that $h$ itself is
convex \citep[Theorem~2.14]{raw} on the set where $\nu > 0$.
We then extend $h$ to the whole of the constraint set (this just adds
the origin to the points already considered) in two steps.
First we define it to have the value $+ \infty$ at all points not
yet considered (those where any component of $\nu$ is nonpositive).
This gives us an extendedrealvalued convex function defined on
all of $\real^2$.
Second we take it to be the
lower semicontinuous (LSC) regularization \citep[p.~14]{raw} of the
function just defined. It is clear that
$$
\liminf_{\substack{b \to \bar{b}\\\nu \searrow 0}} h(b, \nu)
=
\begin{cases}
0, & \bar{b} = 0
\\
+ \infty, & \text{otherwise}
\end{cases}
$$
Thus the LSC regularization is
\begin{equation} \label{eq:hlsc}
h(b, \nu)
=
\begin{cases}
b^2 / \nu, & \nu > 0
\\
0, & \nu = 0 \opand b = 0
\\
+ \infty, & \text{otherwise}
\end{cases}
\end{equation}
The LSC regularization of a convex function is convex
\citep[Proposition~2.32]{raw}, so \eqref{eq:hlsc} defines
an extendedrealvalued convex function.
Note that $h(b, 0)$ considered as a function of $b$
is minimized at $b = 0$ because that is the only point where this function
is finite. Hence this does enforce the constraint that random effects
corresponding to zero variance components must be zero.
Let $k$ denote the map from indices for $b$ to indices for $\nu$ that gives
corresponding components: $\nu_{k(i)}$ is the variance of $b_i$.
Let $\dim(b)$ denote the number of random effects. Then our objective function
can be written
\begin{equation} \label{eq:keylsc}
p(\alpha, b, \nu)
=
 l(a + M \alpha + Z b)
+ \tfrac{1}{2} \sum_{i = 1}^{\dim(b)} h(b_i, \nu_{k(i)})
+ \tfrac{1}{2} \log
\det\bigl[ Z^T \widehat{W} Z D + I \bigr]
\end{equation}
where $h$ is given by \eqref{eq:hlsc}, provided all of the components of
$\nu$ are nonnegative. The proviso is necessary because the third term
on the righthand side is not defined for all values of $\nu$, only those
such that the argument of the determinant is a positive definite matrix.
Hence, we must separately define $p(\alpha, b, \nu) = + \infty$
whenever any component of $\nu$ is negative.
\subsubsection{Subderivatives}
In calculus we learn that the first derivative is zero at a local minimum
and, therefore, to check points where the first derivative is zero.
This is called Fermat's rule. This
rule no longer works for nonsmooth functions, including those that incorporate
constraints, such as \eqref{eq:keylsc}. It does, of course, still work
at points in the interior of the constraint set where \eqref{eq:keylsc}
is differentiable. It does not work to check points on the boundary.
There we need what \citet[Theorem~10.1]{raw}
call \emph{Fermat's rule, generalized:} at a local minimum the
\emph{subderivative function} is nonnegative.
For any extendedrealvalued function $f$ on $\real^d$,
the \emph{subderivative function}, denoted $d f(x)$
is also an extendedrealvalued function on $\real^d$ defined by
$$
d f(x)(\bar{w}) = \liminf_{\substack{\tau \searrow 0 \\ w \to \bar{w}}}
\frac{f(x + \tau w)  f(x)}{\tau}
$$
\citep[Definition~8.1]{raw}. The notation on the lefthand side is read
the subderivative of $f$ at the point $x$ in the direction $\bar{w}$.
Fortunately, we do not have to use this
definition to calculate subderivatives we want, because the calculus
of subderivatives allows us to use simpler formulas in special cases.
Firstly, there is the notion of subdifferential regularity
\citep[Definition~7.25]{raw}, which we can use without knowing the definition.
The sum of regular functions is regular and the subderivative of a sum
is the sum of the subderivatives \citep[Corollary~10.9]{raw}.
A smooth function is regular and the subderivative is given by
\begin{equation} \label{eq:subderivativesmooth}
d f(x)(w) = w^T f'(x),
\end{equation}
where, as in Sections~\ref{sec:complete} and~\ref{sec:laplace} above,
$f'(x)$ denotes the gradient vector (the vector of partial derivatives)
of $f$ at the point $x$ \citep[Exercise~8.20]{raw}. Every LSC convex function
is regular \citep[Example~7.27]{raw}. Thus in computing subderivatives
of \eqref{eq:keylsc} we may compute them term by term, and for the
first and last terms, they are given in terms of the partial derivatives
already computed by \eqref{eq:subderivativesmooth}.
For an LSC convex function $f$, we have the following characterization
of the subderivative \citep[Proposition~8.21]{raw}. At any point $x$
where $f(x)$ is finite, the limit
$$
g(w) = \lim_{\tau \searrow 0} \frac{f(x + \tau w)  f(x)}{\tau}
$$
exists and defines a sublinear function $g$, and then $d f(x)$
is the LSC regularization of $g$.
An extendedrealvalued
function $g$ is \emph{sublinear} if $g(0) = 0$ and
$$
g(a_1 x_1 + a_2 x_2) \le a_1 g(x_1) + a_2 g(x_2)
$$
for all vectors $x_1$ and $x_2$ and positive scalars $a_1$ and $a_2$
\citep[Definition~3.18]{raw}.
The subderivative function of every regular LSC function
is sublinear \citep[Theorem~7.26]{raw}.
So let us proceed to calculate the subderivative of \eqref{eq:hlsc}.
In the interior of the constraint set, where this function is smooth,
we can use the partial derivatives already calculated
$$
d h(b, \nu)(u, v) = \frac{2 b u}{\nu}  \frac{b^2 v}{\nu^2}
$$
where the notation on the lefthand side means the subderivative of $h$
at the point $(b, \nu)$ in the direction $(u, v)$.
On the boundary of the constraint set, which consists of the single point
$(0, 0)$, we take limits. In case $v > 0$, we have
$$
\lim_{\tau \searrow 0}
\frac{h(\tau u, \tau v)  h(0, 0)}{\tau}
=
\lim_{\tau \searrow 0}
\frac{\tau^2 u^2 / (\tau v)}{\tau}
=
\lim_{\tau \searrow 0}
\frac{u^2}{v}
=
\frac{u^2}{v}
$$
In case $v < 0$ or in case $v = 0$ and $u \neq 0$,
\begin{equation} \label{eq:differencequotient}
\frac{h(\tau u, \tau v)  h(0, 0)}{\tau}
\end{equation}
is equal to $+ \infty$ for all $\tau > 0$ so the limit inferior is $+ \infty$.
In case $v = 0$ and $u = 0$, \eqref{eq:differencequotient} is equal to zero
for all $\tau > 0$ so the limit inferior is zero.
Thus we see that the limit inferior already defines an LSC function and
$$
d h(0, 0)(u, v) = h(u, v).
$$
\subsubsection{Applying the Generalization of Fermat's Rule}
The theory of constrained optimization tells us nothing we did not already
know (from Fermat's rule) about smooth functions. The only way we can
have $d f(x)(w) = w^T f'(x) \ge 0$ for all vectors $w$ is if $f'(x) = 0$.
It is only at points where the function is nonsmooth, in the cases of
interest to us, points on the boundary of the constraint set, where
the theory of constrained optimization tells us things we did not know and
need to know.
Even on the boundary, the conclusions of the theory about components
of the state that are not on the boundary agree with what we already knew.
We have
$$
d p(\alpha, b, \nu)(s, u, v)
=
s^T p_\alpha(\alpha, b, \nu)
+
\text{terms not containing $s$}
$$
and the only way this can be nonnegative for all $s$ is if
\begin{equation} \label{eq:descentalpha}
p_\alpha(\alpha, b, \nu) = 0
\end{equation}
in which case
$d p(\alpha, b, \nu)(s, u, v)$ is a constant function of $s$, or, what
is the same thing in other words, the terms of
$d p(\alpha, b, \nu)(s, u, v)$ that appear to involve $s$ are all zero
(and so do not actually involve $s$).
Similarly, $d p(\alpha, b, \nu)(s, u, v) \ge 0$ for all $u_i$ and $v_j$
such that $\nu_j > 0$ and $k(i) = j$ only if
\begin{equation} \label{eq:descentothersmooth}
\begin{split}
p_{\nu_j}(\alpha, b, \nu) & = 0, \qquad \text{$j$ such that $\nu_j > 0$}
\\
p_{b_i}(\alpha, b, \nu) & = 0, \qquad \text{$i$ such that $\nu_{k(i)} > 0$}
\end{split}
\end{equation}
in which case we conclude that
$d p(\alpha, b, \nu)(s, u, v)$ is a constant function of such $u_i$ and $v_j$.
Thus, assuming that we are at a point $(\alpha, b, \nu$)
where \eqref{eq:descentalpha}
and \eqref{eq:descentothersmooth} hold,
and we do assume this throughout the rest of this section,
$d p(\alpha, b, \nu)(s, u, v)$ actually involves only $v_j$ and $u_i$ such
that $\nu_j = 0$ and $k(i) = j$. Define
\begin{equation} \label{eq:keysmooth}
\bar{p}(\alpha, b, \nu)
=
 l(a + M \alpha + Z b)
+ \tfrac{1}{2} \log
\det\bigl[ Z^T \widehat{W} Z D + I \bigr]
\end{equation}
(the part of \eqref{eq:keylsc} consisting of the smooth terms).
Then
\begin{equation} \label{eq:descentothernonsmooth}
\begin{split}
d p(\alpha, b, \nu)(s, u, v)
& =
\sum_{j \in J}
\Biggl[
v_j \bar{p}_{\nu_j}(\alpha, b, \nu)
\\
& \qquad
+
\sum_{i \in k^{ 1}(j)}
\Biggl(
u_i \bar{p}_{b_i}(\alpha, b, \nu)
+
h(u_i, v_j)
\Biggr)
\Biggr]
\end{split}
\end{equation}
where $J$ is the set of $j$ such that $\nu_j = 0$,
where $k^{ 1}(j)$ denotes the set of $i$ such that $k(i) = j$,
and where $h$ is defined by \eqref{eq:hlsc}.
Fermat's rule generalized says we must consider all of the terms
of \eqref{eq:descentothernonsmooth} together.
We cannot consider partial derivatives, because the partial derivatives
do not exist. To check that we are at a local minimum we need to show
that \eqref{eq:descentothernonsmooth} is nonnegative
for all vectors $u$ and $v$.
Conversely, to verify that we are not at a local minimum, we need to find
one pair of vectors $u$ and $v$ such that \eqref{eq:descentothernonsmooth}
is negative. Such a pair $(u, v)$ we call a \emph{descent direction}.
Since Fermat's rule generalized is a necessary but not sufficient condition
(like the ordinary Fermat's rule), the check that we are at a local minimum
is not definitive, but the check that we are not is. If a descent direction
is found, then moving in that direction away from
the current value of $(\alpha, b, \nu)$ will decrease the objective
function \eqref{eq:keylsc}.
So how do we find a descent direction? We want to minimize
\eqref{eq:descentothernonsmooth} considered as a function of $u$ and $v$
for fixed $\alpha$, $b$, and $\nu$.
On further consideration, we can consider the terms of
\eqref{eq:descentothernonsmooth} for each $j$ separately.
If the minimum of
\begin{equation} \label{eq:descentothernonsmoothj}
v_j \bar{p}_{\nu_j}(\alpha, b, \nu)
+
\sum_{i \in k^{ 1}(j)}
\Biggl(
u_i \bar{p}_{b_i}(\alpha, b, \nu)
+
h(u_i, v_j)
\Biggr)
\end{equation}
over all vectors $u$ and $v$ is nonnegative, then the minimum is zero,
because \eqref{eq:descentothernonsmoothj} has the value zero
when $u = 0$ and $v = 0$. Thus we can ignore this $j$ in calculating
the descent direction.
On the other hand, if the minimum is negative,
then the minimum does not occur at $v = 0$ and the minimum is actually
$ \infty$ by the sublinearity of the subderivative,
one consequence of sublinearity being positive homogeneity
$$
d f(x)(\tau w) = \tau d f(x)(w), \qquad \tau \ge 0
$$
which holds
for any subderivative. Thus (as our terminology hints) we are only trying
to find a descent \emph{direction}, the length of the vector $(u, v)$ does
not matter, only its direction. Thus to get a finite minimum we can
do a constrained minimization of \eqref{eq:descentothernonsmoothj},
constraining $(u, v)$ to lie in a ball. This is found by the wellknown
KarushKuhnTucker theory of constrained optimization
(\citealp[Section~9.1]{fletcher}; \citealp[Section~12.2]{nocedalwright})
to be the minimum of the Lagrangian function
\begin{equation} \label{eq:laggard}
L(u, v)
=
\lambda v_j^2
+
v_j \bar{p}_{\nu_j}(\alpha, b, \nu)
+
\sum_{i \in k^{ 1}(j)}
\Biggl(
\lambda
u_i^2
+
u_i \bar{p}_{b_i}(\alpha, b, \nu)
+
\frac{u_i^2}{v_j}
\Biggr)
\end{equation}
where $\lambda > 0$ is the Lagrange multiplier, which would have to be
adjusted if we were interested in constraining $(u, v)$ to lie in a particular
ball. Since we do not care about the length of $(u, v)$ we can use any
$\lambda$. We have replaced $h(u_i, v_i)$ by $u_i^2 / v_j$ because we
know that if we are finding an actual descent direction, then we will
have $v_j > 0$. Now
\begin{align*}
L_{u_i}(u, v)
& =
2 \lambda u_i
+
\bar{p}_{b_i}(\alpha, b, \nu)
+
\frac{2 u_i}{v_j},
\qquad
i \in k^{ 1}(j)
\\
L_{v_j}(u, v)
& =
2 \lambda v_j
+
\bar{p}_{\nu_j}(\alpha, b, \nu)

\sum_{i \in k^{ 1}(j)}
\frac{u_i^2}{v_j^2}
\end{align*}
The minimum occurs where these are zero.
Setting the first equal to zero and solving for $u_i$ gives
$$
\hat{u}_i(v_j)
=
 \frac{\bar{p}_{b_i}(\alpha, b, \nu)}{2 (\lambda + 1 / v_j)}
$$
plugging this back into the second gives
\begin{align*}
L_{v_j}\bigl(\hat{u}(v), v\bigr)
& =
2 \lambda v_j
+
\bar{p}_{\nu_j}(\alpha, b, \nu)

\frac{1}{4 (\lambda v_j + 1)^2}
\sum_{i \in k^{ 1}(j)}
\bar{p}_{b_i}(\alpha, b, \nu)^2
\end{align*}
and we seek zeros of this. The righthand is clearly an increasing
function of $v_j$ so it is negative somewhere only if it is negative when
$v_j = 0$ where it has the value
\begin{equation} \label{eq:descenttest}
\bar{p}_{\nu_j}(\alpha, b, \nu)

\frac{1}{4}
\sum_{i \in k^{ 1}(j)}
\bar{p}_{b_i}(\alpha, b, \nu)^2
\end{equation}
So that gives us a test for a descent direction: we have a descent
direction if and only if \eqref{eq:descenttest} is negative.
Conversely, we appear to have $\hat{\nu}_j = 0$ if \eqref{eq:descenttest}
is nonnegative.
That finishes our treatment of the theory of constrained optimization.
We have to ask is all of this complication really necessary?
It turns out that it is and it isn't. We can partially avoid it
by a change of variables. But the cure is worse than the disease
in some ways. This is presented in the following section.
\subsection{Square Roots} \label{sec:roots}
We can avoid constrained optimization by the following change of parameter.
Introduce new parameter variables by
\begin{align*}
\nu_j & = \sigma_j^2
\\
b & = A c
\end{align*}
where $A$ is diagonal and $A^2 = D$, so the $i$th diagonal component of
$A$ is $\sigma_{k(i)}$.
Then the objective function \eqref{eq:key} becomes
\begin{equation} \label{eq:keyrooted}
\tilde{p}(\alpha, c, \sigma)
=
 l(a + M \alpha + Z A c)
+ \tfrac{1}{2} c^T c
+ \tfrac{1}{2} \log
\det\bigl[ Z^T \widehat{W} Z A^2 + I \bigr]
\end{equation}
There are now no constraints and \eqref{eq:keyrooted} is a continuous
function of all variables.
The drawback is that by symmetry we must have
$\tilde{p}_{\sigma_j}(\alpha, c, \sigma)$ equal to zero when $\sigma_j = 0$.
Thus first derivatives become useless for checking for descent directions,
and second derivative information is necessary.
However, that is not the way unconstrained optimizers like the R functions
\texttt{optim} and \texttt{nlminb} work. They do not expect such pathological
behavior and do not deal with it correctly. If we want to use such optimizers
to find local minima of \eqref{eq:keyrooted}, then we must provide starting
points that have no component of $\nu$ equal to zero,
and hope that the optimizer will never get any component of $\nu$ close to
zero unless zero actually is a solution. But this is only a hope.
The theory that guided the design of these optimizers does not provide
any guarantees for this kind of objective function.
Moreover, optimizer algorithms stop when close to but not exactly
at a solution, a consequence of inexactness of computer arithmetic.
Thus when the optimizer stops and declares convergence with one or more
components of $\nu$ close to zero, how do we know whether the true solution
is exactly zero or not? We don't unless we return to the original
parameterization and apply the theory of the preceding section.
The question of whether the MLE of the variance components are exactly
zero or not is of scientific interest, so it seems that the device of this
section does not entirely avoid the theory of constrained optimization.
We must change back to the original parameters and use \eqref{eq:descenttest}
to determine whether or not we have $\nu_j = 0$.
Finally, there is another issue with this ``square root'' parameterization.
For this new parameterization,
the analogs of the second derivative formulas derived
in Section~\ref{sec:second} above are
extraordinarily illbehaved. The Hessian matrices are badly conditioned
and sometimes turn out to be not positive definite when calculated by
the computer's arithmetic (which is inexact) even though theory says
they must be positive definite. We know this because at one point we
thought that this ``square root'' parameterization was the answer to
everything and tried to use it everywhere. Months of frustration ensued
where it mostly worked, but failed on a few problems. It took us a long
time to see that it is fundamentally wrongheaded. As we said above,
the cure is worse than the disease.
Thus we concluded that, while we may use this ``square root'' parameterization
to do unconstrained rather than constrained minimization,
we should use it only for that.
The test \eqref{eq:descenttest} should be used to determine whether
variance components are exactly zero or not, and the formulas in
Section~\ref{sec:second} should be used to derive Fisher information.
\subsubsection{First Derivatives}
Some of R's optimization routines can use first derivative information,
thus we derive first derivatives in this parameterization.
\begin{align}
\tilde{p}_\alpha(\alpha, c, \sigma)
& =
 M^T [ y  \mu(a + M \alpha + Z A c) ]
\label{eq:keyrootedalpha}
\\
\tilde{p}_c(\alpha, c, \sigma)
& =
 A Z^T [ y  \mu(a + M \alpha + Z A c) ]
+ c
\label{eq:keyrootedc}
\\
\tilde{p}_{\sigma_j}(\alpha, c, \sigma)
& =
 c^T E_j Z^T [ y  \mu(a + M \alpha + Z A c) ]
\nonumber
\\
& \qquad
+
\tr\left( [ Z^T \widehat{W} Z A^2 + I \bigr]^{ 1}
Z^T \widehat{W} Z A E_j
\right)
\label{eq:keyrootedsigma}
\end{align}
where $E_j$ is given by \eqref{eq:eek}.
\subsection{Fisher Information} \label{sec:fisher}
The observed Fisher information matrix is minus the second derivative matrix
of the log likelihood. As we said above, we want to do this in the original
parameterization.
Assembling stuff derived in preceding sections and introducing
\begin{align*}
\mu^*
& =
\mu\bigl(a + M \alpha + Z b^*(\alpha, \nu) \bigr)
\\
W^*
& =
W\bigl(a + M \alpha + Z b^*(\alpha, \nu) \bigr)
\\
H^* & = Z^T W^* Z + D^{ 1}
\\
\widehat{H} & = Z^T \widehat{W} Z D + I
\end{align*}
we obtain
\begin{align*}
q_{\alpha \alpha}(\alpha, \nu)
& =
M^T W^* M  M^T W^* Z (H^*)^{ 1} Z^T W^* M
\\
q_{\alpha \nu_j}(\alpha, \nu)
& =
M^T W^* Z (H^*)^{ 1} D^{ 1} E_j D^{1} b^*
\\
q_{\nu_j \nu_k}(\alpha, \nu)
& =
(b^*)^T D^{ 1} E_j D^{ 1} E_k D^{ 1} b^*
\\
& \qquad

\tfrac{1}{2} \tr \Bigl( \widehat{H}^{ 1} Z^T \widehat{W} Z E_j
\widehat{H}^{ 1} Z^T \widehat{W} Z E_k \Bigr)
\\
& \qquad

(b^*)^T D^{ 1} E_j D^{ 1}
(H^*)^{ 1}
D^{ 1} E_k D^{ 1} b^*
\end{align*}
In all of these $b^*$, $\mu^*$, $W^*$, and $H^*$ are functions of $\alpha$
and $\nu$ even though the notation does not indicate this.
It is tempting to think expected Fisher information simplifies things
because we ``know'' $E(y) = \mu$ and $\var(y) = W$, except we don't know
that! What we do know is
$$
E(y \mid b) = \mu(a + M \alpha + Z b)
$$
but we don't know how to take the expectation of the right hand side
(and similarly for the variance). Rather than introduce further approximations
of dubious validity, it seems best to just use (approximate)
observed Fisher information.
\subsection{REML?}
\citet{b+c} do not maximize the approximate log likelihood
\eqref{eq:logpickle}, but make further
approximations to give estimators motivated by REML (restricted maximum
likelihood) estimators for linear mixed models (LMM).
\citet{b+c} concede that
the argument that justifies REML estimators for LMM does not carry over to
their REMLlike estimators for generalized linear mixed models (GLMM).
Hence these REMLlike estimators have no mathematical justification.
Even in LMM the widely used procedure of following REML estimates of the
variance components with socalled BLUE estimates of fixed effects and BLUP
estimates of random effects, which are actually only BLUE and BLUP if the
variance components are assumed known rather than estimated, is obviously
wrong: ignoring the fact that the variance components are estimated cannot be
justified (and \citeauthor{b+c} say this in their discussion section).
Hence REML is not justified even in LMM when fixed effects
are the parameters of interest. In aster models, because components of the
response vector are dependent and have distributions in different families,
it is very unclear what REMLlike estimators in the style of \citet{b+c}
might be. The analogy just breaks down.
Hence, we do not pursue this REML analogy and stick with what we have
described above.
\section{Practice} \label{sec:reaster}
Our goal is to minimize \eqref{eq:logpickle}.
We replace \eqref{eq:logpickle} with \eqref{eq:logpickletoo}
in some steps because of our inability to differentiate \eqref{eq:logpickle},
but our whole procedure must minimize \eqref{eq:logpickle}.
\subsection{Step 1}
To get close to $(\hat{\alpha}, \hat{c}, \hat{\sigma})$ starting from far away
we minimize
\begin{equation} \label{eq:r}
\begin{split}
r(\sigma)
& =
 l(a + M \tilde{\alpha} + Z A \tilde{c})
+ \tfrac{1}{2} \tilde{c}^T \tilde{c}
\\
& \qquad
+ \tfrac{1}{2} \log
\det\bigl[ Z^T W(a + M \tilde{\alpha} + Z A \tilde{c}) Z A^2 + I \bigr]
\end{split}
\end{equation}
where $\tilde{\alpha}$ and $\tilde{c}$ are the joint minimizers of
\eqref{eq:keyrooted} considered as a function of $\alpha$ and $c$ for
fixed $\sigma$. In \eqref{eq:r}, $\tilde{\alpha}$, $\tilde{c}$, and $A$
are all functions of $\sigma$ although the notation does not indicate
this.
Because we cannot calculate derivatives of \eqref{eq:r}
we minimize using the R function \texttt{optim}
with \texttt{method = "NelderMead"}, the socalled NelderMead simplex
algorithm, a noderivative method nonlinear optimization,
not to be confused with the simplex algorithm
for linear programming.
\subsection{Step 2}
Having found $\alpha$, $c$, and $\sigma$ close to the MLE values
via the preceding step, we then switch to minimization of
\eqref{eq:keyrooted} for which we have the derivative formulas
\eqref{eq:keyrootedalpha}, \eqref{eq:keyrootedc},
and \eqref{eq:keyrootedsigma}.
In this step we can use one of R's optimization functions that uses
first derivative information: \texttt{nlm} or \texttt{nlminb} or
\texttt{optim} with optional argument \texttt{method = "BFGS"}
or \texttt{method = "CG"} or \texttt{method = "LBFGSB"}.
To define \eqref{eq:keyrooted} we also need a $\widehat{W}$, and we take the
value at the current values of $\alpha$, $c$, and $\sigma$.
Because $W$ is typically a very large matrix ($n \times n$, where $n$
is the number of nodes in complete aster graph, the number of nodes
in the subgraph for a single individual times the number of individuals),
we actually store $Z^T \widehat{W} Z$, which is only $r \times r$, where
$r$ is the number of random effects. We set
\begin{equation} \label{eq:zwz}
Z^T \widehat{W} Z
=
Z^T W(a + M \alpha + Z A c) Z
\end{equation}
where $\alpha$, $c$, and $A = A(\sigma)$ are the current values before
we start minimizing $\tilde{p}(\alpha, c, \sigma)$ and this value of
$Z^T \widehat{W} Z$ is fixed throughout the minimization,
as is required by the definition of $\tilde{p}(\alpha, c, \sigma)$.
Having minimized $\tilde{p}(\alpha, c, \sigma)$ we are still not done, because
now \eqref{eq:zwz} is wrong. We held it fixed at the values of
$\alpha$, $c$, and $\sigma$ we had before the minimization, and now
those values have changed. Thus we should reevaluate \eqref{eq:zwz}
and reminimize, and continue doing this until convergence.
When this iteration terminates we are done with this step, and we have
our point estimates $\hat{\alpha}$, $\hat{c}$, and $\hat{\sigma}$.
We also have our point estimates $\hat{b}$ of the random effects
on the original scale given by $A(\hat{\nu}) \hat{c}$
and our point estimates $\nu_j = \sigma_j^2$ of the variance components.
\subsection{Step 3}
Having converted back to the original parameters, if any of the $\nu_j$
are close to zero we use the check \eqref{eq:descenttest} to determine
whether or not they are exactly zero.
\section{R Package Aster}
We use the R statistical computing environment \citep{rcore} in our analysis.
It is free software and can be obtained from
\url{http://cran.rproject.org}. Precompiled binaries
are available for Windows, Macintosh, and popular Linux distributions.
We use the contributed package \verb@aster@ \citep{packageaster}.
If R has been installed, but this package has
not yet been installed, do
\begin{verbatim}
install.packages("aster")
\end{verbatim}
from the R command line
(or do the equivalent using the GUI menus if on Apple Macintosh
or Microsoft Windows). This may require root or administrator privileges.
Assuming the \verb@aster@ package has been installed, we load it
<>=
library(aster)
@
<>=
baz < library(help = "aster")
baz < baz$info[[1]]
baz < baz[grep("Version", baz)]
baz < sub("^Version: *", "", baz)
bazzer < paste(R.version$major, R.version$minor, sep = ".")
@
The version of the package used to make this document
is \Sexpr{baz}.
The version of R used to make this document is \Sexpr{bazzer}.
This entire document and all of the calculations shown were made using
the R command \texttt{Sweave} and hence are exactly reproducible by
anyone who has R and the R noweb (RNW) file from which it was created.
Both the RNW file and and the PDF document produced from it will be made
available at
\url{http://www.stat.umn.edu/geyer/aster}. For further details on
the use of Sweave and R see Chapter~1 of the technical report
\citet{aster2tr} available at the same web site.
Not only can one exactly reproduce the results in the printable document,
one can also modify the parameters of the simulation and get different
results. Anything at all can be changed once one has the RNW file.
In particular, we set the ``seed'' of the random number generator
<>=
set.seed(42)
@
so that every time this RNW file is run it produces the same results.
Changing the argument of \texttt{set.seed} or removing this chunk of
R code will produce different results.
<>=
# options(width = 70)
@
\section{A Digression on Aster Models and Formulas} \label{sec:formula}
\subsection{Observed Fitness}
In an unconditional aster model (and all published examples in the literature
are unconditional aster models except for Example~1 of \citet{aster2} and that
could have also been done using an unconditional aster model) the unconditional
canonical parameter vector $\varphi$ has a multivariate monotone relationship
with the unconditional mean value parameter vector $\mu$
\citep[Appendix]{asterfitness}. The exact relationship between
$\varphi$ and $\mu$ is very complicated.
\citet{asterphilosophic} works through a simple example, and the formulas
become very messy. So the only thing one can have any intuition about is
the multivariate monotone relationship.
The new functions added to the \texttt{aster} package to do random effects
aster models only do unconditional aster models.
Now what unconditional mean values does one want to establish relationships
with? Generally with those that are the best surrogates of overall fitness
(total reproductive success of individuals over their lifetime)
in the data, that is, generally, the last fitness component observed.
In the data in the example in \citet{gws}, that is head count (number of
compound flowers observed). One might think that the earlier components
of fitness are important too, but their effect is already incorporated in
the last component of fitness (you can't have flowers if you are dead,
and similarly for any other component of fitness).
Thus if one has predictors only affecting ``fitness'' nodes of the graph
(those whose sum is the best surrogate of lifetime fitness), an unconditional
aster model will do the right thing by adjusting the unconditional mean values
of those nodes to fit the data. In recent papers we have included an indicator
variable named \texttt{fit} that indicates these ``fitness'' nodes of the graph
(it is one for those nodes and zero for other nodes) that helps in the modeling
and we have done so for the examples in this paper.
The model in \citet{gws} that was deemed the best one for drawing scientific
conclusions (their Model~2) was fit by the formula given in the paper
\begin{verbatim}
resp ~ varb + level:(nsloc + ewloc) + hdct * pop  pop
\end{verbatim}
but that formula can be simplified to
\begin{verbatim}
resp ~ varb + level:(nsloc + ewloc) + hdct:pop
\end{verbatim}
which fits the same model and is easier to understand.
The term \texttt{hdct:pop} is best not thought of as an ``interaction''
although that is the way the R formula minilanguage describes it.
What actually happens is, because \texttt{pop} is categorical with 7 levels
R makes 6 dummy variables (one for each category but throws one away because
all 7 together would be confounded with the intercept dummy variable) then
it multiplies each of these dummy variables by \texttt{hdct} (which, recall,
is zerooronevalued), and this does exactly what is wanted: making those
dummy variables apply to ``fitness'' nodes only.
(It was just noticed that the examples on the help pages for the R functions
\texttt{aster}, \texttt{anova.aster}, and \texttt{predict.aster} had the
oldstyle formulas. Those help pages have now been fixed.)
The only difference between the example just discussed (which now matches
the example on the help page for the \texttt{aster} function) and the ones
in this technical report (other than being random effects models) is that
where it says \texttt{hdct} we would now say \texttt{fit}, taking that for
a conventional name of a dummy (indicator) variable that indicates ``fitness''
nodes of the graph.
Thus the dictum: every variable for which one wants to establish a relationship
with (overall) fitness should enter every formula ``interacted with''
\texttt{fit} (but, as explained above, ``interacted with''
is a bad description, hence the scare quotes). And that ``interacted with''
must be the colon operator (\texttt{:}) in the R formula minilanguage,
not the star operator (\texttt{*}), as in \texttt{pop : fit}.
\subsection{Other Fitness Components}
In general it makes no scientific sense to have terms without interaction
in aster model formulas, except for \texttt{varb} if that is what one is
calling the factor that indicates nodes of the graph (as it does in all of
the examples in this technical report and earlier aster technical reports
and in the aster papers these technical reports accompanied).
The reason why one should not, for example, have \texttt{pop} by itself
is that it makes no sense to have one parameter for the population effect
on all fitness components (survival and fecundity, or survival, flowering
indicator, and number of flowers). That is why the four models tested by
\citet[Table~1]{gws}, which can be simplified (as discussed above) to
\begin{center}
\begin{tabular}{ll}
Model 1 & \verb@resp ~ varb + level:(nsloc + ewloc)@ \\
Model 2 & \verb@resp ~ varb + level:(nsloc + ewloc) + fit:pop@ \\
Model 3 & \verb@resp ~ varb + level:(nsloc + ewloc) + factor(fit)*pop@ \\
Model 4 & \verb@resp ~ varb + level:(nsloc + ewloc) + level*pop@
\end{tabular}
\end{center}
have some ``interaction'' (again with scare quotes because they should
not be interpreted as interactions) with either \texttt{fit} or \texttt{level}
or \texttt{varb}, all of which are indicators for certain nodes of the graph
or groups of nodes of the graph. These ``interactions'' make certain
regression coefficients only apply to certain nodes of the graph.
In model 1 there are no \texttt{pop} effects.
In model 2 there are \texttt{pop} effects only for fitness nodes (the head
count nodes that together are the best surrogate of observed fitness).
In model 3 there are \texttt{pop} effects not only for fitness node
but also for nonfitness nodes. We would now say this model doesn't really
make scientific sense and \citet{gws} said the same thing (``it is difficult
to interpret Model 3 scientifically $\ldots$'') because nonfitness nodes
includes two different fitness components (survival and flowering indicator)
and these should have separate parameters not the same parameters. Thus
we should not use Model 3 \emph{even though it fits best according to purely
statistical criteria}.
In model 4 there are \texttt{pop} effects for each population and each
``level'' (which is shorthand for component of fitness), that is, there
are separate \texttt{pop} effects for survival, for flowering, and for head
count. And that also makes scientific sense.
We could also have a model 5 in which the last term in the formula is
\verb@varb*pop@ but that would be a lot of parameters.
Whenever possible, one wants to have \texttt{fit} interaction with all
predictors with which one wants to establish a relation with fitness.
As \citet{gws} discuss, models other than Model~2 are difficult to interpret.
The fact that \texttt{fit:pop} is statistically significant has the direct
interpretation that individuals having different parental populations have
different fitness. Model~4, in contrast, says that all fitness components
vary with respect to parental population, not necessarily in the same
direction, and the effect on overall fitness is unclear.
\section{Radish Data}
We use data on the invasive California wild radish (\emph{Raphanus sativus})
described by \citet{ridley} and contained in the dataset \texttt{radish}
in the R contributed package \texttt{aster}. For each individual, three
response variables are observed, connected by the following graphical model
$$
\begin{CD}
1
@>\text{Ber}>>
y_1
@>\text{0Poi}>>
y_2
@>\text{Poi}>>
y_3
\end{CD}
$$
$y_1$ being an indicator of whether any flowers were produced,
$y_2$ being the count of the number of flowers produced,
$y_3$ being the count of the number of fruits produced,
the unconditional distribution of $y_1$ being Bernoulli,
the conditional distribution of $y_2$ given $y_1$ being zerotruncated Poisson,
and the conditional distribution of $y_3$ given $y_2$ being Poisson.
We load the data
<>=
data(radish)
names(radish)
levels(radish$varb)
@
This is a ``long format'' data frame produced by the R command reshape
from ``wide format'' data. The variable \texttt{varb} indicates which
components of the response vector (variable \texttt{resp}) corresponded
to which original variables in the ``wide format'' data (components of
fitness). The variable \texttt{fit} is the indicator of the best surrogate
of fitness in these data, which is the last node ($y_3$) of the graph and
the \texttt{"Fruits"} level of \texttt{varb}.
Then set up the graphical model
<>=
pred < c(0,1,2)
fam < c(1,3,2)
sapply(fam.default(), as.character)[fam]
@
These data come from a designed experiment started with seeds collected from
three large wild populations of northern, coastal California wild radish
and three populations of southern, inland California wild radish. Thus we
have populations nested within region.
<>=
reg2pop < split(as.character(radish$Pop), as.character(radish$Region))
reg2pop < lapply(reg2pop, unique)
reg2pop
@
Plants were grown at two experimental sites,
one northern, coastal California field site located
at Point Reyes National Seashore and one
southern, inland site located
at the University of California Riverside Agricultural Experiment Station.
Blocks were nested within site.
<>=
sit2blk < split(as.character(radish$Block), as.character(radish$Site))
sit2blk < lapply(sit2blk, unique)
sit2blk
@
The issue of main scientific interest is the interaction of region and
site, which is indicative of local adaptation when the pattern of mean
values shows that each population does better in its home environment than
in others. Testing significance of this
interaction is complicated by the nesting of populations within region and
blocks within site and the desire of scientists to treat these nested factors
as random effects.
The best surrogate of fitness in these data is the \texttt{Fruits}
component of the response vector. Thus we form the ``interaction'' with
the indicator of this component and all scientifically
interesting predictors (Section~\ref{sec:formula} above).
\subsection{A Fixed Effects Model}
The fixed effects model most closely connected with the random effects model
of interest is
<>=
aout < aster(resp ~ varb + fit : (Site * Region + Block + Pop),
pred, fam, varb, id, root, data = radish)
options(show.signif.stars = FALSE)
summary(aout)
@
<>=
options(width = 70)
@
Note: the variable \texttt{fit} is the same as the dummy variable
\texttt{varbFruits} constructed by the aster model software.
The parameter of main scientific interest is the regression coefficient
named (by R) \texttt{\Sexpr{names(aout$coef)[length(aout$coef)]}}, which is the
region by site interaction. A statistically significantly nonzero value
of this parameter may indicate local adaptation
(more on this in Section~\ref{sec:reanalradish} below).
In the fit above, this parameter is indeed highly statistically significant,
but that $P$value does not take the random effects story into account.
\subsection{Random Effects Model}
The traditional way to deal with a situation like this is to treat the
population effects as random (within region) and the block effects as
random (within site). We now do (an approximation to) this.
To specify a randomeffects aster model using the R function \texttt{reaster}
one supplies either two formulas (when there is only one variance component)
or one formula and a list of formulas (when there is more than one variance
component). The first formula specifies the fixed effects model matrix
($M$ in the notation of Section~\ref{sec:theo}), and the second formula or
list of formulas specifies the random effects model matrix
($Z$ in the notation of Section~\ref{sec:theo}).
When there is a list of formulas, each formula specifies the columns of the
random effects model matrix that go with one variance component.
The components of the list should be named, because the names are taken to name
the variance components.
The first formula (for fixed effects) is just like in an ordinary aster model
(or a linear or generalized linear model). The other formulas (for random
effects) are somewhat different in that they (1) do not need a response
(that is specified by the fixed effects formula) and (2) should not have
an intercept (the way to specify this is to prefix the formula with
\texttt{0 +}). Hence the following. Note that we are following the dictum
of Section~\ref{sec:formula}.
<>=
rout < reaster(resp ~ varb + fit : (Site * Region),
list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop),
pred, fam, varb, id, root, data = radish)
summary(rout)
@
The results are somewhat different from the fixed effects analysis.
One fixed effect \texttt{fit:SitePoint Reyes},
which was statistically significant in the fixed effects model,
is not statistically
significant in the random effects model. The main scientific conclusions,
however, do not change (Section~\ref{sec:reanalradish} below).
We also try out some options of the \texttt{summary.reaster} function,
not because they are particularly interesting here, but just to illustrate
them
<>=
summary(rout, standard.deviation = FALSE)
@
Now the estimates in the \texttt{Variance Components} section of the
printout are variance components $\hat{\nu}_k = \hat{\sigma}_k^2$
rather than their square roots as we had before.
We do not recommend this option
\texttt{standard.deviation = FALSE} because the standard errors,
derived from the delta method go to zero as $\nu_k$ goes to zero,
and this is wrong (it is right ``in asymptopia'' but for sufficiently
small $\nu_k$ whatever sample size one has is too small).
Thus this option seems to provide worse estimates than the default.
\subsection{Does this Reanalysis Change Biological Conclusions?}
\label{sec:reanalradish}
In short, no. \citet{ridley} did not do a random effects aster analysis
because it had not yet been invented. Nevertheless their conclusions hold up.
The main conclusion of interest being local adaptation, which is indicated by
the statistical significance of the fixed effect for regionsite interaction
and the pattern of mean values for different populations and different blocks,
which we do not examine \citep[this was done by][]{ridley}.
The fact that random effects analysis and fixed effects analysis agree
qualitatively on this one example does not, of course, prove that they
would agree on all examples. In these data the regionsite interaction
is very large and almost any sensible statistical analysis would show it.
When the interaction is not so large, the analysis done will make a difference.
\section{Bootstrapping the Radish Data}
In this section we do a parametric bootstrap to check the standard errors
provided by \texttt{reaster} and \texttt{summary}. First we store the
(approximate) maximum likelihood estimates
<>=
names(rout)
alpha.hat < rout$alpha
sigma.hat < rout$sigma
nu.hat < rout$nu
b.hat < rout$b
c.hat < rout$c
sout < summary(rout)
se.alpha.hat < sout$alpha[ , "Std. Error"]
se.sigma.hat < sout$sigma[ , "Std. Error"]
se.nu.hat < sout$nu[ , "Std. Error"]
fixed < rout$fixed
random < rout$random
modmat.tot < cbind(fixed, Reduce(cbind, random))
nfix < ncol(fixed)
nrand < sapply(random, ncol)
a.hat < rep(sigma.hat, times = nrand)
@
To simulate new aster data we first need to change from unconditional
canonical parameters to conditional canonical parameters (because that's
what the R function \texttt{raster} requires).
<>=
c.star < rnorm(sum(nrand))
b.star < a.hat * c.star
eff.star < c(alpha.hat, b.star)
phi.star < as.numeric(as.vector(rout$obj$origin) + modmat.tot %*% eff.star)
theta.star < astertransform(phi.star, rout$obj, to.cond = "conditional",
to.mean = "canonical")
y.star < raster(theta.star, pred, fam, rout$obj$root)
y.star < as.vector(y.star)
@
Now we need to redo the above analysis on the new data. We can take
the simulation truth as starting values.
<>=
rout.star < reaster(y.star ~ varb + fit : (Site * Region),
list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop),
pred, fam, varb, id, root, data = radish,
effects = c(alpha.hat, c.star), sigma = sigma.hat)
sout.star < summary(rout.star)
print(sout.star)
@
Now we want to get estimates and standard errors for this fit
<>=
alpha.star < rout.star$alpha
sigma.star < rout.star$sigma
nu.star < rout.star$nu
as.vector(alpha.star  alpha.hat)
as.vector(sigma.star  sigma.hat)
as.vector(nu.star  nu.hat)
se.alpha.star < sout.star$alpha[ , "Std. Error"]
se.sigma.star < sout.star$sigma[ , "Std. Error"]
se.nu.star < sout.star$nu[ , "Std. Error"]
@
So this is the analysis we bootstrap. All we need to do is put it in a loop.
The bootstrap takes a long time (see below) if done with a reasonable
sample size. Here we use bootstrap sample size
<>=
nboot < 199
@
Thus we save the results and restore them here
<>=
suppressWarnings(foo < try(load("radishboot.rda"), silent = TRUE))
done < (! inherits(foo, "tryerror")) &&
identical(.Random.seed, signature$seed) &&
identical(nboot, signature$nboot) &&
identical(alpha.hat, signature$alpha) &&
identical(sigma.hat, signature$sigma)
done
@
If \texttt{done} is \texttt{TRUE} then the bootstrap is already done
and its results restored. Otherwise we have to do it, either because it
has never been done or because one or more of the values of the R objects
that (partially) determine it has changed. The reason for using both
the \texttt{suppressWarnings} function and the \texttt{silent = TRUE}
argument to the \texttt{try} function is that when the file we are
trying to load does not exist the \texttt{load} function gives both
an error and a warning.
So now we are ready to try it out.
<>=
if (! done) {
signature < list(seed = .Random.seed, nboot = nboot, alpha = alpha.hat,
sigma = sigma.hat)
boot.start.time < proc.time()
alpha.star < matrix(NaN, nboot, length(alpha.hat))
sigma.star < matrix(NaN, nboot, length(sigma.hat))
nu.star < matrix(NaN, nboot, length(nu.hat))
se.alpha.star < alpha.star
se.sigma.star < sigma.star
se.nu.star < nu.star
for (iboot in 1:nboot) {
c.star < rnorm(sum(nrand))
b.star < a.hat * c.star
eff.star < c(alpha.hat, b.star)
phi.star < as.numeric(as.vector(rout$obj$origin) +
modmat.tot %*% eff.star)
theta.star < astertransform(phi.star, rout$obj,
to.cond = "conditional", to.mean = "canonical")
y.star < raster(theta.star, pred, fam, rout$obj$root)
y.star < as.vector(y.star)
rout.star < reaster(y.star ~ varb + fit : (Site * Region),
list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop),
pred, fam, varb, id, root, data = radish,
effects = c(alpha.hat, c.star), sigma = sigma.hat)
sout.star < suppressWarnings(summary(rout.star))
alpha.star[iboot, ] < rout.star$alpha
sigma.star[iboot, ] < rout.star$sigma
nu.star[iboot, ] < rout.star$nu
se.alpha.star[iboot, ] < sout.star$alpha[ , "Std. Error"]
se.sigma.star[iboot, ] < sout.star$sigma[ , "Std. Error"]
se.nu.star[iboot, ] < sout.star$nu[ , "Std. Error"]
}
boot.stop.time < proc.time()
save.random.seed < .Random.seed
save(signature, alpha.star, sigma.star, nu.star, se.alpha.star,
se.sigma.star, se.nu.star,
boot.start.time, boot.stop.time, save.random.seed,
file = "radishboot.rda")
} else {
.Random.seed < save.random.seed
}
@
<>=
ttime < (boot.stop.time  boot.start.time)[1]
tsec < round(ttime %% 60, 1)
tmin < floor((ttime / 60) %% 60)
thr < floor(ttime / (60 * 60))
@
The bootstrap with bootstrap sample size \Sexpr{nboot} took
\Sexpr{thr} hours, \Sexpr{tmin} minutes, and \Sexpr{tsec} seconds.
Occasionally, the \texttt{summary} function fails to calculate standard errors.
The following tells us how many times this happened.
<>=
sum(! is.finite(se.alpha.star[ , 1]))
@
Now we make standardized quantities
<>=
z < cbind(alpha.star, sigma.star)
z < sweep(z, 2, c(alpha.hat, sigma.hat))
se.z < cbind(se.alpha.star, se.sigma.star)
z < z / se.z
apply(z, 2, mean)
apply(z, 2, sd)
@
Not exactly standard normal (mean zero, standard deviation one).
In fact, some seem quite far from standard normal.
If we apply more robust estimators of location and scale
<>=
apply(z, 2, median)
apply(z, 2, mad)
@
we see that some of the nonnormality appears to be due to outliers,
but the distributions are still nowhere near standard normal.
This means if we really want accurate tests and confidence intervals,
they cannot be based on the standard errors printed out
by \texttt{summary.reaster}.
Suppose we want a 95\% confidence interval for the coefficient for
the fixed effect named \texttt{"fit:SiteRiverside:RegionS"}.
Instead of using
$$
\text{point estimate} \pm 1.96 \times \text{standard error}
$$
we should use critical values derived from the bootstrap distribution.
Here's how to do that.
<>=
conf.level < 0.95
n < "fit:SiteRiverside:RegionS"
colnames(z) < c(names(alpha.hat), names(sigma.hat))
myz < z[ , n]
crit < quantile(myz, probs = c((1  conf.level) / 2, (1 + conf.level) / 2))
crit
alpha.hat[n]  rev(crit) * se.alpha.hat[n]
@
In this case, the bootstrap critical values are actually smaller than
1.96 so the bootstrap confidence interval is actually shorter than
the confidence interval based on the standard error printed out by
\texttt{summary.reaster}. But that won't be true in general.
Let us look at what seems to be the worst (most nonstandardnormal)
bootstrap distribution, that for $\nu_1$.
<>=
n < "block"
myz < z[ , n]
myz < myz[! is.na(myz)]
@
makes the bootstrap $z$scores and
<