Using a variable both as fixed and random effect

statistics
mixed-models
Author

Thierry Onkelinx

Published

August 23, 2017

One of the questions to answer when using mixed models is whether to use a variable as a fixed effect or as a random effect. Sometimes it makes sense to use a variable both as fixed and random effect. In this post I will try to make clear in which cases it can make sense and what are the benefits of doing so. I will also handle cases in which it doesn’t make sense. Much will depend on the nature of the variable. Therefore this post is split into three sections: categorical, discrete and continuous. I will only display the most relevant parts of the code. The full code is available on GitHub.

Categorical variable

To make this clear, we start by creating a dummy dataset with 3 categorical covariates. `B` is nested within `A`. The resulting dataset is displayed in Figure 1.

``````library(tidyverse)
set.seed(20170823)
n_a <- 6
n_b <- 2
n_sample <- 3
sd_a <- 2
sd_b <- 1
sd_noise <- 1
dataset <- expand.grid(
B = paste0("b", seq_len(n_a * n_b)),
sample = seq_len(n_sample)
) %>%
mutate(
A = paste0("a", as.integer(B) %% n_a) %>%
factor(),
mu = rnorm(n_a, sd = sd_a)[A] + rnorm(n_a * n_b, sd = sd_b)[B],
Y = mu + rnorm(n(), sd = sd_noise)
)``````

The first model is one that doesn’t make sense. Using a categorical variable both as random and a fixed effect. In this case both effects are competing for the same information. Below is the resulting fit from `lme4` and `INLA`. Note the warning in the `lme4` output, the model failed to converge. Nevertheless, both `lme4` and `INLA` yield the same parameter estimate (Figure 2), albeit the much wider confidence intervals for `lme4`. The estimates for the random effects in both packages are equivalent to zero (Figure 3). Again the `lme4` estimate has more uncertainty.

