emmreml {EMMREML}

Solver for Gaussian mixed model with known covariance structure.

Description

This function estimates the parameters of the model

y=X\beta+Z u+ e

where y is the n vector of response variable, X is a n x q known design matrix of fixed effects, Z is a n x l known design matrix of random effects, \beta is q x 1 vector of fixed effects coefficients and u and e are independent variables with N_l(0, \sigma^2_u K) and N_n(0, \sigma^2_e I_n) correspondingly. It also produces the BLUPs and some other useful statistics like large sample estimates of variances and PEV.

Usage

emmreml(y, X, Z, K,varbetahat=FALSE,varuhat=FALSE, PEVuhat=FALSE, test=FALSE)

Arguments

y

n x 1 numeric vector

X

n x q matrix

Z

n x l matrix

K

l x l matrix of known relationships

varbetahat

TRUE or FALSE

varuhat

TRUE or FALSE

PEVuhat

TRUE or FALSE

test

TRUE or FALSE

Value

Vu

Estimate of \sigma^2_u

Ve

Estimate of \sigma^2_e

betahat

BLUEs for \beta

uhat

BLUPs for u

Xsqtestbeta

\chi^2 test statistics for testing whether the fixed effect coefficients are equal to zero.

pvalbeta

pvalues obtained from large sample theory for the fixed effects. We report the pvalues adjusted by the "padjust" function for all fixed effect coefficients.

Xsqtestu

\chi^2 test statistic values for testing whether the BLUPs are equal to zero.

pvalu

pvalues obtained from large sample theory for the BLUPs. We report the pvalues adjusted by the "padjust" function.

varuhat

Large sample variance for the BLUPs.

varbetahat

Large sample variance for the \beta's.

PEVuhat

Prediction error variance estimates for the BLUPs.

loglik

loglikelihood for the model.

Examples

n=200
M1<-matrix(rnorm(n*300), nrow=n)
K1<-cov(t(M1))
K1=K1/mean(diag(K1))

covY<-2*K1+1*diag(n)

Y<-10+crossprod(chol(covY),rnorm(n))


#training set
Trainset<-sample(1:n, 150)

funout<-emmreml(y=Y[Trainset], X=matrix(rep(1, n)[Trainset], ncol=1),
 Z=diag(n)[Trainset,], K=K1)

cor(Y[-Trainset], funout$uhat[-Trainset])