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)
)Using a variable both as fixed and random effect
statistics
mixed-models
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.
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 gradientWarning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvaluessummary(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 eigenvalueslibrary(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)')
A + (1|A)
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 gradientWarning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvaluessummary(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 eigenvaluesmodel_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.000model_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)
)
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 ' ' 1summary(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.000summary(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')
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 
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
)
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
──────────────────────────────────────────────────────────────────────────────