``````library(lme4)
model_1 <- lmer(Y ~ 0 + A + (1 | A), data = dataset)``````
``````Warning in checkConv(attr(opt, "derivs"), opt\$par, ctrl = control\$checkConv, :
``````Warning in checkConv(attr(opt, "derivs"), opt\$par, ctrl = control\$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues``````
``summary(model_1)``
``````Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ 0 + A + (1 | A)
Data: dataset

REML criterion at convergence: 109.5

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.5304 -0.7336 -0.0650  0.5862  2.2205

Random effects:
Groups   Name        Variance Std.Dev.
A        (Intercept) 6.651    2.579
Residual             1.576    1.255
Number of obs: 36, groups:  A, 6

Fixed effects:
Estimate Std. Error t value
Aa0  -4.0282     2.6294  -1.532
Aa1   1.1165     2.6294   0.425
Aa2  -1.2266     2.6294  -0.466
Aa3   2.6855     2.6294   1.021
Aa4   0.4843     2.6294   0.184
Aa5  -2.9922     2.6294  -1.138

Correlation of Fixed Effects:
Aa0   Aa1   Aa2   Aa3   Aa4
Aa1 0.000
Aa2 0.000 0.000
Aa3 0.000 0.000 0.000
Aa4 0.000 0.000 0.000 0.000
Aa5 0.000 0.000 0.000 0.000 0.000
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues``````
``````library(INLA)
model_2 <- inla(Y ~ 0 + A + f(A, model = "iid"), data = dataset)
summary(model_2)``````
``````
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.992, Running = 0.838, Post = 0.0399, Total = 1.87
Fixed effects:
mean    sd 0.025quant 0.5quant 0.975quant   mode kld
Aa0 -4.027 0.508     -5.029   -4.027     -3.025 -4.027   0
Aa1  1.116 0.508      0.114    1.116      2.118  1.116   0
Aa2 -1.226 0.508     -2.228   -1.226     -0.224 -1.226   0
Aa3  2.685 0.508      1.683    2.685      3.687  2.685   0
Aa4  0.484 0.508     -0.518    0.484      1.486  0.484   0
Aa5 -2.991 0.508     -3.993   -2.991     -1.989 -2.991   0

Random effects:
Name    Model
A IID model

Model hyperparameters:
mean       sd 0.025quant 0.5quant
Precision for the Gaussian observations 6.82e-01 1.71e-01      0.399 6.64e-01
Precision for A                         2.17e+04 2.31e+04   1582.690 1.45e+04
0.975quant     mode
Precision for the Gaussian observations       1.07    0.631
Precision for A                           83100.47 4388.898

Marginal log-Likelihood:  -91.89
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')``````

What if we want to add variable `B` as a nested random effect? We already know that adding `A` to both the fixed and the random effects is nonsense. The correct way of doing this is to use `A` as a fixed effect and `B` as an implicit nested random effect.

``model_1 <- lmer(Y ~ 0 + A + (1 | A / B), data = dataset)``
``````Warning in checkConv(attr(opt, "derivs"), opt\$par, ctrl = control\$checkConv, :
``````Warning in checkConv(attr(opt, "derivs"), opt\$par, ctrl = control\$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues``````
``summary(model_1)``
``````Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ 0 + A + (1 | A/B)
Data: dataset

REML criterion at convergence: 105.6

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.4004 -0.4931 -0.1464  0.4092  2.1372

Random effects:
Groups   Name        Variance Std.Dev.
B:A      (Intercept) 0.7957   0.892
A        (Intercept) 1.8676   1.367
Residual             1.0983   1.048
Number of obs: 36, groups:  B:A, 12; A, 6

Fixed effects:
Estimate Std. Error t value
Aa0  -4.0282     1.5648  -2.574
Aa1   1.1165     1.5648   0.714
Aa2  -1.2266     1.5648  -0.784
Aa3   2.6855     1.5648   1.716
Aa4   0.4843     1.5648   0.310
Aa5  -2.9922     1.5648  -1.912

Correlation of Fixed Effects:
Aa0   Aa1   Aa2   Aa3   Aa4
Aa1 0.000
Aa2 0.000 0.000
Aa3 0.000 0.000 0.000
Aa4 0.000 0.000 0.000 0.000
Aa5 0.000 0.000 0.000 0.000 0.000
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues``````
``````model_1b <- lmer(Y ~ 0 + A + (1 | B), data = dataset)
summary(model_1b)``````
``````Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ 0 + A + (1 | B)
Data: dataset

REML criterion at convergence: 105.6

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.4004 -0.4931 -0.1464  0.4092  2.1372

Random effects:
Groups   Name        Variance Std.Dev.
B        (Intercept) 0.7957   0.892
Residual             1.0983   1.048
Number of obs: 36, groups:  B, 12

Fixed effects:
Estimate Std. Error t value
Aa0  -4.0282     0.7622  -5.285
Aa1   1.1165     0.7622   1.465
Aa2  -1.2266     0.7622  -1.609
Aa3   2.6855     0.7622   3.523
Aa4   0.4843     0.7622   0.635
Aa5  -2.9922     0.7622  -3.926

Correlation of Fixed Effects:
Aa0   Aa1   Aa2   Aa3   Aa4
Aa1 0.000
Aa2 0.000 0.000
Aa3 0.000 0.000 0.000
Aa4 0.000 0.000 0.000 0.000
Aa5 0.000 0.000 0.000 0.000 0.000``````
``````model_2 <- inla(
Y ~ 0 + A + f(A, model = "iid") + f(B, model = "iid"), data = dataset
)
summary(model_2)``````
``````
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.631, Running = 0.864, Post = 0.0183, Total = 1.51
Fixed effects:
mean    sd 0.025quant 0.5quant 0.975quant   mode kld
Aa0 -4.027 0.586     -5.183   -4.027     -2.870 -4.027   0
Aa1  1.116 0.586     -0.040    1.116      2.273  1.116   0
Aa2 -1.226 0.586     -2.383   -1.226     -0.070 -1.226   0
Aa3  2.685 0.586      1.528    2.685      3.841  2.685   0
Aa4  0.484 0.586     -0.672    0.484      1.641  0.484   0
Aa5 -2.991 0.586     -4.147   -2.991     -1.835 -2.991   0

Random effects:
Name    Model
A IID model
B IID model

Model hyperparameters:
mean       sd 0.025quant 0.5quant
Precision for the Gaussian observations 8.52e-01 2.58e-01      0.414 8.25e-01
Precision for A                         2.17e+04 2.31e+04   1581.172 1.46e+04
Precision for B                         1.81e+01 5.96e+01      0.649 5.80e+00
0.975quant     mode
Precision for the Gaussian observations       1.41    0.789
Precision for A                           83018.44 4384.322
Precision for B                             116.79    1.484

Marginal log-Likelihood:  -97.31
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')``````
``````model_2b <- inla(Y ~ 0 + A + f(B, model = "iid"), data = dataset)
summary(model_2b)``````
``````
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.52, Running = 0.863, Post = 0.016, Total = 1.4
Fixed effects:
mean    sd 0.025quant 0.5quant 0.975quant   mode kld
Aa0 -4.027 0.625     -5.261   -4.027     -2.792 -4.027   0
Aa1  1.116 0.625     -0.119    1.116      2.351  1.116   0
Aa2 -1.226 0.625     -2.461   -1.226      0.009 -1.226   0
Aa3  2.684 0.625      1.450    2.684      3.919  2.684   0
Aa4  0.484 0.625     -0.751    0.484      1.719  0.484   0
Aa5 -2.991 0.625     -4.226   -2.991     -1.756 -2.991   0

Random effects:
Name    Model
B IID model

Model hyperparameters:
mean   sd 0.025quant 0.5quant
Precision for the Gaussian observations 453.14 0.00     216.14   453.14
Precision for B                           0.00 0.00       0.00     0.00
0.975quant   mode
Precision for the Gaussian observations    7261.85 453.14
Precision for B                               0.00   0.00

Marginal log-Likelihood:  -94.33
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')``````

Discrete variable

Intro

A discrete variable is a numerical variable but each interval between two values is an integer multiple of a fixed step size. Typical examples are related to time, e.g. the year in steps of 1 year, months expressed in terms of years (step size 1/12), …

We create a new dummy dataset with a discrete variable. The response variable is a third order polynomial of the discrete variable. The `X` variable is rescaled to -1 and 1.

``````n_x <- 25
n_sample <- 10
sd_noise <- 10
dataset <- expand.grid(
X = seq_len(n_x),
sample = seq_len(n_sample)
) %>%
mutate(
mu =  0.045 * X ^ 3 - X ^ 2 + 10,
Y = mu + rnorm(n(), sd = sd_noise),
X = (X - ceiling(n_x / 2)) / floor(n_x / 2)
)``````

Fit with `lme4`

Suppose we fit a simple linear model to the data. We know that this is not accurate because the real pattern is a third order polynomial. And let’s add the variable also as a random effect. We use first `lme4` to illustrate the principle. (@fit-discrete-fit) illustrate how the fit of the fixed part is poor but the random effect of X compensates the fit.

``````model_1 <- lmer(Y ~ X + (1 | X), data = dataset)
summary(model_1)``````
``````Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ X + (1 | X)
Data: dataset

REML criterion at convergence: 1943.1

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.24578 -0.67660 -0.02946  0.60728  3.00837

Random effects:
Groups   Name        Variance Std.Dev.
X        (Intercept) 1442.26  37.977
Residual               88.59   9.412
Number of obs: 250, groups:  X, 25

Fixed effects:
Estimate Std. Error t value
(Intercept)  -20.886      7.619  -2.741
X             11.388     12.678   0.898

Correlation of Fixed Effects:
(Intr)
X 0.000 ``````

The overall model fit improves when we add a second and third polynomial term. And the variance of the random effect decreases. It reduces even to zero once the third polynomial is in the model. (@fit-discrete-fit2) illustrates how the fit of the fixed effect improves when adding the higher order terms. The effect on the fitted values with the random effect is marginal.

``````model_1b <- lmer(Y ~ X + I(X ^ 2) + (1 | X), data = dataset)
model_1c <- lmer(Y ~ X + I(X ^ 2) + I(X ^ 3) + (1 | X), data = dataset)
anova(model_1, model_1b, model_1c)``````
``````Data: dataset
Models:
model_1: Y ~ X + (1 | X)
model_1b: Y ~ X + I(X^2) + (1 | X)
model_1c: Y ~ X + I(X^2) + I(X^3) + (1 | X)
npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
model_1     4 1963.9 1978.0 -977.93   1955.9
model_1b    5 1914.3 1931.9 -952.15   1904.3 51.555  1  6.961e-13 ***
model_1c    6 1836.9 1858.0 -912.46   1824.9 79.387  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1``````
``summary(model_1b)``
``````Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ X + I(X^2) + (1 | X)
Data: dataset

REML criterion at convergence: 1889.6

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.28558 -0.63406 -0.03279  0.61428  3.04912

Random effects:
Groups   Name        Variance Std.Dev.
X        (Intercept) 184.07   13.567
Residual              88.59    9.412
Number of obs: 250, groups:  X, 25

Fixed effects:
Estimate Std. Error t value
(Intercept)  -59.143      4.173 -14.174
X             11.388      4.623   2.463
I(X^2)       105.943      8.622  12.288

Correlation of Fixed Effects:
(Intr) X
X       0.000
I(X^2) -0.746  0.000``````
``summary(model_1c)``
``````Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ X + I(X^2) + I(X^3) + (1 | X)
Data: dataset

REML criterion at convergence: 1814.9

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.4324 -0.6651  0.0087  0.6556  3.4151

Random effects:
Groups   Name        Variance Std.Dev.
X        (Intercept)  0.00    0.000
Residual             88.05    9.384
Number of obs: 250, groups:  X, 25

Fixed effects:
Estimate Std. Error t value
(Intercept) -59.1432     0.8914  -66.35
X           -37.6013     2.4830  -15.14
I(X^2)      105.9426     1.8419   57.52
I(X^3)       75.5290     3.5123   21.50

Correlation of Fixed Effects:
(Intr) X      I(X^2)
X       0.000
I(X^2) -0.746  0.000
I(X^3)  0.000 -0.917  0.000
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')``````

Fit with `INLA`

`INLA` requires that we alter the data to get the same output. First we copy `X` into `X.copy` because `inla` allows the variable to be used only once in the formula. For some reason this wasn’t needed with the categorical variables. The `lme4` syntax `X + (1|X)` translates into the following `INLA` syntax: `X + f(X.copy, model = "iid")`. Then next thing is that `INLA` does the model fitting and prediction in a single step. Getting predictions for new data requires to add the new data to the original data while setting the response to `NA`. If we want predictions for the fixed effect only, then we need to add rows were all random effect covariates are set to `NA`. Hence `X.copy` must be `NA` while `X` must be non `NA`. Note that this would be impossible without creating `X.copy`.

Let’s fit the three same models with `INLA`. The predictions are given in Figure 5. The results are very similar to the `lme4` results.

``````dataset2 <- dataset %>%
mutate(X.copy = X) %>%
bind_rows(
dataset %>%
distinct(X, mu)
)``````
``````model_2 <- inla(
Y ~ X + f(X.copy, model = "iid"), data = dataset2,
control.compute = list(waic = TRUE)
)
model_2b <- inla(
Y ~ X + I(X ^ 2) + f(X.copy, model = "iid"), data = dataset2,
control.compute = list(waic = TRUE)
)
model_2c <- inla(
Y ~ X + I(X ^ 2) + I(X ^ 3) + f(X.copy, model = "iid"), data = dataset2,
control.compute = list(waic = TRUE)
)``````
``````
*** inla.core.safe:  rerun to try to solve negative eigenvalue(s) in the Hessian ``````

Continuous variable

A continuous variable is a numeric variable where there is not fixed step size between two values. In practice the step size will be several magnitudes smaller than the measured values. Again let’s clarify this with an example dataset. For the sake of simplicity we’ll reuse the true model from the example with the discrete variable. Compare Figure 6 with Figure 4 and you’ll see that Figure 6 has no step size while Figure 4 does.

``````n_x <- 25
n_sample <- 10
sd_noise <- 10
dataset <- data.frame(
X = runif(n_x * n_sample, min = 1, max = n_x)
) %>%
mutate(
mu =  0.045 * X ^ 3 - X ^ 2 + 10,
Y = mu + rnorm(n(), sd = sd_noise),
X = (X - ceiling(n_x / 2)) / floor(n_x / 2),
X.copy = X
)``````

The `lmer` model fails because the random effect has as many unique values as observations.

``````tryCatch(
lmer(Y ~ X + (1 | X), data = dataset), error = identity
)``````
``<simpleError: number of levels of each grouping factor must be < number of observations (problems: X)>``
``lmer(Y ~ X + (1 | X), data = dataset)``
``Error: number of levels of each grouping factor must be < number of observations (problems: X)``

The `INLA` model yields output but the variance of the random effect is very high. A good indicator that there is something wrong.

``````model_2 <- inla(
Y ~ X + f(X.copy, model = "iid"), data = dataset,
control.compute = list(waic = TRUE)
)
inla.tmarginal(
function(x) {
x ^ -1
},
model_2\$marginals.hyperpar\$`Precision for X.copy`
) %>%
inla.zmarginal()``````
``````Mean            0.00011261
Stdev           0.000163188
Quantile  0.025 8.88423e-06
Quantile  0.25  3.16364e-05
Quantile  0.5   6.21869e-05
Quantile  0.75  0.000127047
Quantile  0.975 0.000518334 ``````

Conclusion

Using a variable both in the fixed and random part of the model makes only sense in case of a discrete variable.

Session info

These R packages were used to create this post.

``````─ Session info ───────────────────────────────────────────────────────────────
setting  value
version  R version 4.3.1 (2023-06-16)
os       Ubuntu 22.04.3 LTS
system   x86_64, linux-gnu
ui       X11
language nl_BE:nl
collate  nl_BE.UTF-8
ctype    nl_BE.UTF-8
tz       Europe/Brussels
date     2023-08-30
pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
package      * version    date (UTC) lib source
boot           1.3-28.1   2022-11-22 [1] CRAN (R 4.3.0)
class          7.3-22     2023-05-03 [4] CRAN (R 4.3.1)
classInt       0.4-9      2023-02-28 [1] CRAN (R 4.3.1)
cli            3.6.1      2023-03-23 [1] CRAN (R 4.3.0)
colorspace     2.1-0      2023-01-23 [1] CRAN (R 4.3.0)
DBI            1.1.3      2022-06-18 [1] CRAN (R 4.3.0)
digest         0.6.32     2023-06-26 [1] CRAN (R 4.3.1)
dplyr        * 1.1.2      2023-04-20 [1] CRAN (R 4.3.0)
e1071          1.7-13     2023-02-01 [1] CRAN (R 4.3.1)
evaluate       0.21       2023-05-05 [1] CRAN (R 4.3.0)
fansi          1.0.4      2023-01-22 [1] CRAN (R 4.3.0)
farver         2.1.1      2022-07-06 [1] CRAN (R 4.3.0)
fastmap        1.1.1      2023-02-24 [1] CRAN (R 4.3.0)
forcats      * 1.0.0      2023-01-29 [1] CRAN (R 4.3.0)
generics       0.1.3      2022-07-05 [1] CRAN (R 4.3.0)
ggplot2      * 3.4.2      2023-04-03 [1] CRAN (R 4.3.0)
glue           1.6.2      2022-02-24 [1] CRAN (R 4.3.0)
gtable         0.3.3      2023-03-21 [1] CRAN (R 4.3.0)
hms            1.1.3      2023-03-21 [1] CRAN (R 4.3.0)
htmltools      0.5.5      2023-03-23 [1] CRAN (R 4.3.0)
htmlwidgets    1.6.2      2023-03-17 [1] CRAN (R 4.3.0)
INLA         * 23.06.29   2023-06-30 [1] local
inlabru        2.8.0      2023-06-20 [1] CRAN (R 4.3.1)
jsonlite       1.8.7      2023-06-29 [1] CRAN (R 4.3.1)
KernSmooth     2.23-21    2023-05-03 [1] CRAN (R 4.3.0)
knitr          1.43       2023-05-25 [1] CRAN (R 4.3.0)
labeling       0.4.2      2020-10-20 [1] CRAN (R 4.3.0)
lattice        0.21-8     2023-04-05 [4] CRAN (R 4.3.0)
lifecycle      1.0.3      2022-10-07 [1] CRAN (R 4.3.0)
lme4         * 1.1-34     2023-07-04 [1] CRAN (R 4.3.1)
lubridate    * 1.9.2.9000 2023-05-15 [1] https://inbo.r-universe.dev (R 4.3.0)
magrittr       2.0.3      2022-03-30 [1] CRAN (R 4.3.0)
MASS           7.3-60     2023-05-04 [4] CRAN (R 4.3.1)
Matrix       * 1.5-4.1    2023-05-18 [1] CRAN (R 4.3.0)
MatrixModels   0.5-1      2022-09-11 [1] CRAN (R 4.3.0)
minqa          1.2.5      2022-10-19 [1] CRAN (R 4.3.0)
munsell        0.5.0      2018-06-12 [1] CRAN (R 4.3.0)
nlme           3.1-162    2023-01-31 [1] CRAN (R 4.3.0)
nloptr         2.0.3      2022-05-26 [1] CRAN (R 4.3.0)
numDeriv       2016.8-1.1 2019-06-06 [1] CRAN (R 4.3.0)
pillar         1.9.0      2023-03-22 [1] CRAN (R 4.3.0)
pkgconfig      2.0.3      2019-09-22 [1] CRAN (R 4.3.0)
proxy          0.4-27     2022-06-09 [1] CRAN (R 4.3.1)
purrr        * 1.0.1      2023-01-10 [1] CRAN (R 4.3.0)
R6             2.5.1      2021-08-19 [1] CRAN (R 4.3.0)
Rcpp           1.0.10     2023-01-22 [1] CRAN (R 4.3.0)
readr        * 2.1.4      2023-02-10 [1] CRAN (R 4.3.0)
rlang          1.1.1      2023-04-28 [1] CRAN (R 4.3.0)
rmarkdown      2.23       2023-07-01 [1] CRAN (R 4.3.1)
rstudioapi     0.14       2022-08-22 [1] CRAN (R 4.3.0)
scales         1.2.1      2022-08-20 [1] CRAN (R 4.3.0)
sessioninfo    1.2.2      2021-12-06 [1] CRAN (R 4.3.0)
sf             1.0-13     2023-05-24 [1] CRAN (R 4.3.0)
sp           * 2.0-0      2023-06-22 [1] CRAN (R 4.3.1)
stringi        1.7.12     2023-01-11 [1] CRAN (R 4.3.0)
stringr      * 1.5.0      2022-12-02 [1] CRAN (R 4.3.0)
tibble       * 3.2.1      2023-03-20 [1] CRAN (R 4.3.0)
tidyr        * 1.3.0      2023-01-24 [1] CRAN (R 4.3.0)
tidyselect     1.2.0      2022-10-10 [1] CRAN (R 4.3.0)
tidyverse    * 2.0.0      2023-02-22 [1] CRAN (R 4.3.0)
timechange     0.2.0      2023-01-11 [1] CRAN (R 4.3.0)
tzdb           0.4.0      2023-05-12 [1] CRAN (R 4.3.0)
units          0.8-2      2023-04-27 [1] CRAN (R 4.3.0)
utf8           1.2.3      2023-01-31 [1] CRAN (R 4.3.0)
vctrs          0.6.3      2023-06-14 [1] CRAN (R 4.3.0)
withr          2.5.0      2022-03-03 [1] CRAN (R 4.3.0)
xfun           0.39       2023-04-20 [1] CRAN (R 4.3.0)
yaml           2.3.7      2023-01-23 [1] CRAN (R 4.3.0)

[1] /home/thierry/R/x86_64-pc-linux-gnu-library/4.3
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────``````