How to add density plots for x and z axes in a tensor smooth plot?

I created this plot using tensor smooths..

I would like to add a visualization of the data density along the x-axis and z-axis. Ideally, this could be done with either:

A rug plot on the respective axes, or A density function plot for the data on the axes. I have searched extensively but couldn't find a solution.

This is my current code:

gam.model <- gam(death_or_neurodeficit1y ~ 
te(OtherBiomarker, laktat_0, bs = 'cr') + age + isMALE + Diabetes + trop_t_s + shock + byst_cpr + gfr_0,
  data = df,
  family = "binomial")

par(cex = 1)
vis.gam(gam.model,
        view = c("OtherBiomarker", "laktat_0"),
        plot.type = "persp",
        type = "response",
        theta = 320,
        phi = 25,
        color = "heat",
        xlab = "\n\nOtherBiomarker at admission",
        ylab = "\n\nLactate at admission",
        zlab = "\n\nPredicted Probability",
        border = "black",
        r = 100,
        ticktype = "detailed"
        )

I got the density function but could not manage to add it to the plot.

lactate_density <- density(gam.model$model$laktat_0)
OtherBiomarker_density <- density(gam.model$model$OtherBiomarker)

Here is a data sample:

structure(list(death_or_neurodeficit1y = c(0, 1, 1, 1, 0, 1), 
    age = c(53, 67, 76, 69, 81, 51), isMALE = c(1, 1, 1, 0, 1, 
    1), diabetes = c(0, 0, 0, 0, 0, 0), trop_t_s = c(296, 120, 
    24, 41, 71, 258), shock = c(1, 1, 1, 1, 1, 1), byst_cpr = c(1, 
    0, 1, 1, 1, 1), gfr_0 = c(107.1, 64.1, 64.9, 38, 54.2, 65.6
    ), OtherBiomarker = c(14562.291216349, 26308.5822361957, 
    15525.1522892513, 26273.8258989978, 23385.2426802908, 18708.4888976223
    ), laktat_0 = c(16.8, 55.5, 63.9, 40.1, 139.7, 27.6)), terms = death_or_neurodeficit1y ~ 
    age + isMALE + diabetes + trop_t_s + shock + byst_cpr + gfr_0 + 
        OtherBiomarker + laktat_0, na.action = structure(c(`26` = 26L, 
`33` = 33L, `43` = 43L, `46` = 46L, `57` = 57L, `70` = 70L, `91` = 91L, 
`93` = 93L, `105` = 105L, `115` = 115L, `116` = 116L, `118` = 118L, 
`121` = 121L, `127` = 127L, `135` = 135L, `160` = 160L), class = "omit"), row.names = c(NA, 
6L), class = "data.frame")

Thank you in advance.

This might not be exactly what you're looking for, but it might get you on the way. You can use trans3d() to map 3d points into the 2d window on which the perspective plot sits.

Data

library(mgcv)
df <- structure(list(death_or_neurodeficit1y = c(0, 1, 1, 1, 0, 1), 
    age = c(53, 67, 76, 69, 81, 51), isMALE = c(1, 1, 1, 0, 1, 
    1), diabetes = c(0, 0, 0, 0, 0, 0), trop_t_s = c(296, 120, 
    24, 41, 71, 258), shock = c(1, 1, 1, 1, 1, 1), byst_cpr = c(1, 
    0, 1, 1, 1, 1), gfr_0 = c(107.1, 64.1, 64.9, 38, 54.2, 65.6
    ), OtherBiomarker = c(14562.291216349, 26308.5822361957, 
    15525.1522892513, 26273.8258989978, 23385.2426802908, 18708.4888976223
    ), laktat_0 = c(16.8, 55.5, 63.9, 40.1, 139.7, 27.6)), terms = death_or_neurodeficit1y ~ 
    age + isMALE + diabetes + trop_t_s + shock + byst_cpr + gfr_0 + 
        OtherBiomarker + laktat_0, na.action = structure(c(`26` = 26L, 
`33` = 33L, `43` = 43L, `46` = 46L, `57` = 57L, `70` = 70L, `91` = 91L, 
`93` = 93L, `105` = 105L, `115` = 115L, `116` = 116L, `118` = 118L, 
`121` = 121L, `127` = 127L, `135` = 135L, `160` = 160L), class = "omit"), row.names = c(NA, 
6L), class = "data.frame")

Model

gam.model <- gam(death_or_neurodeficit1y ~ 
                   te(OtherBiomarker, laktat_0, bs = 'cr') + age + isMALE + diabetes + trop_t_s + shock + byst_cpr + gfr_0,
                 data = df,
                 family = "binomial")

par(cex = 1)
v <- vis.gam(gam.model,
        view = c("OtherBiomarker", "laktat_0"),
        plot.type = "persp",
        type = "response",
        cond = conds, 
        theta = 320,
        phi = 25,
        color = "heat",
        xlab = "\n\nOtherBiomarker at admission",
        ylab = "\n\nLactate at admission",
        zlab = "\n\nPredicted Probability",
        border = "black",
        r = 100,
        ticktype = "detailed"
)

Essentially, the solution below is like a rug plot, but uses points rather than ticks. You can control the alpha channel of the points to both decrease their prominence and allow the density to be displayed. In the plot you made, you want to know where the points for OtherBiomarker (the x-variable) are when laktat_0 (the y-variable) takes on its smallest value and the predictions take on their smallest value (0 in this case). This you could find as follows (Note, the plot is saved as v above):

mny <- min(df$laktat_0, na.rm=TRUE)
xvals <- trans3d(df$OtherBiomarker, mny, 0, pmat=v)

We can do the same for the laktat_0 points at the minimum of OtherBiormarker

mnx <- min(df$OtherBiomarker, na.rm=TRUE)
yvals <- trans3d(mnx, df$laktat_0, 0, pmat=v)

Finally, you can just add the points to the graph:

points(xvals, col=rgb(1,0,0,.25), pch=16)
points(yvals, col=rgb(1,0,0,.25), pch=16)

Here's the resulting plot.