[R] Error: system is computationally singular: reciprocal conditionnumber

[Ravi Varadhan]

Your covariance matrix Szz is not positive definite.   It is singular.  The
following test that you are doing is neither necessary nor useful:

    zz.ev <- eigen(Szz)$values
    if(min(zz.ev)[1]<0){

        stop("\'Szz\' is not positive definite!\n")
    }


You may want to use Moore-Penrose inverse, also known as generalized inverse
or pseudoinverse to overcome this problem.  This approach uses
singular-value decomposition (SVD).  Take a look at the "corpcor" package.

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Moumita Das
Sent: Thursday, June 25, 2009 10:59 AM
To: r-help at r-project.org
Subject: [R] Error: system is computationally singular: reciprocal
conditionnumber

I get this error while computing partial correlation.


*Error in solve.default(Szz) :
  system is computationally singular: reciprocal condition number =
4.90109e-18*

Why is it?Can anyone give me some idea ,how do i get rid it it?

This is the function i use for calculating partial correlation.


pcor.mat <- function(x,y,z,method="p",na.rm=T){


    x <- c(x)
    y <- c(y)
    z <- as.data.frame(z)



    if(dim(z)[2] == 0){
        stop("There should be given data\n")
    }

    data <- data.frame(x,y,z)

    if(na.rm == T){
        data = na.omit(data)
    }

    xdata <- na.omit(data.frame(data[,c(1,2)]))
    Sxx <- cov(xdata,xdata,m=method)

    xzdata <- na.omit(data)
    xdata <- data.frame(xzdata[,c(1,2)])
    zdata <- data.frame(xzdata[,-c(1,2)])
    Sxz <- cov(xdata,zdata,m=method)

    zdata <- na.omit(data.frame(data[,-c(1,2)]))
    Szz <- cov(zdata,zdata,m=method)


    # is Szz positive definite?
    zz.ev <- eigen(Szz)$values
    if(min(zz.ev)[1]<0){

        stop("\'Szz\' is not positive definite!\n")
    }

    # partial correlation
    Sxx.z <- Sxx - Sxz %*% solve(Szz) %*% t(Sxz)

    print(Sxx.z) # this gets printed

    rxx.z <- cov2cor(Sxx.z)[1,2] #some problem in this function function (V)
    {
       print("cov2cor")
       p <- (d <- dim(V))[1]
        if (!is.numeric(V) || length(d) != 2L || p != d[2L])
            stop("'V' is not a square numeric matrix")
        Is <- sqrt(1/diag(V))
        if (any(!is.finite(Is)))
            warning("diag(.) had 0 or NA entries; non-finite result is
doubtful")
        r <- V
        r[] <- Is * V * rep(Is, each = p)
        r[cbind(1L:p, 1L:p)] <- 1
        r
    }
    return(rxx.z)
}

--
Thanks
Moumita

	[[alternative HTML version deleted]]

______________________________________________
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.