emmremlMultivariate function

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

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.

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

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