[R] partial correlations

[Timothy H. Keitt]

partial.correlation <- function(x, y=NULL, z=NULL, n=NULL,
                                keep.permutations=F) {
  if (!is.null(y)) x <- cbind(x, y)
  if (!is.null(z)) x <- cbind(x, z)
  c <- cor(x)
  out <- list()
  out$cxy <- c[1, 2]
  out$cxz <- c[1, 3]
  out$pcxy <- c[1, 2] - c[1, 3] * c[2, 3] /
    sqrt((1 - c[1, 3]^2) * (1 - c[2, 3]^2))
  out$pcxz <- c[1, 3] - c[1, 2] * c[2, 3] /
    sqrt((1 - c[1, 2]^2) * (1 - c[2, 3]^2))
  if (!is.null(n)) {
    perms <- matrix(nrow=n, ncol=4)
    for (i in 1:n) {
      x[,1] <- sample(x[,1])
      c <- partial.correlation(x)
      perms[i,] <- c(c$cxy, c$pcxy, c$cxz, c$pcxz)
    }
    out$p.cxy <- sum(perms[, 1] >= out$cxy) / n
    out$p.pcxy <- sum(perms[, 2] >= out$pcxy) / n
    out$p.cxz <- sum(perms[, 3] >= out$cxz) / n
    out$p.pcxz <- sum(perms[, 4] >= out$pcxz) / n
    if (keep.permutations) {
      out$cxy.perms <- perms[, 1]
      out$pcxy.perms <- perms[, 2]
      out$cxz.perms <- perms[, 3]
      out$pcxz.perms <- perms[, 4]
    }
  }
  return(out)
}


Ragnar Beer wrote:

> Howdy!
>
> I need to calculate partial correlations and I just can't find out how to
> do that with R. Can anybody help?
>
> Ragnar
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 
>
> r-help mailing list -- Read 
> http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ 
>


-- 
Timothy H. Keitt
Department of Ecology and Evolution
State University of New York at Stony Brook
Stony Brook, New York 11794 USA
Phone: 631-632-1101, FAX: 631-632-7626
http://life.bio.sunysb.edu/ee/keitt/



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._