Search code examples
rfunctionggplot2area

How to shade region between 2 functions ggplot

I would like to shade region between 2 functions on my plot to make my own confidence interval. Indeed, I made nls function with a Monod equation, and I traced this function with the parameters I obtained. I also made a bootstrap to obtain the values of my confidence interval and trace 2 lines to delimit the interval. enter image description here The problem is I didn't manage to do it, even with some answers I find on forums. Also, I think there will be another problem because these lines crossed themselves.

Here is my code :

pr20Y = subset (donnees_tot_g_J4, color == "20Y")

bstart.values <- list(gm=0.7, s=0.05, k=0.2)
nls1 = nls(growth_rate ~ gm*((quantity-s)/(quantity-s+k)), data = pr20Y, start=start.values)  
summary(nls1)
#gm = 0.65 / s = 0.06 / k = 0.11
nls1.boot <- nlsBoot(nls1, niter = 1000)
summary(nls1.boot)
    #Median       2.5%      97.5%
#gm 0.65636222 0.61994560 0.69597337
#s  0.05735049 0.03765189 0.07138767
#k  0.11705112 0.07833971 0.16174788
growthJ4=ggplot(donnees_tot_g_J4, aes(x=quantity, y=growth_rate, col=color))+
  geom_point(size=3,aes(col=color))+
  stat_function(fun = function (quantity) 0.65*((quantity-0.06)/(quantity-0.06+0.11)), color="darkgoldenrod1", size=1)+
  stat_function(fun = function (quantity) 0.61994560*((quantity-0.03765189)/(quantity-0.03765189+0.07833971)), color="darkgoldenrod1", size=1, alpha=0.5)+
  stat_function(fun =function (quantity) 0.69597337*((quantity-0.07138767)/(quantity-0.07138767+0.16605772)), color="darkgoldenrod1", size=1, alpha=0.5)+
  scale_x_continuous(breaks=c(0.1,0.3,0.6,0.9,1.5), limits=c(0.1,1.5))+
  scale_y_continuous(limits=c(0,0.8))+
  scale_colour_manual(values=c("20S" = "aquamarine1","25S" = "aquamarine3","28S" = "aquamarine4","20Y" = "darkgoldenrod1","25Y" = "darkgoldenrod3", "28Y" = "darkgoldenrod4"))+
  scale_fill_manual(values=c("20S" = "aquamarine1","25S" = "aquamarine3","28S" = "aquamarine4","20Y" = "darkgoldenrod1", "25Y" = "darkgoldenrod3","28Y" = "darkgoldenrod4"))+
  theme_minimal()+
  ggtitle("Taux de croissance à J4 en fonction du traitement de nourriture et de la température")+
  xlab("Quantité") + ylab("Taux de croissance µg/j")+
  theme_grey(base_size = 22)
growthJ4

And my data :

structure(list(name = c("J4_S01AC", "J4_S01CC", "J4_S01EC", "J4_S01FC", 
"J4_S03BC", "J4_S03CC", "J4_S03DC", "J4_S03EC", "J4_S03FC", "J4_S06AC", 
"J4_S06DC", "J4_S06EC", "J4_S06FC", "J4_S06KC", "J4_S06MC", "J4_S06NC", 
"J4_S09AC", "J4_S09BC", "J4_S09CC", "J4_S09DC", "J4_S03AM", "J4_S03BM", 
"J4_S03CM", "J4_S15AM", "J4_S15BM", "J4_S15CM", "J4_S15DM", "J4_S01AAC"
), day = c("J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", 
"J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", 
"J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4"), quality = c("S", 
"S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", 
"S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", 
"S"), quantity = c(0.1, 0.1, 0.1, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.9, 0.9, 0.9, 0.9, 0.3, 0.3, 
0.3, 1.5, 1.5, 1.5, 1.5, 0.1), qual_quant = c("S01", "S01", "S01", 
"S01", "S03", "S03", "S03", "S03", "S03", "S06", "S06", "S06", 
"S06", "S06", "S06", "S06", "S09", "S09", "S09", "S09", "S03", 
"S03", "S03", "S15", "S15", "S15", "S15", "S01"), temperature = c(28, 
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28), time = c(101, 101, 
101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 
101, 102, 102, 102, 102, 107.5, 107.5, 107.5, 107.5, 107.5, 107.5, 
107.5, 109), size = c(1344.366, 1291.962, 1304.811, 1128.264, 
1286.151, 1358.421, 1396.671, 1353.076, 1505.565, 1297.17, 0, 
1243.029, 1323.85, 1368.364, 1506.396, 0, 1663.735, 1632.28, 
2115.303, 1921.46, 1506.581, 1501.196, 1370.99, 1870.489, 1941.425, 
1942.186, 1827.395, 1336.588), weight = c(11, 16, 10, 10, 12, 
12, 13, 16, 14, 12, 11, 12, 10, 10, 15, 25, 46, 35, 66, 46, 20, 
16, 15, 49, 73, 63, 60, 11), growth_rate = c(0.0378568479840844, 0.126892202895244, 
0.0152090656484838, 0.0152090656484838, 0.0585326533474545, 0.0585326533474545, 
0.0775525505364441, 0.126892202895244, 0.0951622248242906, 0.0585326533474545, 
0.0378568479840844, 0.0585326533474545, 0.0152090656484838, 0.0152090656484838, 
0.111556439370444, 0.232939774941222, 0.374129156018743, 0.309825356331568, 
0.459072793066665, 0.374129156018743, 0.169037347727828, 0.119219651092275, 
0.104811166292231, 0.369092608582107, 0.458090402952369, 0.425199566970549, 
0.414306966296783, 0.0350783637283718), color = c("28S", "28S", 
"28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", 
"28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", 
"28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S")), row.names = c(NA, 
-28L), class = c("tbl_df", "tbl", "data.frame"))

Or maybe it exists another way to make confidence interval, because with all my other treatment, the 2 lines crossed themselves (enter image description here)and it will not be great for my report, but I don't know how to better ...

Solution

  • Thank you Gregor, I change the way to create my confidence interval with one of the topic you shared and it works !

    I used this code for all my treatment :

    new.data <- data.frame(quantity=seq(0.1, 1.5, by = 0.01))
interval <- as_tibble(predFit(nls1, newdata = new.data, interval = "confidence", level= 0.95)) %>% 
  mutate(quantity = new.data$quantity)

    and I obtain this : plot with confidence interval