emmremlMultiKernel {EMMREML}R Documentation

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β+Z1u1+Z2u2+...Zkuk+ey=X\beta+Z_1 u_1 +Z_2 u_2+...Z_k u_k+ e where yy is the nn vector of response variable, XX is a nxqn x q known design matrix of fixed effects, ZjZ_j is a nxljn x l_j known design matrix of random effects for j=1,2,...,kj=1,2,...,k, β\beta is nx1n x 1 vector of fixed effects coefficients and U=(u1t,u2t,...,ukt)tU=(u^t_1, u^t_2,..., u^t_k)^t and ee are independent variables with NL(0,blockdiag(σu12K1,σu22K2,...,σuk2Kk))N_L(0, blockdiag(\sigma^2_{u_1} K_1,\sigma^2_{u_2} K_2,...,\sigma^2_{u_k} K_k)) and Nn(0,σe2In)N_n(0, \sigma^2_e I_n) correspondingly. The function produces the BLUPs for the L=l1+l2+...+lkL=l_1+l_2+...+l_k dimensional random effect UU. The variance parameters for random effects are estimated as (w^1,w^2,...,w^k)σu2^(\hat{w}_1, \hat{w}_2,...,\hat{w}_k)*\hat{\sigma^2_u} where w=(w1,w2,...,wk)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

nx1n x 1 numeric vector

X

nxqn x q matrix

Zlist

list of random effects design matrices of dimensions nxl1n x l_1,...,nxlkn x l_k

Klist

list of known relationship matrices of dimensions l1xl1l_1 x l_1,...,lkxlkl_k x l_k

varbetahat

TRUE or FALSE

varuhat

TRUE or FALSE

PEVuhat

TRUE or FALSE

test

TRUE or FALSE

Value

Vu

Estimate of σu2\sigma^2_u

Ve

Estimate of σe2\sigma^2_e

betahat

BLUEs for β\beta

uhat

BLUPs for uu

weights

Estimates of kernel weights

Xsqtestbeta

A χ2\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 χ2\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])

[Package EMMREML version 3.1 Index]