Hi,
after reading the archives I found some methods... adopted and
modified one of them to the following. I think it is correct after
checking and comparing the results with other software... but if
possible someone could have a look and spot any mistakes I would be
grateful. Thanks
pcor3 <- function (x, test = T, p = 0.05, alternative="two.sided")
{
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, sig = t > qt(1 - (p/2), df = t.df)); #
original statement
if (alternative == "two.sided") {
pcc <- list(coefs = pcc, sig = t > qt(1 - (p/2), df = t.df),
p.value = 2 * pmin(pt(t, t.df), 1-pt(t, t.df))) # two.sided
} else if (alternative == "greater") {
pcc <- list(coefs = pcc, sig = t > qt(1 - p, df = t.df), p.value
= 1-pt(t, t.df)) # greater
} else if (alternative == "less") {
pcc <- list(coefs = pcc, sig = t > qt(1 - p, df = t.df), p.value
= 2*(1-pt(t, t.df)))
}
}
str <- sprintf("Partial correlation for:"); print(str,
quote=FALSE);
str <- sprintf("%s", colnames(x)); print(str, quote=FALSE);
str <- sprintf("p: %.2f, alternative: %s", p, alternative);
print(str, quote=FALSE);
if (test) {
str <- sprintf("df: %d", t.df); print(str, quote=FALSE);
}
return(pcc)
}
The function was adopted from the following email:
http://tolstoy.newcastle.edu.au/R/help/00a/0518.html
--
yianni
