% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.




\setkeys{Gin}{width=\textwidth}

\begin{frame}\frametitle{Generalized Linear Mixed Models}
  \begin{itemize}
  \item When using linear mixed models (LMMs) we assume that the
    response being modeled is on a continuous scale.
  \item Sometimes we can bend this assumption a bit if the response is
    an ordinal response with a moderate to large number of levels.
    For example, the Scottish secondary school test results were
    integer values on the scale of 1 to 10.
  \item However, an LMM is not suitable for modeling a binary
    response, an ordinal response with few levels or a response that
    represents a count. For these we use generalized linear mixed
    models (GLMMs).
  \item To describe GLMMs we return to the representation of the
    response as an $n$-dimensional, vector-valued, random variable,
    $\bm{\mathcal{Y}}$, and the random effects as a $q$-dimensional,
    vector-valued, random variable, $\bm{\mathcal{B}}$.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Parts of LMMs carried over to GLMMs}
  \begin{itemize}
  \item Random variables
    \begin{compactitem}
    \item $\bm{\mathcal{Y}}$ the response variable
    \item $\bm{\mathcal{B}}$ the (possibly correlated) random effects
    \item $\bm{\mathcal{U}}$ the orthogonal random effects
    \end{compactitem}
  \item Parameters
    \begin{compactitem}
    \item $\bm{\beta}$ - fixed-effects coefficients
    \item $\sigma$ - the common scale parameter (not always used)
    \item $\bm{\theta}$ - parameters that determine
      $\mathrm{Var}(\bm{\mathcal{B}})=\sigma^2(\bm T\bm S)(\bm T\bm S)\trans$ 
    \end{compactitem}
  \item Some matrices
     \begin{compactitem}
     \item $\bm{X}$ the $n\times p$ model matrix for $\bm{\beta}$
     \item $\bm{Z}$ the $n\times q$ model matrix for $\bm b$
     \item $\bm{P}$ fill-reducing $q\times q$ permutation (from $\bm Z$)
     \item $\bm{S}(\bm\theta)$ non-negative $q\times q$ diagonal scale matrix
     \item $\bm{T}(\bm\theta)$ $q\times q$ unit lower-triangular matrix
     \item $\bm{A}(\bm\theta)=(\bm Z\bm T(\bm\theta)\bm S(\bm\theta))\trans$
     \end{compactitem}
   \end{itemize}
 \end{frame}

\begin{frame}\frametitle{The conditional distribution,
    $\bm{\mathcal{Y}}|\bm{\mathcal{U}}$}
  \begin{itemize}
  \item For GLMMs, the marginal distribution,
    $\bm{\mathcal{B}}\sim\mathcal{N}\left(\bm
      0,\bm\Sigma(\bm\theta)\right)$ is the same as in LMMs except
    that $\sigma^2$ is omitted.  We define
    $\bm{\mathcal{U}}\sim(\bm 0,\bm I_q)$ such that
    $\bm{\mathcal{B}}=\bm T(\bm\theta)\bm S(\bm\theta) \bm
    P\trans\bm{\mathcal{U}}$.
  \item For GLMMs we retain some of the properties of the conditional
    distribution
    \begin{displaymath}
      (\bm{\mathcal{Y}}|\bm{\mathcal{U}}=\bm u)\sim
      \mathcal{N}\left(\bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}},\sigma^2\bm
        I\right)\text{ where }
      \bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}(\bm u)=
      \bm X\bm\beta+\bm A\trans\bm P\trans\bm u
    \end{displaymath}
    Specifically
    \begin{itemize}
    \item The distribution $\bm{\mathcal{Y}}|\bm{\mathcal{U}=\bm u}$ depends
      on $\bm u$ only through the conditional mean,
      $\bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}(\bm u)$.
    \item Elements of $\bm{\mathcal{Y}}$ are \Emph{conditionally
        independent}.  That is, the distribution of
      $\bm{\mathcal{Y}}|\bm{\mathcal{U}=\bm u}$  is completely
      specified by the univariate, conditional distributions,
      $\mathcal{Y}_i|\bm{\mathcal{U}},i=1,\dots,n$.
    \item These univariate, conditional distributions all have the same
      form. They differ only in their means.
    \end{itemize}
  \item GLMMs differ from LMMs in the form of the univariate,
    conditional distributions and in how
    $\bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}(\bm u)$ depends on
    $\bm u$.
  \end{itemize}    
