Emy
Emy

Reputation: 977

Differences between plotting contour() function in base R and using geom_contour() or stat_contour() in ggplot2

I have plotted a density function in base R and I would like to replicate the plot in ggplot2.

This is the plot in base R:

library(tidyverse)
library(mvtnorm)
sd <- 1 / 2
# sigma
s1 <- sd^2

# first two vectors
x.points <- seq(-3, 3, length.out = 100)
y.points <- seq(-3, 3, length.out = 100)

# the third vector is a density
z <- matrix(0, nrow = 100, ncol = 100)

mu1 <- c(0, 0)
sigma1 <- matrix(c(s1^2, 0, 0, s1^2), nrow = 2)
for (i in 1:100) {
  for (j in 1:100) {
    z[i, j] <- dmvnorm(c(x.points[i], y.points[j]),
      mean = mu1, sigma = sigma1
    )
  }
}


contour(x.points, y.points, z, xlim = range(-3, 3), ylim = c(-3, 3), nlevels = 5, drawlabels = TRUE)

To obtain the same result in ggplot2, I am following this example:

library(ggplot2)
library(reshape2) # for melt
volcano3d <- melt(volcano)
names(volcano3d) <- c("x", "y", "z")
# Basic plot
v <- ggplot(volcano3d, aes(x, y, z = z))
v + stat_contour()

But in my case vector z has a different length than x.points and y.points. From the errors I get below, it looks like the three vectors should have the same length. How can I transform the dataset presented above so that it can be run through ggplot2?

data1 <- as.data.frame(cbind(x.points, y.points))

p <- ggplot(data = data1, mapping = aes(x.points, y.points, z=z))
p + geom_contour()

#> Error: Aesthetics must be either length 1 or the same as the data (100): z

p + stat_contour()

#> Error: Aesthetics must be either length 1 or the same as the data (100): z

p + stat_function(fun = contour) + xlim(-3,3)

#> Error: Aesthetics must be either length 1 or the same as the data (100): z

Created on 2021-04-08 by the reprex package (v0.3.0)

Upvotes: 1

Views: 293

Answers (1)

teunbrand
teunbrand

Reputation: 37953

The problem is likely that your data isn't in long format: for every value of the z matrix, you need the x and y position, which is different from the base R approach, wherein you just need these positions for every row/column.

We can transform the matrix z to a long format using reshape2::melt and then grab the correct positions from your vectors.

library(tidyverse)
library(mvtnorm)
sd <- 1 / 2
# sigma
s1 <- sd^2

# first two vectors
x.points <- seq(-3, 3, length.out = 100)
y.points <- seq(-3, 3, length.out = 100)

# the third vector is a density
z <- matrix(0, nrow = 100, ncol = 100)

mu1 <- c(0, 0)
sigma1 <- matrix(c(s1^2, 0, 0, s1^2), nrow = 2)
for (i in 1:100) {
  for (j in 1:100) {
    z[i, j] <- dmvnorm(c(x.points[i], y.points[j]),
                       mean = mu1, sigma = sigma1
    )
  }
}

# Here be the reshaping bit
df <- reshape2::melt(z)
df <- transform(
  df,
  x = x.points[Var1],
  y = y.points[Var2]
)

ggplot(df, aes(x, y)) +
  geom_contour(aes(z = value))

Created on 2021-04-08 by the reprex package (v1.0.0)

Upvotes: 1

Related Questions