32
$\begingroup$

I have data collected from an experiment organized as follows:

Two sites, each with 30 trees. 15 are treated, 15 are control at each site. From each tree, we sample three pieces of the stem, and three pieces of the roots, so 6 level 1 samples per tree which is represented by one of two factor levels (root, stem). Then, from those stem / root samples, we take two samples by dissecting different tissues within the sample, which is represented by one of two factor levels for tissue type (tissue type A, tissue type B). These samples are measured as a continuous variable. Total number of observations is 720; 2 sites * 30 trees * (three stem samples + three root samples) * (one tissue A sample + one tissue B sample). Data looks like this...

        ï..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

I am attempting to fit a mixed effects model using R and lme4, but am new to mixed models. I'd like to model the response as the Treatment + Level 1 Factor (stem, root) + Level 2 Factor (tissue A, tissue B), with random effects for the specific samples nested within the two levels.

In R, I am doing this using lmer, as follows

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

From my understanding (...which is not certain, and why I am posting!) the term:

(1|Tree/Organ/Sample)

Specifies that 'Sample' is nested within the organ samples, which is nested within the tree. Is this sort of nesting relevant / valid? Sorry if this question is not clear, if so, please specify where I can elaborate.

$\endgroup$
31
$\begingroup$

I think this is correct.

  • (1|Tree/Organ/Sample) expands to/is equivalent to (1|Tree)+(1|Tree:Organ)+(1|Tree:Organ:Sample) (where : denotes an interaction).
  • The fixed factors Treatment, Organ and Tissue automatically get handled at the correct level.
  • You should probably include Site as a fixed effect (conceptually it's a random effect, but it's not practical to try to estimate among-site variance with only two sites); this will reduce the among-tree variance slightly.
  • You should probably include all the data within a data frame, and pass this explicitly to lmer via a data=my.data.frame argument.

You may find the glmm FAQ helpful (it's focused on GLMMs but does have stuff relevant to linear mixed models as well).

$\endgroup$
  • $\begingroup$ What if Erik wanted to specify a covariance structure for these intercepts? I.e. one might expect a sample with a positive Tree intercept to also have a positive Organ intercept. Does the nesting take care of this issue automatically? If not, how could one specify such a structure? $\endgroup$ – Sheridan Grant Nov 7 '17 at 6:23
  • $\begingroup$ I think if you try to write out the equations for that case you'll find that it's taken care of. $\endgroup$ – Ben Bolker Nov 7 '17 at 13:29
12
$\begingroup$

I read this question and Dr. Bolker's answer, and tried to replicate the data (not caring much, frankly, about what "length" represents in biological terms or units, and then fit the model as above. I'm posting the results here to share and seek feedback as to the probable presence of misunderstandings.

The code I used to generate this fictional data can be found here, and the data set has the inner structure of the 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

The structure is as follows:

 '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 ...

The data set was "rigged" (feedback here would be welcome) as follows:

  1. For treatment, there is a fixed effect with two distinct intercepts for treatment versus controls (100 versus 70), and no random effects.
  2. I set the values for tissue with prominent fixed effects with very different intercepts for phloem versus xylem (3 versus 6), and random effects with a sd = 3.
  3. For organ there are two random intercept "contributions" from a $N(0,3)$ (i.e. sd = 3) with a fixed effect contribution to the intercept of 6 for both root and stem.
  4. For tree we just have random effects with a sd = 7.
  5. For sample I tried to set up only random effects with sd = 5.
  6. For for site also just random eff's with sd = 3.

There were no slopes set up, due to the categorical nature of the variables.

The results of the mixed effects model:

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

were:

 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

How did it work out:

  1. For treatment the intercept without treatment was 79.8623 (I set up a mean of 70), and with treatment it was 79.8623 + 21.4368 = 101.2991 (we set up a mean of 100.
  2. For tissue there was a 3.1820 contribution to the intercept courtesy of xylem, and I had set up a difference between phloem and xylem of 3. The random effects were not part of the model.
  3. For organ, samples from the stem increased the intercept by 0.1856 - I had set up no difference in fixed effects between stem and root. The standard deviation of what I wanted to act as random effects was not reflected.
  4. The tree random effects with a sd of 7 surfaced nicely as 7.027.
  5. As for sample, the initial sd of 5 was underemphasized as 3.088.
  6. site was not part of the model.

So, overall, it seems as though the model matches the structure of the data.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.