# DEViaN-LM: An R Package for Detecting Abnormal Values in the Gaussian Linear Model

Geoffroy Berthelot, Guillaume Saulière, Jérôme Dedecker

<figure><img src="image_1.png" /></figure>

## To cite this version:

Geoffroy Berthelot, Guillaume Saulière, Jérôme Dedecker. DEViaN-LM: An R Package for Detecting Abnormal Values in the Gaussian Linear Model. 2025. hal-05230549v2

HAL Id: hal-05230549

https://hal.science/hal-05230549v2

Preprint submitted on 18 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DEViaN-LM: An R Package for Detecting Abnormal Values in the Gaussian Linear Model

Geoffroy Berthelot 1,2, Guillaume Saulière 1, and Jérôme Dedecker 3

¹Institut de Recherche bioMédicale et d'Epidémiologie du Sport (IRMES, UPR7329), INSEP (Institut National du Sport, de l'Expertise et de la Performance), Paris, France

²Research Laboratory for Interdisciplinary Studies (RELAIS), 89100 Sens, France

³Université Paris Cité, CNRS, UMR 8145, MAP5, F-75006 Paris, France

## Abstract

The **DEViaN-LM** is a **R** package that allows to detect the values poorly explained by a Gaussian linear model. The procedure is based on the maximum of the absolute value of the studentized residuals, which is a free statistic of the parameters of the model. This approach makes it possible to generalize several procedures used to detect abnormal values during longitudinal monitoring of certain biological markers. In this article, we describe the method used, and we show how to implement it on different real datasets.

## 1 Introduction

The detection of outliers or abnormal values in a series of observations is one of the core problems in many applied fields. These abnormal values may represent errors or events of interest such as a novelty. Outliers detection finds application in different areas including, but not limited to: fraud detection, activity monitoring, fault diagnosis, structural defect detection, satellite image analysis, medical condition monitoring, pharmaceutical research, detection of unexpected entries in database or mislabeled data in a training data set [1].

In medicine, detecting abnormal values is typically useful in the biological monitoring of individuals, where blood or urine markers are collected on a temporal basis. Usually, biological values of a given individual are collected and compared to the statistical parameters of a population, which is deemed representative of the biological characteristics of the tested individual. However, [2] and [3] underlined that the inter-individual variability is higher than the intra-individual one, suggesting that the use and design of population data for estimating intra-individual variability is an issue. In line with [4], our approach rather aims at detecting if some observations are abnormal, based on the characteristics of the individual only. This is particularly relevant in the context of personalized/precision medicine which favors an individualized approach: "Precision Medicine refers to the tailoring of medical treatment to the individual characteristics of each patient" [5].

Let us describe the simplest procedure, based on a Z-score, to detect if the observation $x_n$ is abnormal, when $x_n$ is the last observation from a sequence of independent and identically distributed (iid) Gaussian random variables $X_1, \dots, X_n$. In that case, the Z-score statistic is simply defined by

$$ZS_n = \frac{X_n - \bar{X}_{n-1}}{\hat{\sigma}_{n-1} \sqrt{1 + \frac{1}{n-1}}}$$

where

$$\bar{X}_{n-1} = \frac{1}{n-1} \sum_{k=1}^{n-1} X_k \qquad \hat{\sigma}_{n-1}^2 = \frac{1}{n-2} \sum_{k=1}^{n-1} (X_k - \bar{X}_{n-1})^2$$

The Z-score provides an intuitive mechanism to determine the magnitude by which the observation $x_n$ diverges from the others in a sequence. If this divergence is large enough, the observation can be deemed
as an outlier. This procedure works well, because in this Gaussian context, the distribution of $ZS_n$ is the Student($n-2$) distribution; in particular, it does not depends on the unknown quantities $\mu = \mathbb{E}(X_i)$ and $\sigma^2 = \text{Var}(X_i)$.

