4
$\begingroup$

I have a field experiment looking at the effect of a seed-mix treatment on moth abundance and I am struggling to define the correct random effects structure. My experiment is structured like this:

I have 16 Blocks, each split into 3 Sections, with each Section having a different seed mix (Treatment). Each night, I sample insects in 4 Blocks then the next night move on to the next 4 Blocks. This continues on a rotation Monday - Thursday (see image). So the whole experiment is sampled fully once in each week. I repeated this over 16 weeks over 2 years, amounting to 64 sample nights.

The variation in insect abundance from night-to-night is very large (due to weather) but I am not interested in this effect, so accounting for this variation is important.

The data look like this:

str(Moths)

'data.frame':   768 obs. of  8 variables:
 $ Section  : Factor w/ 48 levels "10BC","10GR",..: 22 23 24 25 26 27 28 29 30 31 ...
 $ Week     : Factor w/ 16 levels "1_2018","1_2019",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Treatment: Factor w/ 3 levels "BC","GR","WF": 1 2 3 1 2 3 1 2 3 1 ...
 $ Abundance: int  5 3 5 7 3 16 6 6 14 8 ...
 $ Year     : Factor w/ 2 levels "2018","2019": 1 1 1 1 1 1 1 1 1 1 ...
 $ Big_block: Factor w/ 4 levels "B_1","B_2","B_3",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Night    : Factor w/ 64 levels "1_2019","10_2018",..: 58 58 58 58 58 58 58 58 58 58 ...
 $ Block    : Factor w/ 16 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...

head(Moths, 10)

   Section   Week Treatment Abundance Year Big_block  Night Block
1      1BC 1_2018        BC         5 2018       B_1 6_2018     1
2      1GR 1_2018        GR         3 2018       B_1 6_2018     1
3      1WF 1_2018        WF         5 2018       B_1 6_2018     1
4      2BC 1_2018        BC         7 2018       B_1 6_2018     2
5      2GR 1_2018        GR         3 2018       B_1 6_2018     2
6      2WF 1_2018        WF        16 2018       B_1 6_2018     2
7      3BC 1_2018        BC         6 2018       B_1 6_2018     3
8      3GR 1_2018        GR         6 2018       B_1 6_2018     3
9      3WF 1_2018        WF        14 2018       B_1 6_2018     3
10     4BC 1_2018        BC         8 2018       B_1 6_2018     4
>

Originally, I thought that this was a partially crossed design as each Block is sampled on multiple Nights and each Night is associated with multiple Blocks. I was originally coding my model (in R - lme4) as so:

Mod1 <- glm.nb(Abundance ~ Treatment + (1|Night) + (1|Block), data = Moths)

With Night as a factor (1:64) and Block as a factor (1:16). The response variable is a count with high over-dispersion, hence the negative binomial error structure.

A statistician at my institute agreed with this formulation, but another statistician said that this does not properly account for the fact that the same Block is being visited repeatedly. Statistician No. 2 said that I also need to account for the fact that the same 4 Blocks are always sampled together on the same night (this level I call Big_block, with 4 unique levels each). Statistician No. 2 recommended the following:

Mod2 <- glm.nb(Abundance ~ Treatment + (1|Big_block/Block/Section/Week), data = Moths)

This includes a random intercept for each Week, nested in each Section, nested in each Block, nested in each Big_block. As each Big_block is only sampled once in each week, this implicitly includes a Night effect... I think.

I am still not convinced by this structure though as I feel like it should be partially crossed, not fully nested. As I see it, the Night happens to 4 Blocks all at the same time, so I don't see how the temporal effect can be nested within Section, rather than 'above' it, as I am visualising it. I think It should be more like this:

Mod3 <- glm.nb(Abundance ~ Treatment + (1|Night) + (1|Block/Section), data = Moths)

Which is almost the same as my original formulation. In fact, I don't know whether (1|Block/Section) is any different to (1|Block) considering it's always the same Section in the same Block.

To recap, I have the following variables:

Abundance: Continuous response
Treatment: Factor (3 levels)
Big_block: Factor (4 levels). Each Big_block contains 4 Blocks
Block: Factor (16 levels). Each Block contains 3 Sections
Section: Factor (48 levels). Each Section contains 1 Treatment
Night: Factor with 64 levels
Week: Factor with 16 levels (Each Week contains 4 Nights).

I've been working on this problem for a long time, reading books and forums, and I'm just going round in circles. I hope someone on here can help put me out of my misery!

