emmremlMultivariate: Function to fit multivariate Gaussian mixed model with with...
In EMMREML: Fitting Mixed Models with Known Covariance Structures
Description
Usage
Arguments
Value
Examples
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
1
2
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
χ^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
χ^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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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,])
Example output
Loading required package: Matrix
[,1] [,2] [,3]
[1,] 374.7281 -514.14482 130.95019
[2,] -514.1448 1639.95151 15.75211
[3,] 130.9502 15.75211 273.52187
l 1 l 2 l 3 l 4 l 5
[1,] -6.71388 -22.665209 -8.51306 2.623597 23.55152
[2,] -46.99962 4.558103 14.32785 -59.957458 -86.59351
[3,] -14.74961 0.628672 -19.59707 -40.182521 14.19887
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.0000000 -0.6197034 -0.6602926 -0.5975204 -0.4673006 -0.6084007
[2,] -0.6197034 1.0000000 0.8757381 0.9990086 0.8234167 0.9997550
[3,] -0.6602926 0.8757381 1.0000000 0.8821911 0.7030935 0.8811386
[4,] -0.5975204 0.9990086 0.8821911 1.0000000 0.8096860 0.9996201
[5,] -0.4673006 0.8234167 0.7030935 0.8096860 1.0000000 0.8225948
[6,] -0.6084007 0.9997550 0.8811386 0.9996201 0.8225948 1.0000000
[,1]
[1,] -0.6116819
[,1]
[1,] 0.882494
[,1]
[1,] 0.8209055
EMMREML documentation built on May 2, 2019, 11:01 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value Examples
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
1 2 | 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 |
χ^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 |
χ^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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 | 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,])
|
Example output
Loading required package: Matrix
[,1] [,2] [,3]
[1,] 374.7281 -514.14482 130.95019
[2,] -514.1448 1639.95151 15.75211
[3,] 130.9502 15.75211 273.52187
l 1 l 2 l 3 l 4 l 5
[1,] -6.71388 -22.665209 -8.51306 2.623597 23.55152
[2,] -46.99962 4.558103 14.32785 -59.957458 -86.59351
[3,] -14.74961 0.628672 -19.59707 -40.182521 14.19887
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.0000000 -0.6197034 -0.6602926 -0.5975204 -0.4673006 -0.6084007
[2,] -0.6197034 1.0000000 0.8757381 0.9990086 0.8234167 0.9997550
[3,] -0.6602926 0.8757381 1.0000000 0.8821911 0.7030935 0.8811386
[4,] -0.5975204 0.9990086 0.8821911 1.0000000 0.8096860 0.9996201
[5,] -0.4673006 0.8234167 0.7030935 0.8096860 1.0000000 0.8225948
[6,] -0.6084007 0.9997550 0.8811386 0.9996201 0.8225948 1.0000000
[,1]
[1,] -0.6116819
[,1]
[1,] 0.882494
[,1]
[1,] 0.8209055
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.