The Z-score has found applications in biology ([6]), neuroimaging ([7]), psychology ([8]), sport ([4]) and medicine ([9]) among others. Let us note that $ZS_n$ is in fact the $n$-th studentized residuals in the simple Gaussian linear model $X_i = \mu + \varepsilon_i$, where $\varepsilon_1, \dots, \varepsilon_n$ is iid with $N(0, \sigma^2)$ marginal distribution. Starting from this simple fact, [4] have generalized this approach to the problem of detecting abnormal values in the series $x_1, \dots, x_n$; they proposed to use as a statistic the maximum of the absolute value of the studentized residuals, providing a table of the quantiles of this statistic.

In this paper, we generalize the approach of [4] to the Gaussian linear model. As a consequence of Proposition 1 below, we get that the distribution of the maximum of the absolute value of the studentized residuals does not depends of the unknown parameters of the linear model. However, it depends on the design matrix, so that a universal quantiles table cannot be used. Our second objective in then to propose the **DEViaN-LM** package (Detecting Extremal Values in a Normal Linear Model) that allows, for any given design matrix, to estimate the quantiles of this statistic via Monte-Carlo simulations, and thus to detect all the abnormal observations for a given risk level.

The paper is organized as follows: in Section 2 we describes in details the statistical procedure, an we prove the main result (Proposition 1). In Section 3 we present some applications of the **DEViaN-LM** package to some real data sets. In Section 4 we show that **DEViaN-LM** has been time-optimized, comparing it to direct procedures using **R** or **MATLAB** functions.

## 2 Detecting abnormal values from studentized residuals

We consider the homoscedastic linear model

$$X = M\theta + \varepsilon \quad (1)$$

where $M$ is an $n \times p$ random matrix, $\varepsilon = (\varepsilon_1, \dots, \varepsilon_n)'$ is a Gaussian random vector independent of $M$, whose coordinates are independent and identically distributed (with mean 0 and variance $\sigma^2$), and $\theta = (\theta_1, \dots, \theta_p)'$ is the unknown vector of parameters.

One observes the random vector $X = (X_1, \dots, X_n)'$ and the random design $M$. In this context we want to detect “abnormal observations”, that is observations $x_i$'s that are poorly explained by the model.

Before describing our method in more details, let us mention that we shall work conditionally on the design matrix $M$, which is a standard approach when dealing with Model (1). Consequently, for the sake of simplicity and without loss of generality, we assume from now that $M$ is fixed.

