Introduction Example

Suppose we have \(i = 1,\dots,n\) twin pairs.

For each twin we observe a vector of variables \(\vec x_{ij}\) for \(j = 1,2\).

We have some outcome \(Y_{ij} \in \{0, 1\}\) for \(j = 1,2\).

Introduction Example (cont.)

We suppose that there is an association between \(\vec x_{ij}\) and \(Y_{ij}\). A simple model is:

\[\logit(\Prob(Y_{ij} = 1)) = \vec\beta^\top\vec x_{ij}\]

Exponential Family

This model is in the exponential family. That is, the density is:

\[ \begin{align*} f_{ij}(y) &= \exp(y \eta_{ij}(\vec\beta) - b(\eta_{ij}(\vec\beta)) + c(y)) \\ \eta_{ij}(\vec\beta) &= \vec\beta^\top\vec x_{ij} \end{align*} \]

where \(b\) is a log-partition function. This is only the canonical form.

The rest of the talk equally applies for the general exponential family.

Adding Random effects

There may be genetic effects or environmental effects that change the association between \(\vec x_{ij}\) and \(Y_{ij}\).

Marginal Likelihood

\[ \begin{align*} f_i(\vec y_i) &= \int g_i(\vec u;\vec\theta)\prod_{j = 1}^2 f_{ij}(y_{ij}\mid \vec u) \der\vec u \\ f_{ij}(y_{ij}\mid \vec u) &= \exp\bigg(y_{ij} \eta_{ij}(\vec\beta,\vec u) + b(\eta_{ij}(\vec\beta,\vec u)) + c(y_{ij})\bigg) \\ \eta_{ij}(\vec\beta,\vec u) &= (\vec u + \vec\beta)^\top\vec x_{ij} \end{align*} \]

where \(g_i\) is the density of \(\vec U_i\) given model parameters \(\vec\theta\). Generally Intractable.

Remarks

Typically, one works with:

\[\logit(\Prob(Y_{ij} = 1 \mid \vec U_i = \vec u)) = \vec\beta^\top\vec x_{ij} + \vec u^\top\vec z_{ij}\]

where \(\vec z_{ij}\) contains a subset of \(\vec x_{ij}\).

Talk Overview

Variational approximations.

Estimation methods.

Examples.

Conclusions.

Lower Bound

Define the integrand

\[ \begin{align*} h_i(\vec y_i, \vec u) &= g_i(\vec u;\vec\theta) \prod_{j = 1}^{l_i}f_{ij}(y_{ij}\mid \vec u) \end{align*} \]

such that

\[ \begin{align*} f_i(\vec y_i) &= \int h_i(\vec y_i, \vec u) \der\vec u \\ f_{ij}(y_{ij}\mid \vec u) &= \exp\bigg(y_{ij} \eta_{ij}(\vec\beta,\vec u) + b(\eta_{ij}(\vec\beta,\vec u)) + c(y_{ij})\bigg) \\ \eta_{ij}(\vec\beta,\vec u) &= \vec\beta^\top\vec x_{ij} + \vec u^\top\vec z_{ij} \end{align*} \]

Lower Bound (cont.)

Select some variational distribution with density \(\nu_i\) parameterized by some set \(\Omega_i\). Then

\[ \begin{align*} \log f_i(\vec y_i) &\geq \int \nu_i(\vec u;\vec\omega) \log\left( \frac{h_i(\vec y_i,\vec u)} {\nu_i(\vec u;\vec\omega)}\right)\der\vec u \\ &= \log \tilde f_i(\vec y_i;\vec\omega) \end{align*} \]

for some \(\vec\omega \in \Omega_i\).

Approximate Maximum Likelihood

Use

\[ \argmax_{\vec\beta,\vec\theta,\vec\omega_1,\dots,\vec\omega_n} \sum_{i = 1}^n \log \tilde f_i(\vec y_i;\vec\omega_i) \]

instead of

\[ \argmax_{\vec\beta,\vec\theta} \sum_{i = 1}^n \log f_i(\vec y_i) \]

Remarks

Often very simple (at worst one dimensional integrals):

\[\int \nu_i(\vec u;\vec\omega) \log\left( h_i(\vec y_i,\vec u)\right)\der\vec u\]

Easily computed with e.g. adaptive Gaussian-Hermite quadrature.

Computational Issues

The number of parameters are \(\bigO{n}\).

Ormerod and Wand (2012) show the formulas for GVAs for GLMMs.

Computational Issues (cont.)

The objective functions is:

