# 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. ```r 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 ```r ####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]) ```

emmremlMultiKernel function