Sampling design for Year 1. Four blocks are sampled each night on a rotation.

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ Your description is very confusing. For one thing, you mentioned Section at the start, which never seems to appear in any model and isn't mentioned again. Then you introduce site for the first time in your second model (where was it in the first model?). Big_block seems to be an artificial construct (which may or may not make sense, it's difficult to say with the way you've described the setup). Please edit the question to make it clear how all your variables are related, along with the output of str(data) and head(data,10). In the recap it would be better not to include Big_block $\endgroup$
    – Robert Long
    Commented Feb 14, 2020 at 13:59
  • 1
    $\begingroup$ Along with Robert's suggestions for clarification, I wanted to add that Night does not make sense as a random factor to me. Week potentially does, but only in as much as you see Week having a systematic effect on all Blocks sampled within the same week. $\endgroup$
    – Erik Ruzek
    Commented Feb 14, 2020 at 14:24
  • 1
    $\begingroup$ @RobertLong Thanks for your reply. I've edited my question to correct the naming inconsistency and hopefully made it clearer. I agree that Big_block is an artificial construct, that is partly what I am disagreeing about with the statistician at work! The Big_block accounts for the fact that I always sample 4 Blocks on the same Night. But I think that this grouping effect is accounted for already with a crossed Night by Block effect. $\endgroup$
    – Dan
    Commented Feb 15, 2020 at 15:16
  • 1
    $\begingroup$ @ErikRuzek As the night-to-night variation is so high, this needs to be taken into account. E.g., if Monday night is hot, Blocks 1:4 will have very high counts and if Tuesday night is cold and windy, Blocks 5:8 will have very low counts, both within the same Week. I don't think using Week would account for this as the variation in weather happens night-to-night rather than week-to-week. $\endgroup$
    – Dan
    Commented Feb 15, 2020 at 15:19
  • 1
    $\begingroup$ This is very helpful, Dan. I think the only remaining confusion I have is around Section. You say that it is always the same Section in every Block. But then you say that each Block is divided into 3 Sections. And Section somehow overlaps with Treamtment such that the 3 Sections in each Block correspond to the 3 values of Treatment. Is that correct? If so, then I would vote for your Mod1 as the winner. $\endgroup$
    – Erik Ruzek
    Commented Feb 15, 2020 at 16:18

1 Answer 1

Reset to default
4
$\begingroup$

After clarification in the question comments, I agree with your intuition that this is a crossed design, so you should fit random intercepts for Night; and also Section nested within Block. You said:

I don't know whether (1|Block/Section) is any different to (1|Block) considering it's always the same Section in the same Block.

This doesn't make complete sense. I think you mean that every Section belongs to one and only 1 Block. If so, then that is the actual definition of nesting.

So your Mod3 is the most appropriate model:

Mod3 <- glm.nb(Abundance ~ Treatment + (1|Night) + (1|Block/Section), data = Moths)

I can't quite get my head around what your statistician that recommended mod2 is getting at:

Mod2 <- glm.nb(Abundance ~ Treatment + (1|Big_block/Block/Section/Week), data = Moths)

First, this does not address the night-to-night variability at all. Second, I don't see any way in which Week can be nested within Section. The only point that remains is whether to include Big_block as a higher level grouping variable. I think that including Block itself is sufficient, but I would suggest fitting 2 models:

Mod3 <- glm.nb(Abundance ~ Treatment + (1|Night) + (1|Block/Section), data = Moths)

Mod4 <- glm.nb(Abundance ~ Treatment + (1|Night) + (1|Big_block_Block/Section), data = Moths)

I wouldn't be surprised if mod4 resulted in a singular fit, or doesn't converge at all. But if it does (and without a singular fit) I would be interested to see the results of both models.

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ Robert, the Section variable and treatment might be collinear within a block, although we need more information from Dan on this. If so, then I don't think (1|Block/Section) is necessary. Treatment as a fixed effect in a model with (1|Block) would capture the within-block variation that Dan is interested in. $\endgroup$
    – Erik Ruzek
    Commented Feb 15, 2020 at 19:27
  • 1
    $\begingroup$ @ErikRuzek Good point. It's hard to say. Perhaps Dan can share his data! $\endgroup$
    – Robert Long
    Commented Feb 15, 2020 at 19:40
  • $\begingroup$ Thanks for all the feedback. @ErikRuzek is correct in that Section and Treatment are colinear within Block. Each Block is split into three Sections and each Section is allocated one seed mix (Treatment) which remains in place permanently. I tried running the model with both (1|Block) and (1|Block/Section) and the test statistics (z and p values) come out very similar. The std of the Block random effect goes from 0.2183 to 0.2178 when the nested Section is added and the std of Section:Block is 0.02654. Would you agree that this would be expected if the (1|Block\Section) effect is unnecessary? $\endgroup$
    – Dan
    Commented Feb 17, 2020 at 11:42
  • $\begingroup$ @RobertLong That would be great if you could have a look at my data! How do I share my data here? Sorry, I'm new to this site. $\endgroup$
    – Dan
    Commented Feb 17, 2020 at 11:44
  • $\begingroup$ @Dan, try running a likelihood ratio test comparing the two models. You can use anova for this in R. A significant test indicates that the more complicated model is a better fit to the data as compared to the simpler version whereas a non-significant test would indicate that the more parsimonious model, i.e., with only (1|Night) + (1|Block), is the better fit. $\endgroup$
    – Erik Ruzek
    Commented Feb 17, 2020 at 12:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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