Estimation of variance components, heritability and the ridge penalty in high-dimensional generalized linear models

TL;DR

This paper compares estimators of variance components and heritability in high-dimensional generalized linear models, introducing efficient methods and demonstrating their superior accuracy over existing approaches.

Contribution

It introduces a novel maximum marginal likelihood estimator for ridge penalty in high-dimensional GLMs, extending previous methods and providing computational efficiency.

Findings

MML estimator outperforms CV in accuracy for Poisson and Binomial models.

Robustness of estimators against model departures like sparsity and non-Gaussian errors.

Software implementation enables reproducibility of results.

Abstract

For high-dimensional linear regression models, we review and compare several estimators of variances $\tau^2$ and $\sigma^2$ of the random slopes and errors, respectively. These variances relate directly to ridge regression penalty $\lambda$ and heritability index $h^2$, often used in genetics. Direct and indirect estimators of these, either based on cross-validation (CV) or maximum marginal likelihood (MML), are also discussed. The comparisons include several cases of covariate matrix $\mathbf{X}_{n \times p}$, with $p \gg n$, such as multi-collinear covariates and data-derived ones. In addition, we study robustness against departures from the model such as sparse instead of dense effects and non-Gaussian errors. An example on weight gain data with genomic covariates confirms the good performance of MML compared to CV. Several extensions are presented. First, to the high-dimensional linear mixed effects model, with REML as an alternative to MML. Second, to the conjugate Bayesian setting, which proves to be a good alternative. Third, and most prominently, to generalized linear models for which we derive a computationally efficient MML estimator by re-writing the marginal likelihood as an $n$-dimensional integral. For Poisson and Binomial ridge regression, we demonstrate the superior accuracy of the resulting MML estimator of $\lambda$ as compared to CV. Software is provided to enable reproduction of all results presented here.

Figures

Tables

Table 1. Table 1 : Computing times for hyper-parameter estimation for linear and Poisson ridge regression

Linear	$n=100$			$n=500$
	$p={10}^3$	$p={10}^4$	$p={10}^5$	$p={10}^3$	$p={10}^4$	$p={10}^5$
MML	0.06	0.15	1.12	2.18	6.07	26.64
Bayes	0.04	0.31	4.38	1.10	7.78	93.25
MoM	0.01	0.08	1.03	0.17	2.32	23.70
PCR	0.05	0.39	5.36	1.39	10.31	116.80
Basic	0.05	0.46	6.56	1.44	12.40	145.18
GCV	0.20	0.46	4.56	12.26	26.41	111.38
CV	0.81	6.57	39.95	2.62	21.69	183.50
HiLMM	0.03	0.17	2.01	0.66	3.14	27.99
Poisson	$n=100$			$n=500$
	$p={10}^3$	$p={10}^4$	$p={10}^5$	$p={10}^3$	$p={10}^4$	$p={10}^5$
MML_mgcv	0.32	0.33	0.31	26.21	40.19	48.17
GCV_mgcv	0.39	0.33	0.62	33.48	41.44	54.01

Equations

$${\mathbf{y}}_{n\times 1}={\mathbf{X}}_{n\times p}{\mathbf{\beta }}_{p\times 1}+{\mathbf{ϵ}}_{n\times 1}$$

$${\mathbf{\beta }}_{p\times 1}\sim \mathcal{N}(0,{\tau }^2{\mathbf{I}}_p)$$

$${\mathbf{ϵ}}_{n\times 1}\sim \mathcal{N}(0,{\sigma }^2{\mathbf{I}}_n).$$

$$\begin{array}{rl}\hat{\mathbf{\beta }}& ={\hat{\mathbf{\beta }}}_{\lambda }={\mathbf{C}}_{\lambda }\mathbf{y}=({\mathbf{X}}^T\mathbf{X}+\lambda I_{p\times p}{)}^{-1}{\mathbf{X}}^T\mathbf{y}\text{i.e. the linear ridge estimator}\\ \mathbf{H}& ={\mathbf{H}}_{\lambda }=\mathbf{X}{\mathbf{C}}_{\lambda }=\mathbf{X}({\mathbf{X}}^T\mathbf{X}+\lambda I_{p\times p}{)}^{-1}{\mathbf{X}}^T\text{i.e. the hat matrix.}\end{array}$$

$$\mathbf{C}=({\mathbf{X}}^T\mathbf{X}+{\mathbf{\Lambda }}_p{)}^{-1}{\mathbf{X}}^T=\mathbf{V}({\mathbf{D}}^2+{\mathbf{\Lambda }}_n{)}^{-1}\mathbf{D}{\mathbf{U}}^T.$$

$$\text{tr}({\mathbf{A}}_{p\times n}{\mathbf{B}}_{n\times p})=\sum _{i=1}^n\sum _{j=1}^p[\mathbf{A}\circ {\mathbf{B}}^T{]}_{ij}.$$

$${\lambda }_{cv}=\mathrm{argmin}{.}_{\lambda }\{\sum _{i=1}^n(y_i-{\mathbf{X}}_i{\hat{\beta }}_{\lambda }^{-f(i)}{)}^2\},$$

$$\text{GCV}(\lambda)=\sum_{i=1}^{n}\bigg{(}\frac{y_{i}-\boldsymbol{X}_{i}^{T}\hat{\boldsymbol{\beta}}_{\lambda}}{n-\text{tr}(\boldsymbol{H}_{\lambda}))}\bigg{)}^{2},$$

$$\mathbf{y}\sim \mathcal{N}(0,h^2{\sigma }^{\ast 2}\mathbf{R}+(1-h^2){\sigma }^{\ast 2}{\mathbf{I}}_n),$$

$$h^{2}=\operatorname*{arg\,max}._{h^{2}}\bigg{(}-\text{log}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}\frac{\tilde{y}_{i}^{2}}{h^{2}(\ell_{i}-1)+1}\bigg{)}-\frac{1}{n}\sum_{i=1}^{n}(\text{log}(h^{2}(\ell_{i}-1)+1)\bigg{)},$$

$${\hat{\sigma }}^2=\frac{(\mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }}{)}^T(\mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }})}{\nu },$$

$${\hat{\sigma }}_r^2=\frac{(\mathbf{y}-{\mathbf{Z}}_r{\hat{\mathbf{\alpha }}}_r{)}^T(\mathbf{y}-{\mathbf{Z}}_r{\hat{\mathbf{\alpha }}}_r)}{n-r}.$$

$$\mathrm{argmax}{.}_{{\sigma }^2,{\tau }^2}P(\mathbf{y})=\mathrm{argmax}{.}_{{\sigma }^2,{\tau }^2}{\int }_{\mathbf{\beta }}\mathcal{L}(\mathbf{y};\mathbf{\beta },{\sigma }^2)\pi (\mathbf{\beta };{\tau }^2)d\mathbf{\beta }.$$

$$P(\mathbf{y})$$

$$\begin{array}{rl}{\hat{\tau }}^2& =\frac{{\sum }_{i\ne k}^{n,n}{y}_i{y}_k}{{\sum }_{i\ne k}^{n,n}{\mathbf{\Sigma }}_{ik}^{\mathbf{X}}}\\ {\hat{\sigma }}^2& =n^{-1}\sum _{i=1}^n(y_i^2-{\hat{\tau }}^2{\mathbf{\Sigma }}_{ii}^{\mathbf{X}})\end{array}$$

$${\mathbf{y}}_{n\times 1}={\mathbf{X}}_{n\times p}{\mathbf{\beta }}_{p\times 1}+{\mathbf{ϵ}}_{n\times 1}$$

$$x_{ij}\frac{\text{buildrel}\mathrm{iid}}{\sim }\mathcal{N}(0,1)$$

$${\beta }_j\frac{\text{buildrel}\mathrm{iid}}{\sim }{p}_0{\delta }_0+(1-p_0)\mathcal{N}(0,{\tau }_0^2).$$

$${\beta }_j\frac{\text{buildrel}\mathrm{iid}}{\sim }\text{Laplace}(\mu =0,b=0.0707)\text{and}{\beta }_j\frac{\text{buildrel}\mathrm{iid}}{\sim }\text{Uniform}(a=-0.17,b=0.17),$$

$${ϵ}_i^{\prime }\frac{\text{buildrel}\mathrm{iid}}{\sim }{\text{t}}_4{ϵ}_i=(10/2{)}^{1/2}{ϵ}_i^{\prime },$$