\end{frame}

\begin{frame}
  \frametitle{Some choices of univariate conditional distributions}
  \begin{itemize}
  \item Typical choices of univariate conditional distributions are:
    \begin{itemize}
    \item The \Emph{Bernoulli} distribution for binary (0/1) data, which has
      probability mass function 
      \begin{displaymath}
        p(y|\mu)= \mu^{y}(1-\mu)^{1 - y},\quad 0<\mu< 1,\quad y = 0,1
      \end{displaymath}
    \item Several independent binary responses can be represented as a
      \Emph{binomial} response, but only if all the Bernoulli
      distributions have the same mean.
    \item The \Emph{Poisson} distribution for count ($0,1,\dots$)
      data, which has probability mass function
      \begin{displaymath}
        p(y|\mu)= e^{-y}\frac{\mu^{y}}{y!},\quad 0<\mu,\quad y = 0,1,2,\dots
      \end{displaymath}
    \end{itemize}
  \item All of these distributions are completely specified by the
    conditional mean.  This is different from the conditional normal
    (or Gaussian) distribution, which also requires the common scale
    parameter, $\sigma$.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{The link function, g}
  \begin{itemize}
  \item When the univariate conditional distributions have constraints
    on $\mu$, such as $0<\mu<1$ (Bernoulli) or $0<\mu$ (Poisson), we
    cannot define the conditional mean,
    $\bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}$, to be equal to the
    linear predictor, $\bm X\bm\beta+\bm A\trans\bm P\trans\bm u$,
    which is unbounded.
  \item We choose an invertible, univariate \Emph{link function}, $g$,
    such that $\eta=g(\mu)$ is unconstrained.  The vector-valued link
    function, $\bm g$, is defined by applying $g$ component-wise.
    \begin{displaymath}
      \bm\eta=\bm g(\bm\mu)\quad\mathrm{where}\quad
      \eta_i=g(\mu_i),\quad i=1,\dots,n
    \end{displaymath}
  \item We require that $g$ be invertible so that $\mu=g^{-1}(\eta)$
    is defined for $-\infty<\eta<\infty$ and is in the appropriate
    range ($0<\mu<1$ for the Bernoulli or $0<\mu$ for the
    Poisson).  The vector-valued inverse link, $\bm g^{-1}$, is defined
    component-wise.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{``Natural'' link functions}
  \begin{itemize}
  \item There are many choices of invertible scalar link functions, $g$,
    that we could use for a given set of constraints.
  \item For the Bernoulli and Poisson distributions, however, one link
    function arises naturally from the definition of the probability
    mass function.  (The same is true for a few other, related but
    less frequently used, distributions, such as the gamma distribution.)
  \item To derive the natural link, we consider the logarithm of the
    probability mass function (or, for continuous distributions, the
    probability density function).
  \item For distributions in this ``exponential'' family, the
    logarithm of the probability mass or density can be written as a
    sum of terms, some of which depend on the response, $y$, only and
    some of which depend on the mean, $\mu$, only.  However, only one
    term depends on \textbf{both} $y$ and $\mu$, and this term has the
    form $y\cdot g(\mu)$, where $g$ is the natural link.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{The natural link for the Bernoulli distribution}
  \begin{itemize}
  \item The logarithm of the probability mass function is
    \begin{displaymath}
      \log(p(y|\mu))=\log(1-\mu)+y\log\left(\frac{\mu}{1-\mu}\right),
      \;0<\mu<1,\;y=0,1 .
    \end{displaymath}
  \item Thus, the natural link function is the
    \Emph{logit} link
    \begin{displaymath}
      \eta=g(\mu)=\log\left(\frac{\mu}{1-\mu}\right).
    \end{displaymath}
  \item Because $\mu=P[\mathcal{Y} = 1]$, the quantity
    $\mu/(1-\mu)$ is the odds ratio (in the range $(0, \infty)$)
    and $g$ is the logarithm of the odds ratio, sometimes called ``log
    odds''.
  \item The inverse link is
    \begin{displaymath}
      \mu=g^{-1}(\eta)=\frac{e^\eta}{1+e^\eta}=\frac{1}{1+e^{-\eta}}
    \end{displaymath}
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Plot of natural link for the Bernoulli distribution}
\includegraphics{figs/Modeldef-BernoulliLink}
\end{frame}

