R help - Oct 2011 - Sparse covariance estimation (via glasso) shrinking to a "nonzero" constant

I've only been using R on and off for 9 months and started using the  
glasso package for sparse covariance estimation. I know the concept is  
to shrink some of the elements of the covariance matrix to zero.

However, say I have a dataset that I know has some underlying  
"baseline" covariance/correlation (say, a value of 0.3), how can I  
change or incorporate that into to the code? Basically instead of  
"let's restrict some of these elements are 0 and others have different
nonzero values", I'd like it to be "let's restrict some of
these
elements are 0.3 and others are different or higher/lower than this  
0.3". Does that make any sense?



Here is the glasso code if it helps:
> library(glasso)
> glasso
function (s, rho, zero = NULL, thr = 1e-04, maxit = 10000, approx = FALSE,
     penalize.diagonal = TRUE, start = c("cold", "warm"),
w.init = NULL,
     wi.init = NULL, trace = FALSE)
{
     n = nrow(s)
     BIG = 1e+10
     if (!is.matrix(rho) & length(rho) != 1 & length(rho) != nrow(s)) {
         stop("Wrong number of elements in rho")
     }
     if (length(rho) == 1 & rho == 0) {
         cat("With rho=0, there may be convergence problems if the  
input matrix s\n is not of full rank",
             fill = TRUE)
     }
     if (is.vector(rho)) {
         rho = matrix(sqrt(rho)) %*% sqrt(rho)
     }
     if (length(rho) == 1) {
         rho = matrix(rho, ncol = n, nrow = n)
     }
     if (!is.null(zero)) {
         if (!is.matrix(zero)) {
             zero = matrix(zero, nrow = TRUE)
         }
         for (k in 1:nrow(zero)) {
             i = zero[k, 1]
             j = zero[k, 2]
             rho[i, j] = BIG
             rho[j, i] = BIG
         }
     }
     start.type = match.arg(start)
     if (start.type == "cold") {
         is = 0
         ww = xx = matrix(0, nrow = n, ncol = n)
     }
     if (start.type == "warm") {
         is = 1
         if (is.null(w.init) | is.null(wi.init)) {
             stop("Warm start specified: w.init and wi.init must be
non-null")
         }
         ww = w.init
         xx = wi.init
     }
     itrace = 1 * trace
     ipen = 1 * (penalize.diagonal)
     ia = 1 * approx
     mode(rho) = "double"
     mode(s) = "double"
     mode(ww) = "double"
     mode(xx) = "double"
     mode(n) = "integer"
     mode(maxit) = "integer"
     mode(ia) = "integer"
     mode(itrace) = "integer"
     mode(ipen) = "integer"
     mode(is) = "integer"
     mode(thr) = "double"
     junk <- .Fortran("glasso", n, s, rho, ia, is, itrace, ipen,
         thr, maxit = maxit, ww = ww, xx = xx, niter = integer(1),
         del = double(1), ierr = integer(1), PACKAGE = "glasso")
     ww = matrix(junk$ww, ncol = n)
     xx = matrix(junk$xx, ncol = n)
     if (junk$ierr != 0) {
         stop("memory allocation error")
     }
     critfun = function(Sigmahati, s, rho, penalize.diagonal = TRUE) {
         d = det(Sigmahati)
         temp = Sigmahati
         if (!penalize.diagonal) {
             diag(temp) = 0
         }
         val = -log(d) + sum(diag(s %*% Sigmahati)) + sum(abs(rho *
             temp))
         return(val)
     }
     crit = critfun(xx, s, rho, penalize.diagonal)
     return(list(w = ww, wi = xx, loglik = -(n/2) * crit, errflag = junk$ierr,
         approx = approx, del = junk$del, niter = junk$niter))
}



If anyone can give any suggestions, I'd greatly appreciate it. Thanks!




Vivian

Hi Vivian,

This may be naive given the method (I am unfamiliar with glasso), but
what about simple subtraction?  If it restricts to 0, you believe you
have .3, then just: obs - .3 and restrict to 0 again?  Here is a
little example (assuming .3 correlation, but using glasso with the
covariance matrix so subtraction is marginally more complex):

dat <- lapply(list((dat <- prcomp(matrix(rnorm(1000), 200))$x), dat %*%
 chol(matrix(c(1, rep(c(rep(.3, 5), 1), 4)), 5, 5))), cov)

i <- sqrt(1/diag(dat[[2]]))
v <- dat[[2]] - matrix(c(0, rep(c(rep(.3, 5), 0), 4)), 5, 5)/(i *
rep(i, each = 5))

require(glasso)
Reduce(`-`, res <- list(true = glasso(dat[[1]], .01)$w, fake glasso(v,
.01)$w))

Again this may not make any sense at all given how the methods actually work.

Cheers,

Josh

On Tue, Oct 18, 2011 at 7:24 PM, Vivian Shih <vivs at ucla.edu>
wrote:
> I've only been using R on and off for 9 months and started using the
glasso
> package for sparse covariance estimation. I know the concept is to shrink
> some of the elements of the covariance matrix to zero.
>
> However, say I have a dataset that I know has some underlying
"baseline"
> covariance/correlation (say, a value of 0.3), how can I change or
> incorporate that into to the code? Basically instead of "let's
restrict some
> of these elements are 0 and others have different nonzero values",
I'd like
> it to be "let's restrict some of these elements are 0.3 and others
are
> different or higher/lower than this 0.3". Does that make any sense?
>
>
>
> Here is the glasso code if it helps:
>>
>> library(glasso)
>> glasso
>
> function (s, rho, zero = NULL, thr = 1e-04, maxit = 10000, approx = FALSE,
> ? ?penalize.diagonal = TRUE, start = c("cold", "warm"),
w.init = NULL,
> ? ?wi.init = NULL, trace = FALSE)
> {
> ? ?n = nrow(s)
> ? ?BIG = 1e+10
> ? ?if (!is.matrix(rho) & length(rho) != 1 & length(rho) != nrow(s))
{
> ? ? ? ?stop("Wrong number of elements in rho")
> ? ?}
> ? ?if (length(rho) == 1 & rho == 0) {
> ? ? ? ?cat("With rho=0, there may be convergence problems if the input
> matrix s\n is not of full rank",
> ? ? ? ? ? ?fill = TRUE)
> ? ?}
> ? ?if (is.vector(rho)) {
> ? ? ? ?rho = matrix(sqrt(rho)) %*% sqrt(rho)
> ? ?}
> ? ?if (length(rho) == 1) {
> ? ? ? ?rho = matrix(rho, ncol = n, nrow = n)
> ? ?}
> ? ?if (!is.null(zero)) {
> ? ? ? ?if (!is.matrix(zero)) {
> ? ? ? ? ? ?zero = matrix(zero, nrow = TRUE)
> ? ? ? ?}
> ? ? ? ?for (k in 1:nrow(zero)) {
> ? ? ? ? ? ?i = zero[k, 1]
> ? ? ? ? ? ?j = zero[k, 2]
> ? ? ? ? ? ?rho[i, j] = BIG
> ? ? ? ? ? ?rho[j, i] = BIG
> ? ? ? ?}
> ? ?}
> ? ?start.type = match.arg(start)
> ? ?if (start.type == "cold") {
> ? ? ? ?is = 0
> ? ? ? ?ww = xx = matrix(0, nrow = n, ncol = n)
> ? ?}
> ? ?if (start.type == "warm") {
> ? ? ? ?is = 1
> ? ? ? ?if (is.null(w.init) | is.null(wi.init)) {
> ? ? ? ? ? ?stop("Warm start specified: w.init and wi.init must be
non-null")
> ? ? ? ?}
> ? ? ? ?ww = w.init
> ? ? ? ?xx = wi.init
> ? ?}
> ? ?itrace = 1 * trace
> ? ?ipen = 1 * (penalize.diagonal)
> ? ?ia = 1 * approx
> ? ?mode(rho) = "double"
> ? ?mode(s) = "double"
> ? ?mode(ww) = "double"
> ? ?mode(xx) = "double"
> ? ?mode(n) = "integer"
> ? ?mode(maxit) = "integer"
> ? ?mode(ia) = "integer"
> ? ?mode(itrace) = "integer"
> ? ?mode(ipen) = "integer"
> ? ?mode(is) = "integer"
> ? ?mode(thr) = "double"
> ? ?junk <- .Fortran("glasso", n, s, rho, ia, is, itrace, ipen,
> ? ? ? ?thr, maxit = maxit, ww = ww, xx = xx, niter = integer(1),
> ? ? ? ?del = double(1), ierr = integer(1), PACKAGE = "glasso")
> ? ?ww = matrix(junk$ww, ncol = n)
> ? ?xx = matrix(junk$xx, ncol = n)
> ? ?if (junk$ierr != 0) {
> ? ? ? ?stop("memory allocation error")
> ? ?}
> ? ?critfun = function(Sigmahati, s, rho, penalize.diagonal = TRUE) {
> ? ? ? ?d = det(Sigmahati)
> ? ? ? ?temp = Sigmahati
> ? ? ? ?if (!penalize.diagonal) {
> ? ? ? ? ? ?diag(temp) = 0
> ? ? ? ?}
> ? ? ? ?val = -log(d) + sum(diag(s %*% Sigmahati)) + sum(abs(rho *
> ? ? ? ? ? ?temp))
> ? ? ? ?return(val)
> ? ?}
> ? ?crit = critfun(xx, s, rho, penalize.diagonal)
> ? ?return(list(w = ww, wi = xx, loglik = -(n/2) * crit, errflag =
junk$ierr,
> ? ? ? ?approx = approx, del = junk$del, niter = junk$niter))
> }
>
>
>
> If anyone can give any suggestions, I'd greatly appreciate it. Thanks!
>
>
>
>
> Vivian
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, ATS Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/