$$\begin{array}{rl}{\mathbf{X}}_{n\times p}& \sim \mathcal{N}(0,\mathbf{\Xi }),\end{array}$$

$$\mathbf{y}={\mathbf{X}}_{\text{f}}\mathbf{\alpha }+{\mathbf{X}}_{\text{r}}\mathbf{\beta }+\mathbf{ϵ},$$

$$P({\mathbf{A}}^T\mathbf{y})=\mathcal{N}(\mathbf{y};\mathbf{\mu }=0,\mathbf{\Sigma }={\mathbf{A}}^T{\mathbf{\Sigma }}_r\mathbf{A})$$

$${\mathbf{y}}_{n\times 1}={\mathbf{X}}_{\text{f},n\times m}{\mathbf{\alpha }}_{m\times 1}+{\mathbf{X}}_{\text{r},n\times p}{\mathbf{\beta }}_{p\times 1}+{\mathbf{ϵ}}_{n\times 1}$$

$$x_{f,ik}\frac{\text{buildrel}\mathrm{iid}}{\sim }\mathcal{N}(0,1)$$

$${\alpha }_k\frac{\text{buildrel}\mathrm{iid}}{\sim }{p}_{0,f}{\delta }_0+(1-p_{0,f})\mathcal{N}(0,{\tau }_{0,f}^2)$$

$$\text{ML}(\lambda )={\int }_{\mathbf{\beta }\in {\mathbb{R}}^p}\mathcal{L}(\mathbf{Y};\mathbf{\beta },\mathbf{X}){\pi }_{\lambda }(\mathbf{\beta })d\mathbf{\beta }={\int }_{\mathbf{\beta }\in {\mathbb{R}}^p}\mathcal{L}(\mathbf{Y};\mathbf{\beta },\mathbf{X})\varphi ({\beta }_1;0,1/\lambda )\cdots \varphi ({\beta }_p;0,1/\lambda )d\mathbf{\beta },$$

$$\text{ML}(\lambda )=E_{{\pi }_{\lambda }(\mathbf{\beta })}[\mathcal{L}(\mathbf{Y};\mathbf{\beta },\mathbf{X})]=E_{{\pi }_{\lambda }(\mathbf{\beta })}[\mathcal{L}(\mathbf{Y};\mathbf{X}\mathbf{\beta })]=E_{{\pi }_{\lambda }^{\prime }(\mathbf{X}\mathbf{\beta })}[\mathcal{L}(\mathbf{Y};\mathbf{X}\mathbf{\beta })],$$

$$\begin{array}{rl}\text{ML}(\lambda )& ={\int }_{\mathbf{\beta }\in {\mathbb{R}}^p}{g}_{\mathbf{Y},\lambda }(\mathbf{\beta })d\mathbf{\beta }={\int }_{\mathbf{\beta }\in {\mathbb{R}}^p}\mathcal{L}(\mathbf{Y};\mathbf{\beta },\mathbf{X})\varphi ({\beta }_1;0,1/\lambda )\cdots \varphi ({\beta }_p;0,1/\lambda )d\mathbf{\beta }\\ & ={\int }_{{\mathbf{\beta }}^{\mathbf{X}}\in {\mathbb{R}}^n}{h}_{\mathbf{Y},\lambda }({\mathbf{\beta }}^{\mathbf{X}})d{\mathbf{\beta }}^{\mathbf{X}}={\int }_{{\mathbf{\beta }}^{\mathbf{X}}\in {\mathbb{R}}^n}\mathcal{L}(\mathbf{Y};{\mathbf{\beta }}^{\mathbf{X}},{\mathbf{I}}_n)\mathbf{\varphi }({\mathbf{\beta }}^{\mathbf{X}};0,{\Sigma }_{\lambda })d{\mathbf{\beta }}^{\mathbf{X}}.\end{array}$$

$$\log h_{\mathbf{Y},\lambda }({\mathbf{\beta }}^{\mathbf{X}})\propto \ell (\mathbf{Y};{\mathbf{\beta }}^{\mathbf{X}},{\mathbf{I}}_n)-({\mathbf{\beta }}^{\mathbf{X}}{)}^T{\Sigma }_{\lambda }^{-1}{\mathbf{\beta }}^{\mathbf{X}}.$$

$$\begin{array}{rl}{\mathbf{y}}_{n\times 1}& ={\mathbf{X}}_{n\times p}{\mathbf{\beta }}_{p\times 1}+{\mathbf{ϵ}}_{n\times 1},\\ \mathbf{\beta }& \sim \mathcal{N}(0,{\sigma }^2{\tau }^2{I}_p),\\ \mathbf{ϵ}& \sim \mathcal{N}(0,{\sigma }^2{I}_n),\\ {\sigma }^{-2}& \sim \mathcal{G}(a,b)\end{array}$$

$$\begin{array}{rl}\pi (\mathbf{y}|\mathbf{\beta },{\sigma }^{-2})& =(2\pi {)}^{-n/2}({\sigma }^{-2}{)}^{n/2}\text{exp}\left\{-\frac{1}{2}{\sigma }^{-2}(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^T(\mathbf{y}-\mathbf{X}\mathbf{\beta })\right\}\\ \pi (\mathbf{\beta }|{\sigma }^{-2},{\tau }^{-2})& =(2\pi {)}^{-p/2}({\sigma }^{-2}{\tau }^{-2}{)}^{p/2}\text{exp}\left\{-\frac{1}{2}{\sigma }^{-2}{\tau }^{-2}{\mathbf{\beta }}^T\mathbf{\beta }\right\}\\ \pi ({\sigma }^{-2})& =\frac{b^a}{\Gamma \left(a\right)}({\sigma }^{-2}{)}^{a-1}\text{exp}\left\{-b{\sigma }^{-2}\right\}.\end{array}$$

$$\pi (\mathbf{\beta }|\mathbf{y})=\int \pi (\mathbf{y}|\mathbf{\beta },{\sigma }^{-2})\pi (\mathbf{\beta }|{\sigma }^{-2},{\tau }^{-2})\pi ({\sigma }^{-2})d{\sigma }^{-2}{=}^d{\mathcal{T}}_{2a+n}({\mathbf{\beta }}_{\nu }^{\ast },{\Sigma }_{\nu }^{\ast }).$$

Full text

An example on weight gain data with genomic covariates confirms the good performance of MML compared to CV. Several extensions are presented. First, to the high-dimensional linear mixed effects model, with REML as an alternative to MML. Second, to the conjugate Bayesian setting, which proves to be a good alternative. Third, and most prominently, to generalized linear models for which we derive a computationally efficient MML estimator by re-writing the marginal likelihood as an $n$-dimensional integral. For Poisson and Binomial ridge regression, we demonstrate the superior accuracy of the resulting MML estimator of $\lambda$ as compared to CV. Software is provided to enable reproduction of all results presented here.

††journal: Journal of LaTeX Templates

1 Introduction

Estimation of hyper-parameters is an essential part of fitting high-dimensional Gaussian random effect regression models, also known as ridge regression. These models are widely applied in genomics and genetics applications, where often the number of variables $p$ is much larger than the number of samples $n$, i.e. $p\ggg n$.

We initially focus on the linear model. The goal is to estimate error variance $\sigma^{2}$ and random effects variance $\tau^{2}$ or functions thereof, in particular the ridge penalty parameter, $\lambda=\frac{\sigma^{2}}{\tau^{2}}$, or heritability index, $h^{2}=\frac{p\tau^{2}}{p\tau^{2}+\sigma^{2}}$. Here, the ridge penalty is used in classical ridge regression to shrink the regression coefficients to zero [1], whereas heritability measures the fraction of variation between individuals within a population that is due their genotypes [2]. The estimators of $\sigma^{2}$ and $\tau^{2}$ can be used to estimate $\lambda$ or $h^{2}$, but also for statistical testing [3]. We review several estimators, based on maximum marginal likelihood (MML), moment equations, (generalized) cross-validation, dimension reduction, or degrees-of-freedom adjustment. Some of these estimators are classical, while others have recently been introduced.

We systematically review and compare the estimators in a broad variety of high-dimensional settings. For estimation of $\lambda$ in low-dimensional settings, we refer to [4, 5, 6]. We address the effect of multi-collinearity and robustness against model misspecifications, such as sparsity and non-Gaussian errors. The comparisons are extended to the linear mixed effects model, with $q\ll n$ fixed effects added to the model and to Bayesian linear regression. The linear model part is concluded by a genomics data application to weight gain prediction after kidney transplantation.

The observed good performance of MML in the linear model setting was a stimulus to consider MML for high-dimensional generalized linear models (GLM). MML is more involved here than in the linear model, because of the non-conjugacy of the likelihood and prior. Therefore, approximations are required, such as Laplace ones. While these have been addressed by others [7, 8], we derive an estimator which is computationally efficient for $p\ggg n$ settings. For Poisson and Binomial ridge regression, we demonstrate the superior accuracy of MML estimation of $\lambda$ as compared to cross-validation.

Our software enables reproduction of all results. In addition, it allows comparisons for one’s own high-dimensional data matrix by simulating the response conditional on this matrix, as we do for two cancer genomics examples. Computational shortcuts and considerations are discussed throughout the paper, and detailed at the end, including computing times.

1.1 The Model

We initially focus on high-dimensional linear regression with random effects. Variables are denoted by $j=1,\ldots,p$ and samples by $i=1,\ldots,n$. Then:

[TABLE]

Here, $\boldsymbol{y}=(y_{1},\ldots,y_{n})$ is the vector of responses, $\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{p})^{T}$ corresponds to the random effects and $\boldsymbol{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$ is a vector of Gaussian errors. Furthermore, $\boldsymbol{X}$ is a fixed $n\times p$ matrix: $(\boldsymbol{X}_{1}\cdots\boldsymbol{X}_{n})^{T},$ with $\boldsymbol{X}_{i}=(x_{i1},\ldots,x_{ip})^{T}$.

1.2 Estimation Methods

We distinguish three categories of estimation methods:

Estimation of functions of $(\sigma^{2},\tau^{2})$, in particular $\lambda=\frac{\sigma^{2}}{\tau^{2}}$ [9], used in ridge regression to minimize $||\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}||_{2}^{2}+\lambda||\bm{\beta}||_{2}^{2}$, and heritability $h^{2}=\frac{p\tau^{2}}{p\tau^{2}+\sigma^{2}}$ [10]. 2. 2.

