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)
  )

Figure 1: Dummy dataset with categorical variables.

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, :
unable to evaluate scaled gradient
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)
unable to evaluate scaled gradient
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, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   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)')

Figure 2: Comparison of fixed effects parameters for model A + (1|A)

Figure 3: Comparison of random effects parameters for model A + (1|A)

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, :
unable to evaluate scaled gradient
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)
unable to evaluate scaled gradient
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, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   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, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   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)
  )

Figure 4: Dummy dataset with a discrete variable. The line represents the true model.

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 

Fitted values (line) and observed values (points) from the lme4 model.

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')

Fitted values from the fixed and random part of the lme4 models. Points represent the true model.

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 

Figure 5: Fitted values from the fixed and random part of the INLA models. Points represent the true model.

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
  )

Figure 6: Dummy dataset with a continuous variable. The line represents the true model.

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

──────────────────────────────────────────────────────────────────────────────