[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