Separate estimation of $\sigma^{2}$ [11, 12], possibly followed by plug-in estimation of $\tau^{2}$. 3. 3.

Joint estimation: estimate $\sigma^{2}$ and $\tau^{2}$ jointly

Below, we discuss several methods for each of these categories. They have several matrices and matrix computations in common, which we therefore introduce first.

1.3 Notation and matrix computations

Throughout the paper, we will use the following notation:

[TABLE]

Many of the estimators below require calculations on potentially very large matrices. The following two well-known equalities can highly alleviate the computational burden.

First, $\boldsymbol{C}=\boldsymbol{C}_{\lambda}$, and hence also $\hat{\boldsymbol{\beta}}$ and $\boldsymbol{H}$, can be efficiently computed by using singular value decomposition (SVD). Decompose $\boldsymbol{X}=\boldsymbol{U}_{n\times n}\boldsymbol{D}_{n\times n}(\boldsymbol{V}_{p\times n})^{T}$ by SVD, and denote $\boldsymbol{\Lambda}_{q}=\lambda\boldsymbol{I}_{q}$. Then,

[TABLE]

The latter requires inversion of an $n\times n$ matrix only. Second, the following efficient trace computation for matrix products applies to $\text{tr}(\boldsymbol{H})=\text{tr}(\boldsymbol{X}\boldsymbol{C}_{\lambda}):$

[TABLE]

2 Methods

2.1 Estimating functions of $\sigma^{2}$ and $\tau^{2}$

2.1.1 Estimating $\lambda$ by $K$-fold CV

A benchmark method that is used extensively to estimate $\lambda=\sigma^{2}/\tau^{2}$ is cross-validation. Here, we use $K$-fold CV, as implemented in the popular R-package glmnet [13]. Let $f(i)$ denote the set of samples left out for testing at the same fold as sample $i$. Then, CV-based estimation of $\lambda$ pertains to minimizing the cross-validated prediction error:

[TABLE]

where $\hat{\bm{\beta}}_{\lambda}^{-f(i)}$ denotes the estimate of $\bm{\beta}$ based on training samples $\{1,\ldots,n\}\setminus f(i)$ and penalty $\lambda$. Note that for leave-one-out-cross-validation ($n$-fold CV) the analytical solution of (5) is the PRESS statistic [14].

2.1.2 Estimating $\lambda$ by Generalized Cross Validation

Generalized Cross Validation (GCV) is a rotation-invariant form of the PRESS statistic. It is more robust than the latter to (near-diagonal) hat matrices $\boldsymbol{H}_{\lambda}$ [9]. For the linear model, the criterion is [15]:

[TABLE]

where the trace of $\boldsymbol{H}_{\lambda}$ can be computed efficiently by (4). Then, $\lambda_{\text{gcv}}=\operatorname*{arg\,min}._{\lambda}\text{GCV}(\lambda)$.

2.1.3 Estimating heritability by HiLMM

Heritability is defined by $h^{2}=\frac{p\tau^{2}}{p\tau^{2}+\sigma^{2}}$. A recent method which estimates heritability directly using maximum likelihood is proposed in [10]. Analogously to equation (12), it is based on writing:

[TABLE]

where $\sigma^{*2}=p\tau^{2}+\sigma^{2}$ and $\boldsymbol{R}=\boldsymbol{X}\boldsymbol{X}^{T}/p$. Now, apply an eigen-decomposition to $\boldsymbol{R}$: $\boldsymbol{R}=\boldsymbol{Q}\boldsymbol{L}\boldsymbol{Q}^{T}$. Then, heritability is estimated by [10]:

[TABLE]

with $\ell_{i}$ and $\tilde{y}_{i}$ the $i$th element of $\boldsymbol{L}$ and $\tilde{\boldsymbol{y}}=\boldsymbol{Q}^{T}\boldsymbol{y}$, respectively. The authors provide rigorous consistency results for their estimator, as well as theoretical confidence bounds, also for mixed models and sparse settings.

2.2 Estimation of $\sigma^{2}$

The two methods below rely on an estimate $\hat{\boldsymbol{\beta}}=\hat{\boldsymbol{\beta}}_{\lambda}$, where $\lambda=\sigma^{2}/\tau^{2}$ is estimated by (G)CV. Then $\sigma^{2}$ is estimated conditional on $\hat{\boldsymbol{\beta}}.$ If desired, $\tau^{2}$ may then be estimated by $\hat{\tau}^{2}=\hat{\sigma}^{2}/\hat{\lambda}.$

2.2.1 Basic estimate

A basic estimate of $\sigma^{2}$, and often used in practice, is given by [16]:

[TABLE]

which is the residual mean square error. Here, the residual effective degrees of freedom [16] equals $\nu=n-\text{tr}(2\boldsymbol{H}-\boldsymbol{HH}^{T})$, with $\boldsymbol{H}$ as in (2). We also considered (9) with $\nu=n-\text{tr}(\boldsymbol{H})$, as in [17], which rendered similar, slightly inferior results.

2.2.2 PCR-based estimate

The estimator for $\sigma^{2}$ may also be based on Principal Component Regression (PCR). PCR is based on the eigen-decomposition $\boldsymbol{X}^{T}\boldsymbol{X}=\tilde{\boldsymbol{Q}}\boldsymbol{D}^{2}\tilde{\boldsymbol{Q}}^{T}$. Denoting $\boldsymbol{Z}=\boldsymbol{X}\tilde{\boldsymbol{Q}}$ and $\boldsymbol{\alpha}=\tilde{\boldsymbol{Q}}^{T}\boldsymbol{\beta}$, we have $\boldsymbol{y}=\boldsymbol{Z}\boldsymbol{\alpha}+\boldsymbol{\epsilon}$. Then, $\boldsymbol{Z}$ is reduced from $p$ columns to $r\leq\text{min}(n,p)$ principal components, a crucial step [12]. Using the reduced model, $\sigma^{2}$ is estimated by the residual mean square error [12]:

[TABLE]

2.3 Joint estimation of $\sigma^{2}$ and $\tau^{2}$

2.3.1 MML

An Empirical Bayes estimate of $\sigma^{2}$ and $\tau^{2}$ is obtained by maximizing the marginal likelihood (MML), also referred to as model evidence in machine learning [18]. This corresponds to:

