Gaussian variational approximation with composite likelihood for crossed random effect models

Generalized linear mixed models with crossed random effects are useful for the analysis and inference of cross-classified data, such as arise in educational studies (Chung and Cai, 2021; Menictas et al., 2022), psychometric research studies (Baayen et al., 2008; Jeon et al., 2017), and medical studies (Coull et al., 2001), among others.

Suppose $Y_{ij},i=1,\dots,m;j=1,\dots,n$ are conditionally independent, given random effects $U_{i}$ and $V_{j}$. A generalized linear model for $Y_{ij}$ has density function

and we relate this to the covariates in the usual manner:

where $\mu_{ij}=E(Y_{ij})$. In (2) $X_{ij}=(1,X_{ij})^{{\mathrm{\scriptscriptstyle T}}}$, $\beta$ is a $(p+1)$-vector of fixed-effects parameters, and $U_{i}$ and $V_{j}$ are independent random effects assumed to follow $N(0,\sigma^{2}_{u})$ and $N(0,\sigma^{2}_{v})$ distributions, respectively.

In this paper, we develop a Gaussian variational approximation to a form of composite likelihood for model (1), in order to make the computations more tractable than in a full likelihood approach. We focus on two examples: Poisson regression , where $Y_{ij}\sim\text{Poisson}(\mu_{ij})$, with $\theta_{ij}=\log(\mu_{ij})$, $b(\theta_{ij})=\exp(\theta_{ij})$, $g(\mu_{ij})=\theta_{ij}$, and $\phi=1$, and Gamma regression, where $Y_{ij}$ distributed as a Gamma distribution with shape parameter $\alpha=\phi^{-1}$ and expectation parameter $\mu_{ij}$, so that $\theta_{ij}=-\mu_{ij}^{-1}$, $b(\theta_{ij})=\log(\mu_{ij})$, and $g(\mu_{ij})=\log(\mu_{ij})$.

Other approaches to simplifying the likelihood for crossed random effects have been proposed. Penalized quasi-likelihood is discussed in Breslow and Clayton (1993); Schall (1991), based on Laplace approximation, although for binary data the resulting estimates are not guaranteed to be consistent (Lin and Breslow, 1996). Coull et al. (2001) developed a Monte Carlo Newton-Raphson algorithm based on expectation-maximization (McCulloch, 1997), and Ghosh et al. (2022a) proposed a backfitting algorithm for a Gaussian linear model, where the integral can be evaluated explicitly. This approach was extended to logistic regression in Ghosh et al. (2022b).

Composite likelihood methods have also been proposed to simplify computation in models with random effects. Renard et al. (2004) and Bartolucci and Lupparelli (2016) proposed pairwise likelihood for nested random effects models, and Bellio and Varin (2005) developed a pairwise likelihood approach to crossed random effects.

Our approach combines composite likelihood with variational inference. The theory is developed by extending results of Ormerod and Wand (2012), Hall et al. (2011a), and Hall et al. (2011a) which treat Poisson models with nested random effects. Related variational approaches are developed in Menictas et al. (2022) for linear models, and in Rijmen and Jeon (2013), Jeon et al. (2017), Hui et al. (2017), Hui et al. (2019) for generalized linear mixed models.

We introduce a row-column composite likelihood that requires the computation of only one-dimensional integrals and apply a Gaussian variational approximation to this composite likelihood. We prove consistency and asymptotic normality of the variational estimates for the Poisson and Gamma regression models. Simulations demonstrate that our method is faster than the Gaussian variational approximation to the full likelihood function.

The log-likelihood function for the generalized linear model with crossed random effects is constructed from the marginal distribution of the responses, $Y_{ij}$, given the covariates $X_{ij}$, and depends on the parameters $\Psi=(\beta,\sigma^{2}_{u},\sigma^{2}_{v})$:

