SolveHMME {HMMEsolver}

Solve Henderson's Mixed Model Equation.

Description

Consider a linear mixed model with normal random effects,

Y_{ij} = X_{ij}^T\beta + v_i + \epsilon_{ij}

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

Y = X\beta + Zv + \epsilon

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.

Usage

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)

Value

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.

References

Henderson CR, Kempthorne O, Searle SR, von Krosigk CM (1959). “The Estimation of Environmental and Genetic Trends from Records Subject to Culling.” Biometrics, 15(2), 192. ISSN 0006341X, doi: 10.2307/2527669, http://www.jstor.org/stable/2527669?origin=crossref.

Robinson GK (1991). “That BLUP is a Good Thing: The Estimation of Random Effects.” Statistical Science, 6(1), 15–32. ISSN 0883-4237, doi: 10.1214/ss/1177011926, http://projecteuclid.org/euclid.ss/1177011926.

McLean RA, Sanders WL, Stroup WW (1991). “A Unified Approach to Mixed Linear Models.” The American Statistician, 45(1), 54. ISSN 00031305, doi: 10.2307/2685241, http://www.jstor.org/stable/2685241?origin=crossref.

Kim J (2017). “A Fast Algorithm for Solving Henderson's Mixed Model Equation.” ArXiv e-prints.

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)