[TABLE]

Since $\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}$, $P(\boldsymbol{y})$ is simply derived from the convolution of Gaussian random variables, implying $E[\boldsymbol{y}]=E[\boldsymbol{X}\boldsymbol{\beta}]+E[\boldsymbol{\epsilon}]=\mathbf{0}$, and $V[\boldsymbol{y}]=V[\boldsymbol{X}\boldsymbol{\beta}]+V[\boldsymbol{\epsilon}]=\boldsymbol{X}\boldsymbol{X}^{T}\tau^{2}+\sigma^{2}\boldsymbol{I}_{n}$, so

[TABLE]

This is easily maximized over $\sigma^{2}$ and $\tau^{2}$. Note that after computing $\boldsymbol{X}\boldsymbol{X}^{T}$ (12) requires operations on $n\times n$ matrices only.

2.3.2 Method of Moments (MoM)

An alternative to MML is to match the empirical second moments of $\boldsymbol{y}$ to their theoretical counterparts. From (12) we observe that the covariances depend on $\tau^{2}$ only. Hence, we obtain an estimator of $\tau^{2}$ by equating the sum of $y_{i}y_{k}$ to that of the theoretical covariances, $\boldsymbol{\Sigma}_{ik}=\mathbb{E}[\boldsymbol{y}_{i}\boldsymbol{y}_{k}]$, with $\boldsymbol{\Sigma}$ as in (12). Then, with $\boldsymbol{\Sigma}^{\boldsymbol{X}}=\boldsymbol{X}\boldsymbol{X}^{T}$, an estimator for $\sigma^{2}$ is obtained by substituting $\hat{\tau}^{2}$ and equating the sum of $y_{i}^{2}$ to the sum of theoretical variances, $\boldsymbol{\Sigma}_{ii}=\mathbb{E}[\boldsymbol{y}_{i}^{2}]$:

[TABLE]

These equations also hold for non-Gaussian error terms, which could be an advantage over MML. Moreover, no optimization over $\sigma^{2}$ and $\tau^{2}$ is required, so MoM is computationally very attractive.

3 Comparisons

For the linear random effects model (ridge regression) we study the following settings:

$\boldsymbol{\beta}$ and $\boldsymbol{\epsilon}$ generated from model (1), independent $\boldsymbol{X}$

• 2.

$\boldsymbol{\beta}$ or $\boldsymbol{\epsilon}$ generated from non-Gaussian distributions, independent $\boldsymbol{X}$

• 3.

$\boldsymbol{\beta}$ and $\boldsymbol{\epsilon}$ from model (1), multicollinear $\boldsymbol{X}$

• 4.

$\boldsymbol{\beta}$ and $\boldsymbol{\epsilon}$ from model (1), data-based $\boldsymbol{X}$.

As is common for real data, the variables, i.e. the rows of $\boldsymbol{X}$, were always standardized for the $L_{2}$-penalty to have the same effect on all variables. All the results are based on 100 simulated data sets. Cross-validation is applied on 10 folds. Results from $n$-fold CV (leave-one-out) were generally fairly similar. We focus on the high-dimensional setting with $n=100,p=1000$, with excursions to larger data sets and dimensions of real data. In all visualizations below the red dotted lines indicate true values. Moreover, values larger than 20 times the true value were truncated and slightly jittered. Discussion of all results is postponed to Section 3.4.

3.1 Independent $\boldsymbol{X}$

In correspondence to model (1) we sample :

[TABLE]

Figures 1(c)(a) and (b) display the results for $n=100,p=1000,\tau^{2}=0.01,\sigma^{2}=10$ and for a large data setting $n=1000,p=15000,\tau^{2}=0.01,\sigma^{2}=150$ (which both imply $h^{2}=0.5$).

3.2 Departures from a normal effect size distribution

We study the robustness of the methods against (sparse) non-Gaussian effect size distribution or error distribution. In sparse settings, many variables do not have an effect. To mimic this, we simulated the $\beta$’s from a mixture distribution with a ‘spike’ and a Gaussian ‘slab’:

[TABLE]

Here, we set $p_{0}=0.9,\tau^{2}_{0}=0.1,$ which implies $\tau^{2}=\mathbb{V}(\beta_{j})=\mathbb{E}(\beta_{j}^{2})-\mathbb{E}(\beta_{j})^{2}=(1-p_{0})\tau_{0}^{2}=0.01$, as in the Gaussian $\beta_{j}$ setting. Moreover, we also considered:

[TABLE]

where again the parameters are chosen such that $E(\beta_{j})=0$ and $\tau^{2}=\mathbb{V}(\beta_{j})=0.01.$ Apart from $\bm{\beta}$ all other quantities are simulated as in (LABEL:basicmodel). Results are displayed for $\sigma^{2}=10,\tau^{2}=0.01,n=100,p=1000$ in Figure 1(c)(c) for the Laplace (= lasso) effect size distribution and in Supplementary Figure 3 for the spike-and-slab and uniform effect size distribution.

Moreover, we considered heavy-tailed errors by sampling

[TABLE]

where the scalar is chosen such that $\sigma^{2}=\mathbb{V}(\epsilon_{i})=10$, as in the Gaussian error setting. Apart from $\boldsymbol{\epsilon}$, all other quantities are simulated as in (LABEL:basicmodel). Results are displayed in Supplementary Figure 3(c).

3.3 Multicollinear $\boldsymbol{X}$

3.3.1 Simulated $\boldsymbol{X}$

Next, the design matrix $\boldsymbol{X}$ is sampled using block-wise correlation. We replace the sampling of $\boldsymbol{X}$ in simulation model (LABEL:basicmodel) by:

[TABLE]

where $\boldsymbol{\Xi}$ is a unit variance covariance matrix with blocks of size $p^{*}\ll p$ with correlations $\rho$ on the off-diagonal. Figure 2(c)(a) shows the results for $\rho=0.5,p^{*}=10,n=100,p=1000$.

3.3.2 Real data $\boldsymbol{X}$

Finally, we consider the estimation of $\tau^{2}$ and $\sigma^{2}$ in a high- and medium-dimensional setting where $\boldsymbol{X}$ are real data, with likely collinear columns. The first data set (TCGA KIRC) concerns gene expression data of $p=18,391$ genes for $n=71$ kidney tumors. The second data set (TCPA OV) holds expression data of $p=224$ proteins for $n=408$ ovarian tumor samples. Details on both data sets are supplied in the Supplementary Information. To generate response $\boldsymbol{y}$ we use model (LABEL:basicmodel) with $\boldsymbol{X}$ given by the data. Here, $\tau^{2}=0.01$ and $\sigma^{2}$ is set such that $h^{2}=0.5.$ Figures 2(c)(b) and (c) show the results.

3.4 Discussion of results

3.4.1 MML vs MoM, Basic and PCR

Figures 1(c) and 2(c) and Supplementary Figure 3 clearly show superior performance of MML compared to MoM: both the bias and variability are much smaller for MML. Generally, MML also outperforms the Basic and PCR estimators of $\sigma^{2}$. The PCR estimator approaches the performance of MML for the KIRC and TCPA data (Figures 2(c)(b) and 2(c)(c)), and the Basic estimator performs reasonably well for the latter ($p<n$) data set. For other settings, the Basic estimator performs equally inferior as MoM. The results highlight the importance of joint estimation of $\sigma^{2}$ and $\tau^{2}$ in high-dimensional settings, because of their delicate interplay.

3.4.2 MML vs GCV and CV

For the estimation of $\lambda$ MML seems slightly superior to GCV and CV. GCV shows more estimates that deviate towards too small values of $\lambda$ (e.g. Figures 1(c)(b) and 2(c)(b), i.e. the large $p$ settings), whereas CV tends to render somewhat more skewed results, either to the right (Figures 1(c)(a) and 1(c)(c), 2(c)(a)), or to the left (Figure 2(c)(b)). For the spike-and-slab and uniform effects sizes and the $t_{4}$ errors the right-skewness of the CV-results is more pronounced (Supplementary Figure 3), indicating that minimization of the cross-validated prediction error (5) is more vulnerable to non-Gaussian $\boldsymbol{y}$ than MML and GCV. Note that the Laplace setting (Figure 1(c)(c)) relates directly to the lasso prior with scale parameter $1/\lambda_{1}$ [19]. The results indicate that MML with Gaussian prior could be useful to find the lasso penalty, or serve as a fast initial estimate by simply setting the lasso penalty $\lambda_{1}=\sqrt{(}2)/\hat{\tau}$, which follows from the variance of the lasso prior.

