1 Contexte

Lors de ce TP, nous allons nous intéresser à la part de variabilité héréditaire de la durée de fécondité chez des coléoptères femelles. Nous disposons du fichier beetles_fecondity.csv contenant les données suivantes:

Sire: le numéro d’identification du mâle
dam: le numéro d’identification de la femelle
progeny: le numéro d’identification de la descendance (on ne considère ici que les individus femelles)
fec: la durée de fécondité
Nous savons aussi que chaque femelle a été fécondée par un seul mâle.

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.

mydata <- read.csv("beetles_fecondity.csv")
head(mydata, 15)

   Sire dam progeny fec
1     1   1       1  58
2     1   1       2  64
3     1   1       3  50
4     1   1       4  53
5     1   1       5  50
6     1   2       1  54
7     1   2       2  85
8     1   2       3  48
9     1   2       4  55
10    1   2       5  48
11    1   3       1  51
12    1   3       2  76
13    1   3       3  58
14    1   3       4  48
15    1   3       5  55

Comme on s’apperçoit que l’identifiant des individus est un numéro, il faut indiquer à R que ce sont des facteurs et non des variables continues:

mydata$Sire <- as.factor(mydata$Sire)
mydata$dam <- as.factor(mydata$dam)
mydata$progeny <- as.factor(mydata$progeny)

1. Chargez les données, explorez-les et visualisez-les
Essayez de représenter les descendantes (“progeny”) de chaque père (“Sire”) et chaque mère (“dam”). Utilisez par exemple la fonction facet_wrap() dans un ggplot.

pander(summary(mydata))

Sire	dam	progeny	fec
1 : 25	1:120	1:110	Min. : 4.00
3 : 25	2:120	2:110	1st Qu.: 45.00
4 : 25	3:120	3:110	Median : 57.00
6 : 25	4:119	4:110	Mean : 55.38
7 : 25	5: 70	5:109	3rd Qu.: 68.00
11 : 25			Max. :104.00
(Other):399

ggplot(mydata) + geom_point(aes(x=dam, y=fec)) + facet_wrap(~Sire)

2. Combien de mâles différents y a-t-il dans vos données? Combien de femelles différentes? De descendantes?
Il y a 24 mâles différents. Pour connaitre le nombre de femelles, il faut regarder combien de combinaisons mâle-femelle nous retrouvons dans les données. Le nombre de descendantes devrait être égal au nombre de lignes dans le dataframe. Mais pour être sûrs qu’il n’y ait pas de doublons, c’est toujours mieux de vérifier:

nfemelles = nrow(unique(mydata[,c("Sire", "dam")])) #Unique : same dataframe without any duplicate
nprog = nrow(unique(mydata[,c("Sire", "dam", "progeny")]))

Il y a 110 femelles et 549 descendantes.

3. Vos données sont-elles balancées? Qu’est-ce que cela implique?
Dans la figure que vous avez réalisé au point 1, vous observez tout de suite que les mâles n’ont pas tous été associés à 5 femelles. Les données ne sont donc pas balancées. On peut aussi le voir dans la table ci-dessous:

table(mydata[,c("Sire","dam")])

    dam
Sire 1 2 3 4 5
  1  5 5 5 5 5
  2  5 5 5 5 0
  3  5 5 5 5 5
  4  5 5 5 5 5
  5  5 5 5 5 0
  6  5 5 5 5 5
  7  5 5 5 5 5
  8  5 5 5 5 0
  9  5 5 5 5 0
  10 5 5 5 5 0
  11 5 5 5 5 5
  12 5 5 5 5 5
  13 5 5 5 5 5
  14 5 5 5 5 5
  15 5 5 5 5 5
  16 5 5 5 5 0
  17 5 5 5 5 5
  18 5 5 5 5 0
  19 5 5 5 4 5
  20 5 5 5 5 0
  21 5 5 5 5 0
  22 5 5 5 5 5
  23 5 5 5 5 5
  24 5 5 5 5 0

Cela implique que les estimations des paramètres calculées avec la méthode lm ne seront pas correctes.

4. Quelle est selon vous l’équation du modèle que l’expérimentateur cherche à estimer? Combien de paramètres cherche-t-on à estimer?
\(y = \mu + a_i + b_{j(i)} + \epsilon_{ijk}\)
Nous cherchons à estimer 4 paramètres: \(\mu\), \(\sigma^2_{a}\), \(\sigma^2_{b}\) et \(\sigma^2\)

5. Donnez l’équation matricielle de ce modèle.
\(Y = X\beta + Zu + \epsilon\)
Avec: \(Y \sim N(X\beta,V)\)
\(u \sim N(0,G)\)
\(\epsilon \sim N(0,R)\)
Note: comme \(\beta\) est une constante (car facteurs fixes), il n’a pas de variance.

6. Estimez votre modèle à l’aide de lme4::lmer

mod <- lme4::lmer(fec ~ 1 + (1|Sire) + (1|dam:Sire), mydata)

7. Affichez le résumé de votre modèle. (1) Quelle est la valeur des différents paramètres? (2) Comment est calculée la Std. Error de l’intercept ? (3) Qu’indique le test t de la partie fixe ?