\begin{frame}
  \frametitle{Plot of inverse natural link for the Bernoulli distribution}
\includegraphics{figs/Modeldef-BernoulliinvLink}
\end{frame}

\begin{frame}
  \frametitle{The natural link for the Poisson distribution}
  \begin{itemize}
  \item The logarithm of the probability mass is
    \begin{displaymath}
      \log(p(y|\mu))=\log(y!)-\mu+y\log(\mu)
    \end{displaymath}
  \item Thus, the natural link function for the Poisson is the
    \Emph{log} link
    \begin{displaymath}
      \eta=g(\mu)=\log(\mu)
    \end{displaymath}
  \item The inverse link is
    \begin{displaymath}
      \mu=g^{-1}(\eta)=e^{\eta}
    \end{displaymath}
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{The natural link related to the variance}
  \begin{itemize}
  \item For the natural link function, the derivative of its inverse
    is the variance of the response.
  \item For the Bernoulli, the natural link is the logit and the
    inverse link is $\mu=g^{-1}(\eta)=1/(1+e^{-\eta})$. Then
    \begin{displaymath}
      \frac{d\mu}{d\eta}=\frac{e^{-\eta}}{(1+e^{-\eta})^2}=
      \frac{1}{1+e^{-\eta}}\frac{e^{-\eta}}{1+e^{-\eta}}=\mu(1-\mu)=
      \mathrm{Var}(\mathcal{Y})
    \end{displaymath}
  \item For the Poisson, the natural link is the log and the inverse
    link is $\mu=g^{-1}(\eta)=e^\eta$.  Then
    \begin{displaymath}
      \frac{d\mu}{d\eta}=e^\eta=\mu=\mathrm{Var}(\mathcal{Y})
    \end{displaymath}
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{The unscaled conditional density of
    $\bm{\mathcal{U}}|\bm{\mathcal{Y}}=\bm y$}
  \begin{itemize}
  \item As in LMMs we evaluate the likelihood of the parameters, given
    the data, as
    \begin{displaymath}
      L(\bm\theta,\bm\beta|\bm y)=
      \int_{\mathbb{R}^q}[\bm{\mathcal{Y}}|\bm{\mathcal{U}}](\bm y|\bm u)
      \,[\bm{\mathcal{U}}](\bm u)\, d\bm u ,
    \end{displaymath}
  \item The product $[\bm{\mathcal{Y}}|\bm{\mathcal{U}}](\bm y|\bm u)
    [\bm{\mathcal{U}}](\bm u)$ is the unscaled (or
    \Emph{unnormalized}) density of the conditional distribution
    $\bm{\mathcal{U}}|\bm{\mathcal{Y}}$.
  \item The density $[\bm{\mathcal{U}}](\bm u)$ is a spherical
    Gaussian density $\frac{1}{(2\pi)^{q/2}} e^{-\|\bm u\|^2/2}$.
  \item The expression $[\bm{\mathcal{Y}}|\bm{\mathcal{U}}](\bm y|\bm
    u)$ is the value of a probability mass function or a probability
    density function, depending on whether
    $\mathcal{Y}_i|\bm{\mathcal{U}}$ is discrete or continuous.
  \item The linear predictor is
    $\bm g(\bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}})=\bm\eta=\bm
    X\bm\beta+\bm A(\bm\theta)\trans\bm P\trans\bm u$.
    Alternatively, we can write the conditional mean of
    $\bm{\mathcal{Y}}$, given $\bm{\mathcal{U}}$, as
    \begin{displaymath}
      \bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}(\bm u)=
      \bm g^{-1}\left(\bm X\bm\beta+
        \bm A(\bm\theta)\trans\bm P\trans\bm u\right)
    \end{displaymath}
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{The conditional mode of
    $\bm{\mathcal{U}}|\bm{\mathcal{Y}}=\bm y$}
  \begin{itemize}
  \item In general the likelihood, $L(\bm\theta,\bm\beta|\bm y)$ does
    not have a closed form. To approximate this value, we first determine
    the \Emph{conditional mode}
    \begin{displaymath}
      \tilde{\bm u}(\bm y|\bm\theta,\bm\beta)=\arg\max_{\bm u}
      [\bm{\mathcal{Y}}|\bm{\mathcal{U}}](\bm y|\bm u)\, 
      [\bm{\mathcal{U}}](\bm u)
    \end{displaymath}
    using a quadratic approximation to the logarithm of the
    unscaled conditional density.
  \item This optimization problem is (relatively) easy because the
    quadratic approximation to the logarithm of the unscaled
    conditional density can be written as a penalized, weighted
    residual sum of squares,
    \begin{displaymath}
      \tilde{\bm u}(\bm y|\bm\theta,\bm\beta)=\arg\min_{\bm u}\left\|
        \begin{bmatrix}
          \bm W^{1/2}(\bm\mu)\left(\bm y -
            \bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}(\bm u)\right)\\
          -\bm u
        \end{bmatrix}\right\|^2
    \end{displaymath}
    where $\bm W(\bm\mu)$ is the diagonal weights matrix.  The weights
    are the inverses of the variances of the $\mathcal{Y}_i$.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{The PIRLS algorithm}
  \begin{itemize}
  \item Parameter estimates for generalized linear models (without
    random effects) are usually determined by iteratively reweighted
    least squares (IRLS), an incredibly efficient algorithm.  PIRLS is
    the penalized version.  It is iteratively reweighted in the
    sense that parameter estimates are determined for a fixed weights
    matrix $\bm W$ then the weights are updated to the current
    estimates and the process
    repeated.
  \item For fixed weights we solve
    \begin{displaymath}
      \min_{\bm u}\left\|
        \begin{bmatrix}
          \bm W^{1/2}\left(\bm y -
            \bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}(\bm u)\right)\\
          -\bm u
        \end{bmatrix}\right\|^2
    \end{displaymath}
    as a nonlinear least squares problem with update, $\bm\delta_{\bm
      u}$, given by
    \begin{displaymath}
      \bm P\left(\bm A\bm M\bm W\bm M\bm A\trans+\bm I\right)P\trans
      \bm\delta_{\bm u}=\bm P\bm A\bm M\bm W(\bm y-\bm\mu) - \bm u
    \end{displaymath}
    where $\bm M=d\bm\mu/d\bm\eta$ is the (diagonal) Jacobian matrix.
    Recall that for the natural link, $\bm M=
    \mathrm{Var}(\bm{\mathcal{Y}}|\bm{\mathcal{U}})=\bm W^{-1}$.
  \end{itemize}