3.4.3 MML vs HiLMM

For the estimation of heritability $h^{2}$ Figures 1(c) and 2(c) and Supplementary Figure 3 show very comparable performance of MML and HiLMM. This similar performance is not surprising given that both methods are likelihood-based. Hence, while reparametrizing the likelihood (7) is certainly useful to study it as function of $h^{2}$ [10], the reparametrization seems not beneficial for the purpose of estimating $h^{2}$. In addition, unlike HiLMM, MML also returns estimates of $\tau^{2}$ and $\sigma^{2}$. Finally, comparing Figures 1(c)(a) and 1(c)(b) we observe that both MML and HiLMM clearly benefit from the larger $n$ and $p$.

4 Data example

We re-analyse the weight gain data, recently discussed in [17]. Details on the data are presented there, we provide a summary. The data consists of expression profiles of $n=26$ individuals with kidney transplants, where profiles consists of 28,869 genes as measured by Affymetrix Human Gene 1.0 ST arrays. The data is available in the EMBL-EBI ArrayExpress database (www.ebi.ac.uk/arrayexpress) under accession number E-GEOD-33070. It is known that kidney transplantation may lead to weight gain, and the study [20] investigates whether gene expression can be used to predict this. Such a prediction can be used to decide upon additional measures to prevent excessive weight gains. We reproduced the analysis by [17] as much as possible, including their prior selection of 1000 genes. Details on minor discrepancies, and an alternative analysis that accounts for the gene selection are discussed in the Supplementary Material. These did not affect the comparison qualitatively.

In [17], the authors illustrate their focused ridge (fridge) method and compare it with conventional ridge. In short, fridge estimates sample-specific ridge penalties, based on minimizing a per sample mean squared error (MSE) criterion on the level of the linear predictor $\boldsymbol{X}_{i}\boldsymbol{\beta}$. Since $\boldsymbol{\beta}$ is not known, it is replaced by an initial ridge estimate, $\hat{\boldsymbol{\beta}}_{\lambda}.$ Their sample specific penalty then depends on $\boldsymbol{X}_{i}$, and also on both $\hat{\lambda}$ and $\hat{\sigma}^{2}$. The authors use GCV (6) to obtain $\lambda$, and a slight variation of (9) to estimate $\sigma^{2}$. They show that fridge improves upon GCV-based ridge estimation. We wish to investigate whether i) MML estimation of $\lambda=\sigma^{2}/\tau^{2}$ also improves the performance of GCV-based ridge regression; and ii) whether MML estimation further boosts the performance of the fridge estimator. Here, predictive performance is measured by the mean squared prediction error (MSPE) using leave-one-out cross-validation (loocv).

The estimates of MML differ markedly from those of GCV: $(\hat{\lambda}_{\text{MML}},\hat{\sigma}^{2}_{\text{MML}})=(0.77,0.59)$, while $(\hat{\lambda}_{\text{GCV}},\hat{\sigma}^{2}_{\text{GCV}})=(20.92,8.08)$. Using $\hat{\lambda}_{\text{MML}}$ instead of $\hat{\lambda}_{\text{GCV}}$ for the estimation of $\boldsymbol{\beta}$ substantially reduced the mean squared prediction error: $\text{MSPE}_{\text{MML}}=14.40,$ while $\text{MSPE}_{\text{GCV}}=16.38$, a relative decrease of 12.1%. Using $\hat{\lambda}_{\text{GCV}}$, as in [17], fridge also reduced the MSPE, but to a lesser extent: $\text{MSPE}_{\text{fridge}}=15.80,$ a relative decrease of 3.5% with respect to $\text{MSPE}_{\text{GCV}}.$ Application of fridge using $\hat{\lambda}_{\text{MML}}$ did not further decrease $\text{MSPE}_{\text{MML}}$, nor did it increase it. Possibly, the already fairly small value of $\hat{\lambda}_{\text{MML}}$ left little room for improvement. Figure 3 displays absolute prediction errors per sample and illustrates the improved prediction by ridge using $\lambda_{\text{MML}}$ (and to a lesser extent by fridge) with respect to ridge using $\lambda_{\text{GCV}}$.

5 Extensions

5.1 Extension 1: Mixed effects model

A natural extension of the high-dimensional random effects model (1) is the mixed effects model:

[TABLE]

where we assume that the $n\times m$ design matrix for the fixed effects, $\boldsymbol{X}_{\text{f}}$, is of low-rank, so $m\ll n$, as opposed to the random effects design matrix $\boldsymbol{X}_{\text{r}}$. Restricted maximum likelihood (REML) deals with the fixed effects by contrasting them out. For the error contrast vector $\boldsymbol{y}-\boldsymbol{X}_{\text{f}}\hat{\alpha}^{\text{OLS}}=\boldsymbol{A}^{T}\boldsymbol{y},$ with $\boldsymbol{A}=\boldsymbol{I}_{n}-\boldsymbol{X}_{\text{f}}(\boldsymbol{X}_{\text{f}}^{T}\boldsymbol{X}_{\text{f}})^{-1}\boldsymbol{X}_{\text{f}}^{T}$, the marginal likelihood for the variance components equals (see e.g. [21]):

[TABLE]

with $\boldsymbol{\Sigma}_{r}=\boldsymbol{X}_{r}\boldsymbol{X}_{r}^{T}\tau^{2}+\sigma^{2}\boldsymbol{I}_{n}.$ In addition to maximizing (18) as a function of $(\sigma^{2},\tau^{2})$, we attempted solving the set of two estimation equations suggested by [22], but this rendered instable results inferior to maximizing (18) directly.

Alternatively, MML may be used, but it has to be adjusted to also estimate the fixed effects in the model. This implies replacing $\boldsymbol{0}$ in Gaussian likelihood (11) by $\boldsymbol{X}_{\text{f}}\boldsymbol{\alpha}$, and optimizing (11) with respect to $2+m$ parameters, where $m$ is the number of fixed parameters. The mixed model simulation setting is as follows:

[TABLE]

where $n=100,p=1000,m=10,p_{0}=0.9,\tau_{0}^{2}=0.1$ (implying variance $\tau^{2}=(1-p_{0})\tau^{2}_{0}=0.01$ for generating random effects) and $p_{0,f}=0.5,\tau^{2}_{0,f}=0.20$ (implying variance $\tau^{2}_{f}=0.1$ for generating fixed effects). Note that we focused on a fairly sparse setting for the random effects and larger prior variance of fixed effects than of random effects, which enables a stronger impact of the small number of fixed effects. Figure 4 shows the results of REML, MML and CV (by glmnet, using penalty factor 0 for the fixed effects) for the estimation of $\tau^{2},\sigma^{2},\lambda$ and $h^{2}$.

From Figure 4 we observe that REML indeed improves MML in terms of bias, however at the cost of increased variability. For the estimation of $\lambda$, CV is fairly competitive to REML and MML, although it renders markedly more over-penalization.

5.2 Extension 2: Bayesian linear regression

So far, we focused on classical methods. Bayesian methods may be a good alternative. We applied the standard Bayesian linear regression model, i.e. the conjugate model with i.i.d. priors $\pi(\beta_{j})=N(0,\sigma^{2}\tau^{2})$, with $\tau^{2}$ fixed and $\sigma^{2}$ endowed with a vague inverse-gamma prior (see Supplementary Material for details). For this model the maximum marginal likelihood estimator for $\tau^{2}$ is still analytical [23], and so is the posterior mode estimate of $\sigma^{2}$. Figure 5 shows the results in comparison to MML, i.e. maximization of (12), for the random effects case with multi-collinear $\boldsymbol{X}$, as in Section 3.3.1. Results for other settings were in essence very similar.

From the results we conclude that the conjugate Bayes estimates are very close to those of MML. This is in line with the fact that this conjugate model with prior variance $\tau^{2}=\sigma^{2}/\lambda$ is known to render posterior mean estimates of $\boldsymbol{\beta}$ that equal the $\lambda$-penalized ridge regression estimates.

