Mixed Model
Arthur Leroy
Workshop du LaPsyDE - 22 janvier 2021
Introduction
This document aims at providing a first and basic introduction to the implementation of mixed models in R. We use an applicative example coming from the APEX project, studing the impact of a cognitive program. There might be many ways to study this dataset, and uncountable models that could be defined, keep in mind that the following are just illustrative examples. Although various R packages are available to fit mixed models, which are widely studied objects (still an active research topic), the most famous is probably lme4. The syntaxe is pretty straightforward and the implementation is powerfull. Actually, the full presentation of the package is way beyond the scope of the present workshop. You can find the full documentation here: https://cran.r-project.org/web/packages/lme4/lme4.pdf. Besides, here is an additional (and awesome) reference to go (way) further on mixed models with R: https://m-clark.github.io/mixed-models-with-R/.
As for any statistical analysis, let's first load the packages that we will need later on. The meta-package tidyverse provides a really convenient and unified framework to handle data management such as cleaning, tidying, and visualizing data. I must recommend you to use it as often as possible.
## Uncomment if you don't have installed those packages yet
## install.packages(c('tidyverse', 'lme4', 'BEST'))
library(tidyverse)
library(lme4)Then, we may import the dataset and clean the data as needed. For our simple following examples, we only focus on the variable Stroop_Cog_diff that represents (maybe? probably? Nothing's clear, I'm way out of my safety zone here) a score to a cognitive test performed by children. Since, the variable is originally coded as a string of characters, I rename the variable output and transform it back as actual numbers.
raw_db = read_csv2("documents/APEXAdo_Glucose_dataset.csv")
db = raw_db %>% mutate(output = as.numeric(Stroop_Cog_diff))For information let's display a quick teasing of the raw data. Notice, that we will use the variable groupe, which codes whether a subject was in the control or the intervention group for these tests. There are way more variables, but it's often better to start slowly, and avoid introducing too many cofounding varables and possibly many biases.
head(db)## # A tibble: 6 x 17
## code_sujet groupe sexe age age_T0_mois revenus Time_num Time
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 01ZA001 CA F 16 197 2 1_pre pre
## 2 01LJ004 CA F 16 199 6 1_pre pre
## 3 01FL007 CA F 16 206 6 1_pre pre
## 4 01GI009 CA F 16 203 NA 1_pre pre
## 5 01AS011 CA F 16 202 5 1_pre pre
## 6 01ML012 CA F 16 197 4 1_pre pre
## # ... with 9 more variables: heure_aprox_bat_T0 <chr>, heure_1cat_gly <chr>,
## # glycemie_mg_dl_BAT <dbl>, gly2 <dbl>, gly_cat <chr>, Session <chr>,
## # Stroop_Cog_diff <chr>, Stroop_Emo_diff <chr>, output <dbl>Using the lme4 package, defining a mixed model is pretty straightforward. We can merely use the function lmer() in which we precise the structure of the relationship as first argument, and the dataset as second. First, we chose to study the effect of the intervention alone and thus specify the variable groupe (coding whether we are in control or mindful meditation group) as a fixed effect and the suject as a random effect (since there are several observations, probably dependent, for each subject). Below, we would consider that the random effects only impacts the intercept, meaning that we consider each subject as having its own skills in performing thoses tasks, regarless of the intervention. The intercept is represented by the number 1, between the brackets, which them indicate the random effect.
model_simple = lmer(output ~ groupe + (1|code_sujet), data = db)
summary(model_simple)## Linear mixed model fit by REML ['lmerMod']
## Formula: output ~ groupe + (1 | code_sujet)
## Data: db
##
## REML criterion at convergence: 2536.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.72910 -0.60567 -0.07112 0.48903 3.00458
##
## Random effects:
## Groups Name Variance Std.Dev.
## code_sujet (Intercept) 1595 39.93
## Residual 11902 109.10
## Number of obs: 207, groups: code_sujet, 53
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 82.069 12.908 6.358
## groupePC -9.407 18.771 -0.501
##
## Correlation of Fixed Effects:
## (Intr)
## groupePC -0.688The function summary allows us to dig into the results of the estimation of the model. There are many different thing to look at, but the more important are the estimation of the coefficient and the analysis of there variance. The analysis of those results strongly depends on the context and the expert knowledge, thus we will remain careful in this document and let this aspect to the workshop comments and discussion.
If we acknoledge that the intervention should have an effect, we may introduce the variable Time coding the pre/post timing, to define something slightly more sophisticate. We could imagine that not only the individuals have their own natural skills, but they also present their indivudal effect of the intervention (more or less habits with a tablet for example), and thus specify the random effect on Time coefficient as well.
model = lmer(output ~ groupe + (1 + Time|code_sujet), data = db)
summary(model)## Linear mixed model fit by REML ['lmerMod']
## Formula: output ~ groupe + (1 + Time | code_sujet)
## Data: db
##
## REML criterion at convergence: 2528.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3324 -0.5011 -0.0318 0.4109 2.9337
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## code_sujet (Intercept) 5195 72.08
## Timepre 8090 89.94 -0.76
## Residual 9192 95.88
## Number of obs: 207, groups: code_sujet, 53
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 84.66 12.94 6.543
## groupePC -13.26 18.81 -0.705
##
## Correlation of Fixed Effects:
## (Intr)
## groupePC -0.688As an illustration of the potential to define more (possibly too much) complex models, here is an example that seemed relevant to me, among the many other possibilities.
model = lmer(output ~ sexe + age + glycemie_mg_dl_BAT + Time + Session + (1 |groupe:code_sujet), data = db)
summary(model)## Linear mixed model fit by REML ['lmerMod']
## Formula: output ~ sexe + age + glycemie_mg_dl_BAT + Time + Session + (1 |
## groupe:code_sujet)
## Data: db
##
## REML criterion at convergence: 2082
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5626 -0.4915 -0.1031 0.5219 2.9550
##
## Random effects:
## Groups Name Variance Std.Dev.
## groupe:code_sujet (Intercept) 2358 48.56
## Residual 10889 104.35
## Number of obs: 173, groups: groupe:code_sujet, 52
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -140.0846 359.4520 -0.390
## sexeM -4.5971 22.3617 -0.206
## age 1.1767 21.9824 0.054
## glycemie_mg_dl_BAT 2.0233 0.8403 2.408
## Timepre 19.5856 16.8530 1.162
## SessionEnd -1.6245 15.8747 -0.102
##
## Correlation of Fixed Effects:
## (Intr) sexeM age g___BA Timepr
## sexeM 0.086
## age -0.974 -0.101
## glycm___BAT -0.113 -0.026 -0.110
## Timepre -0.015 -0.030 -0.054 0.224
## SessionEnd -0.020 -0.003 -0.004 0.012 -0.003Once again, not being a specialist in cognitive sciences nor even in applying mixed models, I would restrain myself from over-interpreting such results. However, on another (but really important, I believe) topic, I would like to provide an additional package that would deserve to be way more widespread among the experimental science practitioners, desiring to use simple and meaningful statistical analyses. I am indeed devoted to a endless war, against the decision based on a t-test, a method that no statistician would ever use, altough the entire world seems to enjoy it more than anything (probably because it makes the concept of talking nonsense really straightfoward and appearing scientifically sound. (Am I sarcastic enough here ? Probably)).
Bonus: replace (easily) the ugly t-tests with the BEST approach
For this ambitious task, we will use an awesome package developed by the brilliant John K. Kruschke, Professor of Psychological and Brain Sciences (cheers to your field!). Namely, BEST means 'Bayesian Estimation Superseds T-test' (see ? I'm not the only one who is sarcastic), and propose a really user-friendly framework to conduct a Bayesian analysis of the mean comparisons problem. Whether you want to compare the mean of two different groups (like pre-post test for example) or compare the mean a single group to 0, the syntax remains completly obious. We juste have to use the function BESTmcmc and provides as first (and second if necessary) argument, the vector of observation you want to test. But yet, you will not TEST! Stop testing. Testing is nonsense, it's hiding behind false objectivity and forgetting the experts knowledges of each field. And welcome in the world of estimation, probability, straightforward statements, uncertainty quantification, and correct measure of effect sizes. First of all, let's create two vectors of data, for illustration, to see if we could exhibit a difference in cognitive test based on the gender.
library(BEST)
db_h = raw_db %>%
filter(sexe == 'M') %>%
filter(!is.na(Stroop_Cog_diff)) %>%
pull(Stroop_Cog_diff) %>%
as.numeric
db_f = raw_db %>%
filter(sexe == 'F') %>%
filter(!is.na(Stroop_Cog_diff)) %>%
pull(Stroop_Cog_diff) %>%
as.numericWe could visualize on the boxplot that there seems to be nothing really interesting here. Ok, but how much ? Maybe not a slight effect ? Let's see it more deeply !
ggplot(raw_db) + geom_boxplot(aes(x = sexe, y = as.numeric(Stroop_Cog_diff), color = sexe))
Below, we use the wunerfully simple function BESTmcmc function to analyze the difference between those 2 groups.
res_hf = BESTmcmc(db_h, db_f, numSavedSteps = 10000)
res_hf## MCMC fit results for BEST analysis:
## 10002 simulations saved.
## mean sd median HDIlo HDIup Rhat n.eff
## mu1 78.28 15.051 78.17 48.552 107.67 1.000 5605
## mu2 73.52 9.756 73.57 54.765 93.07 1.002 5720
## nu 29.76 25.380 21.63 3.461 81.66 1.005 1420
## sigma1 120.10 12.398 119.88 95.660 144.52 1.001 3566
## sigma2 104.86 8.223 104.96 88.983 121.30 1.005 2840
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## 'Rhat' is the potential scale reduction factor (at convergence, Rhat=1).
## 'n.eff' is a crude measure of effective sample size.On the results, we have really detailled analysis of the different quantities that have been estimated. However, for a one-look estimation of the gobal effect, let's display the graph of results !
plot(res_hf, 'mean')
As you can see. The difference (or the absence of difference here) is now perfectly quantified. No test needed here, just a nice statistical estimation and your brain that is already expert in making decision by itself! On a second example, we can study the effect of the intervention (this is want we want to know afterall isn't it ?)
db_pre = raw_db %>%
filter(groupe == 'PC') %>%
filter(Time == 'pre') %>%
filter(!is.na(Stroop_Cog_diff)) %>%
pull(Stroop_Cog_diff) %>%
as.numeric()
db_post = raw_db %>%
filter(groupe == 'PC') %>%
filter(Time == 'post') %>%
filter(!is.na(Stroop_Cog_diff)) %>%
pull(Stroop_Cog_diff) %>%
as.numeric()
raw_db %>%
filter(groupe == 'PC') %>%
ggplot() + geom_boxplot(aes(x = Time, y = as.numeric(Stroop_Cog_diff), color = Time))
We can then run the same framework.
res_pre_post = BESTmcmc(db_post, db_pre, numSavedSteps = 10000)
res_pre_post## MCMC fit results for BEST analysis:
## 10002 simulations saved.
## mean sd median HDIlo HDIup Rhat n.eff
## mu1 89.40 15.22 89.50 59.361 118.51 1.001 5678
## mu2 56.55 14.22 56.43 29.985 85.89 1.000 5971
## nu 38.93 30.36 30.53 2.683 98.19 1.004 1977
## sigma1 97.98 11.27 97.04 77.029 120.92 1.002 4827
## sigma2 97.25 11.24 96.57 76.193 119.98 1.000 4612
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## 'Rhat' is the potential scale reduction factor (at convergence, Rhat=1).
## 'n.eff' is a crude measure of effective sample size.And observe the nice results exhibiting a highly probable effect. High enough to mention ? It's up to you to tell now. Enjoy !
plot(res_pre_post, 'mean', showCurve = T)
And more results on the difference of standard errors between groups.
plot(res_pre_post, 'sd', showCurve = T)
And even on the effect size of the intervention! Life's beautiful when it's simple and pretty.
plot(res_pre_post, 'effect', showCurve = T)