Modèle à effets mixtes avec imbrication

34

J'ai des données recueillies à partir d'une expérience organisée comme suit:

Deux sites de 30 arbres chacun. 15 sont traités, 15 sont contrôlés sur chaque site. De chaque arbre, nous échantillonnons trois morceaux de la tige et trois morceaux des racines, soit 6 échantillons de niveau 1 par arbre, représentés par l'un des deux niveaux de facteur (racine, tige). Ensuite, à partir de ces échantillons de tige / racine, nous prélevons deux échantillons en disséquant différents tissus dans l'échantillon, ce qui est représenté par l'un des deux niveaux de facteur pour le type de tissu (type de tissu A, type de tissu B). Ces échantillons sont mesurés en tant que variable continue. Le nombre total d'observations est 720; 2 sites * 30 arbres * (trois échantillons de tige + trois échantillons de racines) * (un échantillon de tissu A + un échantillon de tissu B). Les données ressemblent à ceci ...

        ï..Site Tree Treatment Organ Sample Tissue Total_Length
    1        L  LT1         T     R      1 Phloem           30
    2        L  LT1         T     R      1  Xylem           28
    3        L  LT1         T     R      2 Phloem           46
    4        L  LT1         T     R      2  Xylem           38
    5        L  LT1         T     R      3 Phloem          103
    6        L  LT1         T     R      3  Xylem           53
    7        L  LT1         T     S      1 Phloem           29
    8        L  LT1         T     S      1  Xylem           21
    9        L  LT1         T     S      2 Phloem           56
    10       L  LT1         T     S      2  Xylem           49
    11       L  LT1         T     S      3 Phloem           41
    12       L  LT1         T     S      3  Xylem           30

J'essaie d'adapter un modèle à effets mixtes en utilisant R et lme4, mais je suis novice dans les modèles mixtes. J'aimerais modéliser la réponse sous la forme traitement + facteur de niveau 1 (tige, racine) + facteur de niveau 2 (tissu A, tissu B), avec des effets aléatoires pour les échantillons spécifiques imbriqués dans les deux niveaux.

En R, je le fais en utilisant lmer, comme suit

fit <- lmer(Response ~ Treatment + Organ + Tissue + (1|Tree/Organ/Sample)) 

D'après ma compréhension (... ce qui n'est pas certain, et pourquoi je poste!), Le terme:

(1|Tree/Organ/Sample)

Spécifie que 'Sample' est imbriqué dans les échantillons d'organes, ce qui est imbriqué dans l'arborescence. Ce type d'imbrication est-il pertinent / valide? Désolé si cette question n'est pas claire, dans l'affirmative, veuillez préciser où je peux élaborer.

Erik
la source

Réponses:

33