The conjugate Bayesian model is scale-invariant, because the $\beta$ prior contains the error variance $\sigma^{2}$. Recently, it was criticized for its non-robustness against misspecification of the fixed $\tau^{2}$ when estimating $\sigma^{2}$ [24]. However, in practice one needs to estimate $\tau^{2}$ by either empirical Bayes (e.g. maximum marginal likelihood) or full Bayes. We repeated the simulation by [24] (see Supplementary Material). The results show that the estimates of $\sigma^{2}$ are much better when estimating $\tau^{2}$ by empirical Bayes instead of fixing it, and in fact very competitive to alternatives proposed by [24].

5.3 Extension 3: Generalized linear models

5.3.1 Setting

Motivated by the good results for MML in the linear setting, we wish to extend MML estimation to the high-dimensional generalized linear model (GLM) setting, where the likelihood depends on the regression parameter $\bm{\beta}$ only via the linear predictor, $\boldsymbol{X}\bm{\beta}$. Hence, likelihood $\mathcal{L}(\boldsymbol{Y};\bm{\beta},\boldsymbol{X})$ is defined by a density $f_{\mu}(\boldsymbol{Y})$ (e.g. Poisson), where $\boldsymbol{X}\bm{\beta}$ is mapped to $\mu$ by a link function (e.g. $\log$). As before, we a priori assume i.i.d. $\beta_{j}\sim N(0,\tau^{2})$, here equivalent to an $L_{2}$ penalty $\lambda=1/\tau^{2}$ when estimating $\bm{\beta}$ by penalized likelihood. In [7] an iterative algorithm to estimate $\lambda$ is derived which alternates estimation of $\bm{\beta}$ by maximization w.r.t. $\lambda$, requiring the computation of the trace of a Hessian of a $p\times p$ matrix. Here, the estimation of $\bm{\beta}$ itself is much slower than in the linear case, because it is not analytic and requires iterative weighted least squares approximation. Below we show how to substantially alleviate the computational burden in the $p\ggg n$ setting by re-parameterizing the marginal likelihood implying computations in $\mathbb{R}^{n}$ instead of $\mathbb{R}^{p}$.

5.3.2 Method

We have for the marginal likelihood:

[TABLE]

where $\phi(\beta,\mu,\tau^{2})$ denotes the normal density with mean $\mu$ and variance $\tau^{2}$. Now a crucial observation is that for GLM:

[TABLE]

because the likelihood depends on $\bm{\beta}$ only via the linear predictor $\boldsymbol{X}\bm{\beta}$. Here, $\pi^{\prime}_{\lambda}(\boldsymbol{X}\bm{\beta})$ is the implied $n$-dimensional prior distribution of $\boldsymbol{X}\bm{\beta}$. This is a multivariate normal: $\boldsymbol{\phi}(\bm{\beta}^{\boldsymbol{X}};\boldsymbol{\mu}=\boldsymbol{0},\Sigma_{\lambda}=\boldsymbol{X}\boldsymbol{X}^{T}/\lambda)$. Therefore, we have:

[TABLE]

Hence, the $p$-dimensional integral may be replaced by an $n$-dimensional one, with obvious computational advantages when $p\ggg n$. Moreover, the use of (22) allows applying implemented Laplace approximations, which tend to be more accurate in lower dimensions. The Laplace approximation requires $\hat{\bm{\beta}}^{\boldsymbol{X}}=\operatorname*{arg\,max}._{\bm{\beta}^{\boldsymbol{X}}}\{h_{\boldsymbol{Y},\lambda}(\bm{\beta}^{\boldsymbol{X}})\}$. We emphasize that this does generally not equal $\boldsymbol{X}\hat{\bm{\beta}},$ where $\hat{\bm{\beta}}=\operatorname*{arg\,max}._{\bm{\beta}}\{g_{\boldsymbol{Y},\lambda}(\bm{\beta})\}$: the maximum of the commonly used $L_{2}$ penalized (log)-likelihood. However, $\hat{\bm{\beta}}^{\boldsymbol{X}}$ can be computed by noting that

[TABLE]

In other words, this is the penalized log-likelihood when regressing $\boldsymbol{Y}$ on the identity design matrix $\mathbf{I}_{n}$ using an $L_{2}$ smoothing penalty matrix $(\bm{\beta}^{\boldsymbol{X}})^{T}\Sigma_{\lambda}^{-1}\bm{\beta}^{\boldsymbol{X}}=\lambda(\bm{\beta}^{\boldsymbol{X}})^{T}(\boldsymbol{X}\boldsymbol{X}^{T})^{-1}\bm{\beta}^{\boldsymbol{X}}$. The latter fits conveniently into the set-up of [8], as implemented in the R-package mgcv. This also facilitates MML estimation of $\lambda$ by maximizing $\text{ML}(\lambda)$, with $h_{\boldsymbol{Y},\lambda}(\bm{\beta}^{\boldsymbol{X}})$ as in (23). If the columns of $\boldsymbol{X}$ are standardized (common in high-dimensional studies), $\boldsymbol{X}\boldsymbol{X}^{T}$ has rank $n-1$ instead of $n$, implying that $(\boldsymbol{X}\boldsymbol{X}^{T})^{-1}$ does not exist and should be replaced by a pseudo-inverse $(\boldsymbol{X}\boldsymbol{X}^{T})^{+}$, such as the Moore-Penrose inverse.

In a full Bayesian linear model setting, dimension reduction is also discussed by [25], where $\boldsymbol{X}\bm{\beta}$ is substituted by a $n$-dimensional factor analytic representation, which requires an SVD of $\boldsymbol{X}$. In addition, there it is not used for hyper-parameter estimation by marginal likelihood, but instead for specifying (hierarchical) priors for the factors.

5.3.3 Results

R packages like glmnet [13] and penalized [26] estimate $\lambda$ by cross-validation, and also mgcv allows, next to the MML estimation, (generalized) CV estimation [8]. Figures 6(a)and 6)(b) show the results for Poisson ridge regression, with $\boldsymbol{Y}_{i}\sim\text{Pois}(\lambda_{i}),\lambda_{i}=\exp(\boldsymbol{X}_{i}\bm{\beta})$, $\bm{\beta}$ generated as in (LABEL:basicmodel), and $\boldsymbol{X}$ generated as in (LABEL:basicmodel) and (16), which denote the independent $\boldsymbol{X}$ and multi-collinear $\boldsymbol{X}$ setting, respectively.

Figure 6 clearly shows the superior performance of MML based on (22) over CV. In particular, glmnet and penalized render strongly upward biased values. The mgcv GCV values are still inferior to MML based ones, but much better than the latter two, which may be due to the different regression estimators used (Laplace approximation versus iterative weighted least squares). We should stress that CV does not target for the estimation of $\lambda$ as such, but merely for minimizing prediction error. Nevertheless, the difference is remarkably larger than in the corresponding linear case (see Figures 1(c) and 2(c)).

The Supplementary Material shows the results for Binomial ridge regression. While the differences in performance are less dramatic than for the Poisson setting, MML still renders much better estimates of $\lambda$ than CV-based approaches.

6 Computational aspects and software

All methods and simulations presented here are implemented in a few wrapper R scripts: one for the linear random effects model (which includes the conjugate Bayes estimator), one for the linear mixed effects model, and one for Poisson and Binomial ridge regression. Parallel computations are supported. The scripts allow exact reproduction of the results in this manuscript as well as comparisons for other simulation or user-specific real data $\boldsymbol{X}$ cases. In addition, a script is supplied to produce the box-plots as in this manuscript.

HiLMM, PCR and CV implementations are provided by the R-packages HiLMM, v1.1 [10], ridge, v1.8-16 [12] (code slightly adapted for computational efficiency) and glmnet, v2.0-16 [13]. The methods MML, REML, Bayes, MoM, Basic and GCV were implemented by us for the linear random and mixed effects models. For Poisson and Binomial ridge regression we applied mgcv, v1.8-16 [8] after our re-parametrization (22) to obtain MML and GCV results, while for CV glmnet and penalized, v0.9-50 [26] were applied. For all methods that required optimization the R routine optim was used, with default settings. CV was based on 10 folds.

Computing times of the various methods largely depend on $n$ and $p$, much less so on the exact simulation setting. These are displayed for $n=100,500$ and $p=10^{3},10^{4},10^{5}$ in Table 1, based on computations with one CPU of an Intel® Xeon® CPU E5-2660 v3 @ 2.60GHz server. For Poisson ridge regression, we only report the computing times of MML and GCV, because, as reported in Figure 6, the performance of CV-based methods was very inferior.

