Statistiques au LAPEC – Analyse de données répétées
Marc Manceau
2024-01
Pour les utilisateurs de R, on commence par charger un certain nombre de modules intéressants.
# le package moderne standard pour manipuler ses données
library(tidyverse)
# pour faire des graphes multiples
library(cowplot)
# pour afficher joliment les tableaux dans un document rmarkdown
library(knitr)
# pour les GLMM
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack# pour les LMM avec tests associés
library(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step# quelques couleurs manuelles
bleufonce <- "#3d5468"
bleuclair <- "#5b7c98"
rose <- "#ff5555"On charge ensuite les jeux de données sur lesquels nous ferons des tests. Tout d’abord le jeu de données d’essai clinique sur patients diabétiques auquel nous sommes habitués :
d <- read.csv("etude_fictive_diabetes.csv")
head(d)## X id group sex age height weight_t0 diabetes steps_t0 t_compliance
## 1 1 1 control F 72 1.61 80.25 II 2642 32
## 2 2 2 control M 71 1.81 79.66 I 1675 48
## 3 3 3 control F 32 1.72 67.25 II 2698 8
## 4 4 4 control M 76 1.74 73.59 I 2519 2
## 5 5 5 control F 34 1.53 83.98 I 2764 12
## 6 6 6 control F 70 1.71 69.14 II 2606 24
## age_group bmi bmi_group diff_weight diff_steps is_depressed
## 1 plus de 60 30.95945 obese -3.6168091 -26.66633 1
## 2 plus de 60 24.31550 normal -2.6491665 656.62903 1
## 3 moins de 40 22.73188 normal -1.5288284 17.31293 0
## 4 plus de 60 24.30638 normal -3.0515988 -589.42996 1
## 5 moins de 40 35.87509 obese 0.8933788 761.91017 1
## 6 plus de 60 23.64488 normal -0.3074968 216.47307 1
## feels_sleepy n_meals_glycemia_out_of_target weight_tf steps_tf
## 1 morning 53 76.63319 2615.334
## 2 morning 19 77.01083 2331.629
## 3 morning 4 65.72117 2715.313
## 4 morning 17 70.53840 1929.570
## 5 afternoon 188 84.87338 3525.910
## 6 always 9 68.83250 2822.473Ainsi qu’un second jeu de données qui se prêtera bien à l’analyse de données répétées : une simulation d’expérience sur des souris où chaque souris est soumise à différentes conditions expérimentales.
d2 <- read.csv("etude_fictive_souris.csv")
head(d2)## X id diet exercise weight lineage sex crp fc attitude
## 1 1 1 low sugar I 16 KO M 1054 616.0048 sociable
## 2 2 1 low sugar II 16 KO M 1054 609.8807 sociable
## 3 3 1 low sugar III 16 KO M 1054 572.6794 prostrate
## 4 4 1 low sugar IV 16 KO M 1054 477.5350 sociable
## 5 5 1 high sugar I 16 KO M 1054 644.7328 sociable
## 6 6 1 high sugar II 16 KO M 1054 502.1937 sociable
## n_extrasystole
## 1 68
## 2 51
## 3 68
## 4 65
## 5 149
## 6 187# on transforme le type de la variable "attitude" et "id" en facteurs
d2 <- d2 %>%
mutate(attitude = as.factor(attitude),
id = as.factor(id))Nous sommes prêts à travailler sur ces jeu de données, et, pour ce qui nous intéresse ici, à analyser les données répétées qu’ils renferment.
Tests appariés
Je rajoute une section à propos de cette option d’appariement, qui bien souvent embrouille plus qu’elle n’aide.
On possède des observations appariées, typiquement la même quantité mesurée “avant” \((X_i)\) et “après” \((Y_i)\) intervention sur chaque individu. Par exemple, le poids du patient avant intervention et après intervention. Y a t’il un changement de poids du patient avant et après intervention ?
Cela revient à tester : * H0 : la moyenne des \((Y_i - X_i)\) vaut zero. * H1 : la moyenne des \((Y_i - X_i)\) est différente de zéro.
T test
La façon de faire ça en R est la suivante :
t.test(x=d$weight_tf, y=d$weight_t0, paired=TRUE)##
## Paired t-test
##
## data: d$weight_tf and d$weight_t0
## t = -11.129, df = 499, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -1.509690 -1.056617
## sample estimates:
## mean difference
## -1.283153On vérifie que c’est équivalent au premier test qu’on avait appliqué directement sur la différence de poids :
t.test(d$diff_weight)##
## One Sample t-test
##
## data: d$diff_weight
## t = -11.129, df = 499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.509690 -1.056617
## sample estimates:
## mean of x
## -1.283153Test de Wilcoxon
La même chose s’applique au test de Wilcoxon, la “version non-paramétrique” du T test :
wilcox.test(x=d$weight_tf, y=d$weight_t0, paired=TRUE)##
## Wilcoxon signed rank test with continuity correction
##
## data: d$weight_tf and d$weight_t0
## V = 30482, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0On vérifie que c’est équivalent au test qu’on avait appliqué directement sur la différence de poids :
wilcox.test(d$diff_weight)##
## Wilcoxon signed rank test with continuity correction
##
## data: d$diff_weight
## V = 30482, p-value < 2.2e-16
## alternative hypothesis: true location is not equal to 0On retiendra donc que l’option “test apparié” s’applique à des données mesurées 2 fois sur le même individu. Raisonner directement sur la différence entre les deux mesures permet de s’affranchir de cette option d’appariement.
Dans la suite, on s’intéressera à des données répétées avec plus que deux mesures par individu. On sait que ces mesures sont probablement corrélées entre elles, comment le prendre en compte ?
LMM
Principe
Le LMM signifie Linear Mixed Model. Par rapport à un “simple” LM (Linear Model), on rajoute donc le M de Mixed effect.
Cet effet aléatoire sert à modéliser des effets auxquels on ne s’intéresse pas particulièrement en soi, mais dont on sait qu’ils peuvent corréler artificiellement nos mesures.
Le cas de figure typique est celui de mesures répétées à l’échelle de l’individu, potentiellement dans des conditions différentes. Si \(Y_{ij}\) est la mesure numéro \(j\) au sein du même individu \(i\), associée à un \(n\) variables explicatives \((X_{ij}^{(k)})\), alors on propose le modèle suivant :
\[ Y_{ij} = a_0 + a_1 X_{ij}^{(1)} + a_2 X_{ij}^{(2)} + ... + a_n X_{ij}^{(n)} + \epsilon_i + \epsilon_{ij}\\ \epsilon_i \sim \mathcal{N}(0, \sigma^2)\\ \epsilon_{ij} \sim \mathcal{N}(0, \tau^2) \]
Les \((\epsilon_i)\) sont des variables aléatoires indépendantes et identiquement distribuées entre les individus, qui permettent de “corréler” entre elles les observations d’un même individu. Les \((\epsilon_{ij})\) représentent le bruit associé à chaque mesure, comme dans un LM classique.
Application
Dans notre exemple, chaque souris \(i \in \{1, 2, ..., 50\}\) voit sa fréquence cardiaque fc mesurée après un exercice \(j \in \{\text{I, II, III, IV}\}\) et un régime alimentaire \(k \in \{\text{low sugar, high sugar}\}\).
On peut chercher à explorer les données avec quelques graphes, pour se faire une idée de l’importance de la corrélation des mesures au sein de chaque individu par rapport à l’effet des changements de régime :
d2 %>% group_by(id, diet) %>%
summarize(fc = mean(fc)) %>%
ggplot(aes(x=id, y=fc, shape=diet)) +
geom_point() +
xlab("numéro de la souris") +
ylab("fréquence cardiaque") +
theme_bw()
Et la même chose par rapport aux changements d’exercice :
p_illu_lmm <- d2 %>% group_by(id, exercise) %>%
summarize(fc = mean(fc)) %>%
ggplot(aes(x=id, y=fc, col=exercise)) +
geom_point() +
xlab("numéro de la souris") +
ylab("fréquence cardiaque") +
theme_bw()
p_illu_lmm
On note \(Y_{ijk}\) la mesure de la souris \(i\) sous exercice \(j\) et régime alimentaire \(k\), et on propose le modèle suivant :
\[ Y_{ijk} = a + \sum_{l \in \{\text{II, III, IV}\}} b_l \delta_{j==l} + c_m \delta_{k==\text{"low sugar"}} + \epsilon_i + \epsilon_{ijk} \]
Dans ce modèle, on a
- un effet fixe associé à l’ensemble des mesures,
- un effet fixe associé à chaque catégorie d’exercice,
- un effet fixe additif associé à chaque catérogie de régime,
- un effet aléatoire associé à la souris chez qui la mesure est effectuée,
- un petit bruit aléatoire supplémentaire sur chaque mesure.
En R, on fitte ce modèle avec le code suivant :
res_lmm_1 <- lmer("fc ~ diet + exercise + (1 | id)", data=d2)
summary(res_lmm_1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: "fc ~ diet + exercise + (1 | id)"
## Data: d2
##
## REML criterion at convergence: 4350
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.07098 -0.61984 0.05554 0.68635 2.01146
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2155 46.43
## Residual 2603 51.02
## Number of obs: 400, groups: id, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 669.881 8.697 109.623 77.021 < 2e-16 ***
## dietlow sugar -46.537 5.102 346.000 -9.122 < 2e-16 ***
## exerciseII -3.020 7.215 346.000 -0.419 0.676
## exerciseIII -3.619 7.215 346.000 -0.502 0.616
## exerciseIV -54.998 7.215 346.000 -7.623 2.41e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dtlwsg exrcII exrIII
## dietlowsugr -0.293
## exerciseII -0.415 0.000
## exerciseIII -0.415 0.000 0.500
## exerciseIV -0.415 0.000 0.500 0.500Que concluons-nous ?
Quel modèle plus complexe pourrait-on proposer pour inclure un effet fixe de la lignée de souris et du sexe ?
GLMM
Principe
De même que le LMM rajoutait un effet aléatoire à un LM, le GLMM propose de rajouter un effet aléatoire à un GLM. GLMM signifie donc : Generalized Linear Mixed Model.
Le principe intuitif est en tous points similaire au principe du LMM : on souhaite prendre en compte le fait que les variables mesurées au sein d’un même individu soient corrélées entre elles, mais on ne s’intéresse pas, en soi, à cette corrélation.
On suppose donc que cet effet est bien modélisé par une variable aléatoire indépendante des autres et identiquement distribuée chez tous les individus. Les mesures \(Y_i\) chez un individu \(i\), avec variables explicatives \(X\), ont pour espérance :
\[ \mathbb{E}( Y_i | X ) = g^{-1}( a_0 + a_1 x_1 + a_2 x_2 + ... + a_n x_n + \epsilon_i )\\ \epsilon_i \sim \mathcal{N}(0, \sigma^2) \]
Régression logistique mixte
Dans l’étude sur les souris, on dispose d’une variable binaire attitude, qui est mesurée plusieurs fois chez chaque souris, dans des conditions expérimentales différentes.
On tente de se faire une idée de l’auto-corrélation des données mesurées au sein d’un individu avec un graphe. C’est bien moins parlant que précédemment, puisqu’il n’y a que deux outcomes possibles :
p_glmm_1 <- d2 %>%
ggplot(aes(x=id, y=attitude, shape=diet, col=exercise)) +
geom_point() +
geom_jitter(width=0, height=0.1) +
xlab("numéro de la souris") +
ylab("") +
theme_bw()
p_glmm_1
On imagine qu’une certaine part de la variabilité des données provient du caractère “de base” de la souris, i.e. nos données mesurées au sein d’un individu sont probablement corrélées. On décide de prendre en compte cet effet en rajoutant un effet aléatoire à une régression logistique.
Avec R, la syntaxe pour fitter ce modèle est relativement similaire à ce qu’on a vu précédemment :
res_glmm_1 <- glmer("attitude ~ diet + exercise + crp + weight + (1 | id)",
family = binomial,
data=d2)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0258435 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?summary(res_glmm_1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: attitude ~ diet + exercise + crp + weight + (1 | id)
## Data: d2
##
## AIC BIC logLik deviance df.resid
## 469.8 501.7 -226.9 453.8 392
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6463 -0.6790 0.3268 0.6296 1.9887
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.8307 0.9114
## Number of obs: 400, groups: id, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9982016 1.1113114 -2.698 0.00698 **
## dietlow sugar -1.1798303 0.2496913 -4.725 2.30e-06 ***
## exerciseII -0.4058869 0.3419182 -1.187 0.23519
## exerciseIII -0.2343162 0.3427029 -0.684 0.49415
## exerciseIV -0.4059024 0.3419181 -1.187 0.23518
## crp 0.0019176 0.0004107 4.669 3.02e-06 ***
## weight 0.0265621 0.0417641 0.636 0.52477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dtlwsg exrcII exrIII exrcIV crp
## dietlowsugr -0.015
## exerciseII -0.134 0.037
## exerciseIII -0.143 0.021 0.511
## exerciseIV -0.134 0.037 0.514 0.511
## crp -0.597 -0.141 -0.036 -0.020 -0.036
## weight -0.665 -0.020 -0.005 -0.003 -0.005 -0.132
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0258435 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?Les warnings qui s’affichent ici sont typiques de ce qui vous arrivera dans la vraie vie : ces modèles, bien que très classiques et interprétables, commencent à être compliqués à fitter d’un point de vue numérique. Dans ce cas précis, le warning nous donne une idée de méthode pour que l’algorithme converge mieux : re-scaler des données numériques un peu grandes.
d2bis <- d2 %>%
mutate(crp = (crp - mean(crp))/sd(crp),
weight = (weight - mean(weight))/sd(weight))
res_glmm_2 <- glmer("attitude ~ diet + exercise + crp + weight + (1 | id)",
family = binomial,
data=d2bis)
summary(res_glmm_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: attitude ~ diet + exercise + crp + weight + (1 | id)
## Data: d2bis
##
## AIC BIC logLik deviance df.resid
## 469.8 501.7 -226.9 453.8 392
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6463 -0.6790 0.3268 0.6296 1.9887
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.8306 0.9114
## Number of obs: 400, groups: id, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4057 0.3186 4.413 1.02e-05 ***
## dietlow sugar -1.1798 0.2498 -4.723 2.32e-06 ***
## exerciseII -0.4059 0.3419 -1.187 0.235
## exerciseIII -0.2343 0.3427 -0.684 0.494
## exerciseIV -0.4058 0.3419 -1.187 0.235
## crp 0.9538 0.2055 4.641 3.47e-06 ***
## weight 0.1136 0.1786 0.636 0.525
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dtlwsg exrcII exrIII exrcIV crp
## dietlowsugr -0.471
## exerciseII -0.573 0.037
## exerciseIII -0.560 0.021 0.511
## exerciseIV -0.573 0.037 0.514 0.511
## crp 0.161 -0.144 -0.036 -0.021 -0.036
## weight 0.015 -0.019 -0.005 -0.003 -0.005 -0.133Que concluons-nous ?
Quel autre modèle auriez-vous envie d’envisager ?
Régression de Poisson mixte
Dernier exemple de GLMM : la régression de Poisson avec effets aléatoires !
Les données auxquelles on s’intéresse sont n_extrasystole. Ce sont des données de comptage, qu’on aimerait modéliser comme suivant une loi de Poisson de paramètre \(\lambda\), où ce paramètre \(\lambda\) intègre un effet aléatoire correspondant à la souris chez laquelle la mesure est effectuée.
De cette façon, si certaines souris, de base, sont plus susceptibles de faire des extra-systoles, ce sera pris en compte dans le modèle. On peut encore une fois regarder le nuage de points des mesures de chaque souris pour se faire une idée de cette auto-corrélation des mesures d’un individu :
p_glmm_2 <- d2 %>%
ggplot(aes(x=id, y=n_extrasystole, shape=diet, col=exercise)) +
geom_point() +
xlab("numéro de la souris") +
ylab("nombre d'extra-systoles enregistrées") +
theme_bw()
p_glmm_2
A vous de proposer un modèle qui vous semble intéressant d’un point de vue biologique.
J’en propose un ci-dessous pour que vous ayez la syntaxe permettant de fitter le modèle avec R.
res_glmm_3 <- glmer("n_extrasystole ~ diet + exercise + crp + weight + sex + lineage + (1 | id)",
family = poisson,
data=d2bis)
summary(res_glmm_3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: n_extrasystole ~ diet + exercise + crp + weight + sex + lineage +
## (1 | id)
## Data: d2bis
##
## AIC BIC logLik deviance df.resid
## 1859.7 1899.6 -919.8 1839.7 390
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9427 -0.7381 -0.1427 0.5704 3.3465
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.022 1.011
## Number of obs: 400, groups: id, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.037783 0.252469 12.032 < 2e-16 ***
## dietlow sugar -0.984296 0.028425 -34.627 < 2e-16 ***
## exerciseII 0.006326 0.035515 0.178 0.85863
## exerciseIII -0.020522 0.035755 -0.574 0.56599
## exerciseIV -0.032918 0.035867 -0.918 0.35873
## crp 0.017095 0.149704 0.114 0.90909
## weight 0.516434 0.155256 3.326 0.00088 ***
## sexM 0.050634 0.311217 0.163 0.87076
## lineagewild type -2.040132 0.300829 -6.782 1.19e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dtlwsg exrcII exrIII exrcIV crp weight sexM
## dietlowsugr -0.031
## exerciseII -0.071 0.000
## exerciseIII -0.070 0.000 0.498
## exerciseIV -0.070 0.000 0.497 0.493
## crp -0.049 0.000 0.000 0.000 0.000
## weight -0.078 0.000 0.000 0.000 0.000 -0.118
## sexM -0.461 0.000 0.000 0.000 0.000 0.145 0.237
## linegwldtyp -0.586 0.000 0.000 0.000 0.000 -0.054 -0.102 -0.150Que concluons-nous ?
Conclusion
Les take-home messages de cet item sont les suivants :
un test apparié s’applique à des données mesurées deux fois sur chaque individu. Ca revient à mesurer la différence des deux mesures et à comparer cette différence à zéro.
le LMM rajoute un effet aléatoire à un LM, pour capter la corrélation entre des mesures réalisées au sein d’un même individu. Les autres composantes du modèle sont appelées effets fixes.
le GLMM rajoute un effet aléatoire à un GLM, pour la même raison.
un effet aléatoire est tout indiqué pour modéliser l’effet de corrélations entre les variables dues à
- des mesures répétées à l’échelle d’un individu,
- des mesures répétées par des expérimentateurs différents,
- des mesures répétées dans des centres différents.
nous avons vus ici des effets aléatoires correspondant toujours à un “random intercept”, car c’est ce que vous êtes le plus susceptible de vouloir faire plus tard. On peut cependant également proposer des “random slopes”.