emmremlMultiKernel: Function to fit Gaussian mixed model with multiple mixed...
In EMMREML: Fitting Mixed Models with Known Covariance Structures

Description Usage Arguments Value Examples

This function is a wrapper for the emmreml to fit Gaussian mixed model with multiple mixed effects with known covariances. The model fitted is y=Xβ+Z_1 u_1 +Z_2 u_2+...Z_k u_k+ e where y is the n vector of response variable, X is a n x q known design matrix of fixed effects, Z_j is a n x l_j known design matrix of random effects for j=1,2,...,k, β is n x 1 vector of fixed effects coefficients and U=(u^t_1, u^t_2,..., u^t_k)^t and e are independent variables with N_L(0, blockdiag(σ^2_{u_1} K_1,σ^2_{u_2} K_2,...,σ^2_{u_k} K_k)) and N_n(0, σ^2_e I_n) correspondingly. The function produces the BLUPs for the L=l_1+l_2+...+l_k dimensional random effect U. The variance parameters for random effects are estimated as (\hat{w}_1, \hat{w}_2,...,\hat{w}_k)*\hat{σ^2_u} where w=(w_1,w_2,..., w_k) are the kernel weights. The function also provides some useful statistics like large sample estimates of variances and PEV.

1 2	emmremlMultiKernel(y, X, Zlist, Klist, varbetahat=FALSE,varuhat=FALSE, PEVuhat=FALSE, test=FALSE)

`y`	n x 1 numeric vector
`X`	n x q matrix
`Zlist`	list of random effects design matrices of dimensions n x l_1,...,n x l_k
`Klist`	list of known relationship matrices of dimensions l_1 x l_1,...,l_k x l_k
`varbetahat`	TRUE or FALSE
`varuhat`	TRUE or FALSE
`PEVuhat`	TRUE or FALSE
`test`	TRUE or FALSE

`Vu`	Estimate of σ^2_u
`Ve`	Estimate of σ^2_e
`betahat`	BLUEs for β
`uhat`	BLUPs for u
`weights`	Estimates of kernel weights
`Xsqtestbeta`	A χ^2 test statistic based for testing whether the fixed effect coefficients are equal to zero.
`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 χ^2 test statistic based for testing whether the BLUPs are equal to zero.
`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 β's.
`PEVuhat`	Prediction error variance estimates for the BLUPs.
`loglik`	loglikelihood for the model.

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