[R] Summary: Partial correlation coefficients in R. Thanks everybody!
Kaspar Pflugshaupt
pflugshaupt at geobot.umnw.ethz.ch
Fri Feb 25 17:01:21 CET 2000
Hello all,
here's a collection of answers I got on my question concerning partial
correlation coefficients:
Some people gave a simple formula for the three-variable-case, as did Dave
Lucy:
pcor <- function(v1, v2, v3)
{
c12 <- cor(v1, v2)
c23 <- cor(v2, v3)
c13 <- cor(v1, v3)
partial <- (c12-(c13*c23))/(sqrt(1-(c13^2)) * sqrt(1-(c23^2)))
return(partial)
}
The general algorithm was given by John Logsdon, among others:
The partial correlation coefficients are the negative
scaled off-diagonal inverse variance. So compute the variance-covariance
matrix, invert, scale to the diagonal and negate and you have it.
And an implementation by Martyn Plummer (here, too, I received several):
pcor2 <- function(x){
conc <- solve(var(x))
resid.sd <- 1/sqrt(diag(conc))
pcc <- - sweep(sweep(conc, 1, resid.sd, "*"), 2, resid.sd, "*")
return(pcc)
}
This is the version I'm using now, together with a test for significance of
each coefficient (H0: coeff=0):
f.parcor <-
function (x, test = F, p = 0.05)
{
nvar <- ncol(x)
ndata <- nrow(x)
conc <- solve(cor(x))
resid.sd <- 1/sqrt(diag(conc))
pcc <- -sweep(sweep(conc, 1, resid.sd, "*"), 2, resid.sd,
"*")
colnames(pcc) <- rownames(pcc) <- colnames(x)
if (test) {
t.df <- ndata - nvar
t <- pcc/sqrt((1 - pcc^2)/t.df)
pcc <- list(coefs = pcc, significant = t > qt(1 - (p/2),
df = t.df))
}
return(pcc)
}
Thanks to everybody for your help!
Kaspar
--
Kaspar Pflugshaupt
Geobotanisches Institut
Zuerichbergstr. 38
CH-8044 Zuerich
Tel. ++41 1 632 43 19
Fax ++41 1 632 12 15
mailto:pflugshaupt at geobot.umnw.ethz.ch
privat:pflugshaupt at mails.ch
http://www.geobot.umnw.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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list