Skip to content

partial correlations

3 messages · Ragnar Beer, Ko-Kang Kevin Wang, Tim Keitt

Original
Ragnar Beer
# 
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ko-Kang Kevin Wang
# 
I got this from Associate Professor Brian McArdle's multivariate statistics course.  However I have not tested it yet...



parcor.test <- function( x, y, q=1, alternative = "two.sided" ) {
    if( is.matrix( x ) && ncol( x ) >= 2 ) {
        y <- x[, 2]
        x <- x[, 1]
    }
    if( length( x ) != length( y ) )
        stop( "x and y should be the same length " )
    n <- length( x )
    if( n <= 2 )
        stop( "x and y should effectively be longer than 2" )


# Pearson's product moment correlation:
        coef <- cor( x, y )
        if( is.na( coef ) )
            stop( "too small variance" )
        if( abs( coef ) > 0.999999 )
            stop( paste( "correlation is", coef ) )
        attr( coef, "names" ) <- "cor"
        stat <- coef * sqrt( ( n - 2-q )/( 1 - coef^2 ) )  # t-statistic
        attr( stat, "names" ) <- "t"

    p.value <- if( names( stat ) == "t" ) {
        switch( alternative,
                greater = 1 - pt( stat, n - 2 - q ),
                less = pt( stat, n - 2 - q ),
                two.sided = 2 * pt(  - abs( stat ), n - 2 - q ),
              )
    }

    attr( p.value, "Names" ) <- NULL    #a kludge for nicer output
    if( names( stat ) == "t" ) {
        pars <- n - 2 - q
        attr( pars, "names" ) <- "df"
    }
    else pars <- NULL
    null.value <- 0
    attr( null.value, "names" ) <- c( "coef" )
    z <- list( statistic = stat, parameters = pars, p.value =
        p.value, estimate = coef, null.value = null.value,
        alternative = alternative )
    z$data.name <- paste( deparse( substitute( x ) ), "and", deparse(
        substitute( y ) ) )
    attr( z, "class" ) <- "htest"
    return( z )
}






-------------------------------------------------------------------------------------------- 
Ko-Kang Kevin Wang
Head of Statistical Analysis Division
Software Developers' Klub (SDK)
University of Auckland
New Zealand

-----Original Message-----
From: owner-r-help at stat.math.ethz.ch [mailto:owner-r-help at stat.math.ethz.ch]On Behalf Of Ragnar Beer
Sent: Wednesday, August 01, 2001 9:46 PM
To: r-help at stat.math.ethz.ch
Subject: [R] partial correlations


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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Tim 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: