R:计算并绘制两个密度等高线之间的差异

Question

我有两个包含两个连续变量的数据集:duration 和 waiting。

library("MASS")
data(geyser)

geyser1 <- geyser[1:150,]

geyser2 <- geyser[151:299,]
geyser2$duration <- geyser2$duration - 1
geyser2$waiting <- geyser2$waiting - 20

对于每个数据集，我输出一个二维密度图

ggplot(geyser1, aes(x = duration, y = waiting)) +
  xlim(0.5, 6) + ylim(40, 110) +
  stat_density2d(aes(alpha=..level..),
                 geom="polygon", bins = 10)

ggplot(geyser2, aes(x = duration, y = waiting)) +
  xlim(0.5, 6) + ylim(40, 110) +
  stat_density2d(aes(alpha=..level..),
                 geom="polygon", bins = 10)

我现在想生成一个图，指示两个图具有相同密度(白色)、负差异(从白色到蓝色的渐变，其中 geyser2 比 geyser1 更密集) ) 和正差异(从白色到红色的渐变，其中 geyser1 比 geyser2 更密集)。

如何计算和绘制密度的差异？

Answer 1

您可以先使用 kde2d 计算密度，然后将它们相互相减。然后，您进行一些数据 reshape ，将其转化为可以提供给 ggplot2 的形式。

library(reshape2) # For melt function

# Calculate the common x and y range for geyser1 and geyser2
xrng = range(c(geyser1$duration, geyser2$duration))
yrng = range(c(geyser1$waiting, geyser2$waiting))

# Calculate the 2d density estimate over the common range
d1 = kde2d(geyser1$duration, geyser1$waiting, lims=c(xrng, yrng), n=200)
d2 = kde2d(geyser2$duration, geyser2$waiting, lims=c(xrng, yrng), n=200)

# Confirm that the grid points for each density estimate are identical
identical(d1$x, d2$x) # TRUE
identical(d1$y, d2$y) # TRUE

# Calculate the difference between the 2d density estimates
diff12 = d1 
diff12$z = d2$z - d1$z

## Melt data into long format
# First, add row and column names (x and y grid values) to the z-value matrix
rownames(diff12$z) = diff12$x
colnames(diff12$z) = diff12$y

# Now melt it to long format
diff12.m = melt(diff12$z, id.var=rownames(diff12))
names(diff12.m) = c("Duration","Waiting","z")

# Plot difference between geyser2 and geyser1 density
ggplot(diff12.m, aes(Duration, Waiting, z=z, fill=z)) +
  geom_tile() +
  stat_contour(aes(colour=..level..), binwidth=0.001) +
  scale_fill_gradient2(low="red",mid="white", high="blue", midpoint=0) +
  scale_colour_gradient2(low=muted("red"), mid="white", high=muted("blue"), midpoint=0) +
  coord_cartesian(xlim=xrng, ylim=yrng) +
  guides(colour=FALSE)