ANOVA represenation of mixed model with crossed random effects and uncorrelated variance structure

Question

I was expecting this to be simple when I started, it seems to not be the case.

I am trying to model yield of bean genotypes. There were only 6 genotypes included. 3 trials were carried out in separate years. Each year, genotypes were grown in complete blocks with each genotype replicated 4 times. That gives be 72 observations made in total.

I m now in dilemma as to how I should specify random effect terms in my model. I think that crossed effects of replicationr nested within yeary and genotypeg x yeary interaction form random effects terms. I am specifying genotype as fixed term because, I wish to estimate overall genotype performance and make contrasts among them.

Currently I am fitting the following:

lmer(yield ~ g + (1 | y/r) + (1 | g:y))

A few of my several dilemmas are:

1. Is it the best possible model to have genotypes' effects estimated ?
2. How is the model written in the ANOVA (vector) form ?
3. How is the random effects term read for the above model ? (Should I say, it is the crossing of interaction random effects of slope ?)
4. I have no specification for variance structure of random effects. Is it important to specify an autoregressive variance smoother ? If so, I don't know how with lme4 package.
5. On a related note, I am estimating heritability based all these information and estimated variance components ? I am faced with a situation where components of phenotypic variance are in random effect terms and the only fixed effect term g holds the numerator of heritability estimate ? Is this situation common ? How do I interpret the heritability estimate from this model (if at all it can be calculated)? A relevant literature referral should be helpful for the latter case.

Please note that asreml for model fitting is not an option for me.

Update 1

The output of summary for the model above is:

REML criterion at convergence: 118.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.11440 -0.49144  0.07263  0.41995  3.05095 

Random effects:
 Groups   Name        Variance Std.Dev.
 year:g   (Intercept) 0.27383  0.52329 
 r:year   (Intercept) 0.00755  0.08689 
 year     (Intercept) 0.93938  0.96922 
 Residual             0.17985  0.42409 
Number of obs: 71, groups:  year:g, 18; r:year, 12; year, 3

Fixed effects:
                  Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)        1.09122    0.64809  3.16536   1.684   0.1860  
gChaumae           0.97434    0.46282 10.14350   2.105   0.0612 .
gDhankute Chirrke  0.34906    0.46101  9.98960   0.757   0.4664  
gTrishuli          1.01198    0.46101  9.98960   2.195   0.0529 .
gWhite OP          0.00806    0.46101  9.98960   0.017   0.9864  
gWP Con Bean      -0.13489    0.46101  9.98960  -0.293   0.7758  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gChaum gDhnkC gTrshl gWhtOP
gChaumae    -0.354                            
gDhnktChrrk -0.356  0.498                     
gTrishuli   -0.356  0.498  0.500              
gWhite OP   -0.356  0.498  0.500  0.500       
gWPConBean  -0.356  0.498  0.500  0.500  0.500

