TP1 : Modèle linéaire général (Rappel)
LBIR2222
janvier 25, 2023
1 Contexte et plan d’expérience
Lors d’une expérience en serre, une chercheuse souhaite estimer
l’impact de différents régimes d’irrigation sur la croissance de plants
de maïs de variétés différentes. Sa serre étant bordée d’arbres au
sud-ouest, les plants sont soumis à des expositions différentes
susceptibles d’impacter leur croissance. Elle décide donc de prendre en
compte l’intensité lumineuse à laquelle sont soumis ses plants. Elle
fait donc l’expérience suivante : 
Elle rassemble ensuite ces données dans un fichier CSV (Comma Separated Values) et nous demande de les analyser. Nous avons donc 72 lignes pour les variables suivantes:
M: Biomasse accumulée sur une semaine ( FW
[g/plant])
PPFD: PPFD cumulé sur une semaine (mol/m²)
W: régime d’irrigation (WW=Well Watered, WLD=Water
light deficit, WD=Water deficit)
V: Variété (A, B ou C)
Note 1 : Le PPFD (Photosynthetic Photon Flux Density) est le nombre de photons dans la gamme de fréquence de 400 à 700 nm par unité de temps sur une unité de surface.
Note 2 : Il s’agit d’un plan factoriel complet avec huits répétitions : \(3(lvl_W)*3(lvl_V)*8(répétitions)=72\)
2 Importation et exploration des données
mydata <- read.csv("plant_growth_data.csv")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.
- Déterminez combien de plantes sont testées pour chaque niveau des deux facteurs catégoriels
Commencez par un rapide résumé (à l’aide de str(),
summary() ou table() par exemple):
pander(summary(mydata))| M | PPFD | W | V |
|---|---|---|---|
| Min. :125.2 | Min. :136.0 | Length:72 | Length:72 |
| 1st Qu.:182.9 | 1st Qu.:144.2 | Class :character | Class :character |
| Median :212.6 | Median :169.5 | Mode :character | Mode :character |
| Mean :223.2 | Mean :171.5 | ||
| 3rd Qu.:248.6 | 3rd Qu.:196.8 | ||
| Max. :360.6 | Max. :211.0 |
head(table(mydata[, c("W","V", "PPFD")])), , PPFD = 136
V
W A B C
WD 2 2 2
WLD 2 2 2
WW 2 2 2
, , PPFD = 147
V
W A B C
WD 2 2 2
WLD 2 2 2
WW 2 2 2
, , PPFD = 192
V
W A B C
WD 2 2 2
WLD 2 2 2
WW 2 2 2
, , PPFD = 211
V
W A B C
WD 2 2 2
WLD 2 2 2
WW 2 2 2A l’aide de ce résumé nous pouvons conclure que 24 plantes sont testées pour chacun des niveaux des facteurs W et V.
- Représentez vos données en représentant les différents niveaux des
facteurs à l’aide de couleurs ou de formes (avec un
ggplot())
ggplot(mydata,aes(x=PPFD,y=M,color=W, shape=V)) + geom_point(alpha = 0.5, size=3) +
labs(shape="Maize varieties",
color = "Watering",
x= "Cumulated PPFD [mol/m²]",
y= "Biomass (FW) [g/plant]")
3 Modélisation
Comme vous suspectez des interactions entre les effets des trois variables explicatives, vous aimeriez estimer un modèle linéaire général complet.
3.1 Premier modèle
Testez un premier modèle contenant les trois variables et toutes leurs interactions. Sur base des résultats de l’analyse ANOVA, quels sont les facteurs qui sont significatifs ?
mod <- lm(M ~ PPFD * W * V, mydata)
anova(mod)Analysis of Variance Table
Response: M
Df Sum Sq Mean Sq F value Pr(>F)
PPFD 1 83842 83842 242.0816 < 2.2e-16 ***
W 2 111637 55818 161.1677 < 2.2e-16 ***
V 2 450 225 0.6497 0.5262
PPFD:W 2 10445 5222 15.0792 6.266e-06 ***
PPFD:V 2 711 356 1.0272 0.3649
W:V 4 2375 594 1.7141 0.1603
PPFD:W:V 4 2368 592 1.7095 0.1613
Residuals 54 18702 346
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1On observe que la variété (V) n’est pas un facteur significatif pour expliquer la croissance des plantes.
3.2 Deuxième modèle
Sur base des résultats de l’ANOVA sur le premier modèle, ajustez votre modèle et refaites une nouvelle ANOVA. Que pouvez-vous en conclure ? Nous allons réaliser un nouveau modèle sans le facteur variété :
mod <- lm(M ~ PPFD * W , mydata)
summary(mod)
Call:
lm(formula = M ~ PPFD * W, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-41.48 -11.56 -3.09 11.17 46.80
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 74.7621 22.1664 3.373 0.00125 **
PPFD 0.6110 0.1272 4.804 9.34e-06 ***
WWLD -53.5261 31.3480 -1.707 0.09243 .
WWW -67.5980 31.3480 -2.156 0.03470 *
PPFD:WWLD 0.5201 0.1799 2.891 0.00519 **
PPFD:WWW 0.9507 0.1799 5.285 1.52e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.31 on 66 degrees of freedom
Multiple R-squared: 0.8933, Adjusted R-squared: 0.8852
F-statistic: 110.5 on 5 and 66 DF, p-value: < 2.2e-16anova(mod)Analysis of Variance Table
Response: M
Df Sum Sq Mean Sq F value Pr(>F)
PPFD 1 83842 83842 224.882 < 2.2e-16 ***
W 2 111637 55818 149.717 < 2.2e-16 ***
PPFD:W 2 10445 5222 14.008 8.498e-06 ***
Residuals 66 24607 373
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nous constatons que l’adaptation faite dans le modèle 2 semble plus intéressante puisque tous les effets sont significatifs. Cependant, nous pouvont observer que la \(Sum Sq\) des résidus a augmenté. Effectivement, lorsque nous avons enlevé le facteur variété, une part de la variabilité expliquée par ce facteur (même si non significative) se retrouve maintenant dans la variabilité résiduelle.
Avant d’aller plus loin, vérifiez bien que vous comprenez chacune des valeurs données dans ce tableau ainsi que les tests effectués.
3.3 Vérification des hypothèses sous-jacentes au modèle
Avant de se lancer dans l’interprétation de votre modèle, vérifiez bien les hypothèses sous-jacentes à la régression (Normalité, indépendance des données, constance de la variance). Ces hypothèses sont-elles respectées ?
par(mfrow = c(2,2))
plot(mod)
par(mfrow=c(1,1))Ici, R vous indique plusieurs points suspects. La valeur standardisée
de ces résidus s’éloigne de 0 de plus de 2 et ils semblent s’éloigner de
la distribution normale sur le qqplot. Avant de les retirer, il serait
intéressant de contacter l’équipe de recherche afin de savoir si des
erreurs de mesures auraient pu se produire avec ces individus.
Par contre, la variance des résidus semble constante.
3.4 Visualisation du modèle avec intervalles de confiance
Lorsque vous visualisez votre modèle, soyez attentifs aux labels des axes et au titre de la légende. Les intervales affichés sont les intervales de confiance pour la réponse moyenne (IC au modèle). Quelle information déjà observée dans l’ANOVA vous montrent les pentes des droites ?
ggplot(mydata, aes(x=PPFD,y=M,color=W, fill=W)) + geom_point(size=1) + geom_smooth(method = 'lm') +
labs(x= "Cumulated PPFD [mol/m²]", y= "Biomass (FW) [g/plant]")+
guides(fill=guide_legend("Watering"), color=guide_legend("Watering"))`geom_smooth()` using formula = 'y ~ x'
Des droites de pentes différentes indiquent une interaction entre l’effet du PPFD et du régime d’irrigation.
4 Prédictions
Notre chercheuse aimerait prédire la croissance de ses plants de maïs lorsque le PPFD est de 150 ou 200, avec un régime d’irrigation WW ou WD (bien irrigué ou déficit en eau). Utilisez votre modèle pour effectuer cette prédiction pour toutes les combinaisons possibles de ces valeurs. Calculez aussi l’intervalle de confiance pour la réponse moyenne et l’intervale de prédiction pour chacune de ces situations.
4.1 Intervale de confiance pour la réponse moyenne (IC au modèle)
xpred <- data.frame("W" = c("WW","WD","WW","WD"),
"PPFD" = c(150, 150, 200, 200))
modIC <- predict(mod, newdata = xpred, interval = "confidence")
res <- data.frame(xpred, modIC)
pander(res)| W | PPFD | fit | lwr | upr |
|---|---|---|---|---|
| WW | 150 | 241.4 | 231.8 | 251 |
| WD | 150 | 166.4 | 156.8 | 176 |
| WW | 200 | 319.5 | 308.8 | 330.2 |
| WD | 200 | 197 | 186.3 | 207.6 |
4.2 Intervale de prédiction (IC aux valeurs)
predIC <- predict(mod, newdata = xpred, interval = "prediction")
res <- data.frame(xpred, predIC)
pander(res)| W | PPFD | fit | lwr | upr |
|---|---|---|---|---|
| WW | 150 | 241.4 | 201.7 | 281.1 |
| WD | 150 | 166.4 | 126.7 | 206.1 |
| WW | 200 | 319.5 | 279.5 | 359.5 |
| WD | 200 | 197 | 156.9 | 237 |
5 Matrice X
On peut afficher la matrice X utilisée par R pour
estimer le modèle à l’aide de model.matrix(). Quelle
différence avec la matrice de votre modèle de régression telle que vue
au cours? Pourquoi ?
pander(model.matrix(mod))| (Intercept) | PPFD | WWLD | WWW | PPFD:WWLD | PPFD:WWW |
|---|---|---|---|---|---|
| 1 | 136 | 0 | 1 | 0 | 136 |
| 1 | 147 | 0 | 1 | 0 | 147 |
| 1 | 192 | 0 | 1 | 0 | 192 |
| 1 | 211 | 0 | 1 | 0 | 211 |
| 1 | 136 | 0 | 1 | 0 | 136 |
| 1 | 147 | 0 | 1 | 0 | 147 |
| 1 | 192 | 0 | 1 | 0 | 192 |
| 1 | 211 | 0 | 1 | 0 | 211 |
| 1 | 136 | 0 | 1 | 0 | 136 |
| 1 | 147 | 0 | 1 | 0 | 147 |
| 1 | 192 | 0 | 1 | 0 | 192 |
| 1 | 211 | 0 | 1 | 0 | 211 |
| 1 | 136 | 0 | 1 | 0 | 136 |
| 1 | 147 | 0 | 1 | 0 | 147 |
| 1 | 192 | 0 | 1 | 0 | 192 |
| 1 | 211 | 0 | 1 | 0 | 211 |
| 1 | 136 | 0 | 1 | 0 | 136 |
| 1 | 147 | 0 | 1 | 0 | 147 |
| 1 | 192 | 0 | 1 | 0 | 192 |
| 1 | 211 | 0 | 1 | 0 | 211 |
| 1 | 136 | 0 | 1 | 0 | 136 |
| 1 | 147 | 0 | 1 | 0 | 147 |
| 1 | 192 | 0 | 1 | 0 | 192 |
| 1 | 211 | 0 | 1 | 0 | 211 |
| 1 | 136 | 1 | 0 | 136 | 0 |
| 1 | 147 | 1 | 0 | 147 | 0 |
| 1 | 192 | 1 | 0 | 192 | 0 |
| 1 | 211 | 1 | 0 | 211 | 0 |
| 1 | 136 | 1 | 0 | 136 | 0 |
| 1 | 147 | 1 | 0 | 147 | 0 |
| 1 | 192 | 1 | 0 | 192 | 0 |
| 1 | 211 | 1 | 0 | 211 | 0 |
| 1 | 136 | 1 | 0 | 136 | 0 |
| 1 | 147 | 1 | 0 | 147 | 0 |
| 1 | 192 | 1 | 0 | 192 | 0 |
| 1 | 211 | 1 | 0 | 211 | 0 |
| 1 | 136 | 1 | 0 | 136 | 0 |
| 1 | 147 | 1 | 0 | 147 | 0 |
| 1 | 192 | 1 | 0 | 192 | 0 |
| 1 | 211 | 1 | 0 | 211 | 0 |
| 1 | 136 | 1 | 0 | 136 | 0 |
| 1 | 147 | 1 | 0 | 147 | 0 |
| 1 | 192 | 1 | 0 | 192 | 0 |
| 1 | 211 | 1 | 0 | 211 | 0 |
| 1 | 136 | 1 | 0 | 136 | 0 |
| 1 | 147 | 1 | 0 | 147 | 0 |
| 1 | 192 | 1 | 0 | 192 | 0 |
| 1 | 211 | 1 | 0 | 211 | 0 |
| 1 | 136 | 0 | 0 | 0 | 0 |
| 1 | 147 | 0 | 0 | 0 | 0 |
| 1 | 192 | 0 | 0 | 0 | 0 |
| 1 | 211 | 0 | 0 | 0 | 0 |
| 1 | 136 | 0 | 0 | 0 | 0 |
| 1 | 147 | 0 | 0 | 0 | 0 |
| 1 | 192 | 0 | 0 | 0 | 0 |
| 1 | 211 | 0 | 0 | 0 | 0 |
| 1 | 136 | 0 | 0 | 0 | 0 |
| 1 | 147 | 0 | 0 | 0 | 0 |
| 1 | 192 | 0 | 0 | 0 | 0 |
| 1 | 211 | 0 | 0 | 0 | 0 |
| 1 | 136 | 0 | 0 | 0 | 0 |
| 1 | 147 | 0 | 0 | 0 | 0 |
| 1 | 192 | 0 | 0 | 0 | 0 |
| 1 | 211 | 0 | 0 | 0 | 0 |
| 1 | 136 | 0 | 0 | 0 | 0 |
| 1 | 147 | 0 | 0 | 0 | 0 |
| 1 | 192 | 0 | 0 | 0 | 0 |
| 1 | 211 | 0 | 0 | 0 | 0 |
| 1 | 136 | 0 | 0 | 0 | 0 |
| 1 | 147 | 0 | 0 | 0 | 0 |
| 1 | 192 | 0 | 0 | 0 | 0 |
| 1 | 211 | 0 | 0 | 0 | 0 |
R a enlevé un niveau de chaque facteur afin d’avoir une matrice inversible.
6 Pour rappel, vous vous
intéressez à connaître l’effet du niveau d’irrigation sur la biomasse
des plants de maïs. Sur base des résultats que vous avez, que
pouvez-vous dire ? N’hésitez pas à utiliser les fonctions
emmeans() et ref_grid() pour illustrer vos
explications ?
Lorsque l’on fait un emmeans(), nous observons qu’il y a
des différences significatives entre les biomasses en fonction de si W=
WD, WLD, WW. De manière générale, une plante qui reçoit plus d’eau aura
une biomasse plus élevée.
Cependant, comme l’indique la console de R, ces comparaisons sont effectuées pour un PPFD moyen de 172 mol/m². Hors, nous savons que nos données ne se comportent pas de la même manière si la valeur de PPFD est faible ou élevée (cfr. intéraction W*PPFD). Ce test est donc à prendre avec des pincettes !
emm_s.t <- emmeans(mod, pairwise ~ W | PPFD)
emm_s.t$emmeans
PPFD = 172:
W emmean SE df lower.CL upper.CL
WD 180 3.94 66 172 187
WLD 215 3.94 66 207 223
WW 275 3.94 66 267 283
Confidence level used: 0.95
$contrasts
PPFD = 172:
contrast estimate SE df t.ratio p.value
WD - WLD -35.7 5.57 66 -6.399 <.0001
WD - WW -95.4 5.57 66 -17.123 <.0001
WLD - WW -59.8 5.57 66 -10.724 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates Nous pouvons calculer les valeurs estimées pour chaque valeur de PPFD :
mod.rg1 <- pairs(ref_grid(mod, at = list(PPFD = c(136))))
summary(mod.rg1) contrast estimate SE df t.ratio p.value
PPFD136 WD - PPFD136 WLD -17.2 8.48 66 -2.030 0.1131
PPFD136 WD - PPFD136 WW -61.7 8.48 66 -7.279 <.0001
PPFD136 WLD - PPFD136 WW -44.5 8.48 66 -5.249 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates mod.rg2 <- pairs(ref_grid(mod, at = list(PPFD = c(147))))
summary(mod.rg2) contrast estimate SE df t.ratio p.value
PPFD147 WD - PPFD147 WLD -22.9 7.11 66 -3.226 0.0055
PPFD147 WD - PPFD147 WW -72.2 7.11 66 -10.154 <.0001
PPFD147 WLD - PPFD147 WW -49.2 7.11 66 -6.928 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates mod.rg3 <- pairs(ref_grid(mod, at = list(PPFD = c(192))))
summary(mod.rg3) contrast estimate SE df t.ratio p.value
PPFD192 WD - PPFD192 WLD -46.3 6.68 66 -6.932 <.0001
PPFD192 WD - PPFD192 WW -114.9 6.68 66 -17.197 <.0001
PPFD192 WLD - PPFD192 WW -68.6 6.68 66 -10.265 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates mod.rg4 <- pairs(ref_grid(mod, at = list(PPFD = c(211))))
summary(mod.rg4) contrast estimate SE df t.ratio p.value
PPFD211 WD - PPFD211 WLD -56.2 9.03 66 -6.224 <.0001
PPFD211 WD - PPFD211 WW -133.0 9.03 66 -14.727 <.0001
PPFD211 WLD - PPFD211 WW -76.8 9.03 66 -8.503 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates Avec les résultats observés ci-dessus, nous pouvons conclure que lorsque les valeurs de PPFD sont faibles, les différences de biomasses sont moins marquées que lorsque les valeurs de PPFD sont élevées. Nous remarquons d’ailleurs qu’à des PPFD de 136 mol/m², la différence entre les biomasses des maïs WD et WLD n’est pas significative.
Cela veut dire qu’en fonction de l’expérience et de la question de recherche il peut être plus intéressant de travailler à un PPFD faible (moins de différences dues au niveau d’irrigation) ou à un PPFD élevé (on observe plus de différences dues au niveau d’irrigation).