Under the hood of REGENIE

1 Key idea

• Another fast approach to fit the GWAS-LMM.
• Builds on the developments in BOLT-LMM

1. Project known covariates \(\mathbf{X}_c\) out
2. Take BOLT-LMM-inf idea to use LOCO-residuals under the null model.
3. Exploit link between mixed models and ridge regression for fast approximation of LOCO-BLUP under \(H_0\).
4. Model LOCO residuals using candidate SNP.

\[ \mathbf{y}_\text{residual-LOCO, H_0} = \tilde{\mathbf{x}}_\text{test} \boldsymbol{\beta}_\text{test} + \epsilon \]

5. Use a score test similar to BOLT-LMM-inf

\[ T = \frac{(\tilde{\mathbf{x}}^T_\text{test}\mathbf{y}_\text{residual-LOCO, H_0})^2}{\hat\sigma^2_\epsilon\tilde{\mathbf{x}}^T\tilde{\mathbf{x}}} \]

6. Estimate \(\hat{\sigma}^2_\epsilon\) under the null model:

\[ \hat\sigma^2_\epsilon = \frac{\vert\vert\mathbf{y}_\text{residual-LOCO, H_0}\vert\vert^2_2}{N-C} \]

• In contrast to BOLT-LMM-inf no calibration factor used in denominator
• they “found in applications that the results obtained using this simple form match up closely to those using a calibration factor”

Key contribution efficiently obtain \(\mathbf{y}_\text{residual-LOCO, H_0}\)

• Check method section paper
• As usual in publications: to understand the methods, it is all about the supplement

2 Intermezzo ridge regression

• Penalized regression useful in a high dimensional context when the number of covariates are larger than the number of samples (e.g. M >> N)
• When covariates are highly correlated

The ridge parameter estimator is defined as the parameter \(\mathbf\beta\) that minimises the penalised least squares loss function

\[ \text{SSE}_\text{pen}=\Vert\mathbf{Y} - \mathbf{X\beta}\Vert_2^2 + \lambda \Vert \boldsymbol{\beta} \Vert_2^2 \]

• \(\Vert \boldsymbol{\beta} \Vert_2^2=\sum_{j=1}^p \beta_j^2\) is the \(L_2\) penalty term

• \(\lambda>0\) is the penalty parameter (to be chosen by the user).

Note, that that is equivalent to minimizing \[ \Vert\mathbf{Y} - \mathbf{X\beta}\Vert_2^2 \text{ subject to } \Vert \boldsymbol{\beta}\Vert^2_2\leq s \]

Note, that \(s\) has a one-to-one correspondence with \(\lambda\)

2.1 Graphical interpretation

• Ridge estimator can switch sign as compared to OLS.

2.2 Estimator

The loss function to be minimised is \[ L_\text{ridge}(\mathbf{Y},\boldsymbol{\beta},\lambda) = \text{SSE}_\text{pen}=\Vert\mathbf{Y} - \mathbf{X\beta}\Vert_2^2 + \lambda \Vert \boldsymbol{\beta} \Vert_2^2. \]

First we re-express the loss function in matrix notation: \[ L_\text{ridge}(\mathbf{Y},\boldsymbol{\beta},\lambda) = (\mathbf{Y}-\mathbf{X\beta})^T(\mathbf{Y}-\mathbf{X\beta}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}. \]

Solving \(L_\text{ridge}(\mathbf{Y},\boldsymbol{\beta},\lambda)=0\) gives

[ \[\begin{array}{rcl} \frac{\partial}{\partial \boldsymbol{\beta}}L_\text{ridge}(\mathbf{Y},\boldsymbol{\beta},\lambda) &=& 0 \\\\ -2\mathbf{X}^T(\mathbf{Y}-\mathbf{X\beta})+2\lambda\boldsymbol{\beta} &=& 0\\\\ \hat{\boldsymbol{\beta}} &=& (\mathbf{X^TX}+\lambda \mathbf{I})^{-1} \mathbf{X^T Y} \end{array}\]

]

It can be shown that \((\mathbf{X^TX}+\lambda \mathbf{I})\) is always of rank \(p\) if \(\lambda>0\).

Hence, \((\mathbf{X^TX}+\lambda \mathbf{I})\) is invertible and \(\hat{\boldsymbol{\beta}}\) exists even if \(p>>>n\).

2.3 Link with SVD

Write the SVD of \(\mathbf{X}\) (\(p>N\)) as \[ \mathbf{X} = \sum_{i=1}^N \delta_i \mathbf{u}_i \mathbf{v}_i^T = \sum_{i=1}^p \delta_i \mathbf{u}_i \mathbf{v}_i^T = \mathbf{U}\boldsymbol{\Delta} \mathbf{V}^T , \] with

• \(\delta_{n+1}=\delta_{n+2}= \cdots = \delta_p=0\)

