psolymos / ResourceSelection

Resource Selection (Probability) Functions for Use-Availability Data in R
https://peter.solymos.org/ResourceSelection/
9 stars 4 forks source link

[kdepairs]: Overlay Density to onto Diagonal Histogramme #6

Closed KathyGCY closed 7 years ago

KathyGCY commented 7 years ago

Dear Mr Solymos, Thank you so much for this great package! For almost all the data sets I've ever worked on, I always try the kdepairs as the very first step. However recently I've ran into a bit of a trouble. I was trying to overlay the density (univariate) onto the histogrammes in the diagonal elements.

I've tried to change the source code but I'm a bit stuck. HELP!

P.S. Please reply either on github or to my email: cg2908@columbia.edu

#########################################################################
## Scatter Plot Matrix with 2d Contour
kdepairs.default <- function(x, n=25, density=TRUE, contour=TRUE, ...) {
  y <- data.frame(x)
  fun.lower <- function(x1, x2, ...) {
    if (is.factor(x1)) {
      x1 <- as.integer(x1)
    }
    if (is.factor(x2)) {
      x1 <- as.integer(x2)
    }
    OK <- length(unique(x1))>2 && length(unique(x2))>2
    if (!density && !contour)
      n <- 0
    if (n>0 && OK) {
      if (density || contour)
        d <- MASS::kde2d(x1, x2, n=n)
      if (density)
        image(d, col=terrain.colors(50), add=TRUE)
      if (contour)
        graphics::contour(d,add=TRUE)
    } else points(x1, x2)
  }
  fun.upper <- function(x1, x2, ...) {
    if (is.factor(x1)) {
      x1 <- as.integer(x1)
    }
    if (is.factor(x2)) {
      x1 <- as.integer(x2)
    }
    points(x1,x2, col="lightgrey")
    lines(lowess(x1,x2), col="#81D8D0", lty=1, lwd = 2)
    lines(abline(lm(x2~x1), col = "pink3",lty = 1, lwd = 2))
    COR <- cor(x1, x2)
    text(mean(range(x1,na.rm=TRUE)), mean(range(x2,na.rm=TRUE)),
         round(COR, 3), cex=1+abs(COR))
  }
  panel.hist <- function(x, ...) {
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
      #myden = density(x)
      #lines(myden)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col="#81D8D0", ...)

    box()
  }

  pairs.default(y, lower.panel=fun.lower, upper.panel=fun.upper, diag.panel=panel.histogram)
  invisible(NULL)
}
psolymos commented 7 years ago

Your function was implementing the right thing but I would add the lines on top of the bars drawn by rect like this:

f <-
function(x, n=25, density=TRUE, contour=TRUE, ...) {
    y <- data.frame(x)
    fun.lower <- function(x1, x2, ...) {
        if (is.factor(x1)) {
            x1 <- as.integer(x1)
        }
        if (is.factor(x2)) {
            x1 <- as.integer(x2)
        }
        OK <- length(unique(x1))>2 && length(unique(x2))>2
        if (!density && !contour)
            n <- 0
        if (n>0 && OK) {
            if (density || contour)
            d <- MASS::kde2d(x1, x2, n=n)
            if (density)
                image(d, col=terrain.colors(50), add=TRUE)
            if (contour)
                graphics::contour(d,add=TRUE)
        } else points(x1, x2)
    }
    fun.upper <- function(x1, x2, ...) {
        if (is.factor(x1)) {
            x1 <- as.integer(x1)
        }
        if (is.factor(x2)) {
            x1 <- as.integer(x2)
        }
        points(x1,x2, col="lightgrey")
        lines(lowess(x1,x2), col="darkgreen", lty=1)
        COR <- cor(x1, x2)
        text(mean(range(x1,na.rm=TRUE)), mean(range(x2,na.rm=TRUE)),
            round(COR, 3), cex=1+abs(COR))
    }
    panel.hist <- function(x, ...) {
        usr <- par("usr"); on.exit(par(usr))
        par(usr = c(usr[1:2], 0, 1.5) )
        h <- hist(x, plot=FALSE)
        breaks <- h$breaks; nB <- length(breaks)
        #y <- h$counts
        y <- h$density
        Max <- max(y)
        y <- y/Max
        rect(breaks[-nB], 0, breaks[-1], y, col="gold", border="orange", ...)
        den <- density(x)
        den$y <- den$y/Max
        lines(den)
        box()
    }
    pairs.default(y, lower.panel=fun.lower, upper.panel=fun.upper, diag.panel=panel.hist)
    invisible(NULL)
}

f(iris[1:4])

Added this to the development version by commit 0d813d79c0e465e5eca4196d4c225e312b2e929b .

KathyGCY commented 7 years ago

Thank you so much!