13 Variance & Variance Components

Mixed models provide the advantage of being able to estimate the variance of random variables. Instead of looking at a variable as a collection of specific levels to estimate, random effects view variables as being a random drawn from a normal distribution with a standard deviation. The decision of how to designate a variable as random or fixed depends on

library(nlme); library(emmeans); library(performance)
library(lme4)

13.1 Unequal Variance

13.1.1 Case 1: Unequal Variance Due to a Factor

var_ex1 <- here::here(read.csv("data", "MET_trial_variance.csv"))

var_ex1$block <- as.character(var_ex1$block)
hist(var_ex1$yield)
boxplot(yield ~ site, data = var_ex1)

m1_a <- lme(yield ~ site:variety + variety, 
                random = ~ 1 |site/block, 
                na.action = na.exclude, 
                data = var_ex1)

m1_b <- update(m1_a, weights = varIdent(form = ~1|site))

m1_b <- update(m1_a, weights = varIdent(form = ~1|site))

is equivalent to

m1_b <- lme(yield ~ site:variety + variety, 
                random = ~ 1 |site/block,
                weights = varIdent(form = ~1|site), 
                na.action = na.exclude, 
                data = var_ex1)

13.1.2 Case 2: Variance is related to the fitted values

var_ex2 <- read.csv(here::here("data", "single_trial_variance.csv"))

var_ex1$block <- as.character(var_ex1$block)
hist(var_ex2$yield)

m2_a <- lme(yield ~ variety, 
               random = ~ 1 |block, 
               na.action = na.exclude, 
               data = var_ex2)

m2_b <- update(m2_a, weights = varPower())

13.2 Coefficient of Variation

m2_ave <- fixef(m2_b)[1]
names(m2_b) <- NULL

m2_cv = sigma(m2_b)/m2_ave*100
m2_cv

13.2.1 Looking at Variance Components

var_comps <- read.csv(here::here("data", "potato_tuber_size.csv"))