Greg.em {iDINGO}

Fitting precision regression models

Description

This function fits the covariance regression model by Hoff and Niu (2012) using EM algorithm with the restriction of diagonal matrix for the noise variance

Usage

Greg.em(formula, data = NULL, R = 1, tol = 1e-10, itmax = 1000, verbose = F)

Arguments

formula

an object of class "formula" used in model.frame function

data

a data frame used in model.frame function

R

rank of the model

tol

a stopping criterion

itmax

maximum number of iteration

verbose

If true, estimation results for each iteration are printed

Value

A

MLE of the baseline covariance matrix

B

MLE of the regression coefficients

Author(s)

Min Jin Ha <[e-mail]>

References

Hoff, P. D. and Niu, X. (2012) A covariance regression model. Statistica Sinica, 22, 729-753.

Examples

library(glasso)
data(gbm)
x = gbm[,1]
Y = as.matrix(gbm[,-1])
p = ncol(Y)
# Estimating inverse covariance matrix using GLasso #
S = cov(Y)
w.upper = which(upper.tri(S))

rhoarray = exp(seq(log(0.001),log(1),length=100))
BIC = rep(0,length(rhoarray))
for (rh in 1:length(rhoarray)) {
    fit.gl1 = glasso(S,rho=rhoarray[rh])
    BIC[rh] = extendedBIC(gamma=0,omegahat=fit.gl1$wi,S=S,n=nrow(Y))
}
rho = rhoarray[which.min(BIC)]
fit.gl2 = glasso(S,rho=rho)
Omega = fit.gl2$wi

# Fitting (Covariance Regression on transformed data)
diag.Omega = diag(Omega)
P = -Omega/diag.Omega
diag(P) = 0

tY = Y 
mdat = apply(tY,2,mean)
sdat = apply(tY,2,sd)
std.tY = t((t(tY) - mdat)/sdat)
smat = diag(sdat)

## rank 1 covariance regression
fit.g = Greg.em(std.tY~x,R=1)