with maximum likelihood estimator $\hat{\Psi}=\arg\max_{\Psi}\ell(\Psi)$. To simplify the computation of the $(m+n)$-dimensional integral we apply a Gaussian variational approximation to (3) by introducing pairs of variational parameters $(\mathbf{\mu}_{u_{i}},\lambda_{u_{i}}),i=1,...,m$ and $(\mathbf{\mu}_{v_{j}},\lambda_{v_{j}}),j=1,...,n$. By Jensen’s inequality and concavity of the logarithm function:

where $\phi(U_{i})$ and $\phi(V_{j})$ are Gaussian densities with means $\mu_{u_{i}}$, $\mu_{v_{j}}$, and variances $\lambda_{u_{i}}$, $\lambda_{v_{j}}$, respectively. The function $\underline{\ell}(\Psi,\xi)$ is the variational lower bound on $\ell(\Psi)$; the variational parameters are $\xi=(\mu_{u_{1}},...,\mu_{u_{m}},\mu_{v_{1}},...,\mu_{v_{n}},\lambda_{u_{i}},...,\lambda_{u_{m}},\lambda_{v_{1}},...,\lambda_{v_{n}})^{{\mathrm{\scriptscriptstyle T}}}$. Ignoring some constants, this lower bound simplifies to

The Gaussian variational approximation estimators maximize the parameters of the evidence lower bound; $(\widehat{\underline{\Psi}},\widehat{\underline{\xi}})={\arg\max}_{\Psi,\xi}\underline{\ell}(\Psi,\xi).$ The advantage of using the variational lower bound $\underline{\ell}(\Psi,\xi)$ over $\ell(\Psi)$ is that the former only contains terms $E_{\phi(u_{i})\phi(v_{j})}(\theta_{ij})$ or $E_{\phi(u_{i})\phi(v_{j})}\{b(\theta_{ij})\}$ involving a two-dimensional integral, and in our models these can be evaluated explicitly. For the Poisson regression model,

For the Gamma regression model,

In general models, computation of the variational approximation requires $mn$ two-dimensional integrals; for the Poisson and Gamma regression models this simplifies to $mn$ function evaluations. To reduce computation further we develop a version of composite likelihood by considering the two random effects $U_{i}$ and $V_{j}$ separately.

First, we ignore the dependence among rows, $Y^{(1)}_{i}=(Y_{i1},...,Y_{in}),i=1,...,m$, by ignoring the column random effects $V_{j}$’s; i.e., in (2) we replace $X^{{\mathrm{\scriptscriptstyle T}}}_{ij}\beta+U_{i}+V_{j}$ with $X^{{\mathrm{\scriptscriptstyle T}}}_{ij}\beta^{r}+U_{i}$ , and define the row-composite likelihood function by

The row-composite log-likelihood function is

By a similar operation ignoring row random effects, we can define the column-composite log-likelihood function by

The notation $\beta^{r}=(\beta^{r}_{0},\beta_{1},\dots,\beta_{p})$ and $\beta^{c}=(\beta^{c}_{0},\beta_{1},\dots,\beta_{p})$ is needed, as the misspecification caused by ignoring the column (or row) random effects changes the intercept from $\beta_{0}$, say, to $\beta^{r}_{0}=\beta_{0}+\sigma^{2}_{v}/2$ for the row-composite likelihood and $\beta^{c}_{0}=\beta_{0}+\sigma^{2}_{u}/2$ for the column-composite likelihood, for both the Poisson and the Gamma regression models. The other components, $\beta_{1},\beta_{2},..,\beta_{p},$ are unchanged.

Based on (5) and (6), we propose the misspecified row-column composite log-likelihood function by

where $\Psi^{rc}=(\beta^{r}_{0},\beta^{c}_{0},\beta_{1},\dots,\beta_{p},\sigma^{2}_{u},\sigma^{2}_{v})$. This definition of misspecified composite likelihood reduces the computation of an $(m+n)$-dimensional integral in (3) to $m+n$ one-dimensional integrals. Bartolucci et al. (2017) proposed a different composite likelihood function in which column (or row) random effects are marginalized over instead.

