EMMREML (version 3.1)

emmremlMultivariate: Function to fit multivariate Gaussian mixed model with with known covariance structure.

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.

Usage

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

Value

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