3D Density/Relative Frequency Histogram

Question

I've been tasked with studying the duration and delay of a doctor's appointment book.

My approach so far has been to plot and analyze density charts by duration and by deviation of the appointed time, separately.

Appointment ID	Duration	Deviation	Online/In person	Type of Punctuality/Delay
1	5 min	- 10 min (late)	Online	Patient is Late
2	2 min	- 5 min	In person	Doctor is late
3	10 min	+5 min (ahead of time)	In person	On time
...	...	...	...	...
1487	15 min	0.2 min	Online	On time

Now, I just want to print 3D bivariate kernel density plot, in my Markdown PDF OUTPUT

I can do a 2D density plot.

ggplot(final, aes(x=duration, y=deviation)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  theme_bw() +
  theme(legend.position='none') +
  xlim(0, quantile(as.numeric(final$duracion), 1-0.1)) +
  ylim(quantile(as.numeric(final$atraso), 0.01), quantile(as.numeric(final$atraso), 1-0.01))

However, I would like to extend the plot along the z-axis, in a manner similar to what the following link shows.

(https://www.researchgate.net/figure/Distribution-of-sequence-and-genotype-derived-allele-frequencies-r-067-in-the-SNP_fig1_5555304)

After using kde2d(), rgl() and plot_ly() allow me to do so, but I can't embed the plots in my PDF output. Not to mention that plot_ly messes with one of my axis. And worst of all, I lose my non-numerical variables, which allow me to split the data and do useful comparisons with facet_grid(), or aes(x, y, colour = ).

Is there any way I can plot the density of each pair of durations and deviations, in a 3d graph? I've een researching about wireframe() + outer(), but there does not seem to be any density function for outer()

Thanks beforehand.

Answer 1

You don't give us enough data or code to work with, so I'll show you how to embed a 3D plot of an unrelated function in a pdf document using rgl:

---
title: "Untitled"
date: "2023-02-16"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rgl)
setupKnitr(autoprint = TRUE)
```


```{r}
f <- function(x, y) { 
  r <- sqrt(x^2 + y^2)
  ifelse(r == 0, 10, 10 * sin(r)/r)
}
open3d()
plot3d(f, col = colorRampPalette(c("blue", "white", "red")), 
       xlab = "X", ylab = "Y", zlab = "Sinc( r )", 
       xlim = c(-10, 10), ylim = c(-10, 10),
       aspect = c(1, 1, 0.5))
```

This produces output that looks like this: