# Function to fit multivariate Gaussian mixed model with with known covariance structure. This function estimates the parameters of the model [REMOVE_ME]$$ Y=BX+GZ+ E [REMOVE_ME_2]$$ where $Y$ is the $d x n$ matrix of response variable, $X$ is a $q x n$ known design matrix of fixed effects, $Z$ is a $l x n$ known design matrix of random effects, $B$ is $d x q$ matrix of fixed effects coefficients and $G$ and $E$ are independent matrix variate variables with $N_{d x l}(0, V_G, K)$ and $N_{d x n}(0, V_E, I_n)$ correspondingly. It also produces the BLUPs for the random effects G and some other statistics. ## Description This function estimates the parameters of the model $$ Y=BX+GZ+ E $$ where $Y$ is the $d x n$ matrix of response variable, $X$ is a $q x n$ known design matrix of fixed effects, $Z$ is a $l x n$ known design matrix of random effects, $B$ is $d x q$ matrix of fixed effects coefficients and $G$ and $E$ are independent matrix variate variables with $N_{d x l}(0, V_G, K)$ and $N_{d x n}(0, V_E, I_n)$ correspondingly. It also produces the BLUPs for the random effects G and some other statistics. ```r emmremlMultivariate(Y, X, Z, K,varBhat=FALSE,varGhat=FALSE, PEVGhat=FALSE, test=FALSE,tolpar=1e-06, tolparinv=1e-06) ``` ## Arguments - `Y`: $d x n$ matrix of response variable - `X`: $q x n$ known design matrix of fixed effects - `Z`: $l x n$ known design matrix of random effects - `K`: $l x l$ matrix of known relationships - `varBhat`: TRUE or FALSE - `varGhat`: TRUE or FALSE - `PEVGhat`: TRUE or FALSE - `test`: TRUE or FALSE - `tolpar`: tolerance parameter for convergence - `tolparinv`: tolerance parameter for matrix inverse ## Returns - **Vg**: Estimate of $V_G$ - **Ve**: Estimate of $V_E$ - **Bhat**: BLUEs for $B$ - **Gpred**: BLUPs for $G$ - **XsqtestB**: $\chi^2$ test statistics for testing whether the fixed effect coefficients are equal to zero. - **pvalB**: pvalues obtained from large sample theory for the fixed effects. We report the pvalues adjusted by the "padjust" function for all fixed effect coefficients. - **XsqtestG**: $\chi^2$ test statistic values for testing whether the BLUPs are equal to zero. - **pvalG**: pvalues obtained from large sample theory for the BLUPs. We report the pvalues adjusted by the "padjust" function. - **varGhat**: Large sample variance for BLUPs. - **varBhat**: Large sample variance for the elements of B. - **PEVGhat**: Prediction error variance estimates for the BLUPs. ## Examples ```r l=20 n<-15 m<-40 M<-matrix(rbinom(m*l,2,.2),nrow=l) rownames(M)<-paste("l",1:nrow(M)) beta1<-rnorm(m)*exp(rbinom(m,5,.2)) beta2<-rnorm(m)*exp(rbinom(m,5,.1)) beta3<- rnorm(m)*exp(rbinom(m,5,.1))+beta2 g1<-M%*%beta1 g2<-M%*%beta2 g3<-M%*%beta3 e1<-sd(g1)*rnorm(l) e2<-(-e1*2*sd(g2)/sd(g1)+.25*sd(g2)/sd(g1)*rnorm(l)) e3<-1*(e1*.25*sd(g2)/sd(g1)+.25*sd(g2)/sd(g1)*rnorm(l)) y1<-10+g1+e1 y2<--50+g2+e2 y3<--5+g3+e3 Y<-rbind(t(y1),t(y2), t(y3)) colnames(Y)<-rownames(M) cov(t(Y)) Y[1:3,1:5] K<-cov(t(M)) K<-K/mean(diag(K)) rownames(K)<-colnames(K)<-rownames(M) X<-matrix(1,nrow=1,ncol=l) colnames(X)<-rownames(M) Z<-diag(l) rownames(Z)<-colnames(Z)<-rownames(M) SampleTrain<-sample(rownames(Z),n) Ztrain<-Z[rownames(Z)%in%SampleTrain,] Ztest<-Z[!(rownames(Z)%in%SampleTrain),] ##For a quick answer, tolpar is set to 1e-4. Correct this in practice. outfunc<-emmremlMultivariate(Y=Y%*%t(Ztrain), X=X%*%t(Ztrain), Z=t(Ztrain), K=K,tolpar=1e-4,varBhat=FALSE, varGhat=FALSE, PEVGhat=FALSE, test=FALSE) Yhattest<-outfunc$Gpred%*%t(Ztest) cor(cbind(Ztest%*%Y[1,],Ztest%*%outfunc$Gpred[1,], Ztest%*%Y[2,],Ztest%*%outfunc$Gpred[2,],Ztest%*%Y[3,],Ztest%*%outfunc$Gpred[3,])) outfuncRidgeReg<-emmremlMultivariate(Y=Y%*%t(Ztrain),X=X%*%t(Ztrain), Z=t(Ztrain%*%M), K=diag(m),tolpar=1e-5,varBhat=FALSE,varGhat=FALSE, PEVGhat=FALSE, test=FALSE) Gpred2<-outfuncRidgeReg$Gpred%*%t(M) cor(Ztest%*%Y[1,],Ztest%*%Gpred2[1,]) cor(Ztest%*%Y[2,],Ztest%*%Gpred2[2,]) cor(Ztest%*%Y[3,],Ztest%*%Gpred2[3,]) ```

emmremlMultivariate function