Response surface plots with original experimental points overlaid in R

Question

It's been two weeks and I'm REALLY stuck now. I want to plot the original experimental points on a response surface, How to I do that from the R code https://www.dropbox.com/s/xme9lhkp3s2zwj0/stackoverflow_EXAMPLE.R?dl=0

    library(readxl)
library(rsm)
project.dir <- "/some/User/Documents/R_PhD_Data/"
setwd(project.dir)
data.dir <- paste0(project.dir, "RSM/data_HPLC")
mydata<-"LA_EF.xlsx"
mydata<-read_excel(paste0(data.dir,"/",mydata))
str(mydata)
names(mydata)
# data management
mydata.coded<-coded.data(mydata, x1 ~ (Tr_n - 90.5)/89.5, x2 ~ (Ni_acid - 5)/5, x3 ~ (The - 2.5)/2.5)
# check the well designed data.frame
mydata.coded

# START MY MODEL
# je récupère que les colonnes de paramètres analyser
# I only get the parameter columns analyzed
myvar<-names(mydata)[4 :ncol(mydata)] #starts from propionic acid [4] to the last variable (SBCFA) [7]
#x1 <- mydata$Tr_n
# je lance l'analyse sur toutes ces colonnes, l'une après l'autre et le résultat se met dans une liste
# I run the analysis on all these columns, one after the other, and the result goes into a list
mymodel<-lapply(myvar,
                function(k){
                  rsm(eval(substitute(j~SO(x1,x2,x3),list(j=as.name(k)))),data = mydata.coded)
                }
)
names(mymodel)<-myvar
length(mymodel)

# verification sur quelques objets
# verification on some objects
#summary(mymodel[[7]])
summary(mymodel[[1]])

# j'affiche les résumé des modèles pour chaque variable testé
# I display the summary of the models for each variable tested

# je fais une boucle et je parcours la liste contenant les résultats des modèles
# I loop through the list containing the results of the models

for (i in 1:length(myvar)){
  print(summary(mymodel[[i]]))
}

#### START TO PRINT MY IMAGES
# This piece of code displays teh last response variable
# Ideally I'd like to print all responses in individual images in case I will send these in a report to my supervisors

par(mfrow=c(2,2),cex.axis=0.4,cex.lab=0.65,cex.sub=0.7)
persp(
  mymodel[[i]],~x1+x2+x3,
  at=summary(mymodel[[i]])$canonical$xs,
  zlab=myvar[i],
  col=rainbow(100),
  contours="colors",
  ticktype = "detailed")
# from here I get an error
  lines(trans3d( ), col = 3)  # problem line here 
                              # basically I'm not sure of the syntax to make it use my model
                              # and automatically draw my orginical points from the first variables of my data file
  
  
  # I also tried this line of code
  # with(mymodel[[i]], points(trans3d(mydata.coded$Tr_n,mydata.coded$Ni_acid,mydata.coded$The,pmat), pch=20))
  # still get an error
  
  #### END TO PRINT MY IMAGES
  
  # HERE I PRINT ALL THE IMAGES INTO A SINGLE PDF, ideally with the experimental points displayed 
  # to see the how different the response surface is from the experimental points
  # I also want to display the experimental points on these images.
  
  
  #this code works fine, just that it doesn't have an option to print the experimental points (Trp, Nic_acd, & Thiam)
  pdf("SO_LA_HPLC_BB1_PR.pdf" ,paper="a4r")
  for (i in 1:length(myvar)){
    par(mfrow=c(2,2),cex.axis=0.4,cex.lab=0.65,cex.sub=0.7)
    persp(mymodel[[i]],~x1+x2+x3,at=summary(mymodel[[i]])$canonical$xs,zlab=myvar[i],col=rainbow(100),contours="colors")
    #contour(mymodel[[i]],~x1+x2+x3,at=summary(mymodel[[i]])$canonical$xs,zlab=myvar[i],col=rainbow(100),contours="colors")
    print(i)
  }
  dev.off()
  summary(mymodel)

This is a typical dataset that I'd use https://www.dropbox.com/scl/fi/shmx1t5gzjn06fxshvmxv/LA_EF_stackOverflow.xlsx?dl=0&rlkey=jf1w9isqtvl8eff31p1gs8kac

my data <- structure(list(Tr_n = c(90.5, 90.5, 1, 90.5, 180, 1, 1, 180, 
1, 180, 90.5, 90.5, 90.5, 90.5, 180, 90.5), Ni_acid = c(5, 0.1, 
0.1, 0.1, 5, 5, 5, 5, 10, 0.1, 10, 10, 5, 5, 10, 5), The = c(2.5, 
5, 2.5, 0.015, 0.015, 0.015, 5, 5, 2.5, 2.5, 0.015, 5, 2.5, 2.5, 
2.5, 2.5), `propionic acid` = c(3.3388, 5.1659, 4.7514, 4.1143, 
3.8044, 3.5787, 3.3176, 3.5711, 3.3084, 5.0962, 3.3948, 3.3281, 
3.5362, 3.4132, 3.3848, 3.3734), `2MBA` = c(0.332, 0.1454, 0.0923, 
0.0352, 0.1251, 0.1071, 0.3233, 0.3307, 0.3012, 0.2972, 0.0912, 
0.3407, 0.3391, 0.3232, 0.3251, 0.3216), IBA = c(0.3443, NA, 
0.0012, 0.0103, 0.1046, 0.1297, 0.3049, 0.3345, 0.3328, 0.0093, 
0.1028, 0.4319, 0.347, 0.3361, 0.4084, 0.3453), methCry = c(0.023, 
0.2008, 0.1565, 0.245, 0.0848, 0.2888, 0.0208, 0.0234, 0.0214, 
0.234, 0.0986, 0.0398, 0.0329, 0.0337, 0.0293, 0.0175), `2K3MV` = c(0.332, 
0.1454, 0.0923, 0.0352, 0.1251, 0.1071, 0.3233, 0.3307, 0.3012, 
0.2972, 0.0912, 0.3407, 0.3391, 0.3232, 0.3251, 0.3216), ILE = c(0.77, 
0.4, 0.35, 0.33, 0.47, 0.5, 0.74, 0.73, 0.71, 0.35, 0.48, 0.7, 
0.73, 0.78, 0.75, 0.78), SBCFA = c(76.3479575, 33.9246734, 36.2150499, 
46.8351514, 82.1730982, 70.9092779, 76.6372188, 74.6911517, 76.9682743, 
39.4543268, 84.4243206, 72.2931598, 72.6038991, 73.9702339, 75.4052285, 
77.2390033)), row.names = c(NA, 16L), class = "data.frame")