\[ l(\vec\beta,\vec\theta,\vec\omega_1,\dots,\vec\omega_n) = \sum_{i = 1}^n \log \tilde f_i(\vec y_i;\vec\omega_i) \]

Newton’s Method

Want to minimize some function \(f\) with \(o\) inputs.

Good convergence properties if started close to an optimum. Inversion is \(\bigO {o^3}\) in general.

BFGS

Still a \(\bigO{o^2}\) iteration cost and memory cost.

Limited-memory BFGS

Make an approximation of \(\nabla^2 f(\vec x^{(k)})\) using only the last couple of gradients, \(\nabla f(\vec x^{(k)})\).

\(\bigO o\) iteration cost and memory cost.

Poor convergence properties.

Partial Separability

Suppose that

\[f(\vec x) = \sum_{i = 1}^n f_i(\vec x_{\mathcal I_i})\]

where

\[\begin{align*} \vec x_{\mathcal I_i} &= (\vec e_{j_{i1}}^\top, \dots ,\vec e_{j_{im_i}}^\top)\vec x\\ \mathcal I_i &= (j_{i1}, \dots, \mathcal j_{im_i}) \subseteq \{1, \dots, o\} \end{align*},\]

\(\vec e_k\) is the \(k\)’th column of the \(o\) dimensional identity matrix, and \(m_i \ll o\) for all \(i = 1,\dots,n\).

Partial Separability (cont.)

Make \(n\) BFGS approximations of \(\nabla^2 f_1(\vec x^{(k)}_{\mathcal I_1}),\dots,\nabla^2 f_n(\vec x^{(k)}_{\mathcal I_n})\).

Preserves the sparsity of \(\nabla^2 f(\vec x^{(k)})\) and provides potentially quicker convergence towards a good approximation.

Implementation

Implemented in the psqn package (Christoffersen 2020).

Provides a header-only C++ library and a R interface.

Supports computation in parallel.

Remarks

Closely followed Nocedal and Wright (2006).

1D Random Effects

Three dimensional fixed effects, \(\vec x_{ij}, \vec\beta\in\mathbb R^3\).

One dimensional random effect, \(\vec U_i \in \mathbb R\) and \(\vec z_{ij} = (1)\).

Compare with the Laplace approximation and the adaptive Gauss-Hermite quadrature (AGHQ) from the lme4 package (Bates et al. 2015).

1D Random Effects (small)

\(n = 1000\) clusters. Bias estimates are:

Monte Carlo standard errors are in parentheses. The first three columns are the fixed effects, \(\vec\beta\). Inter.: intercept, cont.: continuous covariate, and bin.: binary covariate. The last column is the random effect standard deviation. The GVA row uses the psqn package.

1D Random Effects (comp. time)

Computation times are in seconds. “GVA LBFGS” uses the limited-memory BFGS implementation in the lbfgs package (Coppola, Stewart, and Okazaki 2020).

1D Random Effects (large)

\(n = 5000\) clusters. Bias estimates are:

Monte Carlo standard errors are in parentheses. The first three columns are the fixed effects, \(\vec\beta\). Inter.: intercept, cont.: continuous covariate, and bin.: binary covariate. The last column is the random effect standard deviation. The GVA row uses the psqn package.

1D Random Effects (comp. time)

Computation times are in seconds. “GVA LBFGS” uses the limited-memory BFGS implementation in the lbfgs package.

3D Random Effects

Three dimensional random effect, \(\vec U_i \in \mathbb R^3\) and \(\vec z_{ij} = \vec x_{ij}\).

\(n = 1000\) clusters. Ten observations per cluster.

3D Random Effects

Bias estimates are:

Monte Carlo standard errors are in parentheses. The first three columns are the fixed effects, \(\vec\beta\). Inter.: intercept, cont.: continuous covariate, and bin.: binary covariate. The last three columns are the random effect standard deviations. The GVA row uses the psqn package.

3D Random Effects

Bias estimates of the correlation coefficients. Monte Carlo standard errors are in parentheses.

3D Random Effects (comp. time)

Computation times are in seconds. “GVA LBFGS” uses the limited-memory BFGS implementation in the lbfgs package.

Conclusions

Introduced variational approximations for generalized linear mixed models.

The psqn package exploits the sparsity of the Hessian and yields very fast estimation times.

Extension

Easy to extend to other members of the exponential family and use other link functions.

Can use other variational distributions.

Extension to the psqn Package

Add Newton conjugate gradient method.

Add constraints.

Add a trust region method.

Add other sparsity patterns of the Hessian.