2

I have a code that generates a density scatter plot like this (found here How to add an ellipse to a scatter plot colored by density in R? ) :

## Data in a data.frame
x1 <- rnorm(n=1000, sd=2)
x2 <- x1*1.2 + rnorm(n=1000, sd=2)
df <- data.frame(x1,x2)

## Use densCols() output to get density at each point
x <- densCols(x1,x2, colramp=colorRampPalette(c("black", "white")))
df$dens <- col2rgb(x)[1,] + 1L

## Map densities to colors
cols <-  colorRampPalette(c("#FF3100", "#FF9400", "#FCFF00", 
                            "#45FE4F", "#00FEFF", "#000099"))(6)
df$col <- ifelse(df$dens >= 250, cols[1], ifelse(df$dens >= 200, cols[2], ifelse(df$dens >= 150, cols[3], ifelse(df$dens >= 100, cols[4], ifelse(df$dens >= 50, cols[5], cols[6])))))

## Plot it, reordering rows so that densest points are plotted on top
plot(x2~x1, data=df[order(df$dens),], pch=20, col=col, cex=1)

enter image description here

I'd like to add contours exactly like in the image below (found here https://fr.mathworks.com/matlabcentral/fileexchange/8577-scatplot) :

enter image description here

I have to do this without ggplot2.

Does anyone have an idea please?

Matt 38
  • 57
  • 1
  • 8

1 Answers1

4

These contour lines show the estimated density of the data at given levels. You may compute the density estimate using MASS::kde2d() and then plot using contour():

biv_kde <- MASS::kde2d(x1, x2)
contour(biv_kde, add = T)

enter image description here

Edit:

Here's another approach to get colored density contours using KernSmooth:

(Note that I used set.seed(123) before generating the data.)

# estimate bandwidths
h <- c(dpik(df$x1), dpik(df$x2))

# obtain density estimatte
f1 <- bkde2D(df[, 1:2], bandwidth = h, gridsize = c(200, 200))

# setup vector of density levels to obtain contours for
contour_levels <- pretty(f1$fhat, 9)

# setup color palette
crp <- colorRampPalette(rev(c("#FF3100", "#FF9400", "#FCFF00", 
                              "#45FE4F", "#00FEFF", "#000099")))

# density based colors 
df$col <- densCols(x1, x2, bandwidth = h, colramp = crp)

# plot
plot(x2 ~ x1, data = df, pch = 20, col = col, cex = 0.5)
contour(f1$x1, f1$x2, f1$fhat, levels = contour_levels, col = crp(9),
        lwd = 2, add = T)

enter image description here

Martin C. Arnold
  • 9,483
  • 1
  • 14
  • 22
  • This is a good start, thank you! However, it does not work in all cases : if I have something like this : https://i.imgur.com/HPNvc09.png , I get something like that : https://i.imgur.com/mMnqfbT.png which does not correspond to the colors. Could this issue be fixed ? – Matt 38 Apr 13 '21 at 09:45
  • Hard to tell without data but I think this is linked to the bandwidth used in estimating the density... note that `densCol()` and `kde2d()` have a bandwidth argument (`bandwidth` / `h`) that determines the smoothness of the estimate. Could you add a reprex to your question? – Martin C. Arnold Apr 13 '21 at 11:29
  • Wouldn't it be because your code doesn't use `df$col`? Maybe I should add it somewhere in the code to match the contour lines to the colors? My real data looks exactly like `df`, with a column for the x-axis, another one for the y-axis and a third one for the colors. By reprex you mean share my real data and my code to reproduce the same plot as in my imgur? I'll try to play with the `bandwidth` argument though! – Matt 38 Apr 13 '21 at 15:33
  • Sorry, I didn't notice that the contours are also colored! I've edited my answer and added another example that yields a smoother-looking result. (: – Martin C. Arnold Apr 13 '21 at 21:01
  • It wasn't about the contour colors sorry, what I needed is to surround the color areas like https://i.imgur.com/GOV32zU.jpg, I thought it didn't work on my real data because you didn't use `df$col`. Your edit looks better on the test data though! I tried it on my real data but got `Error in contour.default(coordPoints$x, coordPoints$y, f1$fhat, levels = contour_levels, : increasing 'x' and 'y' values expected`. I sorted my data by x1, removed dups and NAs but still get this error. I don't understand why because my real data is really like the test data. Maybe I should share them here but how? – Matt 38 Apr 14 '21 at 13:02
  • Not sure about the error message. I think it would be best to open a new question – after all, this one asks how density contours can be added which I have answered in detail (: – Martin C. Arnold Apr 14 '21 at 13:27