From Table 1 we conclude that MML is also computationally very attractive. Its efficiency is explained by the fact that, unlike many of other methods, it does not require an SVD or other matrix decomposition of $\boldsymbol{X}$. Moreover, the only computation that involves dimension $p$ is the product $\boldsymbol{X}\boldsymbol{X}^{T}$.

7 Discussion

We compared several estimators in a large variety of high-dimensional settings. The results showed that plain maximum marginal likelihood works well in many settings. MML is generally superior to methods that aim to separately estimate $\sigma^{2}$ (9, 10). Apparently, the estimates of $\sigma^{2}$ and $\tau^{2}$ are so intrinsically linked in the high-dimensional setting that separate estimation is sub-optimal. The moment estimator (MoM) is generally not competitive to MML. It may, however, be useful in large systems with multiple hyper-parameters to estimate relative penalties, which are less sensitive to scaling issues than the global penalty parameter [27]. MoM may also be a useful initial estimator for more complex estimators that are based on optimization, such as MML.

Possibly somewhat surprising is the good performance of MML for estimating $\lambda$ and $h^{2}$, as these are functions of $\sigma^{2}$ and $\tau^{2}$. For the estimation of $\lambda$ it is generally better than or competitive to (generalized) CV, an observation also made for the low-dimensional setting [8]. The inferior performance of the basic estimator of $\sigma^{2}$ (9) implies that alternative estimators of $\lambda$ that use $\hat{\sigma}^{2}$ as a plug-in are unlikely to perform well in high-dimensional settings. Such estimators, including the original one by Hoerl and Kennard [1], are compared by [4, 6], who show that some do perform well in the low-dimensional setting. For Poisson ridge regression, similar estimators of $\lambda$ are available [5], but these rely on an initial maximum likelihood estimator of $\bm{\beta}$, and hence do not apply to the high-dimensional setting. For estimating $h^{2}$ it should be noticed that HiLMM [10] aims to compute a confidence interval for $h^{2}$ as well. For that purpose their direct estimator (8) is likely more useful than MML on the pair $(\tau^{2},\sigma^{2})$. We also used Esther [28], which precedes HiLMM by sure independence screening. It did not improve HiLMM in our (semi-)sparse settings, and requires manual steps. However, it likely improves HiLMM results in very sparse settings [28].

For mixed effect models with a small number of fixed effects, MML compares fairly well to REML, with a larger bias, but smaller variance. Probably the potential advantage of contrasting out the fixed effects is small when the number of random effects is large. REML may have a larger advantage in very sparse settings [29] or when the number of fixed effects is large with respect to $n$. Estimates from the conjugate Bayes model are very similar to those by MML. We show that estimating $\tau^{2}$ along with $\sigma^{2}$ highly improves the $\sigma^{2}$ estimates presented by [24], where a fixed value of $\tau^{2}$ is used. In the case of many variance components or multiple similar regression equations, Bayesian extensions that shrink the estimates by a common prior are appealing, in particular in combination with efficient posterior approximations such as variational Bayes [30].

Our model (1) implies a dense setting, but we have demonstrated that the MML and REML estimators of $\tau^{2}$ and $\sigma^{2}$ are fairly robust against moderate sparsity, which corroborates the results by [29]. Nevertheless, true sparse models may be preferable when variable selection is desired, which depends on accurate estimation of $\bm{\beta}$. On the other hand, post-hoc selection procedures can be rather competitive [31]. Moreover, the sparsity assumption is questionable for several applications. E.g. in genetics, it was suggested that many complex traits (such as height or cholesterol levels) are not even polygenic, but instead “omnigenic” [32].

The extension of MML to high-dimensional GLM settings (22) is promising given its computational efficiency and performance for Poisson and Binomial regression. A special case of the latter, logistic regression, requires further research, because the Laplace approximations of the marginal likelihood are less accurate here [8]. Extension to survival is a promising avenue, because Cox regression is directly linked to Poisson regression [33]. Alternatively, parametric survival models may be pursued. To what extent the estimates of hyper-parameters impact predictions depends on the sensitivity of the likelihood to these parameters. For the linear setting, a re-analysis of the weight-gain data showed that predictions based on $\hat{\lambda}_{\text{MML}}$ improved those based on $\hat{\lambda}_{\text{CV}}$.

The MML estimator can be extended to estimation of multiple variance components or penalty parameters, which was addressed by iterative likelihood minorization [34] and by parameter-based moment estimation [27]. The latter extends to non-Gaussian response such as survival or binary. Further comparison of these methods with multi-parameter MML, both in terms of performance and computational efficiency, is left for future research. Finally, in particular in genetics applications, extensions of estimation of variance components by MML to non-independent individuals can be implemented by use of a well-structured between-individual covariance matrix $\boldsymbol{\Sigma}$ [3].

Although our simulations cover a fairly broad spectrum of settings, many other variations could be of interest. We therefore supply fully annotated R scripts https://github.com/markvdwiel/Hyperpar that allow i) comparison of all algorithms discussed here, also for one’s ‘own’ real covariate set $\boldsymbol{X}$; and ii) reproduction of all results presented here.

Acknowledgement

Gwenaël Leday was supported by the Medical Research Council, grant number MR/M004421. We thank Jiming Jiang and Can Yang for their correspondence, input and software for the MM algorithm. In addition, we thank Kristoffer Hellton for providing the fridge software and data. Finally, Iuliana Ciocǎnea-Teodorescu is acknowledged for preparing the TCGA KIRC data.

Supplementary Material

7.1 Contents

This Supplementary Material contains:

Details on the estimation of $\sigma^{2}$ and $\tau^{2}$ with the conjugate Bayesian model

• 2.

Estimation results from this Bayesian model for a simulation by [24]

• 3.

Details on the TCGA KIRC and TCPA OV gene and protein expression data

• 4.

An alternative analysis of the Weight gain data

• 5.

Supplementary Figures:

• (a)

Results for $\lambda$ estimation for Binomial ridge regression

• (b)

Robustness of estimates for non-Gaussian $\beta$’s and errors

7.2 Bayesian linear regression

7.2.1 Method

The conjugate Bayesian linear regression model is:

[TABLE]

where $\boldsymbol{X}$ is the design matrix, $\boldsymbol{\beta}$ is the vector of unknown regression parameters, $\epsilon$ is the vector of random errors and the prior shape and rate hyper-parameters $a$ and $b$ are fixed (say $a=1$ and $b=0.001$) so as to induce a non-informative prior. In model (24) the ridge regularization parameter corresponds to $\nu=\tau^{-2}$. In the following we provide the marginal posterior distribution of $\boldsymbol{\beta}$ and $\sigma^{-2}$, as well as a reformulation of the marginal likelihood [23] that allows important computational savings.

The likelihood function and prior densities are

[TABLE]

Therefore, the joint posterior is given by

$\pi(\boldsymbol{y}|\boldsymbol{\beta},\sigma^{-2})\pi(\boldsymbol{\beta}|\sigma^{-2},\tau^{-2})\pi(\sigma^{-2})=$

$\displaystyle{(2\pi)^{-\frac{n+p}{2}}\frac{b^{a}}{\Gamma\left(a\right)}(\sigma^{-2})^{\frac{n+p+2a+2}{2}}(\tau^{-2})^{\frac{p}{2}}\text{exp}\left\{-\frac{1}{2}\sigma^{-2}\left[(\mathbf{y}-\boldsymbol{X}\boldsymbol{\beta})^{T}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})+\tau^{-2}\boldsymbol{\beta}^{T}\boldsymbol{\beta}+2b\right]\right\}}$

The marginal posterior distribution of $\boldsymbol{\beta}$ is recognized to be a Student distribution with $2a+n$ degrees of freedom:

[TABLE]

Here $\boldsymbol{\beta}^{\ast}_{\nu}=V^{\ast}_{\nu}\boldsymbol{X}^{T}y$, $\displaystyle{\Sigma^{\ast}_{\nu}=\left(\frac{b^{\ast}_{\nu}}{a^{\ast}}\right)V^{\ast}_{\nu}}$ and $V^{\ast}_{\nu}=\left(\boldsymbol{X}^{T}\boldsymbol{X}+\nu\boldsymbol{\mathbb{I}}_{p}\right)^{-1}$. The marginal posterior distribution of $\sigma^{-2}$ is a Gamma distribution:

[TABLE]