• \(\boldsymbol{\Delta}\) a \(p\times p\) diagonal matrix of the \(\delta_1,\ldots, \delta_p\)

• \(\mathbf{U}\) an \(n\times p\) matrix and \(\mathbf{V}\) a \(p \times p\) matrix. Note that only the first \(n\) columns of \(\mathbf{U}\) and \(\mathbf{V}\) are informative.

With the SVD of \(\mathbf{X}\) we write \[ \mathbf{X}^T\mathbf{X} = \mathbf{V}\boldsymbol{\Delta }^2\mathbf{V}^T. \] The inverse of \(\mathbf{X}^T\mathbf{X}\) is then given by \[ (\mathbf{X}^T\mathbf{X})^{-1} = \mathbf{V}\boldsymbol{\Delta}^{-2}\mathbf{V}^T. \] Since \(\boldsymbol{\Delta}\) has \(\delta_{n+1}=\delta_{n+2}= \cdots = \delta_p=0\), it is not invertible.

Note that it can be shown that \[ \mathbf{X^TX}+\lambda \mathbf{I} = \mathbf{V} (\boldsymbol{\Delta}^2+\lambda \mathbf{I}) \mathbf{V}^T , \] i.e. adding a constant to the diagonal elements does not affect the eigenvectors, and all eigenvalues are increased by this constant.
\(\longrightarrow\) zero eigenvalues become \(\lambda\).

Hence, \[ (\mathbf{X^TX}+\lambda \mathbf{I})^{-1} = \mathbf{V} (\boldsymbol{\Delta}^2+\lambda \mathbf{I})^{-1} \mathbf{V}^T , \] which can be computed even when some eigenvalues in \(\boldsymbol{\Delta}^2\) are zero.

Note, that for high dimensional data (\(p>>>N\)) many eigenvalues are zero because \(\mathbf{X^TX}\) is a \(p \times p\) matrix and has rank \(N\).

2.4 Properties ridge

• The Ridge estimator is biased! The \(\boldsymbol{\beta}\) are shrunken to zero! \[\begin{eqnarray} \text{E}[\hat{\boldsymbol{\beta}}] &=& (\mathbf{X^TX}+\lambda \mathbf{I})^{-1} \mathbf{X}^T \text{E}[\mathbf{Y}]\\\\ &=& (\mathbf{X}^T\mathbf{X}+\lambda \mathbf{I})^{-1} \mathbf{X}^T \mathbf{X}\boldsymbol{\beta}\\ \end{eqnarray}\]

• Note, that the shrinkage is larger in the direction of the smaller eigenvalues.

\[\begin{eqnarray} \text{E}[\hat{\boldsymbol{\beta}}]&=&\mathbf{V} (\boldsymbol{\Delta}^2+\lambda \mathbf{I})^{-1} \mathbf{V}^T \mathbf{V} \boldsymbol{\Delta}^2 \mathbf{V}^T\boldsymbol{\beta}\\\\ &=&\mathbf{V} (\boldsymbol{\Delta}^2+\lambda \mathbf{I})^{-1} \boldsymbol{\Delta}^2 \mathbf{V}^T\boldsymbol{\beta}\\\\ &=& \mathbf{V} \left[\begin{array}{ccc} \frac{\delta_1^2}{\delta_1^2+\lambda}&\ldots&0 \\ &\vdots&\\ 0&\ldots&\frac{\delta_r^2}{\delta_r^2+\lambda} \end{array}\right] \mathbf{V}^T\boldsymbol{\beta} \end{eqnarray}\]

• Ridge regression can lead to parameters that switch sign if penality increases

2.5 Toxicogenomics example

• N = 30 chemical compounds have been screened for toxicity

• Bioassay data on toxicity screening

• Gene expressions in a liver cell line are profiled for each compound (M = 4000 genes)

• Predict Bioassay score in function of gene expression.

toxData <- read_csv(
  "https://raw.githubusercontent.com/statOmics/HDA2020/data/toxDataCentered.csv",
  col_types = cols()
)
dim(toxData)

## [1]   30 4001

lmFit <- lm (BA ~. , toxData)
lmFit %>%  
  coef %>% 
  head(40)

##   (Intercept)            X1            X2            X3            X4 
## -2.059917e-17 -7.456994e+00  3.571348e-01  1.124923e+01  1.083540e+01 
##            X5            X6            X7            X8            X9 
## -1.374339e+01  5.683387e+00  6.553878e+01  4.340456e+00  7.910392e+00 
##           X10           X11           X12           X13           X14 
##  3.702961e+01 -5.483687e+01 -5.555478e+01  5.792467e+00  2.314280e+01 
##           X15           X16           X17           X18           X19 
## -6.961036e+00 -2.852506e+01 -2.255090e+01 -9.796237e+01 -3.041718e+01 
##           X20           X21           X22           X23           X24 
## -3.269917e+01 -1.428088e+01 -1.614313e+01 -2.274987e+01  7.316352e+01 
##           X25           X26           X27           X28           X29 
## -5.706583e+00  3.747454e+01 -2.019991e+01  1.499068e+01  9.960810e+01 
##           X30           X31           X32           X33           X34 
##            NA            NA            NA            NA            NA 
##           X35           X36           X37           X38           X39 
##            NA            NA            NA            NA            NA

