Analysis of Variance

Data

The table from the slides

tbl_vc <- tibble::tibble(Animal = c(4:7),
                         Sire  = c(2, 1, 3, 2),
                         WWG   = c(2.9, 4.0, 3.5, 3.5))
tbl_vc

Model

Model for a single observation of animal \(i\) with sire \(j\)

\[y_i = \mu + s_{ij} + e_i\]

in matrix-vector notation

\[y = X\beta + Zs + e\]

Components

The vector \(\beta\) of fixed effects contains just one component which is the general mean \(\mu\). As a consequence, the matrix \(X\) has 4 rows and 1 column.

mat_X <- matrix(1, nrow=nrow(tbl_vc))
mat_X

     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1

The matrix \(Z\) relates observations to sire effects. In total there are three sires and therefore matrix \(Z\) has three columns.

mat_Z <- matrix(c(0, 1, 0,
                  1, 0, 0,
                  0, 0, 1,
                  0, 1, 0), nrow = nrow(tbl_vc), byrow = TRUE)
mat_Z

     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    0
[3,]    0    0    1
[4,]    0    1    0

The vector of observations \(y\) is taken from the datatable

vec_y <- tbl_vc$WWG
vec_y

[1] 2.9 4.0 3.5 3.5

Estimates

Estimatesof variance components depend on \(F\), \(S\) and \(R\). From the slides we know that

\[F = y^TX(X^TX)^{−1}X^Ty\]

Computation of \(F\)

(ytX <- crossprod(vec_y, mat_X))

     [,1]
[1,] 13.9

(xtx <- crossprod(mat_X))

     [,1]
[1,]    4

(xtxinv <- solve(xtx))

     [,1]
[1,] 0.25

(xty <- crossprod(mat_X, vec_y))

     [,1]
[1,] 13.9

(comp_F <- ytX %*% xtxinv %*% xty)

        [,1]
[1,] 48.3025

Computation of \(S\)

\[S = y^TZ(Z^TZ)^{−1}Z^Ty−y^TX(X^TX)^{−1}X^Ty = y^TZ(Z^TZ)^{−1}Z^Ty - F\]

(ytz <- crossprod(vec_y, mat_Z))

     [,1] [,2] [,3]
[1,]    4  6.4  3.5

(ztz <- crossprod(mat_Z))

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    2    0
[3,]    0    0    1

(ztzinv <- solve(ztz))

     [,1] [,2] [,3]
[1,]    1  0.0    0
[2,]    0  0.5    0
[3,]    0  0.0    1

(zty <- t(ytz))

     [,1]
[1,]  4.0
[2,]  6.4
[3,]  3.5

(comp_S <- ytz %*% ztzinv %*% zty - comp_F)

       [,1]
[1,] 0.4275

Computation of \(R\)

\[y^Ty − y^TZ(Z^TZ)^{−1}Z^Ty - y^TX(X^TX)^{−1}X^Ty= y^Ty - S - F\]

(yty <- crossprod(vec_y))

      [,1]
[1,] 48.91

(comp_R <- yty - comp_S - comp_F)

     [,1]
[1,] 0.18

Results

Estimates are obtained by replacing expected values of componets with observed values and replacing variance components by their estimates (method of moments)

\[\widehat{\sigma_e^2} = \frac{R}{n-q}\] with \(n\) is the number of observations and \(q\) is the number of sires.

(est_err_var <- comp_R / (4-3))

     [,1]
[1,] 0.18

The estimate for \(\sigma_s^2\) is obtained from

\[M = I − X(X^TX)^{−1}X^T\]

(xxtxinvxt <- mat_X %*% xtxinv %*% t(mat_X))

     [,1] [,2] [,3] [,4]
[1,] 0.25 0.25 0.25 0.25
[2,] 0.25 0.25 0.25 0.25
[3,] 0.25 0.25 0.25 0.25
[4,] 0.25 0.25 0.25 0.25

(mat_M <- diag(1,nrow = nrow(tbl_vc)) - xxtxinvxt)

      [,1]  [,2]  [,3]  [,4]
[1,]  0.75 -0.25 -0.25 -0.25
[2,] -0.25  0.75 -0.25 -0.25
[3,] -0.25 -0.25  0.75 -0.25
[4,] -0.25 -0.25 -0.25  0.75

(ztmz <- t(mat_Z) %*% mat_M %*% mat_Z)

      [,1] [,2]  [,3]
[1,]  0.75 -0.5 -0.25
[2,] -0.50  1.0 -0.50
[3,] -0.25 -0.5  0.75

trztmz <- sum(diag(ztmz))
(est_sire_var <- (comp_S - (3-1)*est_err_var) / trztmz)

      [,1]
[1,] 0.027

In R

R provides the function aov() which can be used to get to these results. For this the sire effect is converted into a fixed effect by changing its data type into a factor.

tbl_vc$Sire <- as.factor(tbl_vc$Sire)
aov_sire <- aov(WWG ~ 1 + Sire, data = tbl_vc)
summary(aov_sire)

            Df Sum Sq Mean Sq F value Pr(>F)
Sire         2 0.4275  0.2137   1.187  0.544
Residuals    1 0.1800  0.1800               

Likelihood

Likelihood-based approach can be done with the pedigreemm package. For our example dataset

n_nr_sire <- nlevels(as.factor(tbl_vc$Sire))
ped_sire <- pedigreemm::pedigree(sire = c(NA, NA, NA), dam = c(NA, NA, NA), label = as.character(1:n_nr_sire))
ped_sire

The pedigreemm package contains a function pedigreemm which computed predicted breeding values and variance components.

lme_sire <- pedigreemm::pedigreemm(formula = WWG ~ 1 + (1|Sire),
                                   pedigree = list(Sire = ped_sire), 
                                   data = tbl_vc)
summary(lme_sire)

Linear mixed model fit by REML ['lmerpedigreemm']
Formula: WWG ~ 1 + (1 | Sire)
   Data: tbl_vc

REML criterion at convergence: 5.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.2199 -0.3086  0.1308  0.4394  0.9582 

Random effects:
 Groups   Name        Variance Std.Dev.
 Sire     (Intercept) 0.04662  0.2159  
 Residual             0.16296  0.4037  
Number of obs: 4, groups:  Sire, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)   3.5025     0.2401   14.59