emmremlMultiKernel {EMMREML}

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