Component-based regularisation of multivariate generalised linear mixed models

TL;DR

This paper introduces a component-based regularisation method for multivariate GLMMs that effectively handles high-dimensional, redundant explanatory variables by constructing orthogonal components, improving model estimation for grouped data.

Contribution

It proposes an extension of SCGLR for regularising multivariate GLMMs, combining component construction with an adapted estimation algorithm, tested on simulated and real data.

Findings

The method outperforms ridge and LASSO regularisations in grouped data scenarios.

Orthogonal components effectively capture relevant structural information.

The approach is validated on both simulated and real datasets.

Abstract

We address the component-based regularisation of a multivariate Generalised Linear Mixed Model (GLMM) in the framework of grouped data. A set Y of random responses is modelled with a multivariate GLMM, based on a set X of explanatory variables, a set A of additional explanatory variables, and random effects to introduce the within-group dependence of observations. Variables in X are assumed many and redundant so that regression demands regularisation. This is not the case for A, which contains few and selected variables. Regularisation is performed building an appropriate number of orthogonal components that both contribute to model Y and capture relevant structural information in X. To estimate the model, we propose to maximise a criterion specific to the Supervised Component-based Generalised Linear Regression (SCGLR) within an adaptation of Schall's algorithm. This extension of SCGLR is tested on both simulated and real grouped data, and compared to ridge and LASSO regularisations. Supplementary material for this article is available online.

Tables

Table 1. Table 1: Optimal regularisation parameter values obtained through cross–validation over 100 100 100 simulations.

	GLMM–LASSO	LMM–ridge	mixed–SCGLR
	shrinkage parameter ${\lambda }_{\text{LASSO}}^{\star }$	shrinkage parameter ${\lambda }_{\text{ridge}}^{\star }$	number of components $K^{\star }$	trade–off parameter $s^{\star }$
$\tau =0.1$	65	24	15	0.50
$\tau =0.3$	92	54	5	0.58
$\tau =0.5$	124	73	3	0.70
$\tau =0.7$	163	78	3	0.73
$\tau =0.9$	175	85	2	0.80

Table 2. Table 2: Mean Upper Relative Squared Error (MURSE) values associated with the optimal parameter values.

	LMM	GLMM–LASSO	LMM–ridge	mixed–SCGLR
	(no regularisation)
$\tau =0.1$	$0.12$	$0.05$	$0.08$	$0.12$
$\tau =0.3$	$0.33$	$0.12$	$0.13$	$0.10$
$\tau =0.5$	$0.61$	$0.20$	$0.16$	$0.07$
$\tau =0.7$	$1.32$	$0.25$	$0.20$	$0.06$
$\tau =0.9$	$4.62$	$0.26$	$0.31$	$0.05$

Table 3. Table 3: Mean Relative Squared Error (MRSE) values obtained with Bernoulli and Poisson responses.

	GLMM		GLMM–LASSO		mixed–SCGLR
	(no regularisation)
	Bernoulli	Poisson	Bernoulli	Poisson	Bernoulli	Poisson
$\tau =0.1$	$316.48$	$0.54$	$8.61$	$0.30$	$14.71$	$0.46$
$\tau =0.3$	$398.78$	$0.64$	$9.23$	$0.36$	$7.21$	$0.21$
$\tau =0.5$	$576.68$	$0.87$	$14.48$	$0.44$	$2.01$	$0.09$
$\tau =0.7$	$886.04$	$1.28$	$17.37$	$0.47$	$1.50$	$0.07$
$\tau =0.9$	$2840.10$	$3.72$	$17.24$	$0.59$	$1.31$	$0.05$

Table 4. Table 4: Cross–validation errors for each response variable.

	$E_{\text{cv}}^1$	$E_{\text{cv}}^2$	$E_{\text{cv}}^3$	$E_{\text{cv}}^4$	$E_{\text{cv}}^5$	$E_{\text{cv}}^6$	$E_{\text{cv}}^7$	$E_{\text{cv}}^8$
SCGLR	$1.32$	$2.46$	$3.27$	$1.43$	$2.56$	$1.28$	$1.54$	$3.44$
mixed–SCGLR	$1.24$	$1.95$	$2.92$	$1.32$	$2.27$	$1.15$	$1.31$	$3.01$

Table 5. Table 5: Mean Relative Squared Error (MRSE) values obtained with binomial and Poisson distributions. The R package glmmLasso does not handle binomial outcomes but only Bernoulli ones, which precludes comparison in this case.

	GLMM		GLMM–LASSO	mixed–SCGLR
	(no regularisation)
	Binomial	Poisson	Poisson	Binomial	Poisson
$\tau =0.1$	$2.31$	$0.50$	$0.31$	$0.51$	$0.45$
$\tau =0.3$	$3.07$	$0.60$	$0.33$	$0.28$	$0.18$
$\tau =0.5$	$3.93$	$0.75$	$0.39$	$0.15$	$0.09$
$\tau =0.7$	$6.50$	$1.07$	$0.40$	$0.10$	$0.07$
$\tau =0.9$	$19.29$	$2.71$	$0.42$	$0.07$	$0.05$

Table 6. Table 6: Summary of cor j subscript cor 𝑗 \text{cor}_{j} and err j subscript err 𝑗 \text{err}_{j} values, and presentation of biases and standard errors of estimated variance components.

	$N=10,R=10\left(n=100\right)$			$N=50,R=10\left(n=500\right)$
${\sigma }_{\text{LV}}^2/{\sigma }_{\text{noise}}^2$	$\frac{1}{3}$	$1$	$3$	$\frac{1}{3}$	$1$	$3$
${\text{cor}}_1$	$0.71$	$0.91$	$0.96$	$0.75$	$0.92$	$0.96$
${\text{cor}}_2$	$0.93$	$0.94$	$0.94$	$0.97$	$0.98$	$0.98$
${\text{err}}_1$	$0.47$	$0.15$	$0.06$	$0.38$	$0.13$	$0.05$
${\text{err}}_2$	$0.12$	$0.12$	$0.12$	$0.05$	$0.04$	$0.04$
${\text{err}}_3$	$0.19$	$0.14$	$0.11$	$0.11$	$0.07$	$0.04$
$\text{bias}\left(\widehat{{\sigma }_1^2}\right)$	$-0.02$	$-0.01$	$0.00$	$0.02$	$0.00$	$-0.02$
$\text{sd}\left(\widehat{{\sigma }_1^2}\right)$	$1.04$	$1.05$	$1.06$	$0.41$	$0.40$	$0.39$
$\text{bias}\left(\widehat{{\sigma }_2^2}\right)$	$-0.11$	$-0.08$	$-0.06$	$-0.06$	$-0.06$	$-0.06$
$\text{sd}\left(\widehat{{\sigma }_2^2}\right)$	$0.50$	$0.51$	$0.52$	$0.21$	$0.21$	$0.21$
$\text{bias}\left(\widehat{{\sigma }_3^2}\right)$	$-0.03$	$-0.04$	$-0.04$	$-0.02$	$-0.02$	$-0.02$
$\text{sd}\left(\widehat{{\sigma }_3^2}\right)$	$0.22$	$0.21$	$0.21$	$0.11$	$0.11$	$0.11$

Table 7. Table 7: Mean Relative Squared Error (MRSE) values for fixed–effect estimates, and biases and standard errors for estimated variance components, obtained with 100 100 100 observations and 150 150 150 explanatory variables.

	${𝜷}_{\mathrm{𝟏}}$	${𝜷}_{\mathrm{𝟐}}$	${𝜷}_{\mathrm{𝟑}}$	${𝜷}_{\mathrm{𝟒}}$	${\sigma }_1^2$	${\sigma }_2^2$	${\sigma }_3^2$	${\sigma }_4^2$
$\tau =0.3$	$0.06$	$0.26$	$0.19$	$0.13$	$-0.01$	$-0.03$	$-0.02$	$0.02$
					$0.09$	$0.09$	$0.03$	$0.06$
$\tau =0.5$	$0.03$	$0.20$	$0.10$	$0.07$	$0.01$	$-0.03$	$0.00$	$0.01$
					$0.11$	$0.08$	$0.07$	$0.07$
$\tau =0.7$	$0.01$	$0.10$	$0.05$	$0.04$	$0.01$	$-0.05$	$0.01$	$0.02$
					$0.07$	$0.09$	$0.10$	$0.07$

Table 8. Table 8: Mean Relative Squared Error (MRSE) values for fixed–effect estimates, and biases and standard errors for estimated variance components, obtained with 100 100 100 observations and 200 200 200 explanatory variables.

	${𝜷}_{\mathrm{𝟏}}$	${𝜷}_{\mathrm{𝟐}}$	${𝜷}_{\mathrm{𝟑}}$	${𝜷}_{\mathrm{𝟒}}$	${\sigma }_1^2$	${\sigma }_2^2$	${\sigma }_3^2$	${\sigma }_4^2$
$\tau =0.3$	$0.06$	$0.15$	$0.18$	$0.10$	$-0.04$	$-0.05$	$0.01$	$-0.02$
					$0.04$	$0.09$	$0.05$	$0.05$
$\tau =0.5$	$0.03$	$0.17$	$0.09$	$0.05$	$-0.05$	$0.00$	$-0.02$	$-0.01$
					$0.06$	$0.19$	$0.04$	$0.04$
$\tau =0.7$	$0.01$	$0.15$	$0.04$	$0.03$	$0.03$	$0.00$	$-0.01$	$-0.02$
					$0.08$	$0.14$	$0.05$	$0.05$

Equations

