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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
There may well be a simple way to do partial correlations, but I can't find it either. Here are three ways: 1. Use the formula. For r(xy,z), xy <- cor(x,y) xz <- cor(x,z) yz <- cor(y,z) (xy - xz*yz)/sqrt((1-xz^2)*(1-yz^2) 2. Correlate residuals. cor(lm(x~z)$resid,lm(y~z)$resid) The second method has the advantage that you can replace z with many variables and partial them all. 3. From the coefficients of regressions. a <- summary(lm(y ~ x + z))$coef[2,1] b <- summary(lm(x ~ y + z))$coef[2,1] sqrt(a*b) The significance of the partial correlation should be the same as the significance of the regression coefficient, e.g., summary(lm(y ~ x + z)) or summary(lm(x ~ y + z)) and look at the first coefficient. Jon Baron -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._