Does the blocking factor confound the fixed effect in my model?

Question

I have data from an experiment with the following design:

11 blocks within a forest along an elevational gradient
Within each block, three treatment plots (control, low and high treatment) with unique elevations
From each plot we assessed the height growth of three trees
Sampling was performed annually during 10 years

Here is an illustration of the design:

And here some data:

          year block treatment tree elevation growth 
    1     2001   b1       C      a      10.3   5.34 
    2     2001   b1       C      b      10.3   7.12  
    3     2001   b1       C      c      10.3   4.62
    4     2001   b1       L      a      12.2   4.33
    5     2001   b1       L      b      12.2   5.01
    6     2001   b1       L      c      12.2   6.35
    7     2001   b1       H      a      14.4   7.11
    8     2001   b1       H      b      14.4   6.54
    9     2001   b1       H      c      14.4   4.24
    10    2001   b2       C      a       1.5   5.98
    11    2001   b2       C      b       1.5   4.45
    12    2001   b2       C      c       1.5   5.64
    13    2001   b2       L      a       3.1   5.96
    14    2001   b2       L      b       3.1   5.44
     .
     .
     .

Goal: I want to estimate the effects of treatment, elevation and year and their interaction. To do so, I'm trying to fit a mixed model to the data.

My approach: I specified the following model using nlme in R:
M1 <- lme(growth ~ treatment * elevation * year, random = ~1|block/treatment/tree)

Problem: I have limited statistical knowledge. I'm afraid that the random effect of block confounds the fixed effect of elevation, as each block covers a distinct range of elevations. Could that be the case?

Is this a valid way of specifying the random effects given the design?

It seems like you are using the lme() function from the nlme package and not lmer() as you indicated. Furthermore your random term specifies site and plot which are not present in your example data. Also it seems based on the data you provided you only have one tree per treatment and block? Could you edit your question accordingly? — Stefan, Jan 05 '17 at 18:02
You are right, thanks. I have corrected the formula. There are three trees (a,b,c) per treatment and block. — Patrick Möhl, Jan 05 '17 at 18:48

score 1 · Accepted Answer · edited Apr 13 '17 at 12:44

Here's what I would do. I expanded your dataset a bit so we have something to work with:

dat <- structure(list(year = c(2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 
2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 
2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 
2001L, 2001L, 2001L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 
2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 
2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 
2002L, 2002L, 2002L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L), plot = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 
3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 
9L, 9L, 9L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 
5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 9L, 9L, 9L, 1L, 1L, 
1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 
7L, 7L, 7L, 8L, 8L, 8L, 9L, 9L, 9L), block = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L
), .Label = c("b1", "b2", "b3"), class = "factor"), treatment = structure(c(1L, 
1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 
2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 
2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 
3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 
1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L
), .Label = c("C", "H", "L"), class = "factor"), tree = structure(c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L
), .Label = c("a", "b", "c"), class = "factor"), elevation = c(10.3, 
10.3, 10.3, 12.2, 12.2, 12.2, 14.4, 14.4, 14.4, 1.5, 1.5, 1.5, 
3.1, 3.1, 3.1, 2.4, 2.4, 2.4, 5.5, 5.5, 5.5, 6.4, 6.4, 6.4, 6.9, 
6.9, 6.9, 10.3, 10.3, 10.3, 12.2, 12.2, 12.2, 14.4, 14.4, 14.4, 
1.5, 1.5, 1.5, 3.1, 3.1, 3.1, 2.4, 2.4, 2.4, 5.5, 5.5, 5.5, 6.4, 
6.4, 6.4, 6.9, 6.9, 6.9, 10.3, 10.3, 10.3, 12.2, 12.2, 12.2, 
14.4, 14.4, 14.4, 1.5, 1.5, 1.5, 3.1, 3.1, 3.1, 2.4, 2.4, 2.4, 
5.5, 5.5, 5.5, 6.4, 6.4, 6.4, 6.9, 6.9, 6.9), growth = c(11L, 
7L, 11L, 4L, 5L, 12L, 7L, 8L, 11L, 6L, 10L, 10L, 5L, 10L, 8L, 
5L, 5L, 8L, 11L, 12L, 10L, 9L, 7L, 6L, 10L, 9L, 10L, 6L, 11L, 
6L, 8L, 5L, 6L, 11L, 9L, 12L, 9L, 8L, 8L, 12L, 5L, 8L, 4L, 12L, 
10L, 6L, 8L, 9L, 6L, 5L, 6L, 4L, 4L, 9L, 10L, 9L, 4L, 11L, 9L, 
11L, 5L, 10L, 11L, 4L, 10L, 4L, 10L, 6L, 10L, 4L, 9L, 10L, 9L, 
6L, 5L, 8L, 7L, 11L, 12L, 4L, 12L)), .Names = c("year", "plot", 
"block", "treatment", "tree", "elevation", "growth"), class = "data.frame", row.names = c(NA, 
-81L))

