R help - Jun 2009 - ERROR: system is computationally singular: reciprocal condition number = 4.90109e-18

Hi All,
This is my R-version information:---
> version
               _
platform       i486-pc-linux-gnu
arch           i486
os             linux-gnu
system         i486, linux-gnu
status
major          2
minor          7.1
year           2008
month          06
day            23
svn rev        45970
language       R
version.string R version 2.7.1 (2008-06-23)

While calculating partial correlation for a dataset ,i keep getting this
error :--
*Error in solve.default(Szz) :
  system is computationally singular: reciprocal condition number 4.90109e-18*



On using the traceback() function i get
this:--------------
> traceback()
10: *.Call("La_dgesv", a, b, tol, PACKAGE = "base")*
9: solve.default(Szz)
8: solve(Szz)
7: pcor.mat(firstvalue, secondvalue, third_var, method, na.rm = T)
6: PartialCorr_Calculation(value1, value2, third_var, method =
"pearson",
       na.rm = T)
5: Partial(contrld_third_var(rowvalues$matrix1, rowvalues$matrix2,
       x <- stringOfItemCategoryToDataFrameOfItemCategory,
item_category_table,
       pcor_type <- "ic"), data1, data2, pcor_thirdvar_type <-
"all")
4: main()
3: eval.with.vis(expr, envir, enclos)
2: eval.with.vis(ei, envir)
1: source("correlationFP.R")

In pcor.mat function :--
# By using var-cov matrix
pcor.mat <- function(x,y,z,method="p",na.rm=T){

    #print("pcor.mat")
    x <- c(x)
    y <- c(y)
    z <- as.data.frame(z)

    m<- identical(all.equal(x, y), TRUE)
    #print(m)
    if(  m=="TRUE")
    {
        print("inside equality-------------------------------")
        return(1)
    }


    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) # this gets printed




    rxx.z <- cov2cor(Sxx.z)[1,2] # probably error in this function cov2cor




    return(rxx.z)
}

