emmremlMultiKernel function

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

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.

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

Returns

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