where $a^{\ast}=a+0.5n$ and $b^{\ast}_{\nu}=\displaystyle{b+0.5\left(\boldsymbol{y}^{T}\boldsymbol{y}-{\boldsymbol{\beta}^{\ast}_{\nu}}^{T}{\Sigma^{\ast}_{\nu}}^{-1}\boldsymbol{\beta}^{\ast}_{\nu}\right)}$ The marginal likelihood of the model is:

[TABLE]

The marginal likelihood in (26) involves the determination of $|\Sigma^{\ast}_{\nu}|$, which can be computationally demanding when the number of variables $p$ is large and when we wish to evaluate the marginal likelihood for various values of $\nu$. To tackle this problem it is helpful to consider the singular value decomposition $UDV^{T}$ of $\boldsymbol{X}$, (where $U$ and $V$ are respectively $n\times q$ and $p\times q$ orthogonal matrices and $D=\text{diag}(d_{1},\ldots,d_{q})$ is the diagonal matrix of singular values $d_{1}>\ldots>d_{q}$ with $q=\text{min}(n,p)$) and focus on the linear model $\boldsymbol{y}=^{d}\mathcal{N}(F\boldsymbol{\theta},\sigma^{2}\mathbb{I}_{n})$, where $F=UD$ and $\boldsymbol{\theta}=V^{T}\boldsymbol{\beta}$, instead of $y=^{d}\mathcal{N}(\boldsymbol{X}\boldsymbol{\beta},\sigma^{2}\mathbb{I}_{n})$. Using simple algebra it can be shown that $\mathbb{E}\left[\boldsymbol{\theta}|\boldsymbol{y}\right]=\boldsymbol{\theta}^{\ast}_{\nu}=V^{T}\boldsymbol{\beta}^{\ast}_{\nu}$ and $\mathbb{V}\left[\boldsymbol{\theta}|\boldsymbol{y}\right]=V^{T}\Sigma^{\ast}_{\nu}V=\left(D^{2}+\nu\mathbb{I}_{q}\right)^{-1}$, which suggest that the marginal likelihood (that is invariant to linear transformations) is more easily determined in the orthogonalized space. Indeed, we now have that

[TABLE]

Additionally, using the fact that

[TABLE]

where $\hat{\boldsymbol{\theta}}=\left(F^{T}F\right)^{-1}F^{T}y=D^{-1}U^{T}y$, we deduce that

[TABLE]

and

[TABLE]

Finally, the log-marginal likelihood reduces to

[TABLE]

with constant $\displaystyle{C=a\log b+\log\Gamma(a^{\ast})-\log\Gamma(a)-\frac{n}{2}\log\pi}$. Expression (31) makes the evaluation of $\log\pi(\boldsymbol{y})$ very efficient for different values of $\nu$. Furthermore, the function $\log\pi(\boldsymbol{y})$ is $\log$-concave in $\nu$, which facilitates the use of numerical methods over grid-search approaches. In practice, it might be good to employ a numerical algorithm on a sub-domain of $\mathbb{R}^{+}$ obtained from a rough grid-search.

7.2.2 Results for simulation by [24]

The $\boldsymbol{\beta}$ prior in the conjugate Bayesian model (24) is scale-invariant, because it contains the error variance $\sigma^{2}$. Recently, it was criticized in [24] for its failure to estimate $\sigma^{2}$ well when$\tau^{2}$ is fixed and misspecified. In reality, $\tau^{2}$ is often estimated as well, which may improve results. To study this, we repeated the linear regression simulation by [24], which has the following specifications: $n=100,p=90,$ with $(\beta_{1},\ldots,\beta_{6})=(-2.5,-2,-1.5,1.5,2,2.5)$, and $\beta_{7}=\beta_{8}=\ldots\beta_{90}=0.$ Moreover $\boldsymbol{X}$ was generated from the independent standard Normal, and error variance $\sigma^{2}=3.$

Figure 7 shows the results for three methods. Here, ML serves as a benchmark and refers the ordinary maximum likelihood estimator of $\sigma^{2}$, referred to by [24] as ‘Least Squares’. In addition, BayesEB and BayesFix estimate $\sigma$ by $(b^{\ast}_{\nu}/(a^{\ast}-1))^{1/2},$ i.e. the square-root of the posterior mean of $\sigma^{2}=1/\sigma^{-2}$ (25), with $\nu=1/\tau^{2}$ estimated by empirical Bayes (EB; maximizing (31)) and with $\nu$ fixed to $1/100$ (as in [24]), respectively.

As in [24], the results show that indeed the conjugate Bayes model is not robust against wrongly fixing $\tau^{2}$. However, it also shows that EB estimation of $\tau^{2}$ strongly improves the $\sigma$ estimate. In fact, these estimates are very competitive to the ones for the scale-independent prior, $\beta\sim\mathcal{N}(0,\tau^{2}I_{p})$, which was advocated by [24].

7.3 Details on real data used for simulations

7.3.1 Kidney tumor gene expression

The TCGA KIRC data is downloaded from the harmonised GDC database http://cancergenome.nih.gov/ using the R-package TCGAbiolinks [35]. It contains RNAseq profiles of $n=71$ kidney renal clear cell carcinomas. Only those genes were retained that had more than two counts per million in at least three samples, rendering $p=18,391$ genes. The data were normalised using the trimmed mean of M-values method of the R-package edgeR [36]. Following common practice, all gene expression values were standardized.

7.3.2 Ovarian tumor protein expression

We use protein expression data from the cancer proteome atlas (TCPA; [37]), which holds 408 ovarian serous cystadenocarcinoma profiles measuring 224 proteins by reverse-phase protein arrays. These are available from the TCPA portal: https://www.tcpaportal.org/tcpa/. Data were normalized by median centering per sample. All protein expression values were standardized.

7.4 Weight gain data example: alternative analysis

In the main document we re-analysed the weight gain data recently discussed in [17]. The authors kindly provided the R-code corresponding to the results in [17], which includes their focused ridge (fridge) methodology. Below we report a small discrepancy between our analysis and theirs, and also present results of an alternative analysis.

The authors of [17] opted to present results on $n=25$ (rather than $n=26$) samples, for technical reasons (personal communication). Results differed very little, so we chose to use all $n=26$ samples of the original study [20]. In addition, the R-code by [17] includes a prior selection of those 1,000 genes with the largest absolute marginal correlation correlation to the response, i.e. weight gain. Such a prior selection can indeed enhance performance of ridge-type predictors. The analysis by [17], however, did not include the gene selection as part of the outer leave-one-out cross-validation loop for assessing predictive performance. This may lead to over-optimism of the prediction errors. We therefore repeated the analysis and comparison as presented in the main document, but with the gene selection as part of the outer CV-loop. Indeed the error estimates increased substantially: $\text{MSPE}_{\text{MML}}=38.65$ (was 14.40), $\text{MSPE}_{\text{GCV}}=44.20\ (16.38)$ and $\text{MSPE}_{\text{fridge}}=41.65\ (15.80).$ Nevertheless, the conclusion remains that the MML and fridge MSPEs are lower that those of GCV, with relative decreases of 12.5% and 5.8%, respectively. Figure 8 displays absolute prediction errors per sample and illustrates the improved prediction by ridge using $\lambda_{\text{MML}}$ with respect to ridge using $\lambda_{\text{GCV}}$.

7.5 Supplementary Figures

7.5.1 Robustness: non-Gaussian $\beta$’s and errors

Here, we show additional results for the simulation settings with either non-Gaussian $\beta$’s (spike-and-slab or Uniform) or errors ($t_{4}$), as presented in the main document. In all cases, parameters of the $\beta$ prior or the error ($\epsilon_{i}$) distribution were set such that $\tau^{2}=\mathbb{V}(\beta_{j})=0.01$ and $\sigma^{2}=\mathbb{V}(\epsilon_{i})=10$.

7.5.2 Binomial ridge regression

Note that unlike Poisson ridge regression Binomial ridge regression is not implemented in penalized. Hence, we compare mgcv-based results only with glmnet. Figure 10 displays the results for estimating $\lambda$, in case binomial $N=5$. Results did not qualitatively differ for $N=3$ and $N=10$. While MML is fairly good on target, the estimates from GCV are roughly 10 times too small, whereas the estimates from CV by glmnet are roughly 2-3 times too large, with outliers towards a 10-100 fold overestimation.