emmreml function

Solver for Gaussian mixed model with known covariance structure.

This function estimates the parameters of the model [REMOVE_ME]$y=X\beta+Z u+ e [REMOVE_ME_2]$ 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.

Description

This function estimates the parameters of the model

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.

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

Returns

• 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])