$${\mathbf{\beta }}^{[t+1]}={\left({\mathbf{X}}^{\text{mbox}T}{\mathbf{W}}^{[t]}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\text{mbox}T}{\mathbf{W}}^{[t]}{\mathbf{z}}^{[t]},$$

$${\mathcal{M}}^{[t]}:|\begin{array}{rl}& {\mathbf{z}}^{[t]}=\mathbf{X}\mathbf{\beta }+{\mathbf{\zeta }}^{[t]}\\ & \text{mbox}with:\mathbb{E}\left({\mathbf{\zeta }}^{[t]}\right)=0\text{mbox}and\mathbb{V}\left({\mathbf{\zeta }}^{[t]}\right)={{\mathbf{W}}^{[t]}}^{-1}.\end{array}$$

$${\mathbf{\eta }}_{\mathbf{k}}=\left(\mathbf{Xu}\right){\gamma }_k+\mathbf{A}{\mathbf{\delta }}_{\mathbf{k}},$$

$$\ell \left(\mathbf{Y}|\mathbf{H}\right)=\sum _{i=1}^n\sum _{k=1}^q{\ell }_k\left(y_{k,i}|{\eta }_{k,i}\right),$$

$${\mathcal{M}}_k^{[t]}:|\begin{array}{rl}& {\mathbf{z}}_{\mathbf{k}}^{[t]}=\left(\mathbf{Xu}\right){\gamma }_k+\mathbf{A}{\mathbf{\delta }}_{\mathbf{k}}+{\mathbf{\zeta }}_{\mathbf{k}}^{[t]}\\ & \text{mbox}with:\mathbb{E}\left({\mathbf{\zeta }}_{\mathbf{k}}^{[t]}\right)=0\text{mbox}and\mathbb{V}\left({\mathbf{\zeta }}_{\mathbf{k}}^{[t]}\right)={{\mathbf{W}}_{\mathbf{k}}^{[t]}}^{-1}.\end{array}$$

$$\underset{\mathbf{u}:{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1}{\min }\sum _{k=1}^q{\Vert {\mathbf{z}}_{\mathbf{k}}-{\Pi }_{\mathrm{span}\left\{\mathbf{Xu},\mathbf{A}\right\}}^{{\mathbf{W}}_{\mathbf{k}}}{\mathbf{z}}_{\mathbf{k}}\Vert }_{{\mathbf{W}}_{\mathbf{k}}}^2⟺\underset{\mathbf{u}:{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1}{\max }\sum _{k=1}^q{\Vert {\Pi }_{\mathrm{span}\left\{\mathbf{Xu},\mathbf{A}\right\}}^{{\mathbf{W}}_{\mathbf{k}}}{\mathbf{z}}_{\mathbf{k}}\Vert }_{{\mathbf{W}}_{\mathbf{k}}}^2.$$

$${\psi }_{\mathbf{A}}\left(\mathbf{u}\right)$$

$$=\sum _{k=1}^q{\Vert {{\mathbf{z}}_{\mathbf{k}}}^{}\Vert }_{{\mathbf{W}}_{\mathbf{k}}}^2{\cos }_{{\mathbf{W}}_{\mathbf{k}}}^2\left({\mathbf{z}}_{\mathbf{k}},{\Pi }_{\mathrm{span}\left\{\mathbf{Xu},\mathbf{A}\right\}}^{{\mathbf{W}}_{\mathbf{k}}}{\mathbf{z}}_{\mathbf{k}}\right).$$

$$\underset{\mathbf{u}:{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1}{\max }\varphi \left(\mathbf{u}\right),$$

$$\phi\left(\boldsymbol{u}\right)=\left[\sum_{j=1}^{p}\omega_{j}\left(\left\langle\,\boldsymbol{Xu}\,\middle|\,\boldsymbol{x_{j}}\,\right\rangle_{\boldsymbol{P}}^{2}\right)^{l}\right]^{\frac{1}{l}}=\left[\sum_{j=1}^{p}\omega_{j}\Big{(}\boldsymbol{u}^{\mbox{\tiny T}}\,\boldsymbol{X}^{\mbox{\tiny T}}\boldsymbol{P}\boldsymbol{x_{j}}\boldsymbol{x}_{\boldsymbol{j}}^{\mbox{\tiny T}}\boldsymbol{PX}\,\boldsymbol{u}\Big{)}^{l}\right]^{\frac{1}{l}},\;l\geqslant 1,$$

$$\underset{\mathbf{u}:{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1}{\max }{\left[\varphi \left(\mathbf{u}\right)\right]}^s{\left[{\psi }_{\mathbf{A}}\left(\mathbf{u}\right)\right]}^{1-s}$$

$$\underset{\mathbf{u}:{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1}{\max }s\log \left[\varphi \left(\mathbf{u}\right)\right]+(1-s)\log \left[{\psi }_{\mathbf{A}}\left(\mathbf{u}\right)\right].$$

$$\forall k\in\left\{1,\ldots,q\right\},\;\;\boldsymbol{\xi_{k}}\,\mathrel{\overset{\text{ind.}}{\scalebox{1.5}[1.0]{$\sim$}}}\,\mathcal{N}_{N}\left(\boldsymbol{0},\boldsymbol{D_{k}}\right),$$

$$\forall k\in \left\{1,\dots ,q\right\},{\mathbf{\eta }}_{\mathbf{k}}^{\mathbf{\xi }}=(\mathbf{Xu}){\gamma }_k+\mathbf{A}{\mathbf{\delta }}_{\mathbf{k}}+\mathbf{U}{\mathbf{\xi }}_{\mathbf{k}},$$

$$z_{k,i}^{\xi }$$

$$={\eta }_{k,i}^{\xi }+e_{k,i},\text{where}{e}_{k,i}=\left(y_{k,i}-{\mu }_{k,i}^{\xi }\right)g_k^{\prime }\left({\mu }_{k,i}^{\xi }\right).$$

$$\text{mbox}Var\left({\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}|{\mathbf{\xi }}_{\mathbf{k}}\right)={{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{-1}=\text{Diag}{\left({\left[g_k^{\prime }\left({\mu }_{k,i}^{\xi }\right)\right]}^2{a}_{k,i}({\varphi }_k)v_k\left({\mu }_{k,i}^{\xi }\right)\right)}_{i=1,\dots ,n},$$

$${{\mathcal{M}}_k^{\xi }}^{[t]}:|\begin{array}{rl}& {{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}=\left(\mathbf{Xu}\right){\gamma }_k+\mathbf{A}{\mathbf{\delta }}_{\mathbf{k}}+\mathbf{U}{\mathbf{\xi }}_{\mathbf{k}}+{\mathbf{e}}_{\mathbf{k}}^{[t]}\\ & \text{mbox}with:\mathbb{E}\left({\mathbf{e}}_{\mathbf{k}}^{[t]}|{\mathbf{\xi }}_{\mathbf{k}}\right)=0\text{mbox}and\mathbb{V}\left({\mathbf{e}}_{\mathbf{k}}^{[t]}|{\mathbf{\xi }}_{\mathbf{k}}\right)={{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{-1}}^{[t]}.\end{array}$$

$${\psi }_{\mathbf{A}}\left(\mathbf{u}\right)=\sum _{k=1}^q{\Vert {\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}\Vert }_{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^2{\cos }_{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^2\left({\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }},\mathrm{span}\left\{\mathbf{Xu},\mathbf{A}\right\}\right).$$

$$\left(\begin{array}{ccc}{\mathbf{f}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{f}& {\mathbf{f}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{A}& {\mathbf{f}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{U}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{A}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{f}& {\mathbf{A}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{A}& {\mathbf{A}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{U}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{U}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{f}& {\mathbf{U}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{A}& {\mathbf{U}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{U}+{\mathbf{D}}_{\mathbf{k}}^{-1}\end{array}\right)\left(\begin{array}{c}{\gamma }_k\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{\delta }}_{\mathbf{k}}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{\xi }}_{\mathbf{k}}\end{array}\right)=\left(\begin{array}{c}{\mathbf{f}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{A}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{U}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}\end{array}\right).$$

$${\sigma }_k^2⟵\frac{{\widehat{{\mathbf{\xi }}_{\mathbf{k}}}}^{\text{mbox}T}\widehat{{\mathbf{\xi }}_{\mathbf{k}}}}{N-\frac{1}{{\sigma }_k^2}\text{Trace}\left[{\left({\mathbf{U}}^{\text{mbox}T}{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}\mathbf{U}+{\mathbf{D}}_{\mathbf{k}}^{-1}\right)}^{-1}\right]}.$$

$${\mathbf{u}}^{[t+1]}=\underset{\mathbf{u}:{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1}{\text{mbox}argmax}{\left[\varphi \left(\mathbf{u}\right)\right]}^s{\left[{\psi }_{\mathbf{A}}^{[t]}\left(\mathbf{u}\right)\right]}^{1-s},\text{where }{\psi }_{\mathbf{A}}\left(\mathbf{u}\right)\text{is given by }(\text{ref}Chauvet::ps{i}_{m}odified)\text{ and }\varphi \left(\mathbf{u}\right)\text{ by }(\text{ref}Chauvet::S{R}_{c}riterio{n}_{p}articulier)$$

$${\mathbf{f}}^{[t+1]}=\mathbf{X}{\mathbf{u}}^{[t+1]}$$

$$\left(\begin{array}{ccc}{{\mathbf{f}}^{[t+1]}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}{\mathbf{f}}^{[t+1]}& {{\mathbf{f}}^{[t+1]}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{A}& {{\mathbf{f}}^{[t+1]}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{U}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{A}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}{\mathbf{f}}^{[t]}& {\mathbf{A}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{A}& {\mathbf{A}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{U}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{U}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}{\mathbf{f}}^{[t]}& {\mathbf{U}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{A}& {\mathbf{U}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{U}+{{\mathbf{D}}_{\mathbf{k}}^{[t]}}^{-1}\end{array}\right)\left(\begin{array}{c}{\gamma }_k\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{\delta }}_{\mathbf{k}}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{\xi }}_{\mathbf{k}}\end{array}\right)=\left(\begin{array}{c}{{\mathbf{f}}^{[t+1]}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}{{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{A}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}{{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ {\mathbf{U}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}{{\mathbf{z}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\end{array}\right)$$

$${{\sigma }_k^2}^{[t+1]}=\frac{{{\mathbf{\xi }}_{\mathbf{k}}^{[t+1]}}^{\text{mbox}T}{\mathbf{\xi }}_{\mathbf{k}}^{[t+1]}}{N-\frac{1}{{{\sigma }_k^2}^{[t]}}\text{Trace}\left[{\left({\mathbf{U}}^{\text{mbox}T}{{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t]}\mathbf{U}+{{\mathbf{D}}_{\mathbf{k}}^{[t]}}^{-1}\right)}^{-1}\right]}\text{and}{\mathbf{D}}_{\mathbf{k}}^{[t+1]}={{\sigma }_k^2}^{[t+1]}{\mathbf{I}}_N$$

$${{\mathbf{\eta }}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t+1]}$$

$${{\mu }_{k,i}^{\xi }}^{[t+1]}$$

$${z_{k,i}^{\xi }}^{[t+1]}$$

$${{\mathbf{W}}_{\mathbf{k}}^{\mathbf{\xi }}}^{[t+1]}$$

$$\{\begin{array}{rl}& \text{max}s\log \left[\varphi \left(\mathbf{u}\right)\right]+(1-s)\log \left[{\psi }_{{\mathbf{A}}_{\mathbf{h}}}\left(\mathbf{u}\right)\right]\\ & \text{subject to:}{\mathbf{u}}^{\text{mbox}T}{\mathbf{M}}^{-1}\mathbf{u}=1\text{and}{\mathbf{D}}_{\mathbf{h}}^{\text{mbox}T}\mathbf{u}=0,\end{array}$$

$$\text{cor}\left({\mathbf{X}}_{\mathbf{j}}\right)=\tau 1_{p_j}{1}_{p_j}^{\text{mbox}T}+(1-\tau ){\mathbf{I}}_{p_j},$$

$$\{\begin{array}{rl}{\mathbf{y}}_1& =\mathbf{X}{\mathbf{\beta }}_1+\mathbf{U}{\mathbf{\xi }}_1+{\mathbf{\epsilon }}_1\\ {\mathbf{y}}_2& =\mathbf{X}{\mathbf{\beta }}_2+\mathbf{U}{\mathbf{\xi }}_2+{\mathbf{\epsilon }}_2,\end{array}$$

$${\mathbf{\beta }}_1$$

Full text

Keywords: generalised linear regression, supervised components, random effects, structural relevance.

1 Introduction

In the framework of regression models on a large number of explanatory variables with redundancies and collinearities, the search for a reduced number of relevant dimensions to model responses has been an ongoing research over the last decades. In particular, the case where the explanatory variables outnumber the observations tends to be a new standard. Generalised Linear Models (GLMs) are the most widely used regression models, because they are easy to interpret and address a very large scope of applications with a variety of response distributions. For instance, Epidemiology, Biology and Social Sciences need to model binary outcomes, count data and survival times. All these fields often have to deal with both grouped data and multivariate responses combining variables of different types (e.g. one binary and another Poisson). In this work, we particularly aim at modelling abundances of several tree genera on plots of land grouped into forest concessions, using multiple redundant explanatory variables.

As far as dimension–reduction is concerned, two main approaches have been developed. The first one is variable–selection, whereas the second one builds components, i.e. linear combinations of the explanatory variables, which synthesise the useful part of their information. As far as variable–selection is concerned, the most popular method is currently the LASSO, introduced by Tibshirani, (1996), which combines the likelihood with a penalty based on the $L_{1}$–norm of the coefficient vector. LASSO is one of the penalty–based regularisation methods, as are also ridge (Hoerl and Kennard,, 1970) and elastic–net (Zou and Hastie,, 2005). This LASSO selection approach has proved efficient to explain the phenomenon of interest when some of the explanatory variables are the “true” ones, surrounded by a high number of irrelevant others. Nevertheless, it may be very unstable and helpless when the true explanatory dimensions are latent and indirectly measured through highly correlated proxies. This is where the component–based approach turns out to be useful. Bry et al., (2013) have developed a new methodology named Supervised Component–based Generalised Linear Regression (SCGLR), later extended and refined in Bry et al., (2014, 2016, 2018). As in any PLS–type method, the construction of components in SCGLR is guided both by the correlation–structure of variables in the explanatory space and by the prediction quality of the responses. Nevertheless, unlike PLS, SCGLR involves a general and flexible criterion allowing to specify the type of structure components are wanted to align with in the explanatory space (e.g. variable bundles, principal components, other subspaces). Moreover, SCGLR searches for explanatory directions common to multiple responses with probability distributions in the exponential family, each response being entitled to their own distribution. The current SCGLR method is implemented in the R package SCGLR (Cornu et al.,, 2018) available at https://CRAN.R-project.org/package=SCGLR and https://github.com/SCnext/SCGLR.

In the present work, we aim at modelling responses with a repeated or grouped design. For this purpose, the use of mixed models with random effects is widespread. Research on variance–component estimation in Generalised Linear Mixed Models (GLMMs) has been very active since the 1980s. For the most general distribution assumptions in such models, parameter estimation faces the intractability of the likelihood expressed as an integral with respect to the random effects. Several numerical approximations of the integral have been proposed: Gaussian quadrature (Anderson and Aitkin,, 1985) or adaptive versions of it (Pinheiro and Bates,, 1995), Laplace approximation leading to the definition of the penalised quasi–likelihood (Breslow and Clayton,, 1993) or modified versions of it (Shun and McCullagh,, 1995). An alternative to this type of analytic approximation is a stochastic approximation of the integral calculation via MCMC techniques. In this approach, Zeger and Karim, (1991) described an approximate Gibbs sampling for GLMMs, which was extended by Clayton, (1996) to more general Metropolis–Hastings algorithms. In parallel, McCulloch, (1997) developed the Monte Carlo EM algorithm where the expectation is computed numerically through a Monte Carlo approximation, after generating random effects with a Metropolis–Hastings sampler. Mention can also be made of the recent work by Knudson, (2016): her strategy is to approximate the entire likelihood function using random effects simulated from a parametrised importance sampling distribution depending on the data. Unfortunately, these different approaches are not necessarily suitable for the same types of random effect designs (one–dimensional random effect, embedded random effects, etc). In the wake of the first type of approximations, we here adopt the “Joint–Maximisation” strategy (McCulloch,, 1997), as introduced for instance by Schall, (1991). The model is iteratively linearised conditional on the random effects and variance components are then estimated using adapted linear mixed models methods. This strategy can be used for any random effect design and is less computationally intensive than Monte Carlo methods. Moreover, it provides us with a linear setting more suitable for the computation of components. Once the components calculated, model parameters can be estimated using any of the aforementioned strategies (see Section 7).

Modelling grouped responses through a GLMM with a large number of explanatory variables is the focus of this paper. The need for dimension–reduction and regularisation has to accommodate the presence of random effects in the model, but our main purpose still remains to investigate the explanatory structure and link it to interpretable dimensions. For Gaussian responses, Eliot et al., (2011) proposed to extend the ridge regression to Linear Mixed Models (LMMs). Based on a penalised complete log–likelihood, the adaptation of the Expectation–Maximisation algorithm they suggest includes a new step to find the best shrinkage parameter using a generalised cross–validation scheme at each iteration. More recently, Schelldorfer et al., (2014) — and also Groll and Tutz, (2014) — proposed an $L_{1}$–penalised algorithm for fitting a high–dimensional GLMM, using Laplace approximation and an efficient coordinate gradient descent. In this work, we combine Schall’s iterative model linearisation with regularisation at each step. However, we do not use a penalty on the coefficient vector’s norm — as proposed by Zhang et al., (2017) within the framework of multivariate count data. We rather propose to combine dimension–reduction and predictor–regularisation using supervised components aligning on the most predictive and interpretable directions in the explanatory space.

The paper is organised as follows. In Section 2, we formalise the model and set the main notations used throughout the paper. In Section 3, we present the key features of SCGLR. Section 4 designs an extension of this methodology to mixed models, and particularly to grouped data. In Section 5, our extended method “mixed–SCGLR” is evaluated on simulations and compared to ridge– and LASSO–based regularisations. Finally, in order to highlight the power of mixed–SCGLR in terms of model interpretation, Section 6 presents an application to real data in the Poisson case.

2 Model definition and notations

In the framework of a multivariate GLMM, we consider $q$ response–vectors $\boldsymbol{y_{1}},\ldots,\boldsymbol{y_{q}}$ forming matrix $\boldsymbol{Y}_{n\times q}$, to be explained by two categories of explanatory variables. The first category consists of few weakly correlated variables \boldsymbol{A}_{n\times r}=\big{[}\,\boldsymbol{a_{1}}\mid\ldots\mid\boldsymbol{a_{r}}\,\big{]}. These variables are assumed to be interesting per se and their marginal effects need to be precisely quantified. The second category consists of abundant and highly correlated variables \boldsymbol{X}_{n\times p}=\big{[}\,\boldsymbol{x_{1}}\mid\ldots\mid\boldsymbol{x_{p}}\,\big{]} considered as proxies to latent dimensions which must be found and interpreted. Since explanatory variables in $\boldsymbol{A}$ are few, non–redundant and of interest, they are kept as such in the model. By contrast, $\boldsymbol{X}$ may contain several unknown structurally relevant dimensions $K<p$ important to model and predict $\boldsymbol{Y}$, how many we do not know. $\boldsymbol{X}$ is thus to be searched for an appropriate number of orthogonal components that both capture relevant structural information in $\boldsymbol{X}$ and contribute to model $\boldsymbol{Y}$.

This work addresses grouped data: the $n$ observations form $N$ groups. Within each group, observations are not assumed independent. For each response $\boldsymbol{y_{k}}$, a $N$–level random effect $\boldsymbol{\xi_{k}}$ is used to model the dependence of observations within each group. Hence, each $\boldsymbol{y_{k}}$ is modelled with a GLMM assuming a conditional distribution from the exponential family.

Notations and conventions

• $\blacktriangleright$

All variables (namely the $\boldsymbol{a_{i}}$’s, $\boldsymbol{x_{j}}$’s and $\boldsymbol{y_{k}}$’s) will be identified with $n$–vectors.

• $\blacktriangleright$

We will use bold lowercase letters for vectors (e.g. $\boldsymbol{u}$) and bold capital letters for matrices (e.g. $\boldsymbol{M}$).

• $\blacktriangleright$

$\boldsymbol{M}$ being any matrix, ${\boldsymbol{M}}^{\mbox{\tiny T}}$ denotes the transpose of $\boldsymbol{M}$.

• $\blacktriangleright$

$\boldsymbol{I}_{n}$ denotes the identity matrix of size $n$.

• $\blacktriangleright$

$\mathbf{1}_{m}$ denotes the all–ones vector of size $m$.

• $\blacktriangleright$

Let $\boldsymbol{u}$ and $\boldsymbol{v}$ be non–zero vectors in $\mathbb{R}^{d}$ and let $\boldsymbol{M}$ be a symmetric positive definite matrix of size $d\times d$. Then $\left\langle\,\boldsymbol{u}\,\middle|\,\boldsymbol{v}\,\right\rangle_{\boldsymbol{M}}={\boldsymbol{u}}^{\mbox{\tiny T}}\boldsymbol{M}\boldsymbol{v}$ refers to the Euclidean scalar product of $\boldsymbol{u}$ and $\boldsymbol{v}$ with respect to metric $\boldsymbol{M}$. The cosine of the angle between $\boldsymbol{u}$ and $\boldsymbol{v}$ with respect to $\boldsymbol{M}$ is given by $\cos_{\boldsymbol{M}}\left(\boldsymbol{u},\boldsymbol{v}\right)=\dfrac{\left\langle\,\boldsymbol{u}\,\middle|\,\boldsymbol{v}\,\right\rangle_{\boldsymbol{M}}}{\left\lVert\boldsymbol{u}\right\rVert_{\boldsymbol{M}}\left\lVert\boldsymbol{v}\right\rVert_{\boldsymbol{M}}}$, where $\left\lVert\boldsymbol{u}\right\rVert_{\boldsymbol{M}}=\sqrt{\left\langle\,\boldsymbol{u}\,\middle|\,\boldsymbol{u}\,\right\rangle_{\boldsymbol{M}}}$.

• $\blacktriangleright$

The space spanned by vectors $\boldsymbol{u_{1}},\ldots,\boldsymbol{u_{h}}$ is denoted by $\operatorname{span}\left\{\boldsymbol{u_{1}},\ldots,\boldsymbol{u_{h}}\right\}$. $\boldsymbol{U}$ being any matrix, $\operatorname{span}\left\{\boldsymbol{U}\right\}$ refers to the space spanned by the column–vectors of $\boldsymbol{U}$.

• $\blacktriangleright$

Let $\mathbb{R}^{n}$ be endowed with metric $\boldsymbol{W}$ and let $\boldsymbol{Z}$ be a matrix of size $n\times p$. Then $\Pi_{\operatorname{span}\left\{\boldsymbol{Z}\right\}}^{\boldsymbol{W}}$ refers to the $\boldsymbol{W}$–orthogonal projector onto $\operatorname{span}\left\{\boldsymbol{Z}\right\}$. Let $\boldsymbol{b}$ be a vector in $\mathbb{R}^{n}$. The cosine of the angle between $\boldsymbol{b}$ and $\operatorname{span}\left\{\boldsymbol{Z}\right\}$ with respect to $\boldsymbol{W}$ is defined by $\cos_{\boldsymbol{W}}\left(\boldsymbol{b},\,\operatorname{span}\left\{\boldsymbol{Z}\right\}\right)=\cos_{\boldsymbol{W}}\left(\boldsymbol{b},\,\Pi_{\operatorname{span}\left\{\boldsymbol{Z}\right\}}^{\boldsymbol{W}}\boldsymbol{b}\right)$.

3 SCGLR with additional explanatory variables

In this section we consider the situation where each $\boldsymbol{y_{k}}$ is modelled with a GLM (without random effect). For the sake of simplicity, we focus on the single–component SCGLR ($K=1$). Section 3.1 briefly recalls some standards for univariate GLMs. Section 3.2 defines the linear predictors considered in the SCGLR methodology, in a multivariate GLM framework with additional explanatory variables. Finally, Section 3.3 introduces the criterion SCGLR maximises to compute the component.

3.1 Notations and main features of univariate GLMs

We refer the reader to McCullagh and Nelder, (1989) for a thorough overview of GLMs. This section is only intended to recall the classical iterative scheme performing maximum likelihood (ML) estimation. Let $\boldsymbol{X}$ denote the $n\times p$ matrix of explanatory variables and $\boldsymbol{\beta}$ the $p$–dimensional parameter vector. At iteration $t+1$, the Fisher Scoring Algorithm (FSA) for ML estimation calculates

[TABLE]

where ${\boldsymbol{z}}^{[t]}$ and ${\boldsymbol{W}}^{[t]}$ respectively denote the classical working variable and the associated weight matrix at iteration $t$. As pointed out by Nelder and Wedderburn, (1972), update (1) may be interpreted as a weighted least squares step in the linearised model $\mathcal{M}^{[t]}$ defined by

[TABLE]

3.2 Linear predictors for SCGLR with multiple responses

We are now considering a multivariate GLM (Fahrmeir and Tutz,, 1994). In this context, SCGLR searches for a component common to all the $\boldsymbol{y_{k}}$’s. This component will be denoted $\boldsymbol{f}$ and its $p$–dimensional loading–vector will be denoted $\boldsymbol{u}$, so that $\boldsymbol{f}=\boldsymbol{Xu}$. The linear predictor associated with response–vector $\boldsymbol{y_{k}}$ then writes

[TABLE]

where $\gamma_{k}$ and $\boldsymbol{\delta_{k}}$ are the regression parameters associated respectively with component $\boldsymbol{f}$ and additional explanatory variables $\boldsymbol{A}$. $\boldsymbol{f}$ being common to all the $\boldsymbol{y_{k}}$’s, predictors are collinear in their $\boldsymbol{X}$–part. For identification purposes, we impose ${\boldsymbol{u}}^{\mbox{\tiny T}}\boldsymbol{M}^{-1}\boldsymbol{u}=1$, where $\boldsymbol{M}$ may so far be any $p\times p$ symmetric positive definite matrix. Let us note $y_{k,i}$ the $i$–th observation of the $k$–th response–vector and $\boldsymbol{H}=\left\{\eta_{k,i}\mid 1\leqslant k\leqslant q,1\leqslant i\leqslant n\right\}$ the predictor set. We assume that the $q$ responses are independent conditional on $\boldsymbol{f}$, and that the $n$ observations are independent. The log–density then writes

[TABLE]

where $\ell_{k}$ is the log–density of the $k$–th response, conditional on its linear predictor. As a result, $\boldsymbol{z_{k}}$ being the working variable associated with $\boldsymbol{y_{k}}$ and $\boldsymbol{W}_{\boldsymbol{k}}^{-1}$ its variance matrix, the corresponding linearised model derived from the FSA at iteration $t$ is

[TABLE]

Although linearised models (2) and (4) seem very similar, (4) is no longer linear, owing to the product $\boldsymbol{u}\gamma_{k}$. An alternate version of the FSA must therefore be used:

• (i)

Given current values of all the $\gamma_{k}$’s and $\boldsymbol{\delta_{k}}$’s, a new loading–vector $\boldsymbol{u}$ is obtained by solving an SCGLR–specific program (see Section 3.3 for details).

• (ii)

Given a current value of $\boldsymbol{u}$, each $\boldsymbol{z}_{\boldsymbol{k}}$ is regressed independently on \big{[}\,\boldsymbol{Xu}\mid\boldsymbol{A}\,\big{]} with respect to weight matrix $\boldsymbol{W_{k}}$, yielding new regression parameters $\gamma_{k}$ and $\boldsymbol{\delta_{k}}$.

3.3 Calculating the component maximising an SCGLR–specific criterion

For an easier reading of this part, we omit the $[t]$ index. For each $k\in\left\{1,\ldots,q\right\}$, consider model $\mathcal{M}_{k}$ endowed with weight matrix $\boldsymbol{W_{k}}$. As suggested in Bry et al., (2013), the best loading–vector in the weighted least–squares sense would be the solution of

[TABLE]

The maximisation program also writes   $\underset{\boldsymbol{u}:\,\boldsymbol{u}^{\mbox{\tiny T}}\boldsymbol{M}^{-1}\boldsymbol{u}=1}{\max}\;\,\psi_{\boldsymbol{A}}\left(\boldsymbol{u}\right)$, where

[TABLE]

Now, $\psi_{\boldsymbol{A}}$ is a mere goodness–of–fit (GoF) measure that does not take into account the closeness of component $\boldsymbol{f}=\boldsymbol{Xu}$ to interpretable directions in $\boldsymbol{X}$. The GoF measure, $\psi_{\boldsymbol{A}}$, must therefore be combined with a measure $\phi$ of structural relevance (SR).

Assume matrix $\boldsymbol{X}$ consists of $p$ standardised numeric variables. Consider a weight system $\boldsymbol{\omega}=\left\{\omega_{1},\ldots,\omega_{p}\right\}$ — e.g. $\omega_{j}=\frac{1}{p}\;\forall j\in\left\{1,\ldots,p\right\}$ — reflecting the a priori relative importance of variables. Also consider a weight matrix $\boldsymbol{P}$ — e.g. $\boldsymbol{P}=\frac{1}{n}\boldsymbol{I}_{n}$ — reflecting the a priori relative importance of observations. We define the most structurally relevant loading–vector as the solution of

[TABLE]

where

[TABLE]

for the scalar product is commutative. Formula (6) is in fact a particular case of the SR criterion proposed by Bry and Verron, (2015); Bry et al., (2016). It can be viewed as a generalised average version of the usual dual PCA criterion: $\sum_{j=1}^{p}\cos^{2}_{\boldsymbol{P}}\left(\boldsymbol{Xu},\boldsymbol{x_{j}}\right)=\sum_{j=1}^{p}\left\langle\,\boldsymbol{Xu}\,\middle|\,\boldsymbol{x_{j}}\,\right\rangle_{\boldsymbol{P}}^{2}$. For $\boldsymbol{M}=\left(\boldsymbol{X}^{\mbox{\tiny T}}\boldsymbol{PX}\right)^{-1}$, (6) is called “Variable–Powered Inertia” (VPI). It should be stressed that for $\boldsymbol{X}^{\mbox{\tiny T}}\boldsymbol{PX}$ to be invertible, $\boldsymbol{X}$ must be a column full rank matrix. In case of strict collinearities within $\boldsymbol{X}$, as it always happens in high–dimensional settings, we replace $\boldsymbol{X}$ with the matrix $\boldsymbol{C}$ of its principal components associated with non–zero eigenvalues. The component is then sought as $\boldsymbol{f}=\boldsymbol{Cu}$. We have $\boldsymbol{C}=\boldsymbol{XV}$, where $\boldsymbol{V}$ is the matrix of corresponding unit-eigenvectors. Then, $\boldsymbol{f}=\boldsymbol{Cu}=\boldsymbol{X\widetilde{u}}$ with $\boldsymbol{\widetilde{u}}=\boldsymbol{Vu}$. Bry et al., (2018) show that among all loading–vectors $\boldsymbol{t}$ such that $\boldsymbol{Xt}=\boldsymbol{f}$, $\boldsymbol{\widetilde{u}}$ is that which has the minimum $L_{2}$–norm.

Tuning parameter $l$ allows to draw component towards more (greater $l$) or less (smaller $l$) local bundles of correlated variables, as depicted on Figure 1 in the particular instance of four coplanar variables. Informally, a bundle is a set of variables correlated “enough” to be viewed as proxies to the same latent dimension. The notion of bundle is flexible, and parameter $l$ tunes the level of within–bundle correlation to be considered: the higher the correlation, the more local the bundle. Overall, taking $l=1$ draws the components towards global structural directions (namely the principal components) while taking $l$ higher leads to more local ones (ultimately, the variables themselves). The goal is to focus on the most interpretable directions.

Finally, let $s\in[0,1]$ be a parameter tuning the importance of the SR relative to the GoF. SCGLR attempts a trade–off between (5) and (6) by solving

[TABLE]

or equivalently

[TABLE]

More detail can be found in Bry et al., (2018).

4 Extension to mixed models

We now propose to extend SCGLR to mixed models. This extension will be called “mixed–SCGLR”. A particular focus is placed on grouped data, for which the independence assumption of observations is no longer valid. The within–group dependence of each response is modelled with a random group–effect. Consequently, each $\boldsymbol{y_{k}}$ is modelled with a GLMM. As in SCGLR, the responses are assumed to be independent conditional on the components. Section 4.1 presents the single–component mixed–SCGLR method. The underlying algorithm is given in Section 4.2. Considering only one component is generally not enough to explain the responses making it necessary to search for $K$ explanatory components, with $1\leqslant K\leqslant\text{rank}\left({\boldsymbol{X}}\right)$. The way in which we extract higher rank components is explained in Section 4.3.

4.1 First component

The random group–effect is assumed different across responses. This leads to $q$ random–effect vectors $\boldsymbol{\xi_{1}},\ldots,\boldsymbol{\xi_{q}}$, which are assumed independent and normally distributed:

[TABLE]

where $N$ denotes the number of groups. In this paper, variance components models will be considered. We assume $\boldsymbol{D_{k}}=\sigma_{k}^{2}\,\boldsymbol{I}_{N}$, where $\sigma_{k}^{2}$ is the group variance component associated with response $\boldsymbol{y_{k}}$. Linear predictors involved in mixed–SCGLR are expressed as

[TABLE]

where $\boldsymbol{U}$ is the known random effects’ design matrix. Predictor $\boldsymbol{\eta^{\xi}_{k}}$ epitomises the way we capture the dependence between outcomes. Indeed, as component $\boldsymbol{f}=\boldsymbol{Xu}$ does not depend on $k$, it captures a structural dependence between the various $\boldsymbol{y_{k}}$’s. By contrast, the random effect $\boldsymbol{\xi_{k}}$ models the within–group stochastic dependence of outcomes forming response–vector $\boldsymbol{y_{k}}$.

Recall that the distribution of the data conditional on the random effects is supposed to belong to the exponential family. The FSA was adapted by Schall, (1991) to the GLMM dependence structure. The key idea is to extend Schall’s algorithm to the component–based predictors in (8).

4.1.1 Linearisation step

Let $g_{k}$ denote the link function for response $\boldsymbol{y_{k}}$, $g^{\prime}_{k}$ its first derivative and $\boldsymbol{\mu^{\xi}_{k}}$ the conditional expectation (i.e. $\boldsymbol{\mu^{\xi}_{k}}\vcentcolon=\mathbb{E}\left(\boldsymbol{y_{k}}\,\lvert\,\boldsymbol{\xi_{k}}\right)$). The working variable associated with $y_{k,i}$ is calculated through

[TABLE]

In view of the conditional independence assumption, the conditional variance matrix for $\boldsymbol{z^{\xi}_{k}}$ is

[TABLE]

where $a_{k,i}$ and $v_{k}$ are known functions, and $\phi_{k}$ is the dispersion parameter related to $\boldsymbol{y_{k}}$. At iteration $t$, the conditional linearised model for working vector $\boldsymbol{z_{k}^{\xi}}$ is then defined by

[TABLE]

Besides the variance component estimation, an alternated estimation step has to be developed (as aforementioned in Section 3.2) to deal with the non–linearity of (9).

4.1.2 Estimation step

Calculating the component:

Given current values of all the $\gamma_{k}$’s, $\boldsymbol{\delta_{k}}$’s, $\boldsymbol{\xi_{k}}$’s and $\sigma_{k}^{2}$’s, a new component $\boldsymbol{f}=\boldsymbol{Xu}$ is calculated by solving a (7)–type program. However, (5) has to be adapted to conditional linearised models $\mathcal{M}^{\xi}_{k}$’s, involving weight matrices $\boldsymbol{W^{\xi}_{k}}$’s. The appropriate goodness–of–fit measure is

[TABLE]

Estimating the regression parameters and variance–components:

Given a current value of component $\boldsymbol{f}$, we apply Schall’s method with the linear predictors given in (8). New values of parameters $\gamma_{k}$ and $\boldsymbol{\delta_{k}}$ as well as new prediction $\boldsymbol{\xi_{k}}$ are obtained by solving the following Henderson system (Henderson,, 1975):

[TABLE]

Finally, as mentioned by Schall, (1991), given prediction $\widehat{\boldsymbol{\xi_{k}}}$ for $\boldsymbol{\xi_{k}}$, the update of the ML estimation of variance component $\sigma_{k}^{2}$ is

[TABLE]

4.2 The algorithm

The conditional linearised models considered at iteration $t$ are given by (9). Algorithm 1 describes the $(t+1)$–th iteration of the single–component mixed–SCGLR. It is repeated until stability of parameters is reached.

Step 1:

Computing the component. Set

$\displaystyle\boldsymbol{u}^{[t+1]}=\underset{\boldsymbol{u}:\,\boldsymbol{u}^{\mbox{\tiny T}}\boldsymbol{M}^{-1}\boldsymbol{u}=1}{\mbox{arg max}}\;\left[\phi\left(\boldsymbol{u}\right)\right]^{s}\>\left[\psi_{\boldsymbol{A}}^{[t]}\left(\boldsymbol{u}\right)\right]^{1-s},\;\text{where }\psi_{\boldsymbol{A}}\left(\boldsymbol{u}\right)\,\text{is given by }(\ref{Chauvet::psi_modified})\text{ and }\phi\left(\boldsymbol{u}\right)\text{ by }(\ref{Chauvet::SR_criterion_particulier})$

$\displaystyle\boldsymbol{f}^{[t+1]}=\boldsymbol{X}\boldsymbol{u}^{[t+1]}$

Step 2:

Henderson systems. For each $k\in\left\{1,\ldots,q\right\}$, solve the system

$\begin{pmatrix}{\boldsymbol{f}^{[t+1]}}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}{\boldsymbol{f}^{[t+1]}}&{\boldsymbol{f}^{[t+1]}}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\boldsymbol{A}&{\boldsymbol{f}^{[t+1]}}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\boldsymbol{U}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \boldsymbol{A}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}{\boldsymbol{f}^{[t]}}&\boldsymbol{A}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\boldsymbol{A}&\boldsymbol{A}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\boldsymbol{U}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \boldsymbol{U}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}{\boldsymbol{f}^{[t]}}&\boldsymbol{U}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\boldsymbol{A}&\boldsymbol{U}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\boldsymbol{U}+{\boldsymbol{D}_{\boldsymbol{k}}^{[t]}}^{-1}\end{pmatrix}\begin{pmatrix}\gamma_{k}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \boldsymbol{\delta_{k}}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \boldsymbol{\xi_{k}}\end{pmatrix}=\begin{pmatrix}{\boldsymbol{f}^{[t+1]}}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\,\boldsymbol{z^{\xi}_{k}}^{[t]}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \boldsymbol{A}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\,\boldsymbol{z^{\xi}_{k}}^{[t]}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \boldsymbol{U}^{\mbox{\tiny T}}{\boldsymbol{W^{\xi}_{k}}}^{[t]}\,\boldsymbol{z^{\xi}_{k}}^{[t]}\end{pmatrix}$

Call $\gamma_{k}^{[t+1]}$, $\boldsymbol{\delta}_{\boldsymbol{k}}^{[t+1]}$ and $\boldsymbol{\xi}_{\boldsymbol{k}}^{[t+1]}$ the solutions.

Step 3:

Updating variance–component estimates. For each $k\in\left\{1,\ldots,q\right\}$, compute

${\sigma_{k}^{2}}^{[t+1]}=\dfrac{{\boldsymbol{\xi}_{\boldsymbol{k}}^{[t+1]}}^{\mbox{\tiny T}}\;\boldsymbol{\xi}_{\boldsymbol{k}}^{[t+1]}}{N-\frac{1}{{\sigma_{k}^{2}}^{[t]}}\,\text{Trace}\left[\left(\boldsymbol{U}^{\mbox{\tiny T}}\,{\boldsymbol{W^{\xi}_{k}}}^{[t]}\,\boldsymbol{U}+{\boldsymbol{D}_{\boldsymbol{k}}^{[t]}}^{-1}\right)^{-1}\right]}\quad\text{and}\quad\boldsymbol{D}_{\boldsymbol{k}}^{[t+1]}={\sigma_{k}^{2}}^{[t+1]}\boldsymbol{I}_{N}$

Step 4:

Updating working variables and weighting matrices.

For each $k\in\left\{1,\ldots,q\right\}$, compute

$\displaystyle{\boldsymbol{\eta^{\xi}_{k}}}^{[t+1]}$ $\displaystyle=\boldsymbol{f}^{[t+1]}\gamma_{k}^{[t+1]}+\boldsymbol{A}\boldsymbol{\delta}_{\boldsymbol{k}}^{[t+1]}+\boldsymbol{U}\boldsymbol{\xi}_{\boldsymbol{k}}^{[t+1]}$

$\displaystyle{\mu_{k,i}^{\xi}}^{[t+1]}$ $\displaystyle=g_{k}^{-1}\left({\eta_{k,i}^{\xi}}^{[t+1]}\right),\;i=1,\ldots,n$

$\displaystyle{z_{k,i}^{\xi}}^{[t+1]}$ $\displaystyle={\eta_{k,i}^{\xi}}^{[t+1]}+\left(y_{i}^{k}-{\mu_{k,i}^{\xi}}^{[t+1]}\right)g^{\prime}_{k}\left({\mu_{k,i}^{\xi}}^{[t+1]}\right),\;i=1,\ldots,n$

$\displaystyle{\boldsymbol{W^{\xi}_{k}}}^{[t+1]}$ $\displaystyle=\textbf{Diag}\left(\left\{\left[g^{\prime}_{k}\left({\mu_{k,i}^{\xi}}^{[t+1]}\right)\right]^{2}a_{k,i}(\phi_{k})\,v_{k}\left({\mu_{k,i}^{\xi}}^{[t+1]}\right)\right\}^{-1}\right)_{i=1,\ldots,n}$

Incrementing:

$t\longleftarrow t+1$

Algorithm 1 Iteration of the single–component mixed–SCGLR

4.3 Extracting higher rank components

Let {\boldsymbol{F_{h}}}=\big{[}\,\boldsymbol{f_{1}}\mid\ldots\mid\boldsymbol{f_{h}}\,\big{]} be the matrix of the first $h$ components, where $h<K$. An extra component $\boldsymbol{f_{h+1}}$ must best complement the existing ones plus $\boldsymbol{A}$, i.e. \boldsymbol{A_{h}}=\big{[}\,\boldsymbol{F_{h}}\mid\boldsymbol{A}\,\big{]}. So $\boldsymbol{f_{h+1}}$ must be calculated using $\boldsymbol{A_{h}}$ as additional explanatory variables. Moreover, we must impose that $\boldsymbol{f_{h+1}}$ be orthogonal to $\boldsymbol{F_{h}}$, i.e. $\boldsymbol{F}_{\boldsymbol{h}}^{\mbox{\tiny T}}\boldsymbol{P}\boldsymbol{f_{h+1}}=\boldsymbol{0}$. Component $\boldsymbol{f_{h+1}}=\boldsymbol{X}\boldsymbol{u_{h+1}}$ is thus obtained by solving

[TABLE]

where $\psi_{\boldsymbol{A_{h}}}\left(\boldsymbol{u}\right)=\displaystyle{\sum_{k=1}^{q}}\left\lVert\boldsymbol{z^{\xi}_{k}}\right\rVert_{\boldsymbol{W^{\xi}_{k}}}^{2}\cos_{\boldsymbol{W^{\xi}_{k}}}^{2}\left(\boldsymbol{z^{\xi}_{k}},\,\operatorname{span}\left\{\boldsymbol{Xu},\boldsymbol{A_{h}}\right\}\right)$ and $\boldsymbol{D_{h}}=\boldsymbol{X}^{\mbox{\tiny T}}\boldsymbol{PF_{h}}$.

In the online Supplementary Material, we give a simple tool to maximise, at least locally, any criterion on the unit sphere: the Projected Iterated Normed Gradient (PING) algorithm. In particular, PING solves (11)–type programs, which give all components of rank $h>1$. The rank–one component is computed using the same program with $\boldsymbol{A_{0}}=\boldsymbol{A}$ and $\boldsymbol{D_{0}}=\boldsymbol{0}$.

5 Comparative results on simulated data

Five simulation studies have been implemented to assess our method. The first one (discussed in Section 5.1 – 5.4) focuses on LMMs. It compares the performances of mixed–SCGLR, LMM–ridge (Eliot et al.,, 2011) and GLMM–LASSO (Groll and Tutz,, 2014; Schelldorfer et al.,, 2014). The second simulation (Section 5.5) extends the first one to binary and Poisson outcomes. All simulation studies have been performed using R (R Core Team,, 2017). To compute LASSO regressions, we have used the R package glmmLasso (Groll,, 2017). The extension of SCGLR to mixed models is available at https://github.com/SCnext/mixedSCGLR. Three additional simulations are presented in the online Supplementary Material. The first one reproduces the simulation scheme of Section 5.5 with binomial and Poisson outcomes. The second one assesses the performance of mixed–SCGLR on a different bundle structure and presents results concerning variance component estimates. The third one deals with high dimensional data.

5.1 Data generation

To generate grouped data, we consider $N=10$ groups and $R=10$ observations per group (i.e. a total of $n=100$ observations). The random effects’ design matrix is then $\boldsymbol{U}=\boldsymbol{I}_{N}\otimes\mathbf{1}_{R}$. Explanatory variables $\boldsymbol{X}$ consist of three independent bundles: $\boldsymbol{X_{0}}$ (15 variables), $\boldsymbol{X_{1}}$ (10 variables) and $\boldsymbol{X_{2}}$ (5 variables). Each explanatory variable is normally simulated with mean [math] and variance $1$. Parameter $\tau\in\left\{0.1,0.3,0.5,0.7,0.9\right\}$ tunes the level of redundancy within each bundle: the correlation matrix of bundle $\boldsymbol{X_{j}}$ is

[TABLE]

where $p_{j}$ is the number of variables in $\boldsymbol{X_{j}}$. In order to enable comparison with LASSO and ridge and to focus on regularisation, our simulations do not involve additional explanatory variables ($\boldsymbol{A}=\boldsymbol{0}$). Two random responses \boldsymbol{Y}=\big{[}\,\boldsymbol{y_{1}}\mid\boldsymbol{y_{2}}\,\big{]} are generated as

[TABLE]

such that $\boldsymbol{y_{1}}$ is predicted only by bundle $\boldsymbol{X_{1}}$, $\boldsymbol{y_{2}}$ only by bundle $\boldsymbol{X_{2}}$, and bundle $\boldsymbol{X_{0}}$ plays no explanatory role. Our choice for the fixed–effect parameters is

[TABLE]

Finally, for each $k\in\left\{1,2\right\}$, random effect and noise vectors are simulated respectively from

[TABLE]

where $\sigma_{k}^{2}=\omega_{k}^{2}=1$. For each value of $\tau$, $B=100$ samples are generated according to Model (12).

5.2 Parameter calibration

In order to compare mixed–SCGLR with the ridge and LASSO regressions, we recall the regularisation parameters required by each method. For both LMM–ridge and GLMM–LASSO methods, a unique shrinkage parameter has to be calibrated: $\lambda_{\text{ridge}}$ and $\lambda_{\text{LASSO}}$ respectively. For mixed–SCGLR, three tuning parameters need to be calibrated: the number of components $K$ and the trade–off parameter $s$, which are both regularisation parameters, and the bundle–locality parameter $l$. For greater clarity, the simulation focuses on the behaviour of $K$ and $s$. As recommended by Bry et al., (2013), we set $l=4$. In case–studies, $l$ has to be tuned to maximise the interpretability of components.

For both mixed–SCGLR and GLMM–LASSO, optimal regularisation parameters are obtained through a $5$–fold cross–validation, withdrawing $2$ observations from each group every time. This could be termed “leave–two–observations–out per group.” The data are thus divided into five parts $\mathcal{P}_{1},\ldots,\mathcal{P}_{5}$, each $\mathcal{P}_{j}$ containing $20$ observations, $2$ for each of the $10$ groups. Let $y_{k,i}^{(b)}$ be the $i$–th observation of the $k$–th response vector in the $b$–th sample. Let also $\widehat{y^{(b)}_{k,i(-j)}}$ be the fit for $y_{k,i}^{(b)}$ with part $\mathcal{P}_{j}$ removed. The cross–validation error in the $b$–th sample, $E^{(b)}$, is defined as

[TABLE]

where

[TABLE]

In the $b$–th sample, the optimal number of components ${K^{\star}}^{(b)}$, the trade–off parameter ${s^{\star}}^{(b)}$, and the shrinkage parameter $\lambda_{\text{LASSO}}^{\star^{{(b)}}}$ are selected to minimise the cross–validation error (13). We then define

[TABLE]

By contrast, Eliot et al., (2011) suggest to calibrate the ridge parameter at each step of their EM implementation, using the generalised cross–validation. We thus define

[TABLE]

where $\lambda_{\text{ridge}}^{\star^{{(b)}}}$ denotes the average of the ridge parameter values obtained over all the iterations of the EM algorithm in the $b$–th sample.

Table 1 summarises the optimal regularisation parameters selected through cross– validation. In both ridge and LASSO, the shrinkage parameter value increases with the level of redundancy $\tau$. Whereas for mixed-SCGLR, when $\tau$ increases, $K^{\star}$ decreases towards the true number of predictive variable–bundles: the greater the value of $\tau$, the better mixed–SCGLR focuses on the structures in $\boldsymbol{X}$ that contribute to model $\boldsymbol{Y}$. Moreover, when $\tau$ increases, the trade–off parameter $s^{\star}$ increases, meaning that regularisation requires a greater importance of the structural relevance relative to the goodness–of–fit.

5.3 Comparison of the estimate accuracies

Once tuning parameters are obtained, we focus on the fixed–effect estimates’ accuracy. Since the response–vectors $\boldsymbol{y_{1}}$ and $\boldsymbol{y_{2}}$ are normally distributed and have comparable orders of magnitude, the fixed–effect relative errors are on the same scale. Then we consider a risk–averse comparison criterion called “Mean Upper Relative Squared Error” (MURSE) defined as

[TABLE]

where $\widehat{\boldsymbol{\beta}}_{\boldsymbol{k}}^{(b)}$ is the estimate of $\boldsymbol{\beta_{k}}$ associated with sample $b$. The MURSE values for mixed–SCGLR, LMM–ridge and GLMM–LASSO are presented in Table 2. The LMM results obtained without regularisation are also presented. They were computed using the R package lme4 (Bates et al.,, 2015). In the latter case, relative errors increase dramatically with $\tau$. Those of ridge and LASSO increase less drastically (but increase anyway) because these methods suffer from the high correlations among the explanatory variables. Except for $\tau=0.1$, mixed–SCGLR provides the most accurate fixed effect estimates. Indeed, if there are no real bundles in $\boldsymbol{X}$ ($\tau\simeq 0$), searching for structures in $\boldsymbol{X}$ may lead mixed–SCGLR to be slightly less accurate. Conversely, mixed–SCGLR takes advantage of the high correlations among the explanatory variables: the stronger the structures (high $\tau$), the more efficient the method.

5.4 Model interpretation

This section aims at highlighting the power of mixed–SCGLR for model interpretation. Figure 2 presents an example of the first component planes obtained for $\tau=0.5$, with associated optimal parameter values $s^{\star}=0.7$ and $K^{\star}=3$. We still impose $l=4$. The first two components obtained are the ones which explain the responses. It clearly appears that $\boldsymbol{y_{1}}$ is explained by bundle $\boldsymbol{X_{1}}$ and $\boldsymbol{y_{2}}$ by $\boldsymbol{X_{2}}$. Interestingly, although bundle $\boldsymbol{X_{0}}$ is the one with maximum inertia ($26.83\%$), it appears only along the third component, for having no explanatory part.

5.5 Additional simulations involving non–Gaussian outcomes

This section aims at assessing our method in the case of Bernoulli ($\mathcal{B}$) and Poisson ($\mathcal{P}$) distributions of responses. We still consider $N=10$ groups and $R=10$ observations per group. We keep design matrices $\boldsymbol{X}$ and $\boldsymbol{U}$ defined in Section 5.1, as well as the values of $\boldsymbol{\beta_{1}}$, $\boldsymbol{\beta_{2}}$, $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$. The group variance components are given by $\varsigma_{1}^{2}=0.1\sigma_{1}^{2}$ and $\varsigma_{2}^{2}=\sigma_{2}^{2}$ so that for each $k\in\left\{1,2\right\}$, \widetilde{\boldsymbol{\xi_{k}}}\mathrel{\overset{}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{N}\left(\boldsymbol{0},\;\varsigma_{k}^{2}\,\boldsymbol{I}_{N}\right). Then given $\widetilde{\boldsymbol{\xi_{1}}}$ and $\widetilde{\boldsymbol{\xi_{2}}}$, we simulate \boldsymbol{Y}=\big{[}\,\boldsymbol{y_{1}}\mid\boldsymbol{y_{2}}\,\big{]} as

[TABLE]

where $\boldsymbol{\theta_{1}}=0.1\boldsymbol{\beta_{1}}$ and $\boldsymbol{\theta_{2}}=\boldsymbol{\beta_{2}}$. Again, for each value of $\tau$, $B=100$ samples are generated according to Model (14). As in Section 5.2, tuning parameters are calibrated so as to minimise the cross–validation error (13). However, since $\boldsymbol{y_{1}}$ and $\boldsymbol{y_{2}}$ do not have the same range of values, the prediction errors have to be standardised. The cross–validation error for response $\boldsymbol{y_{k}}$ in the $b$–th sample is now given by

[TABLE]

Unlike in Section 5.3, the response–vectors do not come from the same distribution and have different orders of magnitude. The fixed–effect relative errors are thus not comparable. To compare mixed–SCGLR with GLMM–LASSO and classical GLMM (without regularisation), we thus use the Mean Relative Squared Error (MRSE) defined as

[TABLE]

where $\widehat{\boldsymbol{\theta}}_{\boldsymbol{k}}^{(b)}$ is the estimate of $\boldsymbol{\theta_{k}}$ from the $b$–th sample. MRSE values for the GLMM, mixed–SCGLR and GLMM–LASSO are presented in Table 3. For all methods, estimating a Bernoulli model is obviously a more challenging task than estimating a Poisson model. Regardless of the level of redundancy $\tau$, both mixed–SCGLR and GLMM–LASSO outperform classical GLMM estimation. Compared with the Gaussian case (Section 5.3), the results deteriorate but (overall) the same behaviours are observed.

• $\blacktriangleright$

For $\tau=0.1$, fixed–effect estimates provided by mixed–SCGLR are less accurate than those provided by GLMM–LASSO. In this case, GLMM–LASSO has indeed a double advantage. First, many $\theta_{k,j}$’s are true zeros. Unlike mixed–SCGLR, GLMM–LASSO often shrinks their estimates to exactly zero. Second, since the level of redundancy is low, GLMM–LASSO also provides accurate coefficient estimates of active variables.

• $\blacktriangleright$

By contrast, for $\tau\geqslant 0.3$, mixed–SCGLR takes advantage of redundancies within the explanatory variables. Thus, mixed-SCGLR outperforms GLMM–LASSO in this case, despite the sparse structure of the $\boldsymbol{\theta_{k}}$’s.

Even if the response variables are not Gaussian, the power of mixed–SCGLR for model interpretation is preserved. Graphical diagnoses similar to those provided in Section 5.4 are available in the Supplementary Material.

6 An application to forest ecology data

6.1 Data description

The present study is based on the Genus dataset of the CoForChange project (see http://www.coforchange.eu). The subsample we consider gives the abundance of $8$ common tree genera on $2615$ Congo Basin land plots. These plots are grouped into $22$ forest concessions. To predict abundances, we have $56$ environmental variables, plus $2$ explanatory variables which code geology and anthropogenic interference. $\boldsymbol{X}$ consists of all environmental variables which are:

• $\blacktriangleright$

$29$ physical factors linked to topography, rainfall or soil moisture,

• $\blacktriangleright$

$25$ photosynthesis activity indicators (the Enhanced Vegetation Indices, EVI, the Near–InfraRed indices, NIR, and the Mid–InfraRed indices, MIR),

• $\blacktriangleright$

$2$ indicators which describe the tree height.

Physical factors are many and redundant: monthly rainfalls are highly correlated, and so are photosynthesis activity indicators. By contrast, geology and anthropogenic interference are weakly correlated and interesting per se. These variables are then considered as additional explanatory variables and included in matrix $\boldsymbol{A}$.

6.2 Model and parameter calibration

Abundances of species given in Genus are count data. For each $k\in\left\{1,\ldots,8\right\}$, we consider a Poisson regression with $\log$ link

[TABLE]

where $\boldsymbol{\xi_{k}}$ is the $22$–level random–effect vector used to model the dependence between the observations of $\boldsymbol{y_{k}}$ within concessions. The first cross–validations we performed — with different fixed values of parameters $s$ and $l$ — indicated that four components were sufficient to capture most of the information in $\boldsymbol{X}$ needed to model and predict responses. We therefore keep $K^{\star}=4$. The optimal values of trade–off and locality parameter $s^{\star}$ and $l^{\star}$ are then determined through another cross–validation. Using the same procedure and notations as in Section 5.2, the data are divided into five parts $\mathcal{P}_{1},\ldots,\mathcal{P}_{5}$. Let $n_{j}$ be the size of $\mathcal{P}_{j}$.

[TABLE]

where

[TABLE]

On Figure 3, we plot the errors $E$ for parameter pairs $(s,l)\in\mathcal{E}_{s}\times\mathcal{E}_{l}$, where

[TABLE]

Parameter grid $\mathcal{E}_{s}\times\mathcal{E}_{l}$ therefore contains $264$ pair values. Selecting the best parameter pair from $\mathcal{E}_{s}\times\mathcal{E}_{l}$ through a $5$–fold cross–validation requires a computation time of about $65$ minutes (parallel computing on $6$ CPU cores, Intel Core i7–6700HQ, 2.6GHz). It should be noted that there is a risk of non–convergence when the trade–off parameter $s$ is too close to [math]. Indeed, if we consider no structural information ($s$ exactly equal to [math]) in $\boldsymbol{X}$, mixed–SCGLR merely performs classical GLMM estimation and does not converge with this data. When $s=0.025$, our algorithm converges but leads to fairly unstable estimates and high cross–validation errors because regularisation is then very weak. By contrast, the components calculated with $s\in\left\{0.5,0.6,\ldots,1\right\}$ are close to principal components. The associated errors are therefore stable in most cases, but rather high. Finally, $s\in\left\{0.1,\ldots,0.4\right\}$ leads to the lowest cross–validation errors, but only for $l\leqslant 10$. Indeed, when $s$ is not too high, mixed–SCGLR may focus on the most predictive structures of $\boldsymbol{X}$. However, parameter $l$ must not exceed a certain value, in order to avoid being drawn towards too local variable–bundles. As can be seen, choosing $\left(s^{\star},l^{\star}\right)=\left(0.1,10\right)$ minimises the cross–validation error.

6.3 Prediction quality and interpretation results

This part evaluates the benefits obtained by taking within–group dependence into account. The predictions we get with mixed–SCGLR and with initial version of SCGLR are compared with respect to the cross–validation criterion given by (15). Table 4 summarises the $E_{k}$’s for both SCGLR and mixed–SCGLR methods. Optimal parameter value triplet $\left(K^{\star},s^{\star},l^{\star}\right)=\left(4,0.1,10\right)$ is selected for both methods. For each $k\in\left\{1,\ldots,8\right\}$, mixed–SCGLR gives a lower cross–validation error than SCGLR: taking into account the within–group dependence has clearly improved prediction performances.

Moreover, mixed–SCGLR enables to correctly reconstitute observed abundance maps, as illustrated on Figure 4.

As has been seen in Section 5.4, mixed–SCGLR allows an easy interpretation of the model through the decomposition of linear predictors on interpretable components. Figure 5 shows the first two component planes resulting from mixed–SCGLR on real data Genus. Component plane $(1,2)$ reveals two patterns. The first one is a global rain–wind pattern driven by the pluvio’s and wd’s variables which explain the abundances of Species $1,2,5,6$. The second is a rather local pattern driven by variables altitude, wetness and annual pluviometry (pluvio$\_$an) which prove important to model and predict responses $\boldsymbol{y_{3}}$ and $\boldsymbol{y_{7}}$. Lastly, Component $3$ reveals a photosynthesis pattern driven by a part of the Evi’s, which seems useful to predict $\boldsymbol{y_{4}}$ and $\boldsymbol{y_{8}}$.

The decomposition of linear predictors on interpretable components allows to detect the species that tend to share common explanatory dimensions and those which are more idiosyncratic. We can then identify the variable–bundles these dimensions are related to. The underlying goal is a better understanding of the bio– and ecosystem diversity with a view to preserve them. Species 1, 2, 5 and 6 are sensitive to the same rain–wind regime, and Species 4 and 8 are explained by the same photosynthetic pattern. On the contrary, Species 3 and 7 are clearly separated. Species 7 grows at high altitudes where the atmosphere is rather dry while the abundance of Species 3 is favoured by regular rainfall and high humidity.

7 Discussion and Conclusions

Like Sufficient Dimension Reduction (SDR) methods, mixed–SCGLR is based on the construction of a reduction function of dimension less than $p$ which tries to capture all the relevant information that $\boldsymbol{X}$ contains about $\boldsymbol{Y}$. However, the two approaches do not exactly pursue the same objectives. Indeed, SDR methods look for the “central subspace” containing the predictive information irrespective of the structures within $\boldsymbol{X}$ (e.g. dimensions capturing a large part of $\boldsymbol{X}$’s variance, or bundles of correlated variables). Mixed–SCGLR rather aims at basing the explanatory subspace on such structural dimensions so as to both gain interpretability and stabilise prediction. We think that extracting a hierarchy of strong and interpretable dimensions, and decomposing the linear predictor on them, is an essential asset in model–building. The difference in goals entails a difference in means: SDR is based on the sufficiency principle, which is enough to identify a subspace but not to track strong predictive dimensions in it. By contrast, in the wake of PLS regression, mixed–SCGLR uses a criterion combining goodness–of–fit and structural relevance of components.

The supervised–component paradigm has proved effective in situations where regularisation is necessary but where variable selection is inappropriate — for instance when the true explanatory dimensions are latent and indirectly measured through highly correlated proxies.

• $\blacktriangleright$

When $l=1$, trade–off parameter $s$ allows to continuously tune the attraction of components towards the principal components of explanatory variables. This results in a continuum between classical GLMM estimation ($s=0$ is associated with no regularisation) and principal component generalised linear mixed regression (with $s=1$).

• $\blacktriangleright$

When $l>1$, we take better advantage of local predictive structures in $\boldsymbol{X}$. The components we build are then usually closer to local gatherings of variables, thus easier to interpret.

Mixed–SCGLR is able to identify more or less local predictive structures common to all the $\boldsymbol{y_{k}}$’s and performs well on grouped data with Gaussian, Bernoulli, binomial and Poisson outcomes. Compared to penalty–based approaches as ridge or LASSO, the orthogonal components built by mixed–SCGLR reveal the multidimensional explanatory and predictive dimensions, and greatly facilitate the interpretation of the model.

However, a natural question arises as to the accuracy of our methodology under significant deviations from normality. With binary data for instance, variance component estimates are prone to some bias towards zero (McCulloch,, 1997). That is why other estimation strategies might be considered, especially Monte Carlo integration methods which have the advantage of being based on direct approximations of the likelihood. Some examples are the MCMC methods developed by Hadfield, (2010) in the GLMM framework, and the Monte Carlo Likelihood Approximation (MCLA) proposed by Knudson, (2016). Indirect maximisations of the likelihood are also available such as Monte Carlo Expectation–Maximisation (MCEM) and Monte Carlo Newton–Raphson (McCulloch,, 1997). We think that these methodologies and Schall’s could be combined sequentially. Indeed we could first take advantage of the linear approximation of the model in order to build the components, and then use MC–based methods to estimate both fixed–effect parameters and variance components. This would lead to replacing the current iteration of mixed–SCGLR given by Algorithm 1 with the following steps (to keep things simple, we take the canonical link):

Compute components \boldsymbol{F}=\big{[}\,\boldsymbol{f_{1}}\mid\ldots\mid\boldsymbol{f_{K}}\,\big{]} via the PING algorithm on Schall’s linearised models.

• 2.

For each $k\in\left\{1,\ldots,q\right\}$, consider the hierarchy

[TABLE]

where $\boldsymbol{\eta^{\xi}_{k}}=\boldsymbol{F\gamma_{k}}+\boldsymbol{A\delta_{k}}+\boldsymbol{U\xi_{k}}$, and $c$, $d$ are the functions associated with the natural parametrisation of the GLM. For example, for the Bernoulli–logistic regression, we have: $c(x)=\log(1+e^{x})$ and $d(x)=0$.

• 3.

Apply MC–based methods such as MCMC, MCLA, MCEM or MCNR to update $\boldsymbol{\gamma_{k}}$, $\boldsymbol{\delta_{k}}$, $\boldsymbol{\xi_{k}}$ and $\boldsymbol{D_{k}}$, $k\in\left\{1,\ldots,q\right\}$.

• 4.

Update working variables and weight matrices to define the new Schall’s linearised models.

Even though such MC–based methods are computationally much more intensive than the “Joint–Maximisation” and have intrinsic disadvantages (particularly in the assessment of convergence and in the choice of prior distributions), they could give better results in case of binary data.

SUPPLEMENTARY MATERIAL

Additional simulations:

The first simulation reproduces that of Section 5.5 in the case of binomial and Poisson outcomes. The second simulation explores a different structure of variable–bundles, considers Gaussian, binomial and Poisson outcomes, and presents results concerning variance component estimates. The third one involves high dimensional data. (pdf file)

Projected Iterated Normed Gradient (PING) algorithm:

We give some technical details about the PING algorithm, which maximises, at least locally, any criterion on the unit sphere. (pdf file)

R package mixedSCGLR:

We provide an R package to perform mixed–SCGLR, also available at https://github.com/SCnext/mixedSCGLR. It contains the dataset Genus used in Section 6. The package also provides demo codes, in particular for visualising the component planes (mixedSCGLR.tar.gz).

Code for running simulations:

We also provide the R codes required to reproduce most of the simulation results (R and Rdata files).

ACKNOWLEDGMENTS

The extended data Genus required the arrangement and the inventory of 140.000 developed plots across four countries : Central African Republic, Gabon, Cameroon and Democratic Republic of Congo. The authors thank the members of the CoForTips project for allowing the use of this data. We are also grateful to the editor, the associate editor and to the referees for their thorough and constructive review of this work.

**Supplementary file to “Component–based regularisation of multivariate generalised linear mixed models”:

THE PING ALGORITHM**

Jocelyn Chauvet1  Catherine Trottier1,2  Xavier Bry1

1 IMAG, Univ Montpellier, CNRS, Montpellier, France.

jocelyn.chauvet@umontpellier.fr ; xavier.bry@umontpellier.fr

2 Univ Paul-Valéry Montpellier 3, Montpellier, France.

catherine.trottier@univ-montp3.fr

The Projected Iterated Normed Gradient (PING) is a basic extension of the iterated power algorithm, for solving any program of the form

[TABLE]

Note that putting $\boldsymbol{v}=\boldsymbol{M}^{-1/2}\boldsymbol{u}$, $\mathcal{G}_{h}\left(\boldsymbol{v}\right)=\mathcal{J}_{h}\left(\boldsymbol{M}^{1/2}\boldsymbol{v}\right)$ and $\boldsymbol{B_{h}}=\boldsymbol{M}^{1/2}\boldsymbol{\Delta_{h}}$, Program (16) is strictly equivalent to Program (17):

[TABLE]

Denoting

[TABLE]

a Lagrange multiplier-based reasoning gives the basic iteration of the PING algorithm:

[TABLE]

Although Iteration (18) follows a direction of ascent, it does not guarantee that $\mathcal{G}_{h}$ actually increases on every step. We therefore propose a generic iteration of PING (Algorithm 2) and an alternative one (Algorithm 3), which both ensure that the criterion increases.

while convergence of $\boldsymbol{v}$ non reached do

       $\boldsymbol{\kappa}^{[t]}=\dfrac{\Pi_{\operatorname{span}\left\{\boldsymbol{B_{h}}\right\}^{\bot}}\,\Gamma_{h}\left(\boldsymbol{v}^{[t]}\right)}{\left\lVert\Pi_{\operatorname{span}\left\{\boldsymbol{B_{h}}\right\}^{\bot}}\,\Gamma_{h}\left(\boldsymbol{v}^{[t]}\right)\right\rVert}$

      

A unidimensional Newton–Raphson maximisation procedure is used to find

      

the maximum of $\mathcal{G}_{h}\left(\boldsymbol{v}\right)$ on the arc $\left(\boldsymbol{v}^{[t]},\boldsymbol{\kappa}^{[t]}\right)$ and take it as $\boldsymbol{v}^{[t+1]}$.

       $t\longleftarrow t+1$

end while

Algorithm 2 Generic iteration of the PING algorithm

while convergence of $\boldsymbol{v}$ non reached do

       $\boldsymbol{m}\longleftarrow\dfrac{\Pi_{\operatorname{span}\left\{\boldsymbol{B_{h}}\right\}^{\bot}}\,\Gamma_{h}\left(\boldsymbol{v}^{[t]}\right)}{\left\lVert\Pi_{\operatorname{span}\left\{\boldsymbol{B_{h}}\right\}^{\bot}}\,\Gamma_{h}\left(\boldsymbol{v}^{[t]}\right)\right\rVert}$

       while $\mathcal{G}_{h}\left(\boldsymbol{m}\right)<\mathcal{G}_{h}\left(\boldsymbol{v}^{[t]}\right)$ do

             $\boldsymbol{m}\longleftarrow\dfrac{\boldsymbol{v}^{[t]}+\boldsymbol{m}}{\left\lVert\boldsymbol{v}^{[t]}+\boldsymbol{m}\right\rVert}$

       end while

      $\boldsymbol{v}^{[t+1]}=\boldsymbol{m}$

       $t\longleftarrow t+1$

end while

Algorithm 3 Alternative generic iteration of the PING algorithm

First rank component.

Component $\boldsymbol{f_{1}}=\boldsymbol{X}\boldsymbol{u_{1}}$ is obtained by solving

[TABLE]

This corresponds to Program (16) with $h=0$, where

• $\blacktriangleright$

$\mathcal{J}_{0}\left(\boldsymbol{u}\right)=s\log\left[\phi\left(\boldsymbol{u}\right)\right]+(1-s)\log\left[\psi_{\boldsymbol{A_{0}}}\left(\boldsymbol{u}\right)\right]$,

• $\blacktriangleright$

$\boldsymbol{A_{0}}=\boldsymbol{A}$ (the matrix of additional explanatory variables), and

• $\blacktriangleright$

$\boldsymbol{\Delta_{0}}=\boldsymbol{0}$.

In this particular case, we have $\boldsymbol{B_{0}}=\boldsymbol{M}^{1/2}\boldsymbol{\Delta_{0}}=\boldsymbol{0}$, and so:

[TABLE]

Higher rank components.

Let {\boldsymbol{F_{h}}}=\big{[}\,\boldsymbol{f_{1}}\mid\ldots\mid\boldsymbol{f_{h}}\,\big{]} be the matrix of the first $h$ components and \boldsymbol{A_{h}}=\big{[}\,\boldsymbol{F_{h}}\mid\boldsymbol{A}\,\big{]}. Let $\boldsymbol{P}$ denote the weight matrix reflecting the a priori relative importance of observations ($\boldsymbol{P}=\frac{1}{n}\boldsymbol{I}_{n}$ if all observations are of equal importance). Component $\boldsymbol{f_{h+1}}=\boldsymbol{X}\boldsymbol{u_{h+1}}$ is obtained by solving

[TABLE]

This corresponds to Program (16), where

• $\blacktriangleright$

$\mathcal{J}_{h}\left(\boldsymbol{u}\right)=s\log\left[\phi\left(\boldsymbol{u}\right)\right]+(1-s)\log\left[\psi_{\boldsymbol{A_{h}}}\left(\boldsymbol{u}\right)\right]$,

• $\blacktriangleright$

\boldsymbol{A_{h}}=\big{[}\,\boldsymbol{F_{h}}\mid\boldsymbol{A}\,\big{]}, and

• $\blacktriangleright$

$\boldsymbol{\Delta_{h}}=\boldsymbol{X}^{\mbox{\tiny T}}\boldsymbol{P}\boldsymbol{F_{h}}$.

Initialisation.

To quickly find $\boldsymbol{f_{1}}$, algorithm PING is initialised with the first PLS component of the responses on $\boldsymbol{X}$. In like manner, for $h\geqslant 2$, PING is initialised with the first PLS component of the responses on $\boldsymbol{X}$ deflated on components \boldsymbol{F_{h-1}}=\big{[}\,\boldsymbol{f_{1}}\mid\ldots\mid\boldsymbol{f_{h-1}}\,\big{]}.

**Supplementary file to “Component–based regularisation of multivariate generalised linear mixed models”:

ADDITIONAL SIMULATIONS**

Jocelyn Chauvet1  Catherine Trottier1,2  Xavier Bry1

1 IMAG, Univ Montpellier, CNRS, Montpellier, France.

jocelyn.chauvet@umontpellier.fr ; xavier.bry@umontpellier.fr

2 Univ Paul-Valéry Montpellier 3, Montpellier, France.

catherine.trottier@univ-montp3.fr

8 Comparative results with binomial and Poisson outcomes

In this section, we simply extend the simulation scheme presented in Section 5.5 to binomial and Poisson outcomes. We maintain design matrices $\boldsymbol{X}$ and $\boldsymbol{U}$ as defined in Section 5.1. Fixed–effect parameters $\boldsymbol{\theta_{k}}$’s and random–effect vectors $\boldsymbol{\widetilde{\xi_{k}}}$’s are defined in Section 5.5. Given $\boldsymbol{\widetilde{\xi_{1}}}$ and $\boldsymbol{\widetilde{\xi_{2}}}$, we then simulate \boldsymbol{Y}=\big{[}\,\boldsymbol{y_{1}}\mid\boldsymbol{y_{2}}\,\big{]} as

[TABLE]

Table 5 gives the Mean Relative Squared Error (MRSE) values for $\boldsymbol{\theta_{1}}$ and $\boldsymbol{\theta_{2}}$ obtained on 100 samples for each value of $\tau$.

For the Poisson distribution, the results in Table 5 are essentially identical to those in the article: mixed–SCGLR outperforms GLMM–LASSO (Groll,, 2017, R package glmmLasso) except for $\tau=0.1$. As for the binomial distribution, the regularisation provided by mixed–SCGLR improves the results obtained without regularisation (Bates et al.,, 2015, R package lme4), regardless of the level of redundancy within the explanatory variables. Unsurprisingly, the errors are much smaller than in the binary case.

The power of mixed–SCGLR in terms of model interpretation remains the same for non–Gaussian outcomes. Figure 6 (respectively Figure 7) presents an example of the first component planes output by mixed–SCGLR in the binomial/Poisson (respectively Bernoulli/Poisson) case. As for Gaussian outcomes, the component planes reveal that $\boldsymbol{y_{1}}$ is explained by bundle $\boldsymbol{X_{1}}$ and $\boldsymbol{y_{2}}$ by $\boldsymbol{X_{2}}$. In the binomial/Poisson case with $\tau=0.3$ (Figure 6), predictive bundles $\boldsymbol{X_{1}}$ and $\boldsymbol{X_{2}}$ are captured respectively by the first and the second components. The third component aligns on nuisance bundle $\boldsymbol{X_{0}}$, despite its high inertia. Figure 7 illustrates what may happen when the level of redundancy is very high ($\tau=0.7$ here). Since the explanatory variables are highly correlated, mixed–SCGLR regularisation requires that the structural relevance be given a heavy weight with respect to the goodness–of–fit, which leads to a trade–off parameter $s$ close to 1 ($s=0.9$ here). Having the greatest structural strength, the nuisance bundle is captured by the second component despite its lack of explanatory power. This is sometimes the price to be paid for the trade–off. In our example, the second explanatory bundle is captured by the third component, so that the predictive dimensions are accurately represented in component plane $(1,3)$.

9 A new structure for the bundles — Gaussian, Poisson and binomial distributions

This simulation study tests mixed-SCGLR on a slightly more complex bundle structure. Results concerning variance component estimates are also presented.

We consider a fixed–effect design matrix $\boldsymbol{X}_{n\times p}$ partitioned into 3 blocks $\boldsymbol{\mathcal{G}_{1}}$, $\boldsymbol{\mathcal{G}_{2}}$ and $\boldsymbol{\mathcal{G}_{3}}$. Block $\boldsymbol{\mathcal{G}_{1}}$ contains 10 predictive explanatory variables structured about a latent variable \boldsymbol{\varphi_{1}}\mathrel{\overset{}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{n}\left(\boldsymbol{0},\sigma_{\text{LV}}^{2}\boldsymbol{I}_{n}\right). Thus for each $j\in\left\{1,\ldots,10\right\}$, $\boldsymbol{x_{j}}=\boldsymbol{\varphi_{1}}+\boldsymbol{\varepsilon_{j}}$, where \boldsymbol{\varepsilon_{j}}\mathrel{\overset{}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{n}\left(\boldsymbol{0},\sigma_{\text{noise}}^{2}\boldsymbol{I}_{n}\right) such that $\sigma_{\text{LV}}^{2}+\sigma_{\text{noise}}^{2}=1$. The correlation within $\boldsymbol{\mathcal{G}_{1}}$ is tuned by signal to noise (StN) ratio $\sigma_{\text{LV}}^{2}/\sigma_{\text{noise}}^{2}$ (chosen in $\left\{\frac{1}{3},1,3\right\}$ in practice). $\boldsymbol{\mathcal{G}_{2}}$ contains a single predictive variable \boldsymbol{\varphi_{2}}=\boldsymbol{x_{11}}\mathrel{\overset{}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{n}\left(\boldsymbol{0},\boldsymbol{I}_{n}\right). $\boldsymbol{\mathcal{G}_{3}}$ contains 20 unstructured noise variables: for each $j\in\left\{12,\ldots,31\right\}$, \boldsymbol{x_{j}}\mathrel{\overset{}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{n}\left(\boldsymbol{0},\boldsymbol{I}_{n}\right). For each $k\in\left\{1,2,3\right\}$, random–effect vectors are simulated as \boldsymbol{\xi_{k}}\mathrel{\overset{\text{ind.}}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{N}\left(\boldsymbol{0},\sigma_{k}^{2}\,\boldsymbol{I}_{N}\right). Given $\boldsymbol{\xi_{1}}$, $\boldsymbol{\xi_{2}}$, $\boldsymbol{\xi_{3}}$, we simulate 3 responses having different distributions, \boldsymbol{Y}=\big{[}\,\boldsymbol{y_{1}}\mid\boldsymbol{y_{2}}\mid\boldsymbol{y_{3}}\,\big{]}, as

[TABLE]

In our simulations, we set $\alpha_{1}=\sigma_{1}^{2}=2$, $\alpha_{2}=\sigma_{2}^{2}=1$ and $\alpha_{3}=\sigma_{3}^{2}=0.5$.

We consider in turn $N=10$ and $N=50$ groups, and $R=10$ observations per group ($n=100$ and $n=500$ observations in total). $B=100$ samples are generated for each pair of values $\left(N,\text{StN}\right)$. The main goal of the study is to assess the ability of mixed–SCGLR to track down both latent variable $\boldsymbol{\varphi_{1}}$ and predictive variable $\boldsymbol{\varphi_{2}}$. For $j=1$ and $2$, we then define

[TABLE]

where $\boldsymbol{f}_{\boldsymbol{j}}^{(b)}$ is the component most correlated with $\boldsymbol{\varphi_{j}}$ issued from mixed–SCGLR in the $b$–th sample. Consistency of fixed–effect estimates is assessed through criteria $\text{err}_{1}$, $\text{err}_{2}$ and $\text{err}_{3}$ defined by

[TABLE]

where $\widehat{\boldsymbol{\beta}}_{\boldsymbol{j}}^{(b)}$ is the fixed–effect estimate related to response $\boldsymbol{y_{j}}$ associated with sample $b$.

Table 6 summarises the values of the afore–defined criteria and presents biases and standard errors of variance components estimates. For a given value of $N$, $\text{cor}_{1}$ increases towards $1$ with ratio $\sigma_{\text{LV}}^{2}/\sigma_{\text{noise}}^{2}$: the tighter the block $\boldsymbol{\mathcal{G}_{1}}$ is structured about its latent variable, the better mixed–SCGLR can reconstruct it. The associated criterion $\text{err}_{1}$ then naturally decreases towards [math]. On the other hand, $\text{cor}_{2}$ and $\text{err}_{2}$ are very stable, which proves that mixed–SCGLR is able to detect an isolated predictive variable among a large number of irrelevant others. As $\text{err}_{3}$ depends on how accurately mixed–SCGLR recovers $\boldsymbol{\varphi_{1}}$ and $\boldsymbol{\varphi_{2}}$, it slightly decreases when the StN ratio increases. Both variance component biases and standard errors seem rather stable regardless of the value of StN. Finally, when $N$ increases, all the $\text{cor}_{j}$’s increase towards 1 and all the $\text{err}_{j}$’s decrease towards 0. As far as variance component estimates are concerned, the biases are getting slightly closer to 0 and the standard errors decrease significantly.

Model interpretation is revealed by Figure 8 in the case of $N=10$ groups and $R=10$ observations per group. The first component aligns with block $\boldsymbol{\mathcal{G}_{1}}$ which alone explains response $\boldsymbol{y_{1}}$. The second aligns with $\boldsymbol{\mathcal{G}_{2}}$ (containing single explanatory variable $\boldsymbol{x_{11}}$) which alone explains $\boldsymbol{y_{2}}$. Finally, note that the projection of the $\boldsymbol{X}$–part of the linear predictor related to $\boldsymbol{y_{3}}$ is well represented on component plane $(1,2)$. This indicates that $\boldsymbol{y_{3}}$ is explained jointly by $\boldsymbol{\mathcal{G}_{1}}$ and $\boldsymbol{\mathcal{G}_{2}}$.

10 High dimensional data

10.1 Key idea

To cope with high dimensional data, the key idea is to replace the fixed–effect design matrix, $\boldsymbol{X}$, with the matrix $\boldsymbol{C}$ of its principal components associated with non–zero eigenvalues. More precisely, $\lambda_{j}$ being the eigenvalue associated with the $j$–th eigenvector $\boldsymbol{v_{j}}$, the last eigenvector we consider, $\boldsymbol{v_{r}}$, is such that

[TABLE]

where $p$ is the number of columns of matrix $\boldsymbol{X}$. The matrix of the corresponding unit-eigenvectors is denoted \boldsymbol{V}=\big{[}\,\boldsymbol{v_{1}}\mid\ldots\mid\boldsymbol{v_{r}}\,\big{]}, and $\boldsymbol{C}=\boldsymbol{XV}$. The component $\boldsymbol{f}$ is then sought as a combination of the principal components: $\boldsymbol{f}=\boldsymbol{Cu}=\boldsymbol{X\widetilde{u}}$, where $\boldsymbol{\widetilde{u}}=\boldsymbol{Vu}$. Mixed–SCGLR then solves

[TABLE]

where the goodness–of–fit measure, $\psi_{\boldsymbol{A}}$, is given by

[TABLE]

and the structural relevance by

[TABLE]

This idea is tested on simulated data where the number of explanatory variables $p$ exceeds the number of observations $n$.

10.2 Data generation

To generate grouped data, we consider $N=10$ groups and $R=10$ observations per group (i.e. a total of $n=100$ observations). The random effects’ design matrix is then $\boldsymbol{U}=\boldsymbol{I}_{N}\otimes\mathbf{1}_{R}$. Explanatory variables consist of four independent bundles $\boldsymbol{X_{j}},j\in\left\{1,2,3,4\right\}$, such as \boldsymbol{X}=\big{[}\,\boldsymbol{X_{0}}\mid\boldsymbol{X_{1}}\mid\boldsymbol{X_{2}}\mid\boldsymbol{X_{3}}\,\big{]}. Each explanatory variable is normally simulated with mean [math] and variance $1$. Parameter $\tau\in\left\{0.3,0.5,0.7\right\}$ tunes the level of redundancy within each bundle: the correlation matrix of bundle $\boldsymbol{X_{j}}$ is

[TABLE]

where $p_{j}$ is the number of variables in $\boldsymbol{X_{j}}$. For each $k\in\left\{1,2,3,4\right\}$, random–effect vectors are simulated as \boldsymbol{\xi_{k}}\mathrel{\overset{\text{ind.}}{\scalebox{1.5}[1.0]{\sim}}}\mathcal{N}_{N}\left(\boldsymbol{0},\sigma_{k}^{2}\,\boldsymbol{I}_{N}\right). Given $\boldsymbol{\xi_{1}}$, $\boldsymbol{\xi_{2}}$, $\boldsymbol{\xi_{3}}$, $\boldsymbol{\xi_{4}}$, we simulate 4 responses having different distributions, \boldsymbol{Y}=\big{[}\,\boldsymbol{y_{1}}\mid\boldsymbol{y_{2}}\mid\boldsymbol{y_{3}}\mid\boldsymbol{y_{4}}\,\big{]}, as

[TABLE]

Response $\boldsymbol{y_{1}}$ is predicted only by bundle $\boldsymbol{X_{1}}$, $\boldsymbol{y_{2}}$ only by bundle $\boldsymbol{X_{2}}$, $\boldsymbol{y_{3}}$ only by bundle $\boldsymbol{X_{3}}$, $\boldsymbol{y_{4}}$ by both bundles $\boldsymbol{X_{2}}$ and $\boldsymbol{X_{3}}$, and bundle $\boldsymbol{X_{0}}$ plays no explanatory role. Our choice for the fixed–effect parameters is

[TABLE]

We consider in turn $p=150$ ($p_{0}=60$, $p_{1}=45$, $p_{2}=30$, $p_{3}=15$) and $p=200$ ($p_{0}=80$, $p_{1}=60$, $p_{2}=40$, $p_{3}=20$) explanatory variables. Variance components are set to $\sigma_{1}^{2}=\sigma_{2}^{2}=\sigma_{3}^{2}=0.1$, and $\sigma_{4}^{2}=0.05$. For each value of $p$ and for each value of $\tau$, $B=20$ samples are generated according to Model (19).

10.3 Results

Table 7 and Table 8 present the results for respectively $150$ and $200$ explanatory variables. They give the Mean Relative Squared Error (MRSE) values for $\boldsymbol{\beta_{k}},k\in\left\{1,\ldots,4\right\}$, as well as biases and standard errors of estimated variance components, obtained on 20 samples for each value of $\tau$.

Some component planes are given on Figure 9 ($150$ explanatory variables) and Figure 10 ($200$ explanatory variables).