We construct a Gaussian variational approximation to the misspecified row-column composite likelihood (7) using the same variational distributions as in the previous section, leading to

For the Poisson regression model,

For the Gamma regression model,

We define the estimators based on this Gaussian variational approximation by

To get the estimator $\widehat{\underline{\Psi}}$ from $\widehat{\underline{\Psi}}^{rc}$, we only need to convert $\widehat{\underline{\beta}^{r}_{0}}$ and $\widehat{\underline{\beta}^{c}_{0}}$ to $\widehat{\underline{\beta}}_{0}$. This is discussed in §4 Remark 1. The advantage of $\underline{c\ell}(\Psi^{rc},\xi)$ over $\underline{\ell}(\Psi,\xi)$ is that the former only involves $m+n$ one-dimensional integrals, and especially, for Poisson and Gamma regression models with explicit expressions, even one-dimensional integrals are no longer involved; only $m+n$ calculations are needed. The misspecified row-column composite likelihood-based variational approximation enhances the efficiency of calculations compared to the variational approximation to the log-likelihood function.

In this section we present our main results on the consistency and convergence rates of the parameter estimates based on the Gaussian variational approximation introduced above. We denote the true value of the parameter with a superscript $0$.

• •

  Poisson regression model:

  	$\displaystyle Y_{ij}\mid X_{ij},U_{i},V_{j}~\text{follows~independent~Poisson~with~mean}~\exp(\beta^{0}_{0}+\beta^{0}_{1}X_{ij}+U_{i}+V_{j}),$
  	$\displaystyle U_{i}\sim N(0,(\sigma^{2}_{u})^{0}),V_{j}\sim N(0,(\sigma^{2}_{v})^{0}),U_{i}~\text{and}~V_{j}~\text{are independent}.$		(9)

• •

  Gamma regression model:

  	$\displaystyle Y_{ij}\mid X_{ij},U_{i},V_{j}~\text{independent~Gamma~with~$\alpha$ and mean}~\exp(\beta^{0}_{0}+\beta^{0}_{1}X_{ij}+U_{i}+V_{j}),$
  	$\displaystyle U_{i}\sim N(0,(\sigma^{2}_{u})^{0}),V_{j}\sim N(0,(\sigma^{2}_{v})^{0}),U_{i}~\text{and}~V_{j}~\text{are independent},$		(10)

  where the shape parameter $\alpha$ is assumed known.

Under the assumptions (A1)-(A9) in the §S1 of the Supplementary Material, we have the following two theorems.

In addition, we have the following asymptotic normality results.

In this section, we perform simulations to evaluate the effectiveness of the proposed approach. We assess the performance of both Poisson and Gamma regression models with crossed random effects as specified in equations (• ‣ 4) and (• ‣ 4). For both models, we independently generate $X_{ij}$ from $N(1,1)$ for $i=1,\ldots,m$ and $j=1,\ldots,n$. We set $\beta^{0}_{0}=-2$, $\beta^{0}_{1}=-2$, $(\sigma_{u})^{0}=0.5$, and $(\sigma_{v})^{0}=0.5$. For the Gamma regression model, there is an additional shape parameter set as $\alpha=0.8$. We consider two different scenarios: $(m,n)=(50,50)$ and $(m,n)=(100,100)$.

We compare the proposed approach with the Gaussian variational inference based on the full log-likelihood function. We report the mean and standard deviation of the parameter estimates based on 1000 repetitions of Monte Carlo studies. We also calculate the average computational times of both methods, using R version 4.1.1 on a laptop equipped with a 2.8 gigahertz Intel Core i7-1165G7 processor and 16 gigabytes of random access memory. The results are presented in Table 1.

