SolveHMME function

Solve Henderson's Mixed Model Equation.

Consider a linear mixed model with normal random effects, [REMOVE_ME]$Y_{ij} = X_{ij}^T\beta + v_i + \epsilon_{ij} [REMOVE_ME_2]$

where $i=1,\ldots,n,\quad j=1,\ldots,m$, or it can be equivalently expressed using matrix notation, [REMOVE_ME]$Y = X\beta + Zv + \epsilon [REMOVE_ME_2]$

where $Y\in \mathrm{R}^{nm}$ is a known vector of observations, $X \in \mathrm{R}^{nm\times p}$ and $Z \in \mathrm{R}^{nm\times n}$ design matrices for $\beta$ and $v$ respectively, $\beta \in \mathrm{R}^p$ and $v\in \mathrm{R}^n$ unknown vectors of fixed effects and random effects where $v_i \sim N(0,\lambda_i)$, and $\epsilon \in \mathrm{R}^{nm}$ an unknown vector random errors independent of random effects. Note that $Z$

does not need to be provided by a user since it is automatically created accordingly to the problem specification.

SolveHMME(X, Y, Mu, Lambda)

Arguments

• X: an $(nm\times p)$ design matrix for $\beta$.
• Y: a length-$nm$ vector of observations.
• Mu: a length-$nm$ vector of initial values for $\mu_i = E(Y_i)$.
• Lambda: a length-$n$ vector of initial values for $\lambda$, variance of $v_i \sim N(0,\lambda_i)$

Returns

a named list containing

• beta: a length-$p$ vector of BLUE $\hat{beta}$.
• v: a length-$n$ vector of BLUP $\hat{v}$.
• leverage: a length-$(mn+n)$ vector of leverages.

Description

Consider a linear mixed model with normal random effects,

where $i=1,\ldots,n,\quad j=1,\ldots,m$, or it can be equivalently expressed using matrix notation,

where $Y\in \mathrm{R}^{nm}$ is a known vector of observations, $X \in \mathrm{R}^{nm\times p}$ and $Z \in \mathrm{R}^{nm\times n}$ design matrices for $\beta$ and $v$ respectively, $\beta \in \mathrm{R}^p$ and $v\in \mathrm{R}^n$ unknown vectors of fixed effects and random effects where $v_i \sim N(0,\lambda_i)$, and $\epsilon \in \mathrm{R}^{nm}$ an unknown vector random errors independent of random effects. Note that $Z$

does not need to be provided by a user since it is automatically created accordingly to the problem specification.

Examples

## small setting for data generation
n = 100; m = 2; p = 2
nm = n*m;   nmp = n*m*p

## generate artifical data
X = matrix(rnorm(nmp, 2,1), nm,p) # design matrix
Y = rnorm(nm, 2,1)                # observation

Mu = rep(1, times=nm)
Lambda = rep(1, times=n)

## solve
ans = SolveHMME(X, Y, Mu, Lambda)

References

Rdpack::insert_ref(key="henderson_estimation_1959",package="HMMEsolver")

Rdpack::insert_ref(key="robinson_that_1991",package="HMMEsolver")

Rdpack::insert_ref(key="mclean_unified_1991",package="HMMEsolver")

Rdpack::insert_ref(key="kim_fast_2017",package="HMMEsolver")