lmFit %>% 
  coef %>% 
  is.na %>% 
  sum

## [1] 3971

lmFit %>% 
  coef %>% 
  is.na %>% 
  `!` %>% 
  sum 

## [1] 30

library(glmnet)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loaded glmnet 4.1-8

mRidge <- glmnet(
  x = toxData[,-1] %>%
    as.matrix,
  y = toxData %>%
    pull(BA),
  alpha = 0) # ridge: alpha = 0

plot(mRidge, xvar="lambda")

2.6 No penalisation of some parameters

Note, that we typically do not penalise the intercept. We can do this by introducing matrix

\[ \mathbf{D} = \left[\begin{array}{ccc}0& 0 \\ 0&\mathbf{I}\end{array}\right] \] And let \[ \mathbf{C} = \left[ \begin{array}{cc} 1 & \mathbf{1}\\\mathbf{1}&\mathbf{X}_{1\ldots p}\end{array}\right] \text{ and } \boldsymbol{\theta} =\left[\begin{array}{c}\beta_0\\ \boldsymbol{\beta}_{1\ldots p}\end{array}\right] \]

with \(\mathbf{X}_{1\ldots p}\) the matrix of the predictors for which the slope terms \(\boldsymbol {\beta}_{1\ldots p}\) are estimated.

The loss function then becomes

\[ L_\text{ridge}(\mathbf{Y},\boldsymbol{\beta},\lambda) = (\mathbf{Y}-\mathbf{C\theta})^T(\mathbf{Y}-\mathbf{C\theta}) + \lambda \boldsymbol{\beta}^T\mathbf{D}\boldsymbol{\theta}. \]

and the ridge estimator then becomes

\[ \hat{\boldsymbol{\beta}} = (\mathbf{C^TC}+\lambda \mathbf{D})^{-1} \mathbf{C^T Y} \]

2.7 Link LMM

Add Gaussian prior on the model parameters/specify the fixed effects as random effects

\[ \begin{array}{ccl} \mathbf{Y} &=& \mathbf{1}\beta_0 + \mathbf{X}_{1\ldots p} \boldsymbol{\beta}_{1\ldots p} + \boldsymbol{\epsilon}\\ \boldsymbol{\beta}_{1\ldots p} &\sim& \text{MVN}(0,\mathbf{I}\sigma^2_\beta) \\ \boldsymbol{\epsilon} &\sim& \text{MVN}(0,\mathbf{I}\sigma^2_\epsilon) \\ \end{array} \]

Note that we know from LMM theory that the BLUP is

\[ \hat{\boldsymbol{\theta}} = (\mathbf{C}^T\mathbf{C} + \sigma^2_\epsilon H)^{-1}\mathbf{C}^T\mathbf{Y} \] with

\[ \mathbf{H}=\left[\begin{array}{cc} \mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{G}^{-1} \end{array}\right] = \left[\begin{array}{cc} \mathbf{0}&\mathbf{0}\\ \mathbf{0}&\sigma^{-2}_{\beta}\mathbf{I} \end{array}\right] = \sigma^{-2}_\beta \mathbf{D} \]

So the BLUP reduces to

\[ \hat{\boldsymbol{\theta}} = (\mathbf{C}^T\mathbf{C} + \frac{\sigma^2_\epsilon}{\sigma^2_\beta} \mathbf{D})^{-1}\mathbf{C}^T\mathbf{Y} \]

which is equivalent to the ridge solution!

• \(\frac{\sigma^2_\epsilon}{\sigma^2_\beta}\) plays the role of ridge penalty \(\lambda\) in ridge regression.

• we can exploit the mixed model machinery to perform ridge regression and to estimate the penalty parameter using variance components.

\(\rightarrow\) Regenie is exploiting this link to avoid computational complexity of fitting LMMs.
\(\rightarrow\) Revisit to supplement

3 Remarks

• Use of R different penalties \(\lambda_1 \ldots \lambda_R\) as a proxy to allow for different random effect variances
• Construct R predictions for phenotype with different penalties per block
• Combine the predictors in a stacked predictor
• I feel that this also allows deviations of the infinitesimal model
• Danger that it is not well calibrated: indeed nominator

\[ \hat\sigma^2_\epsilon\tilde{\mathbf{x}}^T\tilde{\mathbf{x}} \leftrightarrow \tilde{\mathbf{x}}^T\hat{\mathbf{V}}^{-1}_{\text{LOCO},H_0}\tilde{\mathbf{x}} \]

---

This work is licensed under the CC BY-NC-SA 4.0 licence.