4

Let's say I fit this model to the Oats data set

lmer(yield ~ Variety + (1|Block/Variety), data=Oats)

(1|Block/Variety) expands to (1|Block) + (1|Block:Variety), which means it's a random effect of Block and a random effect Block:Variety interaction. Can someone please explain to me what this interaction means? I thought it meant that the relationship between intercepts for the varieties is different for the different blocks. So Victory might be higher than Golden Rain in Block 1 but Golden Rain might be higher than Victory in Block 2, and so on. If I am understanding this correctly, why would I ever want to remove variance related to differences between varieties (my fixed factor)?

locus
  • 743
  • 4
  • 17
  • 2
    Please look at [this page](https://stats.stackexchange.com/q/228800/28500) and see if that answers your question. If not, please revise this question to specify what still remains unclear to you. As I understand it, this coding is for nested designs and is not appropriate for a situation where a Variety could be in more than one Block. – EdM Feb 18 '19 at 23:02
  • Thanks for your comment @EdM. I'm not sure how to revise my question because that is exactly the problem. I'd read the post you mention above and also thought that the syntax `(1|Block/Variety)` only applied to nested designs. But apparently it doesn't. As I said above this expands to `(1|Block) + (1|Block:Variety)` which can also apply to crossed effects (for example Ben's answer in this [post](https://stats.stackexchange.com/questions/70556/have-i-correctly-specified-my-model-in-lmer)). – locus Feb 19 '19 at 20:59
  • Thanks for setting me straight on the notation; I've always had a mental block about the nesting I associate with "/". Your formula actually expands to provide a term of `(1|Variety:Block)` as shown in the results within my answer, but that only affects the order of presentation of the random effects, not their values. – EdM Feb 19 '19 at 23:19

1 Answers1

4

This interaction between a fixed and a random factor allows for differences in behavior of the fixed factor among the random factors. Let's run that code on the data set, available in the MASS package in R. (I kept the short variable names provided in that copy of the data.)

BVmodel <- lmer(Y ~ V + (1|B/V), data=oats)

> summary(BVmodel)
Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ V + (1 | B/V)
   Data: oats

REML criterion at convergence: 647.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.66511 -0.67545 -0.00126  0.74643  2.11366 

Random effects:
 Groups   Name        Variance Std.Dev.
 V:B      (Intercept)  19.26    4.389  
 B        (Intercept) 214.48   14.645  
 Residual             524.28   22.897  
Number of obs: 72, groups:  V:B, 18; B, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)  104.500      7.798  13.402
VMarvellous    5.292      7.079   0.748
VVictory      -6.875      7.079  -0.971

Correlation of Fixed Effects:
            (Intr) VMrvll
VMarvellous -0.454       
VVictory    -0.454  0.500

This gives fixed effects for two Varieties (expressed relative to the Intercept that represents the Yield of the "Golden rain" Variety) and sets of random effects for Blocks and for the Variety:Block interaction.

Now let's look at the random effects themselves; for this purpose the ranef() function provides the clearest display, as it shows the random effects with respect to the immediately higher level in the hierarchy. (I omit some of the interaction effects, as they aren't needed to make the point.)

> ranef(BVmodel)
$`V:B`
                (Intercept)
Golden.rain:I     0.4264964
Golden.rain:II    0.7807406
Golden.rain:III  -1.4377120
Golden.rain:IV    1.0514971
Golden.rain:V     0.2028329
Golden.rain:VI   -1.0238550
Marvellous:I     -0.7000427
Marvellous:II     1.1277787
...    
$B
    (Intercept)
I     25.421563
II     2.656992
III   -6.529897
IV    -4.706029
V    -10.582936
VI    -6.259694

Notice that the 6 Block random effects ($B) all add up to 0. These represent how Yield differs (randomly) among the Blocks. The 6 random effects representing the interaction between Block and the "Golden rain" Variety also sum to 0, as do those for the 6 interactions of each of the other 2 Varieties with Block (not shown).

The interaction between each Variety and Block allows the Yield of a Variety to differ (randomly) among Blocks, around the overall Yield for the Variety and the overall random effect for the Block. It doesn't alter the fixed effect associated with each Variety, whether shown as the Intercept for "Golden rain" or as differences from "Golden rain" for the other 2 varieties.

Different yields of the same Variety among different Blocks, even after accounting for between-Block random effects, might be expected in practice and might need to be accounted for. So rather than thinking about this as "remov[ing] variance related to differences between varieties," think about this as allowing for a source of variance within each variety that is related to potentially different behaviors among Blocks.

EdM
  • 57,766
  • 7
  • 66
  • 187
  • thank you very much! That was an excellent explanation, and it makes perfect sense. With your last sentence, do you mean that the model can then account for uninteresting variance related to differences among blocks, so that this variance doesn't "spread" to fixed effects? Is that the purpose of the `V:B` random effect interaction? – locus Feb 21 '19 at 21:08
  • @locus it allows for a more complete description of the structure of the experimental design and lower residual error in principle. That doesn't necessarily mean that it prevents "spread" of variance to fixed effects. Note that with this balanced data set the models with the `B` and `B/V` random effects are almost indistinguishable, with numerically a bit larger standard errors of fixed-effect estimates and numerically lower residual variance in the `B/V` model. In unbalanced designs you might find something a good deal different. – EdM Feb 21 '19 at 21:44
  • Thanks @EdM that makes sense. Do you mind if I ask another question? Based on your answer I think the `$V:B` values within each `Variety` represent the difference from each variety's group intercept, is that right? But what is the reference for `$B` values, is it the mean within each `Block`? – locus Feb 25 '19 at 22:43
  • @locus I find it easiest to start with the top of the hierarchy and move down, the opposite of how the `ranef()` results are displayed. Random effects shown for `$B` represent Block-related _differences from the overall intercept of the model_. To that sum of the fixed intercept and a `$B` value, add the fixed effect for a Variety to get a first-order estimate for that Block and Variety. The `$V:B` random effects then represent the differences for each Block:Variety combination from that first-order estimate for the Block and Variety. – EdM Feb 26 '19 at 02:20