Now I am not super familiar with lme() so I am using lmer() from the lme4 package. Given your design, I would remove year as a fixed effect. It has 10 levels, and since it is a growth experiment, it will probably be significant because trees grow over time. Making sense of a model output with 3 main effects and their interactions when year has 10 levels might be a bit too much. That reduces your model to treatment and block (I wouldn't use elevation unless these small differences in elevation within the plots were purposely chosen).

As for the random part, I would include year (to account for multiple measurements on the same trees) and block:treatment, which represents the plots, because the 3 trees in the plots may not be spatially independent(?!), which should be accounted for.

library(lme4)
m1 <- lmer(growth ~ treatment * block + (1|year/block:treatment), data = dat)
summary(m1)

Linear mixed model fit by REML ['lmerMod']
Formula: growth ~ treatment * block + (1 | year/block:treatment)
   Data: dat

REML criterion at convergence: 364

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6991 -0.8580  0.2137  0.6845  1.6875 

Random effects:
 Groups               Name        Variance Std.Dev.
 block:treatment:year (Intercept) 0.4115   0.6415  
 year                 (Intercept) 0.0000   0.0000  
 Residual                         6.6914   2.5868  
Number of obs: 81, groups:  block:treatment:year, 27; year, 3

Fixed effects:
                   Estimate Std. Error t value
(Intercept)          8.3333     0.9384   8.880
treatmentH           1.0000     1.3271   0.753
treatmentL          -0.4444     1.3271  -0.335
blockb2             -0.6667     1.3271  -0.502
blockb3              0.1111     1.3271   0.084
treatmentH:blockb2  -1.2222     1.8769  -0.651
treatmentL:blockb2   1.0000     1.8769   0.533
treatmentH:blockb3  -1.2222     1.8769  -0.651
treatmentL:blockb3  -0.7778     1.8769  -0.414

To get p-values for the fixed effects you could use the mixed() function from the afex package:

library(afex)
mixed(growth~treatment*block+(1|year/block:treatment), data=dat)

Contrasts set to contr.sum for the following variables: treatment, block
Fitting 4 (g)lmer() models:
[....]
Obtaining 3 p-values:
[...]
           Effect    df F.scaling    F p.value
1       treatment 2, 16      1.00 0.27     .76
2           block 2, 16      1.00 0.51     .61
3 treatment:block 4, 16      1.00 0.51     .73

However, note that this may not be the only way to analyze your dataset. Depending on the actual questions you have, you can build multiple models and compare them. To generally get a better idea about fixed and random effects, you could start reading here.

Thank you very much for input! I need to include treatment * year, because we are interested in the effect sizes for each year. I found that the model with 'block' in addition to 'block:treatment' as random effects in the model has a significantly better AIC. I think there are probably also other (unknown) factors beside elevation that differ between the blocks. The model including 'block' as random effect lowers the standard errors of estimates for the treatment (and increases standard errors for the estimate for elevation). Should I maybe use two dfferent models? — Patrick Möhl, Jan 09 '17 at 12:37
@mountainous I am glad I could help. I am not sure about using two different models... I would use one model that best captures your objectives and goals. There might be other opinions too. You could ask this as a new question perhaps. Here's a few more info when working with lmer() and glmer(). See also under REML for GLMMs? for how you can specify those random terms and what they do. — Stefan, Jan 09 '17 at 16:55

score 0 · Answer 2 · answered Jan 05 '17 at 17:58

0

In your experiment you are blocking by elevation, so yes they are one in the same. You should run your analysis with the blocking factor and not include the elevation. If your blocking is not significant you can drop the blocking factor which will simplify your model. Alternatively you could look at running your model as an ANCOVA with elevation as a co-variate.

Are you measuring the same tree each year? If so, you might need to look at a repeated measures ANOVA.

answered Jan 05 '17 at 17:58

D. Bogner

21

Thank your for your answer. I need to include elevation as fixed term (it is part of our hypothesis). The blocking factor is significant, so I don't want to drop it. Yes, the same trees are measured - shouldn't block:treatment:tree as random effect account for differences between individual trees over the years (different baseline per tree)? – Patrick Möhl Jan 09 '17 at 12:42

Does the blocking factor confound the fixed effect in my model?

2 Answers2