```

Answer 1

I'm going to use the following terminology (adapted from my comments to the question).

• Three years of experiments
• Each year there were four fields available
• Each field was broken into 6 sub-fields
• Each sub-field had a genotype applied to it

Each field serves as a randomized complete block (RCB), so each year had 4 RCB replicates within it.

I'll walk through the options for analysis, talking my way through why I would or wouldn't use a particular analysis.

Response:

• yield

Factors:

• year (in theory can be random effect)
• field (in theory can be random effect)
• genotype (fixed effect)

First thing is to generate some fake data to analyze. Each genotype will get an effect equal to 0.35 times it's number (genotype 2 gets a 0.70 effect, etc). I'm going to give some random field-specific noise to each field, as well as an overall random noise component to represent subfield-to-subfield variation.

Additionally, I'm creating an extra field column called field_unique. I'll explain why soon.

# create the grid of variables
years <- c("yr1", "yr2", "yr3")
fields_per_year <- c("f1", "f2", "f3", "f4")
genotypes <- c("g1", "g2", "g3", "g4", "g5", "g6")
df <- expand.grid(genotype=genotypes, field=fields_per_year, year=years)
df$field_unique <- factor(paste(df$year,df$field, sep=""))

# create a field-specific error
between_field_sd <- 1
x1 <- length(years)*length(fields_per_year)
x2 <- length(genotypes)
set.seed(1)
df$noise_field <- rep(rnorm(x1, mean = 0, sd = between_field_sd), each = x2)

# create noise for every experimental measurement
subfield_noise_sd <- 1
df$noise <- rnorm(n = 72, mean = 0, sd = subfield_noise_sd)

# create the yield with effect from genotype and field, but none from year
df$yield <- 0.35*as.numeric(df$genotype) + df$noise_field + df$noise

str(df)
head(df, 15)
#outputs not included here

Options for fixed and random effects

The simplest, though not necessarily the best, approach will be to consider all three factors as fixed effects. In concept I would consider field to be a random effect, since I would think of each field as drawn from a random distribution of fields. Year could also be a random effect for the same reason. However, having just a few units in a stratum sometimes poses problems for treating something as a random effect. Year should probably be treated as a fixed effect because there's only 3 of them. I would probably favor treating field as a random effect, but we'll do an analysis as a fixed effect as well as a random effect and see if that gives problems. Your data is different, so you'll have to decide what you want to do.

Treating everything as a fixed effect

Even here there are several options for the analysis. The first two models below use field as a factor, and the second two use field_unique. The second and fourth use explicit nesting.

# mod1 is wrong because it doesn't nest field
# mod2 does nest field in year
# mod3 same as 2 for anova, different confints
# mod4 same as 3

mod_aov_1 <- aov(yield ~ genotype + year + field, data=df)
mod_aov_2 <- aov(yield ~ genotype + year/field, data=df)
mod_aov_3 <- aov(yield ~ genotype + year + field_unique, data=df)
mod_aov_4 <- aov(yield ~ genotype + year/field_unique, data=df)

> anova(mod_aov_1) # wrong analysis
Analysis of Variance Table

Response: yield
          Df Sum Sq Mean Sq F value  Pr(>F)  
genotype   5 11.283  2.2566  1.6244 0.16710  
year       2  3.592  1.7962  1.2930 0.28186  
field      3 10.429  3.4764  2.5024 0.06765 .
Residuals 61 84.740  1.3892                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> anova(mod_aov_2)
Analysis of Variance Table

Response: yield
           Df Sum Sq Mean Sq F value    Pr(>F)    
genotype    5 11.283  2.2566  2.3412   0.05351 .  
year        2  3.592  1.7962  1.8635   0.16478    
year:field  9 42.156  4.6841  4.8596 8.585e-05 ***
Residuals  55 53.013  0.9639                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Models 2-4 give the same ANOVA table, while mod_aov_1 gives a different (incorrect) one. The reason that mod_aov_1 is incorrect is because it doesn't know to treat f1 from yr2 as a different entity than f1 from yr1. This problem is taken care of by mod_aov_2, where field is nested in year. mod_aov_3 does not suffer from the same problem as mod_aov_1 because there are no longer non-unique entries for fields, so there's no possibility for the model to confuse two different fields as being the same. It's sometimes recommended to create unique names for units that are unique for the very reason that it protects one from accidentally treating terms as crossed instead of nested, as is found in mod_aov_1. Even without the explicit nesting for mod_aov_3, field_unique ends up nested (implicitly) in year anyway.

I said that the ANOVA tables were the same for models 2-4. While that is true, the confidence intervals are not the same for the terms in all three models. The two using field_unique generate the same confidence intervals regardless of whether the nesting is implicit or explicit, but mod_aov_2, for some reason unknown to me, generates different intervals for the years and fields. Genotype confidence intervals are all the same.

> confint(mod_aov_2)
                      2.5 %      97.5 %
(Intercept)     -1.00550355  0.90657021
genotypeg2      -0.75337539  0.85308790
genotypeg3      -0.38197619  1.22448709
genotypeg4      -0.17944142  1.42702186
genotypeg5       0.15489591  1.76135919
genotypeg6       0.20351884  1.80998213
yearyr2          0.12335545  2.39523762
yearyr3          0.62528334  2.89716551
yearyr1:fieldf2  0.01244717  2.28432933
yearyr2:fieldf2 -2.28612089 -0.01423872
yearyr3:fieldf2 -2.14604649  0.12583567
yearyr1:fieldf3 -1.39424291  0.87763925
yearyr2:fieldf3 -1.18220562  1.08967655
yearyr3:fieldf3 -0.90883849  1.36304368
yearyr1:fieldf4  1.19012711  3.46200928
yearyr2:fieldf4 -0.48628748  1.78559469
yearyr3:fieldf4 -2.00357006  0.26831211

> confint(mod_aov_3)
                        2.5 %     97.5 %
(Intercept)       -1.00550355  0.9065702
genotypeg2        -0.75337539  0.8530879
genotypeg3        -0.38197619  1.2244871
genotypeg4        -0.17944142  1.4270219
genotypeg5         0.15489591  1.7613592
genotypeg6         0.20351884  1.8099821
yearyr2            0.77300906  3.0448912
yearyr3           -0.24234563  2.0295365
field_uniqueyr1f2  0.01244717  2.2843293
field_uniqueyr1f3 -1.39424291  0.8776393
field_uniqueyr1f4  1.19012711  3.4620093
field_uniqueyr2f1 -1.78559469  0.4862875
field_uniqueyr2f2 -2.93577449 -0.6638923
field_uniqueyr2f3 -1.83185922  0.4400229
field_uniqueyr3f1 -0.26831211  2.0035701
field_uniqueyr3f2 -1.27841752  0.9934646
field_uniqueyr3f3 -0.04120952  2.2306726

Note that the same analyses could be done with lm() and the same results gotten, except now there are some additional lines with NAs. Also note the differences in the point estimates between mod_lm_2 and mod_lm_3 (likewise mod_lm_4) for the years and fields coefficients. (anyone know why?)

mod_lm_2 <- lm(yield ~ genotype + year/field, data=df) # same as aov2
mod_lm_3 <- lm(yield ~ genotype + year + field_unique, data=df)
mod_lm_4 <- lm(yield ~ genotype + year/field_unique, data=df)

> summary(mod_lm_2)

Call:
lm(formula = yield ~ genotype + year/field, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.46200 -0.51523  0.05362  0.55835  1.94178 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -0.04947    0.47705  -0.104 0.917791    
genotypeg2       0.04986    0.40081   0.124 0.901460    
genotypeg3       0.42126    0.40081   1.051 0.297844    
genotypeg4       0.62379    0.40081   1.556 0.125363    
genotypeg5       0.95813    0.40081   2.391 0.020276 *  
genotypeg6       1.00675    0.40081   2.512 0.014976 *  
yearyr2          1.25930    0.56682   2.222 0.030440 *  
yearyr3          1.76122    0.56682   3.107 0.002987 ** 
yearyr1:fieldf2  1.14839    0.56682   2.026 0.047626 *  
yearyr2:fieldf2 -1.15018    0.56682  -2.029 0.047292 *  
yearyr3:fieldf2 -1.01011    0.56682  -1.782 0.080263 .  
yearyr1:fieldf3 -0.25830    0.56682  -0.456 0.650400    
yearyr2:fieldf3 -0.04626    0.56682  -0.082 0.935245    
yearyr3:fieldf3  0.22710    0.56682   0.401 0.690224    
yearyr1:fieldf4  2.32607    0.56682   4.104 0.000136 ***
yearyr2:fieldf4  0.64965    0.56682   1.146 0.256704    
yearyr3:fieldf4 -0.86763    0.56682  -1.531 0.131579    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9818 on 55 degrees of freedom
Multiple R-squared:  0.5183,    Adjusted R-squared:  0.3781 
F-statistic: 3.698 on 16 and 55 DF,  p-value: 0.0001464

> summary(mod_lm_3)

Call:
lm(formula = yield ~ genotype + year + field_unique, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.46200 -0.51523  0.05362  0.55835  1.94178 

Coefficients: (2 not defined because of singularities)
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.04947    0.47705  -0.104 0.917791    
genotypeg2         0.04986    0.40081   0.124 0.901460    
genotypeg3         0.42126    0.40081   1.051 0.297844    
genotypeg4         0.62379    0.40081   1.556 0.125363    
genotypeg5         0.95813    0.40081   2.391 0.020276 *  
genotypeg6         1.00675    0.40081   2.512 0.014976 *  
yearyr2            1.90895    0.56682   3.368 0.001389 ** 
yearyr3            0.89360    0.56682   1.576 0.120647    
field_uniqueyr1f2  1.14839    0.56682   2.026 0.047626 *  
field_uniqueyr1f3 -0.25830    0.56682  -0.456 0.650400    
field_uniqueyr1f4  2.32607    0.56682   4.104 0.000136 ***
field_uniqueyr2f1 -0.64965    0.56682  -1.146 0.256704    
field_uniqueyr2f2 -1.79983    0.56682  -3.175 0.002453 ** 
field_uniqueyr2f3 -0.69592    0.56682  -1.228 0.224770    
field_uniqueyr2f4       NA         NA      NA       NA    
field_uniqueyr3f1  0.86763    0.56682   1.531 0.131579    
field_uniqueyr3f2 -0.14248    0.56682  -0.251 0.802473    
field_uniqueyr3f3  1.09473    0.56682   1.931 0.058600 .  
field_uniqueyr3f4       NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9818 on 55 degrees of freedom
Multiple R-squared:  0.5183,    Adjusted R-squared:  0.3781 
F-statistic: 3.698 on 16 and 55 DF,  p-value: 0.0001464

I don't know why that is, and I'm going to post a question on it, probably tomorrow, but in the meantime perhaps someone with greater expertise could explain why the difference in confidence intervals for field vs field_unique in the aov() models, and for why the lm() models contain some NAs. Once I create the question, I'll link to it from here.

Treating field as a random effect

Now, instead of treating field and field_unique as fixed effects, we'll do the analyses using lmer() and treating them as random effects. The primary conceptual problem is still how to structure the model to account for the experimental structure. Fortunately, it's not really much different for lmer() than it is for the aov() and lm() implementations.

Here are the four parallel lmer() models to the mod_aov_X models, though this time treating field and field_unique as random:

library(lme4)
library(lmerTest)

mod_lmer_1 <- lmer(yield ~ genotype + year + (1|field), data=df)
mod_lmer_2 <- lmer(yield ~ genotype + year + (1|year:field), data=df)
mod_lmer_u1 <- lmer(yield ~ genotype + year + (1|field_unique), data=df)
mod_lmer_u2 <- lmer(yield ~ genotype + year + (1|year:field_unique), data=df)

# mod_lmer_1 still gets the nesting wrong
# mod_lmer_2 nests field within year, but makes field random
# mod_lmer_u1 doesn't explicitly nest field_unique, but since it's unique, it's implicitly nested
# mod_lmer_u2 explicitly nests field_unique 

> print(anova(mod_lmer_1, ddf="Kenward-Roger"))
Type III Analysis of Variance Table with Kenward-Roger's method
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
genotype 11.2828  2.2566     5    61  1.6244 0.1671
year      3.5924  1.7962     2    61  1.2930 0.2819

> print(anova(mod_lmer_2, ddf="Kenward-Roger"))
Type III Analysis of Variance Table with Kenward-Roger's method
          Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
genotype 11.2828 2.25656     5    55  2.3412 0.05351 .
year      0.7392 0.36962     2     9  0.3835 0.69211  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The mod_lmer_1 p-value for genotype (0.16710) matches mod_aov_1, again both being incorrect in structure The other three lmer models have genotype p-values (0.05351) that match the other three aov models.

Note, you could also have used these to do include the structure and random effect of field:

mod_aov_6 <- aov(yield ~ genotype + year + Error(year/field), data=df)
mod_aov_7 <- aov(yield ~ genotype + year + Error(field_unique), data=df)

I'm more comfortable with lmer() than aov() for this, though, so I chose to go with lmer().

So, what does this all say with respect to the findings for genotype? In this case it doesn't matter whether you use aov/lm or lmer, as long as you get the nesting structure correct. I think that's the case for your data because it's balanced, but such would not be the case if you had unbalanced data. You'd need to use lmer if that were the case.

Adding year*genotype interaction

Ok, what about the year*genotype interaction? I'm pretty sure you could just add the interaction term without any problems, as such:

mod_aov_int_2 <- aov(yield ~ genotype + year:genotype + year/field, data=df)
mod_lmer_int_2 <- lmer(yield ~ genotype + year + year:genotype + (1|year:field), data=df)

> anova(mod_aov_int_2)
Analysis of Variance Table

Response: yield
              Df Sum Sq Mean Sq F value    Pr(>F)    
genotype       5 11.283  2.2566  2.2502 0.0654918 .  
year           2  3.592  1.7962  1.7911 0.1784566    
genotype:year 10  7.885  0.7885  0.7863 0.6416215    
year:field     9 42.156  4.6841  4.6708 0.0002136 ***
Residuals     45 45.128  1.0028                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> print(anova(mod_lmer_int_2, ddf="Kenward-Roger"))
Type III Analysis of Variance Table with Kenward-Roger's method
               Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
genotype      11.2828 2.25656     5    45  2.2502 0.06549 .
year           0.7691 0.38456     2     9  0.3835 0.69211  
genotype:year  7.8852 0.78852    10    45  0.7863 0.64162  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Bayesian analysis using rstanarm functions

Finally, and I won't go into detail as this is already too long, but personally I'd look to go with a Bayesian analysis using stan_glmer as such:

library(rstanarm)
mod_sglmer_2 <- stan_glmer(yield ~ genotype + year + (1|year:field), data=df)
print(summary(mod_sglmer_2), digits=4)
posterior_interval(mod_sglmer_2, prob=0.95)

If you're not familiar with it, there's tons of good information, and here's an excellent starting point:

http://mc-stan.org/rstanarm/

Complete code

# create the grid of variables
years <- c("yr1", "yr2", "yr3")
fields_per_year <- c("f1", "f2", "f3", "f4")
genotypes <- c("g1", "g2", "g3", "g4", "g5", "g6")
df <- expand.grid(genotype=genotypes, field=fields_per_year, year=years)
df$field_unique <- factor(paste(df$year,df$field, sep=""))

# create a field-specific error
between_field_sd <- 1
x1 <- length(years)*length(fields_per_year)
x2 <- length(genotypes)
set.seed(1)
df$noise_field <- rep(rnorm(x1, mean = 0, sd = between_field_sd), each = x2)

# create noise for every experimental measurement
subfield_noise_sd <- 1
df$noise <- rnorm(n = 72, mean = 0, sd = subfield_noise_sd)

# create the yield with effect from genotype and field, but none from year
df$yield <- 0.35*as.numeric(df$genotype) + df$noise_field + df$noise

str(df)
head(df, 15)
#outputs not included here

# aov analysis

# mod1 is wrong because it doesn't nest field
# mod2 does nest field in year
# mod3 same as 2 for anova, different confints
# mod4 same as 3

mod_aov_1 <- aov(yield ~ genotype + year + field, data=df)
mod_aov_2 <- aov(yield ~ genotype + year/field, data=df)
mod_aov_3 <- aov(yield ~ genotype + year + field_unique, data=df)
mod_aov_4 <- aov(yield ~ genotype + year/field_unique, data=df)

anova(mod_aov_1) # wrong analysis
anova(mod_aov_2)
confint(mod_aov_2)
confint(mod_aov_3)

# lm analysis

mod_lm_2 <- lm(yield ~ genotype + year/field, data=df) # same as aov2
mod_lm_3 <- lm(yield ~ genotype + year + field_unique, data=df)
mod_lm_4 <- lm(yield ~ genotype + year/field_unique, data=df)

summary(mod_lm_2)
summary(mod_lm_3)

# lmer analysis

library(lme4)
library(lmerTest)

mod_lmer_1 <- lmer(yield ~ genotype + year + (1|field), data=df)
mod_lmer_2 <- lmer(yield ~ genotype + year + (1|year:field), data=df)
mod_lmer_u1 <- lmer(yield ~ genotype + year + (1|field_unique), data=df)
mod_lmer_u2 <- lmer(yield ~ genotype + year + (1|year:field_unique), data=df)

# mod_lmer_1 still gets the nesting wrong
# mod_lmer_2 nests field within year, but makes field random
# mod_lmer_u1 doesn't explicitly nest field_unique, but since it's unique, it's implicitly nested
# mod_lmer_u2 explicitly nests field_unique 

print(anova(mod_lmer_1, ddf="Kenward-Roger"))
print(anova(mod_lmer_2, ddf="Kenward-Roger"))

# could also have used aov with Error()
mod_aov_6 <- aov(yield ~ genotype + year + Error(year/field), data=df)
mod_aov_7 <- aov(yield ~ genotype + year + Error(field_unique), data=df)

# add the interaction

mod_aov_int_2 <- aov(yield ~ genotype + year:genotype + year/field, data=df)
mod_lmer_int_2 <- lmer(yield ~ genotype + year + year:genotype + (1|year:field), data=df)
anova(mod_aov_int_2)
print(anova(mod_lmer_int_2, ddf="Kenward-Roger"))

# Bayesian analysis using stan_glm from the rstanarm package
library(rstanarm)
mod_sglmer_2 <- stan_glmer(yield ~ genotype + year + (1|year:field), data=df)
print(summary(mod_sglmer_2), digits=4)
posterior_interval(mod_sglmer_2, prob=0.95)

I think I've gotten things right, but hopefully other members will help clarify or correct anything that's needed.