Our results demonstrate that both methods yield similar parameter estimates. However, our proposed approach exhibits a notable advantage in terms of computational efficiency compared to the Gaussian variational approximation to the log-likelihood function. We further investigate additional settings to assess the performance of both methods. The results are included in the supplementary §S6 and Table S1, and the findings are similar.

We illustrate our methods with data concerning motor vehicle insurance in Indonesia, in 2014, obtained from the Financial Services Authority (Adam et al., 2021). The dataset consists of 175,381 claim events from 35 distinct areas over a period of 12 months. Each claim event includes information on the claim amount and deductible. However, 17 out of the 35 areas did not have any claim events during the entire 12-month period. Consequently, our analysis focuses on the remaining 18 areas. Table S3 includes the counts of claim events for each area and month.

We randomly selected one claim event from each area and month. The response variable, denoted as $Y_{ij}$, represents the claim amount, while the explanatory variable, denoted as $X_{ij}$, represents the deductible, where $i=1,\dots,18$ and $j=1,\ldots,12$. A portion of the dataset can be found in Table S4 of the supplementary material. From that table, we can see both $X_{ij}$ and $Y_{ij}$ have large magnitudes, therefore we divide them by $10^{7}$ to avoid singular Hessian matrices in computation (Adam et al., 2021). We then fit a Gamma regression model, introducing random effects for both area and month of occurrence.

To evaluate the proposed method and the Gaussian variational approximation to the log-likelihood function, we regenerated $1000$ datasets using different random seeds at each step of the random event selection process. We then applied both methods and reported the mean and standard deviation of both estimators across these $1000$ datasets. In addition, we record the average time it takes to fit the model. Table 2 presents the results of this analysis.

From Table 2, it is evident that the two methods yield comparable model parameter estimates. However, our proposed approach “GVACL” is much faster in terms of computation time.

In this article, we cover two examples: the Poisson regression and Gamma regression. An interesting future direction is to extend the results to other generalized linear mixed models. For example, the logistic regression model with crossed random effects has

Unlike the Poisson and Gamma regression models, there is no analytic solution to the two-dimensional integral

appearing in (4). One solution is to evaluate the integrals using adaptive Gauss-Hermite quadrature (Liu and Pierce, 1994) with a product of $N_{1}\times N_{2}$ quadrature points over the two-dimensional integral in (25). The one-dimensional integrals arising in the row-column composite likelihood of (8) can also be computed by adaptive Gauss-Hermite quadrature. Compared to the Poisson and Gamma cases, the computational time of the proposed method is longer, especially when $m$ and $n$ are large.

We carried out some limited simulations for the logistic model and report the results in the supplementary §S6 and Table S2. We found that the computational time of the proposed method is much shorter than for the Gaussian variational approximation to the log-likelihood function. Based on the simulation results, we conjecture that the intercept and slope estimates and the sum of random effect variance estimates have the following relationship:

The rigorous derivation of the relationship is beyond the scope of the paper and we leave it for future research.

Variational approximations to the composite log-likelihood function of crossed random effects models establish a connection between two distinct topics in statistics. We introduce the row-composite likelihood function (5) to eliminate the row dependence of responses by disregarding the column random effects. Subsequently, we utilize the Gaussian variational approximation to eliminate the column dependence of responses through the variational distributions of row random effects. Similarly, we construct the column-composite likelihood function (6), to eliminate the column dependence of responses by disregarding the random row effects. We then apply the Gaussian variational approximation to eliminate the row dependence of responses through the variational distributions of column random effects. Comparing the row-column composite likelihood-based Gaussian variational approximation and likelihood-based Gaussian variational approximation, both approaches share the common goal of breaking the dependence of responses through manipulating random effects for the generalized linear mixed models with crossed random effects. For any other links between composite likelihood and variational approximations, we will leave them for future researc

The authors are grateful to Dr. Fia Fridayanti Adam for sharing the motor vehicle insurance data. This research was partially supported by the Natural Sciences and Engineering Council of Canada.