EMMREML (version 3.1)

emmremlMultiKernel: Function to fit Gaussian mixed model with multiple mixed effects with known covariances.

Description

This function is a wrapper for the emmreml to fit Gaussian mixed model with multiple mixed effects with known covariances. The model fitted is $y=X\beta+Z_1 u_1 +Z_2 u_2+...Z_k u_k+ e$ where $y$ is the $n$ vector of response variable, $X$ is a $n x q$ known design matrix of fixed effects, $Z_j$ is a $n x l_j$ known design matrix of random effects for $j=1,2,...,k$, $\beta$ is $n x 1$ vector of fixed effects coefficients and $U=(u^t_1, u^t_2,..., u^t_k)^t$ and $e$ are independent variables with $N_L(0, blockdiag(\sigma^2_{u_1} K_1,\sigma^2_{u_2} K_2,...,\sigma^2_{u_k} K_k))$ and $N_n(0, \sigma^2_e I_n)$ correspondingly. The function produces the BLUPs for the $L=l_1+l_2+...+l_k$ dimensional random effect $U$. The variance parameters for random effects are estimated as $(\hat{w}_1, \hat{w}_2,...,\hat{w}_k)*\hat{\sigma^2_u}$ where $w=(w_1,w_2,..., w_k)$ are the kernel weights. The function also provides some useful statistics like large sample estimates of variances and PEV.

Usage

emmremlMultiKernel(y, X, Zlist, Klist,
varbetahat=FALSE,varuhat=FALSE, PEVuhat=FALSE, test=FALSE)

Arguments

y

$n x 1$ numeric vector

X

$n x q$ matrix

Zlist

list of random effects design matrices of dimensions $n x l_1$,...,$n x l_k$

Klist

list of known relationship matrices of dimensions $l_1 x l_1$,...,$l_k x l_k$

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$

weights

Estimates of kernel weights

Xsqtestbeta

A $\chi^2$ test statistic based 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

A $\chi^2$ test statistic based 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

####example
#Data from Gaussian process with three 
#(total four, including residuals) independent 
#sources of variation 

n=80
M1<-matrix(rnorm(n*10), nrow=n)

M2<-matrix(rnorm(n*20), nrow=n)

M3<-matrix(rnorm(n*5), nrow=n)


#Relationship matrices
K1<-cov(t(M1))
K2<-cov(t(M2))
K3<-cov(t(M3))
K1=K1/mean(diag(K1))
K2=K2/mean(diag(K2))
K3=K3/mean(diag(K3))


#Generate data
covY<-2*(.2*K1+.7*K2+.1*K3)+diag(n)

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

#training set
Trainsamp<-sample(1:80, 60)

funout<-emmremlMultiKernel(y=Y[Trainsamp], X=matrix(rep(1, n)[Trainsamp], ncol=1), 
Zlist=list(diag(n)[Trainsamp,], diag(n)[Trainsamp,], diag(n)[Trainsamp,]),
Klist=list(K1,K2, K3),
varbetahat=FALSE,varuhat=FALSE, PEVuhat=FALSE, test=FALSE)
#weights
funout$weights
#Correlation of predictions with true values in test set
uhatmat<-matrix(funout$uhat, ncol=3)
uhatvec<-rowSums(uhatmat)

cor(Y[-Trainsamp], uhatvec[-Trainsamp])