In the previous chapter, we introduced how linear models are fitted and how coefficients are estimated. However, these coefficients depend on how categorical factors are encoded in the design matrix (for example, using the first level as a reference). As a result, they are not directly expressed in the units of the response variable. In this section, we show in more detail how to recover estimated marginal means (EMMs)—also known as least-squares means—from the model coefficients, along with their variances.
1) Estimating coefficients
data <-read.csv("../data/example_1.csv") |>mutate(gen =as.factor(gen), block =as.factor(block))n <-12n_b <-3n_g <-4X <-model.matrix(yield ~1+ block + gen, data = data)y <-matrix(data[, "yield"])print(X)
Now let’s explore the matrix form. Specifically, we will make use of the \(\boldsymbol{L}\) matrix to compute all desired values simultaneously, so that:
The \(\boldsymbol{L}\) matrix represents the levels we want to compute. For example, if we want to extract the EMMs of the genotypes, we should include the intercept, the fraction of each block level that we have to add to the mean, and the levels that we want to compute from the genotypes. Notice how \(\boldsymbol{L}\) has as many rows as levels are in the factor to be solved.
We will discuss more in detail about the Coefficient Matrix once we start digging into Linear Mixed Models. For now, we will just say that \(\boldsymbol{C_{11}}\) corresponds to \(Var(\boldsymbol{\hat\beta})\), so
# Compute the errors to compute sigma^2y_hat <- X %*% betaerrors <- y - y_hatSSE <-sum(errors^2)sigma_2 <- SSE / (n -6)C_11 <-solve(t(X) %*% X) * sigma_2C_11
We can also use lm to fit the model and extract the coefficients and emmeans to retrieve the estimated marginal means. The function emmeans does the job automatically but it also returns the \(\boldsymbol{L}\) and \(\boldsymbol{C_{11}}\) matrices used, so we can explore them.
mod <-lm(formula = yield ~1+ block + gen, data = data)beta_lm <-coef(mod)beta_lm
# Average standard error of the differences between two genotype means.SE <-sqrt(diag(var_diff))avg_diff <-mean(SE)avg_diff
[1] 0.5163978
Importantly, since \(\boldsymbol{L}_{diff}\) only contains \(1\) and \(-1\), and as mentioned before, the intercept and the block don’t contribute to the estimation, we can calculate the variance of the differences by simply using the subset of \(\boldsymbol{C_{11}}\) that involves the genotypes (or the factor of interest).