Bill Anderegg
2014-May-13 21:33 UTC
[R] Correlelogram with partial correlation coefficients
Hi, I am trying to construct a correlelogram plot but plot partial correlation
coefficients, rather than normal coefficients. I've been using the
corrgram() and ggm() libraries to do correlelograms and partial correlation
analysis respectively, but I can't figure out how to combine them. Here is
my code (my dataset is "climvar4" and is a matrix with 11 climate
variables with 51
observations per variable. so in each cell of the correlelogram I want to
display the partial correlation between variable x and variable y, while
accounting for all
other variables):
library(ggm)
library(corrgram)
panel.shadeNtextP <- function (x, y, corr = NULL, col.regions, ...)
{
corr <- pcor(x, y, var(climvar4))
results <- pcor.test(corr,10,51)
est <- results$p.value
stars <- ifelse(est < 5e-4, "***",
ifelse(est < 5e-3, "**",
ifelse(est < 5e-2, "*", "")))
ncol <- 14
pal <- col.regions(ncol)
col.ind <- as.numeric(cut(corr, breaks = seq(from = -1, to = 1,
length = ncol + 1),
include.lowest = TRUE))
usr <- par("usr")
rect(usr[1], usr[3], usr[2], usr[4], col = pal[col.ind],
border = NA)
box(col = "lightgray")
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- formatC(corr, digits = 2, format = "f")
cex.cor <- .8/strwidth("-X.xx")
fonts <- ifelse(stars != "", 2,1)
# option 1: stars:
text(0.5, 0.4, paste0(r,"\n", stars), cex = cex.cor)
# option 2: bolding:
#text(0.5, 0.5, r, cex = cex.cor, font=fonts)
}
corrgram(climvar4, type="data", order=F,
lower.panel=panel.shadeNtextP,
upper.panel=NULL, cex=1.3)
Thanks so much in advance!
Bill
Kevin Wright
2014-May-14 19:10 UTC
[R] Correlelogram with partial correlation coefficients
Bill,
For package-specific questions, you should email the package maintainer
(for corrgram, that's me).
Michael Friendly's paper on corrgrams includes an example of partial
correlations. The code to reproduce his example is hiding in the
"tests"
directory of the corrgram sourcecode zip file. Here it is below.
Kevin
# Figure 8
require(Matrix) # For block diagonal function
partial <- function(r, xvar){
# r is a correlation matrix
# Calculate partial correlation of y|x
yvar <- setdiff(colnames(r), xvar)
ri <- r[yvar,yvar] - r[yvar,xvar] %*% solve(r[xvar,xvar]) %*% r[xvar,yvar]
s <- diag(ri)
s <- diag(sqrt(1/s))
ri <- s %*% ri %*% s
ri <- as.matrix(bdiag(ri, r[xvar, xvar]))
diag(ri) <- 1 # Should already be 1, but could be 1 + epsilon
colnames(ri) <- rownames(ri) <- c(yvar, xvar)
return(ri)
}
vars8a <- c("Gratio", "Rep78", "Rep77",
"Hroom", "Trunk", "Rseat",
"Length", "Weight", "Displa",
"Turn")
vars8b <- c("MPG", "Price")
vars8 <- c(vars8a, vars8b)
auto.cor <- cor(auto[, vars8], use="pair")
auto.par <- partial(auto.cor, vars8b)
corrgram(auto.par, lower.panel=panel.pie, upper.panel=panel.pie,
main="Auto data, partialing out Price,MPG")
On Tue, May 13, 2014 at 4:33 PM, Bill Anderegg
<anderegg@princeton.edu>wrote:
> Hi, I am trying to construct a correlelogram plot but plot partial
> correlation
> coefficients, rather than normal coefficients. I've been using the
> corrgram() and ggm() libraries to do correlelograms and partial correlation
> analysis respectively, but I can't figure out how to combine them. Here
is
> my code (my dataset is "climvar4" and is a matrix with 11 climate
> variables with 51
> observations per variable. so in each cell of the correlelogram I want to
> display the partial correlation between variable x and variable y, while
> accounting for all
> other variables):
>
> library(ggm)
> library(corrgram)
> panel.shadeNtextP <- function (x, y, corr = NULL, col.regions, ...)
> {
> corr <- pcor(x, y, var(climvar4))
> results <- pcor.test(corr,10,51)
> est <- results$p.value
> stars <- ifelse(est < 5e-4, "***",
> ifelse(est < 5e-3, "**",
> ifelse(est < 5e-2, "*",
"")))
> ncol <- 14
> pal <- col.regions(ncol)
> col.ind <- as.numeric(cut(corr, breaks = seq(from = -1, to = 1,
> length = ncol + 1),
> include.lowest = TRUE))
> usr <- par("usr")
> rect(usr[1], usr[3], usr[2], usr[4], col = pal[col.ind],
> border = NA)
> box(col = "lightgray")
> on.exit(par(usr))
> par(usr = c(0, 1, 0, 1))
> r <- formatC(corr, digits = 2, format = "f")
> cex.cor <- .8/strwidth("-X.xx")
> fonts <- ifelse(stars != "", 2,1)
> # option 1: stars:
> text(0.5, 0.4, paste0(r,"\n", stars), cex = cex.cor)
> # option 2: bolding:
> #text(0.5, 0.5, r, cex = cex.cor, font=fonts)
> }
> corrgram(climvar4, type="data", order=F,
lower.panel=panel.shadeNtextP,
> upper.panel=NULL, cex=1.3)
>
> Thanks so much in advance!
>
> Bill
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/
> posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Kevin Wright
[[alternative HTML version deleted]]