mvcovest: Estimated Covariance Matrix for Multivariate Longitudinal...

In priyakohli5/MLGM: Multivariate Longitudinal Graphical Models

Description

mvcovest returns estimated covariance matrix estimate using the estimates for the elements of modified Cholesky block decomposition (MCBD).

Usage

mvcovest(time,j,gamma.mat,alpha.mat,U,V)

Arguments

time	a vector with T equally or unequally spaced time points.
j	a positive integer for the number of outcomes.
gamma.mat	matrix with rj rows and j columns for estimates for parameters in regressogram models. Here r is the order of time-lag polynomial used for regressogram plots.
alpha.mat	Vector with sj(j+1)/2 estimates for parameters in log innovariogram models. Here s is the order of the time polynomial used for log innovariogram plots.

Value

• bigT is the estimated regression coefficients matrix in the MCBD.

• D.mat is the estimated D matrix in the MCBD.

• Sigma is the estimated covariance matrix.

• logDt.list block of log innovation matrix in a list format.

• logDt elements of log innovation matrix in a matrix format.

References

Kohli, P. Garcia, T. and Pourahmadi, M. 2016 Modeling the Cholesky Factors of Covariance Matrices of Multivariate Longitudinal Data, Journal of Multivariate Analysis, 145, 87-100.

Kohli, P. Du, X. and Shen, H. 2020+ Multivariate Longitudinal Graphical Models (MLGM): An R Package for Visualizing and Modeling Mean \& Dependence Patterns in Multivariate Longitudinal Data, submitted.

Examples

data(Tcells)
time <- c(0, 2, 4, 6, 8, 18, 24, 32, 48, 72)
t <- length(time)
j <- 4
n <- 44
r <- 3
s <- 2
design_mat <- mvrmat(time,j,r,s)
gene.names <- c("FYB", "CD69", "IL2RG", "CDC2")
MVR <- mvr(Tcells,time,j,n,inno=FALSE,inverse=FALSE,loginno=TRUE,plot=FALSE,pch.plot=19,par1.r = 4,par2.r = 4,par1.d=4,par2.d=4)
j2 <- j^2
lags <- 0
for(i in 1:(t-1)){
time.lags <- diff(time,i)
lags <- c(lags,rep(time.lags))
}
lags <- lags[-1]
gamma.init <- matrix(0,nrow=j2,ncol=r+1)
# par(mfrow=c(4,4))
for(i in 1:j2){
 data.phi <- MVR$Phi.plot[[i]]
 best.model <- polyfit(data.phi,lags,degree=r,plot=FALSE,title=gene.names[i],xlabel="Time points",ylabel="Expression Response",indicator=FALSE,index=1,intercept=TRUE)
 gamma.init[i,] <- best.model$output$coefficients[1:(r+1)]
 # if wanted to see initial fits
 # best.model <- polyfit(data.phi,lags,degree=r,plot=TRUE,title=gene.names[i],xlabel="Time points",ylabel="Expression Response",indicator=FALSE,index=1,intercept=TRUE)
 # legend("bottomright",legend=c("observed","fitted"),col=c("black","red"),lty=c(1,2),lwd=c(2,2),cex=0.7)
}
gamma.mat <- matrix(0,((r+1)*j),j)
for(i in 1:j){
row <- ((i-1)*j)+1
col <- ((i-1)*(r+1)) +1
gamma.mat[col:(i*(r+1)),1:j] <- t(gamma.init[row:(i*j),(1:(r+1))])
}
# par(mfrow=c(4,3))
j3 <- j*(j+1)/2
alpha.init <- 0
for(i in 1:j3){
 data.logD <- MVR$logD.elements[i,]
 model.best <- polyfit(data.logD,time,degree=s,plot=FALSE,title="",xlabel="Time points",ylabel="Expression Response",indicator=FALSE,index=4,intercept=TRUE)
 alpha.init <- c(alpha.init,model.best$output$coefficients[1:(s+1)])
 # if wanted to see initial fits
 # model.best <- polyfit(data.logD,time,degree=s,plot=TRUE,title="",xlabel="Time points",ylabel="Expression Response",indicator=FALSE,index=4,intercept=TRUE)
 # legend("topright",legend=c("observed","fitted"),col=c("black","red"),lty=c(1,2),lwd=c(2,2),cex=0.8)
}
alpha.init <- alpha.init[-1]
alpha.mat <- matrix(alpha.init,nrow=(s+1),ncol=j3)
Cov_est <- mvcovest(time,j,gamma.mat,alpha.mat,design_mat$U,design_mat$V)$Sigma

priyakohli5/MLGM documentation built on April 24, 2021, 4:22 p.m.