for some other datasets this pcor.mat function works fine.
The solve.default function is this:---------------
solve.default<-function(a, b, tol = ifelse(LINPACK, 1e-07,
.Machine$double.eps),LINPACK = FALSE, ...)
{

    if (is.complex(a) || (!missing(b) && is.complex(b))) {
        a <- as.matrix(a)
        if (missing(b)) {
            if (nrow(a) != ncol(a))
                stop("only square matrices can be inverted")
            b <- diag(1 + (0+0i), nrow(a))
            colnames(b) <- rownames(a)
        }
        else if (!is.complex(b))
            b[] <- as.complex(b)
        if (!is.complex(a))
            a[] <- as.complex(a)
        return(if (is.matrix(b)) {
            if (ncol(a) != nrow(b)) stop("'b' must be compatible
with 'a'")
            rownames(b) <- colnames(a)
            .Call("La_zgesv", a, b, PACKAGE = "base")
        } else drop(.Call("La_zgesv", a, as.matrix(b), PACKAGE =
"base")))
    }
    if (is.qr(a)) {
        warning("solve.default called with a \"qr\" object: use
'qr.solve'")
        return(solve.qr(a, b, tol))
    }
    if (!LINPACK) {
        a <- as.matrix(a)
        if (missing(b)) {
            if (nrow(a) != ncol(a))
                stop("only square matrices can be inverted")
            b <- diag(1, nrow(a))
            colnames(b) <- rownames(a)
        }
        else storage.mode(b) <- "double"
        storage.mode(a) <- "double"
        return(if (is.matrix(b)) {
            if (ncol(a) != nrow(b)) stop("'b' must be compatible
with 'a'")
            rownames(b) <- colnames(a)
            .Call("La_dgesv", a, b, tol, PACKAGE = "base")
        } else drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE
"base")))
    }
    a <- qr(a, tol = tol)
    nc <- ncol(a$qr)
    if (a$rank != nc)
        stop("singular matrix 'a' in 'solve'")
    if (missing(b)) {
        if (nc != nrow(a$qr))
            stop("only square matrices can be inverted")
        b <- diag(1, nc)
        colnames(b) <- rownames(a$qr)
    }
    qr.coef(a, b)
}

So what has to be done in the "base" package to get rid of this error.


-- 
Thanks
Moumita

	[[alternative HTML version deleted]]

Hi,

You are being rude to those who are trying to help you by ignoring their
suggestions and by sending useless, irrelevant information.

Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Moumita Das <das.moumita.online at gmail.com>
Date: Sunday, June 28, 2009 11:20 am
Subject: [R] ERROR: system is computationally singular: reciprocal condition
number = 4.90109e-18
To: r-help at r-project.org

> Hi All,
>  This is my R-version information:---
>  > version
>                 _
>  platform       i486-pc-linux-gnu
>  arch           i486
>  os             linux-gnu
>  system         i486, linux-gnu
>  status
>  major          2
>  minor          7.1
>  year           2008
>  month          06
>  day            23
>  svn rev        45970
>  language       R
>  version.string R version 2.7.1 (2008-06-23)
>  
>  While calculating partial correlation for a dataset ,i keep getting this
>  error :--
>  *Error in solve.default(Szz) :
>    system is computationally singular: reciprocal condition number > 
4.90109e-18*
>  
>  
>  
>  On using the traceback() function i get this:--------------
>  > traceback()
>  10: *.Call("La_dgesv", a, b, tol, PACKAGE = "base")*
>  9: solve.default(Szz)
>  8: solve(Szz)
>  7: pcor.mat(firstvalue, secondvalue, third_var, method, na.rm = T)
>  6: PartialCorr_Calculation(value1, value2, third_var, method =
"pearson",
>         na.rm = T)
>  5: Partial(contrld_third_var(rowvalues$matrix1, rowvalues$matrix2,
>         x <- stringOfItemCategoryToDataFrameOfItemCategory,
>  item_category_table,
>         pcor_type <- "ic"), data1, data2, pcor_thirdvar_type
<- "all")
>  4: main()
>  3: eval.with.vis(expr, envir, enclos)
>  2: eval.with.vis(ei, envir)
>  1: source("correlationFP.R")
>  
>  In pcor.mat function :--
>  # By using var-cov matrix
>  pcor.mat <- function(x,y,z,method="p",na.rm=T){
>  
>      #print("pcor.mat")
>      x <- c(x)
>      y <- c(y)
>      z <- as.data.frame(z)
>  
>      m<- identical(all.equal(x, y), TRUE)
>      #print(m)
>      if(  m=="TRUE")
>      {
>          print("inside equality-------------------------------")
>          return(1)
>      }
>  
>  
>      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) # this gets printed
>  
>  
>  
>  
>      rxx.z <- cov2cor(Sxx.z)[1,2] # probably error in this function
cov2cor
>  
>  
>  
>  
>      return(rxx.z)
>  }
>  
>  for some other datasets this pcor.mat function works fine.
>  The solve.default function is this:---------------
>  solve.default<-function(a, b, tol = ifelse(LINPACK, 1e-07,
>  .Machine$double.eps),LINPACK = FALSE, ...)
>  {
>  
>      if (is.complex(a) || (!missing(b) && is.complex(b))) {
>          a <- as.matrix(a)
>          if (missing(b)) {
>              if (nrow(a) != ncol(a))
>                  stop("only square matrices can be inverted")
>              b <- diag(1 + (0+0i), nrow(a))
>              colnames(b) <- rownames(a)
>          }
>          else if (!is.complex(b))
>              b[] <- as.complex(b)
>          if (!is.complex(a))
>              a[] <- as.complex(a)
>          return(if (is.matrix(b)) {
>              if (ncol(a) != nrow(b)) stop("'b' must be
compatible with
> 'a'")
>              rownames(b) <- colnames(a)
>              .Call("La_zgesv", a, b, PACKAGE = "base")
>          } else drop(.Call("La_zgesv", a, as.matrix(b), PACKAGE =
"base")))
>      }
>      if (is.qr(a)) {
>          warning("solve.default called with a \"qr\" object:
use 'qr.solve'")
>          return(solve.qr(a, b, tol))
>      }
>      if (!LINPACK) {
>          a <- as.matrix(a)
>          if (missing(b)) {
>              if (nrow(a) != ncol(a))
>                  stop("only square matrices can be inverted")
>              b <- diag(1, nrow(a))
>              colnames(b) <- rownames(a)
>          }
>          else storage.mode(b) <- "double"
>          storage.mode(a) <- "double"
>          return(if (is.matrix(b)) {
>              if (ncol(a) != nrow(b)) stop("'b' must be
compatible with
> 'a'")
>              rownames(b) <- colnames(a)
>              .Call("La_dgesv", a, b, tol, PACKAGE =
"base")
>          } else drop(.Call("La_dgesv", a, as.matrix(b), tol,
PACKAGE >  "base")))
>      }
>      a <- qr(a, tol = tol)
>      nc <- ncol(a$qr)
>      if (a$rank != nc)
>          stop("singular matrix 'a' in 'solve'")
>      if (missing(b)) {
>          if (nc != nrow(a$qr))
>              stop("only square matrices can be inverted")
>          b <- diag(1, nc)
>          colnames(b) <- rownames(a$qr)
>      }
>      qr.coef(a, b)
>  }
>  
>  So what has to be done in the "base" package to get rid of this
error.
>  
>  
>  -- 
>  Thanks
>  Moumita
>  
>  	[[alternative HTML version deleted]]
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.