emmremlMultiKernel {EMMREML} | R Documentation |
Function to fit Gaussian mixed model with multiple mixed effects with known covariances.
Description
This function is a wrapper for the emmreml to fit Gaussian mixed model with multiple mixed effects with known covariances. The model fitted is where
is the
vector of response variable,
is a
known design matrix of fixed effects,
is a
known design matrix of random effects for
,
is
vector of fixed effects coefficients and
and
are independent variables with
and
correspondingly. The function produces the BLUPs for the
dimensional random effect
.
The variance parameters for random effects are estimated as
where
are the kernel weights. The function also provides some useful statistics like large sample estimates of variances and PEV.
Usage
emmremlMultiKernel(y, X, Zlist, Klist,
varbetahat=FALSE,varuhat=FALSE, PEVuhat=FALSE, test=FALSE)
Arguments
y |
|
X |
|
Zlist |
list of random effects design matrices of dimensions |
Klist |
list of known relationship matrices of dimensions |
varbetahat |
TRUE or FALSE |
varuhat |
TRUE or FALSE |
PEVuhat |
TRUE or FALSE |
test |
TRUE or FALSE |
Value
Vu |
Estimate of |
Ve |
Estimate of |
betahat |
BLUEs for |
uhat |
BLUPs for |
weights |
Estimates of kernel weights |
Xsqtestbeta |
A |
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 |
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 |
PEVuhat |
Prediction error variance estimates for the BLUPs. |
loglik |
loglikelihood for the model. |
Examples
####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])