# components to which each variance covariance matrix is estimated print(X) # variance covariance matrices for the components print(Ds) # model matrix for each component Zs=lapply(data.frame(X),function(x){model.matrix(~x-1)}) print(Zs) # heteroscedastic errors (sigma's) sigmas=c(1,2) Sigma=diag(c(Zs[[1]]%*%sigmas)) print(Sigma) # expected variance covariance matirx V=array(0,c(8,8)) for(i in 1:3){ V=V+Zs[[i]]%*%Ds[[i]]%*%t(Zs[[i]]) } V=V+Sigma print(V) # Scaling parameter for gamma distribution Nu=10 # observed data matrix Y0=matrix(rnorm(8000),1000)%*%chol(V)/sqrt(rgamma(1000,Nu/2,Nu/2)) # model fitting # Covariance matrix type: # 1: full matrix # 2: diagonal different variances # 3: diagonal same variance res=lmm(Y0,Zs[[1]],Zs,covType=c(1,2,3)) # estimated sigma's print(res[[1]]) # estimated variance covariance matrices for the components print(res[[2]]) # estimated Nu print(res[[7]])