The univariate mixed-effects model

In univariate mixed-effects model: \(\mathbf{y} = \mathbf{Xb} + \mathbf{Zu} + \mathbf{e}\) (1), we estimate the variance component \(\sigma_u^2\) and \(\sigma_e^2\) for \(\mathbf{u} \sim N(\mathbf{0}, \sigma_u^2\mathbf{I})\) and \(\mathbf{e} \sim N(\mathbf{0}, \sigma_e^2\mathbf{I})\).

In formula (1): \(\mathbf{y}\) is n \(\times\) 1 response vector; \(\mathbf{X}\) and \(\mathbf{Z}\) are n \(\times\) p and n \(\times\) s; \(\mathbf{b}\) is p \(\times\) 1 coefficients vector for the fixed effects; \(\mathbf{u}\) is s \(\times\) 1 matrix for the random effects, \(\mathbf{e}\) is n \(\times\) 1 vector of random errors.

The multivariate mixed-effects model

In multivariate mixed-effects model: \(\mathbf{y} = (\mathbf{I} \otimes \mathbf{X})\mathbf{b} + (\mathbf{I} \otimes \mathbf{Z} )\mathbf{u} + \mathbf{e}\) (2), in which \(\mathbf{y} = \left\{\mathbf{y}_i\right\}_c, \mathbf{b} = \left\{\mathbf{b}_i\right\}_c, \mathbf{u} = \left\{\mathbf{u}_i\right\}_c, \mathbf{e} = \left\{\mathbf{e}_i\right\}_c, i =1, \dots, q\), we estimate the variance-covariance matrix of random effects \(\mathbf{u}\) and \(\mathbf{e}\) on \(q\) response variates, namely \(\mathbf{T}\) and \(\mathbf{E}\), for \(var(\mathbf{u}) = \mathbf{G} = \mathbf{T}\otimes \mathbf{I}_s, var(\mathbf{e}) = \mathbf{R} = \mathbf{E} \otimes \mathbf{I}_n, var(\mathbf{y}) = \mathbf{V} = \mathbf{T}\otimes \mathbf{ZZ}' + \mathbf{E} \otimes \mathbf{I}_n\).

In formula (2): \(\mathbf{y}\) is nq \(\times\) 1 response vector; \(\mathbf{X}\) is n \(\times\) p design matrix for the fixed effects; \(\mathbf{b}\) is pq \(\times\) 1 coefficients vector for the fixed effects; \(\mathbf{Z}\) is n \(\times\) s design matrix for the random effects; \(\mathbf{u}\) is sq \(\times\) 1 vector of the random effects; \(\mathbf{e}\) is nq \(\times\) 1 vector of random errors.

# install.packages("MMeM")
library(MMeM)
data(simdata)

Initialization

# initialize with a positive-definiate var-cov
T.start = matrix(c(10, 5, 5, 15), 2, 2)
E.start = matrix(c(10, 1, 1, 3), 2, 2)

Estimation

# using the Henderson3 estimation mothod
results_henderson = MMeM_henderson3(fml = c(V1,V2) ~ X_vec + (1|Z_vec), data = simdata, factor_X = TRUE)

## Bivariate response: V1 and V2

print(results_henderson)

## $T.estimates
##                 T: V1   T: V1 V2     T: V2
## T.estimates  65.47395   9.969188  9.766204
## T.variance  807.56303 128.824836 20.552160
## 
## $E.estimates
##                E: V1  E: V1 V2     E: V2
## E.estimates 55.74506 11.477502 41.171927
## E.variance  11.09826  8.196899  6.054027