Je pense que c'est correct.

  • (1|Tree/Organ/Sample)étend à / est équivalent à (1|Tree)+(1|Tree:Organ)+(1|Tree:Organ:Sample)(où :désigne une interaction).
  • Les facteurs fixes Treatment, Organet Tissueautomatiquement être manipulé au bon niveau.
  • Vous devriez probablement inclure Siteun effet fixe (conceptuellement, il s'agit d'un effet aléatoire, mais il n'est pas pratique d'essayer d'estimer la variance entre sites avec seulement deux sites); cela réduira légèrement la variance parmi les arbres.
  • Vous devriez probablement inclure toutes les données dans un cadre de données et le transmettre explicitement lmervia un data=my.data.frameargument.

Vous trouverez peut-être que la FAQ de glmm est utile (elle est axée sur les GLMM mais contient également des éléments relatifs aux modèles mixtes linéaires).

Ben Bolker
la source
Et si Erik voulait spécifier une structure de covariance pour ces interceptions? Par exemple, on peut s'attendre à ce qu'un échantillon avec une interception d'arbre positive ait également une interception d'organe positive. Est-ce que la nidification prend en charge ce problème automatiquement? Sinon, comment pourrait-on spécifier une telle structure?
Sheridan Grant
Je pense que si vous essayez d'écrire les équations de ce cas, vous constaterez que le problème est réglé.
Ben Bolker
13

J'ai lu cette question et la réponse de M. Bolker, et j'ai essayé de reproduire les données (ne me souciant pas beaucoup, franchement, de ce que représente la "longueur" en termes biologiques ou en unités, puis de l'adapter au modèle comme ci-dessus. Je publie les résultats ici. pour partager et obtenir des commentaires sur la présence probable de malentendus.

Le code que j'ai utilisé pour générer ces données fictives peut être trouvé ici , et le jeu de données a la structure interne de l'OP:

     site     tree treatment organ sample tissue    length
1    L       LT01         T  root      1  phloem  108.21230
2    L       LT01         T  root      1  xylem   138.54267
3    L       LT01         T  root      2  phloem   68.88804
4    L       LT01         T  root      2  xylem   107.91239
5    L       LT01         T  root      3  phloem   96.78523
6    L       LT01         T  root      3  xylem    88.93194
7    L       LT01         T  stem      1  phloem  101.84103
8    L       LT01         T  stem      1  xylem   118.30319

La structure est la suivante:

 'data.frame':  360 obs. of  7 variables:
     $ site     : Factor w/ 2 levels "L","R": 1 1 1 1 1 1 1 1 1 1 ...
 $ tree     : Factor w/ 30 levels "LT01","LT02",..: 1 1 1 1 1 1 1 1 1 1 ...
     $ treatment: Factor w/ 2 levels "C","T": 2 2 2 2 2 2 2 2 2 2 ...
 $ organ    : Factor w/ 2 levels "root","stem": 1 1 1 1 1 1 2 2 2 2 ...
     $ sample   : num  1 1 2 2 3 3 1 1 2 2 ...
 $ tissue   : Factor w/ 2 levels "phloem","xylem": 1 2 1 2 1 2 1 2 1 2 ...
     $ length   : num  108.2 138.5 68.9 107.9 96.8 ...

L'ensemble de données a été "truqué" (les commentaires ici seraient les bienvenus) comme suit:

  1. En effet treatment, il existe un effet fixe avec deux interceptions distinctes pour le traitement versus contrôle ( 100versus 70), et aucun effet aléatoire.
  2. Je règle les valeurs pour tissueavec des effets fixes importants avec des interceptions très différentes pour phloemversus xylem( 3versus 6) et des effets aléatoires avec a sd = 3.
  3. Pour organil y a deux « contributions » d'interception au hasard à partir d' un (c. -à ) avec une contribution de l' effet fixé à l'intersection de la fois et .N(0,3)sd = 36rootstem
  4. Car treenous avons juste des effets aléatoires avec un sd = 7.
  5. Car samplej'ai essayé de mettre en place uniquement des effets aléatoires avec sd = 5.
  6. Pour siteaussi juste eff eff aléatoire avec sd = 3.

En raison de la nature catégorique des variables, aucune pente n’a été configurée.

Les résultats du modèle à effets mixtes:

fit <- lmer(length ~ treatment + organ + tissue + (1|tree/organ/sample), data = trees) 

étaient:

 Random effects:
 Groups              Name        Variance  Std.Dev. 
 sample:(organ:tree) (Intercept) 9.534e-14 3.088e-07
 organ:tree          (Intercept) 0.000e+00 0.000e+00
 tree                (Intercept) 4.939e+01 7.027e+00
 Residual                        3.603e+02 1.898e+01
Number of obs: 360, groups:  sample:(organ:tree), 180; organ:tree, 60; tree, 30

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  79.8623     2.7011  52.5000  29.567  < 2e-16 ***
treatmentT   21.4368     3.2539  28.0000   6.588 3.82e-07 ***
organstem     0.1856     2.0008 328.0000   0.093    0.926    
tissuexylem   3.1820     2.0008 328.0000   1.590    0.113    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Comment cela s'est-il passé?

  1. Pour treatmentl'interception sans traitement était 79.8623(je mis en place une moyenne 70), et avec le traitement était 79.8623 + 21.4368 = 101.2991(nous avons mis en place une moyenne de 100.
  2. Car tissueil y avait une 3.1820contribution à la politesse d'interception de xylem, et j'avais mis en place une différence entre phloemet xylemde 3. Les effets aléatoires ne faisaient pas partie du modèle.
  3. En effet organ, les échantillons de la stemaugmentation de l'interception de 0.1856- j'avais mis en place aucune différence dans les effets fixes entre stemet root. L'écart type de ce que je voulais faire comme effets aléatoires n'était pas reflété.
  4. Les treeeffets aléatoires avec un SD de 7surface ont bien fait surface 7.027.
  5. En ce qui concerne sample, l'initiale sdde a 5été sous-estimée comme 3.088.
  6. site ne faisait pas partie du modèle.

Donc, dans l’ensemble, il semble que le modèle corresponde à la structure des données.

Antoni Parellada
la source