\end{frame}
\begin{frame}
  \frametitle{The Laplace approximation to the deviance}
  \begin{itemize}
  \item At convergence, the sparse Cholesky factor, $\bm L$, used to
    evaluate the update is
    \begin{displaymath}
      \bm L\bm L\trans =
      \bm P\left(\bm A\bm M\bm W\bm M\bm A\trans+\bm I\right)\bm P\trans
    \end{displaymath}
    or
    \begin{displaymath}
      \bm L\bm L\trans =
      \bm P\left(\bm A\bm M\bm A\trans+\bm I\right)\bm P\trans
    \end{displaymath}
    if we are using the natural link.
  \item The integrand of the likelihood is approximately a constant
    times the density of the $\mathcal{N}(\tilde{\bm u},\bm L\bm
    L\trans)$ distribution.
  \item On the deviance scale (negative twice the log-likelihood) this
    corresponds to
    \begin{displaymath}
      d(\bm\beta,\bm\theta|\bm y)=d_g(\bm y,\bm\mu(\tilde{\bm u}))+
      \|\tilde{\bm u}\|^2+\log(|\bm L|^2)
    \end{displaymath}
    where $d_g(\bm y,\bm\mu(\tilde{\bm u}))$ is the GLM deviance for
    $\bm y$ and $\bm\mu$.
  \end{itemize}