Assume that $n > p + 1$ and that the columns of $M$ are linearly independent (so that $M'M$ is invertible). Let $M_{(i)}$ be the matrix $M$ deprived of its $i$-th row, and $X_{(i)}$ be the vector $X$ deprived of its $i$-th coordinate, and assume that for any $i \in \{1, \dots, n\}$ the matrix $M'_{(i)}M_{(i)}$ is invertible. Let now $\hat{\theta}_{n,i}$ be the least squares estimator of $\theta$ based on the variables $X_j$'s, for $j \neq i$:

$$\hat{\theta}_{n,i} = (M'_{(i)} M_{(i)})^{-1} M'_{(i)} X_{(i)}$$

Let $L_i$ be the $i$-th row of the matrix $M$, and let $\hat{X}_{n,i}$ be the prediction of $X_i$ obtained from $\hat{\theta}_{n,i}$, that is $\hat{X}_{n,i} = L_i \hat{\theta}_{n,i}$. Standard computations show that

$$\text{Var}(X_i - \hat{X}_{n,i}) = \sigma^2 (1 + L_i (M'_{(i)} M_{(i)})^{-1} L'_{i}).$$

Let

$$\hat{\sigma}_{n,i}^2 = \frac{1}{n-p-1} \sum_{k \in \{1, \dots, n\}, k \neq i} (X_k - L_k \hat{\theta}_{n,i})^2$$

be the estimator of $\sigma^2$ based on the variables $X_j$'s, for $j \neq i$. Define then

$$\hat{e}_i(X) = \frac{X_i - \hat{X}_{n,i}}{\hat{\sigma}_{n,i} \sqrt{1 + L_i (M'_{(i)} M_{(i)})^{-1} L'_{i}}}.$$
It is easy to see that the $\hat{e}_i(X)$'s are identically distributed (but not independent), with common distribution $St(n-p-1)$. In fact, one can check that the $\hat{e}_i(X)$'s are the so-called "Studentized residuals" of Model (1) (see for instance [10], Exercise 4.4 p. 87).

We see that the distribution of the random variable $\hat{e}_i(X)$ does not depend on the unknown parameters $(\theta, \sigma^2)$. In fact, we shall prove in Proposition 1 a much stronger result: the distribution of the random vector $(\hat{e}_1(X), \dots, \hat{e}_n(X))'$ does not depend on $(\theta, \sigma^2)$ (but it depends on the design $M$). This important fact justifies the procedure below.

To see if there is at least one abnormal observation in the sequence $x_1, \dots, x_n$, we consider the statistic

$$T_n = \max_{i \in \{1, \dots, n\}} \left| \frac{X_i - \hat{X}_{n,i}}{\hat{\sigma}_{n,i} \sqrt{1 + L_i (M'_{(i)} M_{(i)})^{-1} L'_{i}}} \right| = \max_{i \in \{1, \dots, n\}} |\hat{e}_i(X)|.$$

From the discussion of the previous paragraph, we deduce that the distribution of the statistic $T_n$ does not depend on $(\theta, \sigma^2)$. Hence, for each design matrix $M$, it can be tabulated via a basic Monte-Carlo procedure. For a given risk level $\alpha \in (0, 1)$, we say that there is at least one abnormal observation in the sequence $x_1, \dots, x_n$ if the observed value $t_n$ of $T_n$ is such that $t_n > c_{\alpha,n}$, where $c_{\alpha,n}$ is the quantile of order $1 - \alpha$ of the distribution of $T_n$. As usual, this is equivalent to the following procedure: for a given risk level $\alpha \in (0, 1)$, we say that there is at least one abnormal observation in the sequence $x_1, \dots, x_n$ if the $p$-value $p_n = \mathbb{P}(T_n > t_n)$ is less than the risk $\alpha$.

It is not possible to provide a table of quantiles of $T_n$ for each possible design matrix $M$. The solution is then, for a given risk level $\alpha$ and a given design matrix $M$, to estimate via Monte-Carlo the quantile $c_{\alpha,n}$ and/or the $p$-value $p_n$. This is exactly what the package **DEViaN-LM** does.

Let us now give the main result of this section. Consider Model (1), and define the vector $Y$ by

$$Y = \frac{X - M\theta}{\sigma},$$

in such a way that the coordinates $(Y_i)_{1 \le i \le n}$ are independent and identically distributed, with $\mathcal{N}(0, 1)$ marginal distribution. Let also

$$\hat{\theta}_{n,i}(Y) = (M'_{(i)} M_{(i)})^{-1} M'_{(i)} Y_{(i)}, \quad \hat{\sigma}_{n,i}^2(Y) = \frac{1}{n-p-1} \sum_{k \in \{1, \dots, n\}, k \ne i} (Y_k - L_k \hat{\theta}_{n,i}(Y))^2,$$

and

$$\hat{Y}_{n,i} = L_i \hat{\theta}_{n,i}(Y), \quad \hat{e}_i(Y) = \frac{Y_i - \hat{Y}_{n,i}}{\hat{\sigma}_{n,i}(Y) \sqrt{1 + L_i (M'_{(i)} M_{(i)})^{-1} L'_{i}}}.$$

**Proposition 1** In Model (1), the following equality holds: for any $i \in \{1, \dots, n\}$

$$\hat{e}_i(X) = \hat{e}_i(Y).$$

Consequently, the distribution of $(\hat{e}_1(X), \dots, \hat{e}_n(X))'$ does not depend on $(\theta, \sigma^2)$.

**Remark 2** If we do not assume that the design matrix $M$ is fixed, then the last assertion of Proposition 1 writes as follows: the conditional distribution of $(\hat{e}_1(X), \dots, \hat{e}_n(X))'$ given $M$ does not depend on $(\theta, \sigma^2)$.

**Remark 3** The statistics $T_n$ is the sup-norm of the random vector $(\hat{e}_1(X), \dots, \hat{e}_n(X))'$. From Proposition 1, we see that the $\ell_p$-norm ($p \ge 1$) of $(\hat{e}_1(X), \dots, \hat{e}_n(X))'$ is also distribution-free, and could be used to detect if the series $x_1, \dots, x_n$ is abnormal. The advantage of $T_n$ is the following: once we have computed the quantile $c_{\alpha,n}$ from the distribution of $T_n$, then one can exhibit all the abnormal observations (if any) at level $\alpha$, that is the observations $x_i$ whose residuals $|\hat{e}_i(x)|$ are above the threshold $c_{\alpha,n}$.

**Proof of Proposition 1.** By definition of $Y_i$ and $\hat{Y}_{n,i}$, one can write

$$\begin{align}
Y_i - \hat{Y}_{n,i} &= \frac{X_i}{\sigma} - \frac{L_i \theta}{\sigma} - L_i (M'_{(i)} M_{(i)})^{-1} M'_{(i)} \frac{X_i}{\sigma} + L_i (M'_{(i)} M_{(i)})^{-1} M'_{(i)} M_{(i)} \frac{\theta}{\sigma} \nonumber \\
&= \frac{X_i}{\sigma} - L_i (M'_{(i)} M_{(i)})^{-1} M'_{(i)} \frac{X_i}{\sigma} \nonumber \\
&= \frac{X_i - \hat{X}_{n,i}}{\sigma} \tag{2}
\end{align}$$
In the same way

$$Y_k - L_k \hat{\theta}_{n,i}(Y) = \frac{X_k - L_k \hat{\theta}_{n,i}}{\sigma},$$

and consequently

$$\hat{\sigma}_{n,i}^2(Y) = \frac{\hat{\sigma}_{n,i}^2}{\sigma^2}. \qquad (3)$$

The result follows from (2), (3) and the definition of $\hat{e}_i(X)$ and $\hat{e}_i(Y)$. □

# 3 Applications to biological and sociodemographic data

The package **DEViaN-LM** provides a convenient way to detect the values poorly explained by a Gaussian linear model or, in other words, to detect abnormal values in the Gaussian framework. As described above, the package performs a series of simple operations in order to estimate the quantiles of $T_n$ via Monte-Carlo simulations for a given model design and level $\alpha$. The abnormal observations at level $\alpha$ can then be found. The package takes 5 input parameters: the variable of interest and the explanatory variables, the quantile of interest, the total number of simulations for estimating the quantiles and the number of CPU to use for the simulations. For example, in order to estimate the grip strength *df$grip\_strength* of mouse lemurs with regards to their age *df$age* and body mass *df$body\_mass* using the following model in R: *df$grip\_strength* ~ *df$age* + *df$body\_mass*, the package **DEViaN-LM** is called using the following code:

```R
R> y = df$grip\_strength
R> x = cbind(1, df$age, df$body\_mass)
R> model = devianlm_stats(y, x, n_sims = 1e5, quant = 0.95)
```

where the default input values are set to $10^5$ simulations and one processor. The user needs to build the model design and send it to the package. For instance, for a 2nd order polynomial regression, one can typically use:

```R
R> x <- rnorm(10)
R> e <- rnorm(10, 0, 0.5)
R> y <- 0.5 + 2 * x + x^2 + e
R> model = devianlm_stats(y, cbind(1, x, x^2))
```

Note that `devianlm_stats()` does not add an intercept in the design. If an intercept is needed (as in the preceding examples) the user has to include a column of 1 in the design.

The package returns 4 objects: a (double) vector containing the studentized residuals $\hat{e}_i(X)$, a (integer) vector with the indices (positions in the original data) of observations identified as outliers based on the threshold, the estimated threshold value $c_{\alpha,n}$ and a (integer) binary vector (0 or 1) indicating whether each observation is considered an outlier (1) or not (0). We illustrate the use of the **DEViaN-LM** in 4 different datasets: longitudinal sequences of biological markers in Elite soccer players, wage of a sample of the US population, grip-strength in a mouse lemur cohort and wheel activity in a cohort of mice. In each of these datasets, we set the number of simulations to $10^5$ simulations.

## 3.1 Elite soccer players

The description of the original dataset is provided in [4]. New samples were added and the new dataset now consists in a database of 5 typical biological markers from 3936 male soccer players from the French elite leagues 1 & 2. These biological markers include concentrations of ferritin (mol/L), serum iron ($\mu$mol/L), hemoglobin (g/L), erythrocytes (T/L) and hematocrit levels (%). They were collected every 6 months in July/August and in January/March from 2005 to 2017. The smallest interval between two measures is 6 months, allowing for independent sampling ([11]). The procedure for selecting and testing the data sample is identical as in [4]. We use **DEViaN-LM** for assessing the number of abnormal values of serum iron in the studied population.
## 3.2 US wage sample

The Current Population Survey (CPS) is used to supplement census information between census years. This data consist of a random sample of 599 persons from the 2012 CPS, with information on wages and other characteristics of the workers. We use **DEViaN-LM** for assessing abnormal values of the logarithm of the salary / wage, adjusted to the age, the educational level and the number of children. The linear model is: $log(salary) = \beta_0 + \beta_1age + \beta_2educational level + \beta_3nb$ of children

## 3.3 Grip strength of mouse lemurs

The maximal grip strength of mouse lemurs (*Microcebus murinus*) was measured using a small iron bar mounted on a piezo-electric force platform (Kistler squirrel force plate, ±0.1N; Winterthur, Switzerland), connected to a charge amplifier (Kistler charge amplifier type 9865). The pull strength was recorded during 60s at 1kHz. Animals performed the pull task several times during this interval. The task consisted of letting the animal grab the iron bar and pull them horizontally until they let it go. The maximum force was then extracted using the Bioware software (Kistler). The repeatability of this task was previously verified ([12]). A sample of 359 individuals was measured: 170 females and 189 males, with 39 individuals tested twice. A total of 397 grip-strengths (N) were measured. We use the **DEViaN-LM** for assessing the number of abnormal (log) values in the grip-strengths, in regards to age and body mass. The dataset consists in 340 grip-strengths, with corresponding age and body mass (57 body mass values are missing). See above for the equation used.

## 3.4 Wheel activity in mice

The maximal distance per day in the wheel activity is recorded for mice (*Musdomesticus*) based on their voluntary behavior to practice wheel activity ([13]). A total of 159 mice were genetically selected for high locomotor activity and 14241 daily wheel running performances (7078 for males and 7163 for females) were gathered ([14]). The distance run was converted to km/week and we use the **DEViaN-LM** for assessing abnormal distances with regards to the age and body mass. We select the first lineage, as the wheel running distances of this lineage appear more normally distributed compared to the distances of other lineages. The dataset consists in 1673 wheel running distances from 19 mice (9 females and 10 males) with corresponding age and body mass.

# 4 Software and benchmark

The package performs a series of simple operations in order to estimate the quantiles of $T_n$ via Monte-Carlo simulations for a given model design and level $\alpha$. The precision of these quantiles estimates increases with the number of simulations performed in the Monte-Carlo process. However, increasing the number of simulations can be both time- and resource-intensive. In this section we detail the typical methods for computing the quantiles of $T_n$ using Monte-Carlo simulations in **R** and **MATLAB**. We compare their running times with a custom C++ routine, developed for the proposed package.

## 4.1 R and MATLAB implementations

The quantiles of the statistic $T_n$ can be computed using the following straightforward **R** code example:

```R
R> for(i in 1:nsimul)
R> {
R>     x <- rnorm(n)
R>     reg2 <- lm(x ~ X)
R>     Max[i] <- max(abs(rstudent(reg2)))
R> }
R> Qu <- quantile(Max, 0.95)
```

where $X$ is the variable in the model design (such as $X = AGE$ in WAGE AGE, see section 3.2) and nsimul is the number of simulations to perform for estimating the quantiles. Similarly, the following MATLAB code can be used to compute the quantiles of $T_n$:
<figure><img src="image_3.png" /><figcaption>Figure 1: Application of **DEViaN-LM** to biological data in 4 examples. Black dots represent the data points, red dots show abnormal values, and red lines represent the quantiles for the requested α level (for all pictures we took α = 5%). The upper left panel illustrates the detection of abnormal values in the serum iron samples for one elite soccer player (**A**). In the upper right panel (**B**); the abnormal values of distance per week, adjusted by age and bodymass (weight) are given for mice. The lower left panel (**C**) reveals abnormal salaries adjusted by age, educational level and the number of children. The lower right panel (**D**) shows the mouse lemurs which present abnormal grip strength values, when adjusting for age and weight.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Benchmark results. In the left panel (A), the median runtime is given for different sample-sizes and different number of CPU cores. Each measure (node) is repeated 200 times. Similarly, in the right panel (B), the median runtime is given for an increasing number of simulations *nsimul* using 100% of the WAGE dataset. Each measure is repeated 100 times.</figcaption></figure>

MATLAB> for i=1:nsimul
MATLAB> x = randn(s, 1);
MATLAB> mdl = fitlm(X, x);
MATLAB> Max(i) = max(abs(mdl.Residuals.Studentized));
MATLAB> end
MATLAB> Qu = quantile(Max, 0.95);

However, in both cases, the usage of *lm()* or *fitlm()* in the *for()* loop is time consuming when the number of simulations *nsimul* is large. The *lm()* and *fitlm()* functions compute many additional statistics that are not required for computing *T*<sub>*n*</sub>. In contrast, **DEViaN-LM** only computes the (externally) studentized residuals, resulting in lesser operations overall.

## 4.2 Benchmark

The efficiency or performance of the **DEViaN-LM** package in terms of runtime was assessed using a simple benchmark procedure. It was carried out on a desktop computer equipped with 6 INTEL Xeon E-2186G CPU at 3.80GHz, 16GB of RAM (2666MHz) and running **R** 4.3.2 and **MATLAB®** 9.3.0.713579 (R2017b). For the **R** environment, runtime measurements were obtained using the *microbenchmark* package [15], while the *tic-toc* functionality was employed for **MATLAB®**. The performance was compared to a typical **R** and **MATLAB®** code implementation (see previous section 4.1). Two effects were measured: the sample-size effect and the number of Monte-Carlo simulations. In the first case, a dataset with a large number of observations was simulated. In the second case, the number of simulations was increased. In order to have a consistent measure of the time taken to run *nsimul* simulations, we repeated each measure and took the median value of the runtime (in seconds). For the sampling-size effect, we repeated each measure 200 times and we used a simulated linear model $y_i = 25.0 + 3.4x_i + \epsilon_i$ with $x \sim N(0, 1)$, $\epsilon \sim N(0, 2)$ and $i = \{100, 500, 10^3, 5 \times 10^3, 10^4, 10^5, 10^6\}$. For the Monte-Carlo simulations we repeated each measure 100 times and used the WAGE dataset (section 3.2), with *nsimul* taking the following values: $nsimul = \{100, 500, 10^3, 5 \times 10^3, 10^4, 1.5 \times 10^4, 2.0 \times 10^4, 2.5 \times 10^4\}$.

Results are provided in Figure (Fig. 2) and show that the median runtime increases with the sampling-size and the number of simulations *nsimul* (Fig. 2). The increase in the runtime can be estimated using a linear regression of the form $y = rx + b$ where *r* denotes the rate of increase (Table 1). The package **DEViaN-LM** has the smallest *r* in both cases (Table 1).
<table><caption>Table 1: Estimates of growth rate *r* using linear regression for different samples-sizes and number of simulations for 1 CPU. The value of *r* is given in seconds per sample-size (for samples-sizes, Fig. 2, left panel) and in seconds per number of simulations (for different number of simulations, Fig. 2, right panel). The t-stat is the t-statistic for *r* which is the slope of the regression line and *p*-value is the associated *p* value. The RMSE and adjusted R² are given to assess the quality of the linear regression.</caption><thead><tr><th>type</th><th>language</th><th>method</th><th>r</th><th>t-stat</th><th>p-value</th><th>RMSE</th><th>R2a</th></tr></thead><tbody><tr><td>sample-size</td><td>MATLAB</td><td>fitlm()</td><td>1.80×10<sup>-4</sup></td><td>212.91</td><td>0.77</td><td>0.18</td><td>0.99</td></tr><tr><td>sample-size</td><td>R</td><td>lm()</td><td>2.58×10<sup>-4</sup></td><td>86.84</td><td>2.71</td><td>0.048</td><td>0.99</td></tr><tr><td>sample-size</td><td>R</td><td>DEViaN-LM</td><td>7.72×10<sup>-5</sup></td><td>86.45</td><td>0.82</td><td>0.00087</td><td>0.99</td></tr><tr><td>sim. number</td><td>MATLAB</td><td>fitlm()</td><td>4.43×10<sup>-3</sup></td><td>64.11</td><td>9.68×10<sup>-10</sup></td><td>1.75</td><td>0.99</td></tr><tr><td>sim. number</td><td>R</td><td>lm()</td><td>6.06×10<sup>-4</sup></td><td>493.16</td><td>4.69×10<sup>-15</sup></td><td>0.031</td><td>0.99</td></tr><tr><td>sim. number</td><td>R</td><td>DEViaN-LM</td><td>3.99×10<sup>-5</sup></td><td>47.06</td><td>6.17×10<sup>-9</sup></td><td>0.022</td><td>0.99</td></tr></tbody></table>

# 5 Discussion and concluding remarks

One statistical approach when searching for abnormal values or outliers in a linear model relies on the analyse of the studentized residuals. Studentized residuals quantify how large the residuals are in standard deviation units, and therefore can be used to identify outliers: An observation with an externally studentized residual that is larger than 2 (in absolute value) is generally deemed an outlier (see for instance Section 4.1 and paragraph 4.1.2 in [10]). They are typically used to assess the potential outliers that influence a regression model to such an extent that the estimated regression function is attracted towards the potential outlier. The usual way to proceed is to estimate the studentized residuals then look at which ones seem "off the chart" so to speak.

Note that looking at whether each studentized residual is greater than the threshold (calculated from Student's distribution) gives rise to a multiple testing problem. To be convinced of this, suppose that we observe *n* values $x_1, \dots, x_n$ obtained from a sample $X_1, \dots, X_n$ distributed according to the $N(0, 1)$ law. If we say that a value is abnormally large if its absolute value exceeds the threshold of 1.96, then we know that, if $n = 100$, we should observe between 4 and 6 "abnormal" values. These "abnormal" values are in fact perfectly normal, and simply come from the repetition of draws according to the $N(0, 1)$ distribution. We are here exactly in the issue of multiple tests. One way to correct this issue is to identify the observations whose absolute values are above the threshold calculated from the law of $\max_{1 \le i \le n} |X_i|$. These values can really be considered abnormally large.

The method we propose follows exactly this approach, in the context of the Gaussian linear model. It computes the exact distribution of the maximum values of the studentized residuals, based on the design of the model. The estimate of the quantiles can only be computed using numerical simulations as the quantiles depends on each design. The **DEViaN-LM** package automates this process and provides users with a unique and practical tool for measuring abnormal values in the Gaussian Framework. Its implementation is efficient and compares very well with different build-in functions (Fig. 2). We focused on the free **R** software environment, but the proposed **C++** implementation can be easily adapted to other languages.

# 6 Acknowledgments

We thank Brigitte Gelein, Thibaud Lefèvre, Kenza Charifi and Aliya Mohamed Soultane for providing additional help in coding the package.

# References

[1] Victoria Hodge and Jim Austin. A survey of outlier detection methodologies. *Artificial intelligence review*, 22(2):85–126, 2004.

[2] Luca Malcovati, Cristiana Pascutto, and Mario Cazzola. Hematologic passport for athletes competing in endurance sports: a feasibility study. *Haematologica*, 88(5):570–581, 2003.
[3] Florian Egger, Tim Meyer, and Anne Hecksteden. Interindividual variation in the relationship of different intensity markers—a challenge for targeted training prescriptions. *PloS one*, 11(10):e0165010, 2016.

[4] Guillaume Saulière, Jérôme Dedecker, Laurie-Anne Marquet, Pierre Rochcongar, Jean-François Toussaint, and Geoffroy Berthelot. Z-scores-based methods and their application to biological monitoring: an example in professional soccer players. *Biostatistics*, 20(1):48–64, 2019.

[5] National Research Council et al. Toward precision medicine: building a knowledge network for biomedical research and a new taxonomy of disease. 2011.

[6] Chris Cheadle, Marquis P Vawter, William J Freed, and Kevin G Becker. Analysis of microarray data using z score transformation. *The Journal of molecular diagnostics*, 5(2):73–81, 2003.

[7] K Ishii, M Sasaki, M Matsui, S Sakamoto, S Yamaji, N Hayashi, T Mori, H Kitagaki, N Hirono, and E Mori. A diagnostic method for suspected alzheimer's disease using h 2 15 o positron emission tomography perfusion z score. *Neuroradiology*, 42(11):787–794, 2000.

[8] JP Guilford. Structure-of-intellect abilities in preliterate children and the mentally retarded. *Monograph of the American Association of Mental Deficiency*, (1):46–58, 1973.

[9] Frederic Dallaire, Jean-Luc Bigras, Milan Prsa, and Nagib Dahdah. Bias related to body mass index in pediatric echocardiographic z scores. *Pediatric cardiology*, 36(3):667–676, 2015.

[10] Pierre-André Cornillon and Eric Matzner-Løber. Régression avec R. *Pratique R*, 2011.

[11] Ken Sharpe, Michael J Ashenden, and Yorck O Schumacher. A third generation approach to detect erythropoietin abuse in athletes. *Haematologica*, 91(3):356–363, 2006.

[12] P Thomas, E Pouydebat, ML Brazidec, F Aujard, and Anthony Herrel. Determinants of pull strength in captive grey mouse lemurs. *Journal of Zoology*, 298(2):77–81, 2016.

[13] Theodore J Morgan, Theodore Garland Jr, and Patrick A Carter. Ontogenies in mice selected for high voluntary wheel-running activity. i. mean ontogenies. *Evolution*, 57(3):646–657, 2003.

[14] AM Bronikowski, TJ Morgan, T Garland Jr, and PA Carter. The evolution of aging and age-related physical decline in mice selectively bred for high voluntary exercise. *Evolution*, 60(7):1494–1508, 2006.

[15] Olaf Mersmann. *microbenchmark: Accurate Timing Functions*, 2019. R package version 1.4-7.