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

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

$$$$

• How many years are there? Why is replication nested within year ? Why us g a fixed effect but also included as a grouping variable in an interaction ? Please can you explain the design and your research question is more detail ? Oct 7, 2020 at 17:06
• There are 3 years. Randomized complete blocks were formed in each years' trial. Block/replication effects are not considered to be same for all years as the sites for each year where trial was set differered. Doesn't nesting them take that into account? Further, I m also interested in interaction of genotype and year just because that might explain significant variability if yield. Oct 7, 2020 at 23:56
• Why does lm(y ~ g + y + r) plus any interactions you are interested in not answer your question(s). Please add the output of summary() for your model into your post. Oct 8, 2020 at 8:59
• I am starting to understand the problem, but please post the output of summary(mymodel) into the the post Oct 9, 2020 at 14:28
• @dd_rookie Did you figure this out? Is your structure like this? Three years / Within each year there were 4 RCBs (randomized complete blocks) / Each RCB was a patch of land, broken into 6 sub-patches, and to each of those sub-patches a genotype was applied. // To say it another way, each year you had 4 fields available, and those 4 fields were each broken down into 6 separate sections, and each section had one of the 6 genotypes applied to it. Is that a good description of how you ran the experiment? Oct 23, 2020 at 2:15

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)
#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)

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.

• Thank you for nice exposition. You have discussed some very important ideas here including analogy between fixed component estimates of mixed model and linear model which I hadn't noticed earlier. Furthermore, I would also be interested to know what makes certain effects to be unestimable (with NA coefficients) in the linear model. I suspect the differing coefficient estimates for nested and uniquely level coded mdoel is due to different estimates for year in mod_aov_2 and mod_aov_4. But coefficients for mod_aov_2 and mod_aov_3 should be same. I'd surprised if not. Oct 29, 2020 at 23:24