\end{frame}
% \begin{frame}
%   \frametitle{The Laplace approximation to the deviance}
%   \begin{itemize}
%   \item The variance-covariance matrix
%     of the conditional distribution $\bm U|\bm Y=\bm y$ is
%     approximated by
%     \begin{displaymath}
%       \bm L\bm L\trans = \bm P\left(\bm A(\bm\theta)\frac{d\bm\mu}{d\bm\eta}
%       \bm W(\bm\mu) \frac{d\bm\mu}{d\bm\eta}
%       \bm A(\bm\theta)\trans + \bm I_q\right)\bm P\trans
%     \end{displaymath}
%     where $\bm W(\bm\mu)$ and $\frac{d\bm\mu}{d\bm\eta}$ are evaluated
%     at the conditional mode, $\tilde{\bm u}(\bm y|\bm\theta,\bm\beta)$,
%     and $\bm L$ is the (sparse) Cholesky factor.
%   \item Another characteristic of the natural link functions is that
%     the derivative
%     $\frac{d\bm\mu}{d\bm\eta}=\mathrm{Var}(\bm\mathcal{Y}|\bm\mathcal{B})$.
%     Because the weight matrix, $\bm W$, is the inverse of the
%     variance-covariance matrix, this leads to a stable algorithm that
%     converges quickly.
%   \item The Laplace approximation to the deviance is exactly analogous
%     to the evaluation of the deviance for linear mixed models
%     \begin{displaymath}  %% check this
%       d(\bm\theta,\bm\beta|\bm y)=\log(|\bm L|^2)+
%       h(\tilde{u}(\bm\beta,\bm\theta))
%     \end{displaymath}
%     where $h(\bm u)$ is the unscaled conditional density
%     $\bm{\mathcal{U}}|\bm{\mathcal{Y}}=\bm y$.
%   \end{itemize}
% \end{frame}

\begin{frame}
  \frametitle{Modifications to the algorithm}
  \begin{itemize}
  \item Notice that this deviance depends on the fixed-effects
    parameters, $\bm\beta$, as well as the variance-component
    parameters, $\bm\theta$.  This is because $\log(|\bm L|^2)$
    depends on $\bm\mu_{\bm{\mathcal{Y}}|\bm{\mathcal{U}}}$ and,
    hence, on $\bm\beta$.  For LMMs $\log(|\bm L|^2)$ depends only on
    $\bm\theta$.
  \item It is likely that modifying the PIRLS algorithm to optimize
    simultaneously on $\bm u$ and $\bm\beta$ would result in a value
    that is very close to the deviance profiled over $\bm\beta$.
  \item Another approach, which is being implemented as a Google
    Summer of Code project, is adaptive Gauss-Hermite quadrature
    (AGQ).  This has a similar structure to the Laplace approximation
    but is based on more evaluations of the unscaled conditional
    density near the conditional modes.  It is only appropriate for
    models in which the random effects are associated with only one
    grouping factor
  \end{itemize}
\end{frame}