\(\hat{\mu} = 55.467\), \(\hat{\sigma}^2_{b}=114.6\), \(\hat{\sigma}^2_{a}=31.8\) et \(\hat{\sigma}^2=177.5\)
Notez que dans le résumé du modèle, Std.Dev. correspond aux écarts-types des facteurs aléatoires, alors que Std. Error correspond aux écarts-types des estimateurs des facteurs fixes.

La Std. Error de l’intercept correspond à l’erreur standard de la moyenne. Elle est calculée à partir de \(\sqrt{\frac{MS_{Sire}}{N}}\) (NB : la \(MS_{Sire}\) peut être calculée à partir d’une ANOVA classique avec la fonction \(lm\). Elle vaut \(1477.51\), et N vaut \(549\)).

Le test t n’est pas très important dans notre cas. Effectivement, celui-ci permet de vérifier si la moyenne des fec est égale à zéro.

summary(mod)

Linear mixed model fit by REML ['lmerMod']
Formula: fec ~ 1 + (1 | Sire) + (1 | dam:Sire)
   Data: mydata

REML criterion at convergence: 4572.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6929 -0.6312  0.0394  0.6400  2.9368 

Random effects:
 Groups   Name        Variance Std.Dev.
 dam:Sire (Intercept) 114.63   10.706  
 Sire     (Intercept)  31.79    5.639  
 Residual             177.48   13.322  
Number of obs: 549, groups:  dam:Sire, 110; Sire, 24

Fixed effects:
            Estimate Std. Error t value
(Intercept)   55.467      1.643   33.77

8. Donnez un intervalle de confiance à \(\mu\). Bonus: donnez cet intervalle pour un niveau de confiance de 0.99.
Ici, Std. Error correspond à l’écart-type de \(\hat{\mu}\) (c’est la même valeur que le Std. Error du résumé du modèle).

emmeans(mod, ~ 1)

 1       emmean   SE   df lower.CL upper.CL
 overall   55.5 1.64 22.9     52.1     58.9

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95

emmeans(mod, ~ 1, level=.99)

 1       emmean   SE   df lower.CL upper.CL
 overall   55.5 1.64 22.9     50.9     60.1

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.99

9. Affichez l’analyse de la variance de votre modèle à l’aide de la fonction ranova(). (1) Si l’on considère un niveau de confiance de 0.99, peut-on considérer que l’effet du facteur mâle est significatif ? (2) Quel est le principe du test effectué ?

Non, car la p-valeur associée à la variance estimée du facteur mâle est plus grande que \(0.01\). Nous ne pouvons donc pas dire que la variabilité due au facteur mâle impacte significativement la variation totale des résultats.
Le test calcul pour chaque facteur le LRT à partir des fonctions de vraissemblances du modèle restreint (sans le facteur) et du modèle complet (\(LRT = -2*log(\frac{L(R)}{L(C)})\)). Les hypothèses du test \(\chi^2\) sont les suivantes :

HO : \(\sigma_{Sire} = 0\)
H1 : \(\sigma_{Sire} \ne 0\)

ranova(mod)

ANOVA-like table for random-effects: Single term deletions

Model:
fec ~ (1 | Sire) + (1 | dam:Sire)
               npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>            4 -2286.0 4580.1                         
(1 | Sire)        3 -2288.4 4582.7  4.656  1    0.03094 *  
(1 | dam:Sire)    3 -2335.4 4676.9 98.834  1    < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La variance de l’effet du facteur femelle est-elle significative ?
Oui
Et si l’on considère un niveau de confiance de 0.95? Les variances des effets des facteurs mâle et femelle sont elles significatives ?
Oui pour les deux.

10. Quelles devraient être les dimensions des matrices G, Z, R et V? (lignesxcolonnes)
G: 134x134 (134=le nombre de femelles différentes + le nombre de mâles différents)
Z: 549x134
R: 549x549
V: 549x549

11. Calculez les matrices G, Z, R et V estimées par le modèle sachant que \(R=\sigma^2I\) et \(V=ZGZ'+R\)

G <- as.matrix(getME(mod, "Lambda") * sigma(mod)) ** 2
Z <- as.matrix(getME(mod,"Z"))
R <-  sigma(mod)**2 * diag(nrow(mydata))
V <- Z %*% G %*% t(Z) + R

12. Vérifiez la dimension de chacune de ces matrices à l’aide de dim(). Est-ce que cela vous semble cohérent?
Oui

dim(G);dim(Z);dim(R);dim(V)

[1] 134 134

[1] 549 134

[1] 549 549

[1] 549 549

13. Visualisez les 100 premières lignes des 100 premières colonnes de la matrie V. Comprenez-vous la structure de cette matrice? Pouvez-vous décomposer chaque valeur (représentée par une couleur différente sur le grahique) sur base des variances estimées par le modèle?
Rouge: \(\hat{\sigma}_a^2\)
Orange: \(\hat{\sigma}_a^2 +\hat{\sigma}_b^2\)
Jaune: \(\hat{\sigma}_a^2 + \hat{\sigma}_b^2 + \hat{\sigma}_{\epsilon}^2\)

mybreaks = sort(unique(c(V)) +1)
plot(V[1:100,1:100], border=NA, breaks=mybreaks)

TP2 : Analyse de la variance à deux critères de classification hiérarchisés

LBRAI2222

janvier 25, 2023

1 Contexte