TP3 : Modèles mixtes
LBRAI2222
janvier 25, 2023
1 Introduction
Vous êtes un mémorant de la faculté d’agro de l’UCLouvain, et vous avez lu l’article publié par Fanchon de Serres (chercheuse du TP1) sur l’impact de l’irrigation sur la croissance de plants de maïs. Vous voudriez utiliser le même protocole pour étudier l’impact de l’irrigation sur la croissance d’une autre variété de maïs. Vous rencontrez cependant un problème; le budget qui vous est alloué par la faculté ne vous permet pas d’installer des capteurs de lumière dans la serre que vous avez à votre disposition. Pourtant, vous savez que l’exposition des plants dans la serre est très hétérogène et que cela risque d’influencer vos résultats.
Quelle solution pouvez-vous mettre en place pour palier à ce
problème?
(Répondez à cette question avant de passer à la suite)
Heureusement, vous avez assisté au cours LBRAI2222. Vous faites donc
le plan suivant:
Vous réalisez l’expérience et rassemblez les données dans un fichier CSV data_TP3.csv. Ce fichier contient les variables suivantes :
M: Biomasse accumulée sur une semaine ( poids frais,
[g/plant])
W: Régime d’irrigation (WW=Well Watered, WLD=Water
light deficit, WD=Water deficit)
Loc: Localisation de la table (North, South, West,
East)
1. Avant toute chose, il s’agit d’explorer et de visualiser vos données afin d’avoir un premier apperçu des variables et éventuels problèmes que vous pourriez rencontrer.
Nous remarquons que nos données sont parfaitement balancées. Nous avons un plan factoriel complet à deux facteurs avec 6 répétitions.
mydata <- read.csv("TP3_data.csv")
pander(table(mydata[,c("W","Loc")]))| East | North | South | West | |
|---|---|---|---|---|
| WD | 6 | 6 | 6 | 6 |
| WLD | 6 | 6 | 6 | 6 |
| WW | 6 | 6 | 6 | 6 |
ggplot(mydata,aes(x=W,y=M, color = Loc)) + geom_point(alpha = 0.5, size=3) +
labs(x= "Watering",
y= "Biomass (FW) [g/plant]")
2. Vous pouvez inclure le facteur ‘Loc’ de deux manières
différentes dans votre modèle. Quelles sont-elles? Qu’est-ce qui
justifierait le choix de l’une ou de l’autre?
Il peut être pris en compte comme
- Facteur fixe:
- Si le chercheur veut connaître l’effet des différentes expositions
dans sa serre sur la croissance des plants
- Ou s’il considère que toutes les expositions possibles sont inclues dans ses données
- Si le chercheur veut connaître l’effet des différentes expositions
dans sa serre sur la croissance des plants
- Facteur aléatoire:
- Si le chercheur ne s’intéresse pas à l’effet des différentes
expositions sur la croissance des plants mais souhaite prendre en compte
la variation due à ce facteur dans son analyse
- Si tous les niveaux d’exposition possible ne sont pas présents dans l’expérience mais que le chercheur souhaite pouvoir généraliser ses conclusions à toutes les situations d’exposition.
- Si le chercheur ne s’intéresse pas à l’effet des différentes
expositions sur la croissance des plants mais souhaite prendre en compte
la variation due à ce facteur dans son analyse
3. Quel sera l’impact de ce choix sur la méthode d’estimation
du modèle utilisé?
Pour un modèle contenant uniquement des facteurs fixes, c’est la méthode
des moindres carrés qui sera utilisée pour estimer les paramètres du
modèle. Si ‘Loc’ est considéré comme facteur aléatoire, c’est la méthode
du maximum de vraisemblance restreint qui sera utilisée.
2 Modèle 1: ANOVA à un facteur
4. Estimez un modèle à un facteur sans prendre en compte le
facteur ‘Loc’ avec la fonction lm. À quoi correspondent les
colonnes estimates dans le summary ?
\((intercept) = \hat{\mu} +
\hat{\alpha}_{WD}\)
\(WWLD = \hat{\alpha}_{WLD} -
\hat{\alpha}_{WD}\)
\(WWW = \hat{\alpha}_{WW} -
\hat{\alpha}_{WD}\)
mod1 <- lm(M ~ W, mydata)
summary(mod1)
Call:
lm(formula = M ~ W, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-133.855 -27.738 3.572 33.975 105.792
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 180.58 10.34 17.45 < 2e-16 ***
WWLD 32.91 14.63 2.25 0.0277 *
WWW 92.02 14.63 6.29 2.51e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 50.68 on 69 degrees of freedom
Multiple R-squared: 0.3706, Adjusted R-squared: 0.3524
F-statistic: 20.32 on 2 and 69 DF, p-value: 1.155e-07anova(mod1) # comme c'est balancé, même résultat que type=IIIAnalysis of Variance Table
Response: M
Df Sum Sq Mean Sq F value Pr(>F)
W 2 104360 52180 20.316 1.155e-07 ***
Residuals 69 177220 2568
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15. Quelles sont les moyennes estimées de la biomasse accumulée pour les trois régimes d’irrigation?
\(\hat{\mu}_{WD}=180.58\)
g/plant
\(\hat{\mu}_{WLD}=180.58+32.91=213.49\)
g/plant
\(\hat{\mu}_{WW}=180.58+92.02=272.6\)
g/plant
6. Quel est le test fait dans l’anova?
\(H_0\): \(\alpha_{WD}=\alpha_{WLD}=\alpha_{WW}=0\)
\(H_1\): Au moins un des \(\alpha_{i} \ne 0\)
7. Pourquoi la Std. Error de l’intercept est-elle différente que celle de WWLD et WWW ?
La Std. Error de l’intercept est celle associée à la moyenne alors que celle de WWLD et WWW est associée à une différence entre deux niveaux des effets du facteur watering.
Nous pouvons les retrouver à partir du tableau de l’ANOVA :
\[Std. Error_{intercept} =
\sqrt{\frac{S^2}{n}}= \sqrt{\frac{2568}{24}}=10.34 \]
\[Std. Error_{WWLD} = Std. Error_{WWW} =
\sqrt{\frac{2*S^2}{n}}= \sqrt{\frac{2*2568}{24}}=14.63 \]
3 Modèle 2: ANOVA à 2 facteurs
8. Réalisez maintenant un modèle en intégrant également le facteur ‘Loc’. Comparez le résultat de l’anova pour le facteur W avec l’anova du premier modèle. Pourquoi la valeur F est-elle beaucoup plus grande dans le deuxième modèle?
Pour rappel, \(F_{obs}=\frac{MSA}{MSE}\). Or, dans le premier modèle la variabilité liée à ‘Loc’ est inclue dans le \(MSE\). La valeur de \(MSE\) étant plus grande, le \(F_{obs}\) sera plus petit et la p-valeur plus grande. Le fait d’inclure le facteur ‘Loc’ dans le deuxième modèle permet donc de diminuer la valeur du \(MSE\) et de mettre en évidence avec plus de certitude un effet de \(W\).
mod2 <- lm(M ~ W + Loc, mydata)
summary(mod2)
Call:
lm(formula = M ~ W + Loc, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-86.882 -21.797 1.577 20.931 76.548
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 207.996 9.245 22.499 < 2e-16 ***
WWLD 32.914 9.245 3.560 0.000693 ***
WWW 92.021 9.245 9.954 9.07e-15 ***
LocNorth -74.393 10.675 -6.969 1.84e-09 ***
LocSouth 20.669 10.675 1.936 0.057118 .
LocWest -55.957 10.675 -5.242 1.79e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.02 on 66 degrees of freedom
Multiple R-squared: 0.7596, Adjusted R-squared: 0.7414
F-statistic: 41.71 on 5 and 66 DF, p-value: < 2.2e-16anova(mod2)# comme c'est balancé, même résultat que type=IIIAnalysis of Variance Table
Response: M
Df Sum Sq Mean Sq F value Pr(>F)
W 2 104360 52180 50.881 4.266e-14 ***
Loc 3 109535 36512 35.603 8.343e-14 ***
Residuals 66 67685 1026
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14 Modèle 3: Modèle mixte
9. Réalisez un troisième modèle en utilisant le facteur ‘Loc’
comme aléatoire. A l’aide de la fonction
lmerTest::show_teststrouvez quelle est la matrice \(L\) utilisée par R pour réaliser le test de
l’anova du troisième modèle sous forme de constrast ?
mod3 <- lmer(M ~ W + (1|Loc), mydata)
summary(mod3)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: M ~ W + (1 | Loc)
Data: mydata
REML criterion at convergence: 694.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.75425 -0.65526 0.04876 0.63049 2.34916
Random effects:
Groups Name Variance Std.Dev.
Loc (Intercept) 1971 44.40
Residual 1026 32.02
Number of obs: 72, groups: Loc, 4
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 180.575 23.143 3.346 7.803 0.002938 **
WWLD 32.914 9.245 66.000 3.560 0.000693 ***
WWW 92.021 9.245 66.000 9.954 9.07e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) WWLD
WWLD -0.200
WWW -0.200 0.500anova(mod3) #=anova de type 3Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
W 104360 52180 2 66 50.881 4.266e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ranova(mod3)ANOVA-like table for random-effects: Single term deletions
Model:
M ~ W + (1 | Loc)
npar logLik AIC LRT Df Pr(>Chisq)
<none> 5 -347.22 704.44
(1 | Loc) 4 -373.53 755.07 52.63 1 4.027e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1test <- lmerTest::show_tests(anova(mod3))
L <- test$W
pander(L)| (Intercept) | WWLD | WWW | |
|---|---|---|---|
| WWLD | 0 | 1 | -0.5 |
| WWW | 0 | 0 | 1 |
10. Donnez l’hypothèse \(H_0\) du test effectué.
\(H_0:\) \(\alpha_{WLD} = \frac{\alpha_{WD} +
\alpha_{WW}}{2}\) et \(\alpha_{WD}=\alpha_{WW}\)
11. Donnez la matrice L telle qu’elle serait utilisée avec
votre matrice X complète.
\(L = \left( \begin{array}{cccc}
0&1&-0.5&-0.5\\ 0&0&1&-1 \end{array}
\right)\)
12. Comment pouvez-vous effectuer le même test à l’aide de la
fonction lmerTest::contest?
lmerTest::contest(mod3, L, joint=T) #joint = T signifie que les deux lignes de la matrice sont un seul test et pas deux test différents Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
1 104360 52179.99 2 66 50.88097 4.266449e-145 Comparaison des modèles
5.1 Intervalles de confiance
13. Calculez les intervales de confiances pour chaque niveau du facteur w pour les trois modèles
emmeans(mod1, ~W) W emmean SE df lower.CL upper.CL
WD 181 10.3 69 160 201
WLD 213 10.3 69 193 234
WW 273 10.3 69 252 293
Confidence level used: 0.95 emmeans(mod2, ~W) W emmean SE df lower.CL upper.CL
WD 181 6.54 66 168 194
WLD 213 6.54 66 200 227
WW 273 6.54 66 260 286
Results are averaged over the levels of: Loc
Confidence level used: 0.95 emmeans(mod3, ~W) W emmean SE df lower.CL upper.CL
WD 181 23.1 3.35 111 250
WLD 213 23.1 3.35 144 283
WW 273 23.1 3.35 203 342
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95 14. Comment sont calculés les SE et les df pour les modèles 1
et 2?
Modèle 1:
Tous les \(Y_i\) sont indépendants et
ont la même variance \(S^2=2568\).
Comme nous avons 24 observations par niveau d’irrigation, l’écart-type
de la moyenne pour chaque niveau est \(SE=\sqrt{\frac{2568}{24}}=10.34\)
\(Df=69=72-df_{W}-1=72-2-1=69\)
Modèle 2:
Tous les \(Y_{ij}\) sont indépendants
et ont la même variance \(S^2=1026\).
Comme nous avons 24 observations par niveau d’irrigation, l’écart-type
de la moyenne pour chaque niveau est \(SE=\sqrt{\frac{1026}{24}}=6.54\)
\(Df=66=72-df_{W}- df_{Loc}-1=72-2-3-1=66\)
15. Pourquoi y a-t-il une grande différence entre les intervalles de confiance des différents modèles?
Ici la différence entre les IC est importante parce que la variance liée au facteur ‘Loc’ est grande.
5.2 Tests
16. Trouvez deux manières différentes de faire un test sur la différence attendue de croissance des plants entre les régimes d’irrigation \(WLD\) et \(WW\) (sur base du modèle 3).
pairs(emmeans(mod3, ~W)) contrast estimate SE df t.ratio p.value
WD - WLD -32.9 9.24 66 -3.560 0.0020
WD - WW -92.0 9.24 66 -9.954 <.0001
WLD - WW -59.1 9.24 66 -6.394 <.0001
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates lmerTest::contest(mod3, c(0,1,-1), joint=F) Estimate Std. Error df t value lower upper Pr(>|t|)
1 -59.10791 9.244506 66 -6.393842 -77.56516 -40.65066 1.908268e-0817. Dans la fonction lmerTest::contest, à quoi
sert l’option joint? Quel est le test qui est effectué selon
que l’on utilise joint=T ou joint=F?
Pour répondre à cette question, nous reprenons le code utilisé précédamment avec la matrice L en testant joint=T et joint=F.
Lorsque joint=T est effectué : si la matrice L comporte plusieurs lignes, la fonction réalisera un test F en prenant en compte toutes les lignes ensembles.
Lorsque joint=F est effectué : si la matrice L comporte plusieurs lignes, la fonction réalisera un test t pour chacune des lignes.
lmerTest::contest(mod3, L, joint=T) Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
1 104360 52179.99 2 66 50.88097 4.266449e-14lmerTest::contest(mod3, L, joint=F) Estimate Std. Error df t value lower upper Pr(>|t|)
WWLD -13.09718 8.005977 66 -1.635925 -29.08162 2.887272 1.066158e-01
WWW 92.02147 9.244506 66 9.954179 73.56422 110.478720 9.073053e-15Dans le cas où la matrice L ne contient qu’une seule ligne, les deux tests donnent le même résultat:
lmerTest::contest(mod3, c(0,1,-1), joint=T) Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
1 41924.94 41924.94 1 66 40.88122 1.908268e-08lmerTest::contest(mod3, c(0,1,-1), joint=F) Estimate Std. Error df t value lower upper Pr(>|t|)
1 -59.10791 9.244506 66 -6.393842 -77.56516 -40.65066 1.908268e-08