The QuantiSlakeTest, dynamic weighting of soil under water to measure soil structural stability
Openscience & Reproducible Research of Vanwin & Hardy (2022)

Geek zone
rm(list=ls())
## load(file= "./data-output/qst-openscience.rda")

Geek zone
library(dplyr)
library(ggplot2) 
require(ggpubr) ## for one graph with ggarrange
require(lme4) ## for models
require(gt) # for beautiful table
require(car)
require(emmeans)
library(conflicted)

Geek zone
## ggplot2
conflict_prefer("annotate", "ggplot2")

Geek zone
illustr.qst.curve.early.nb.nd <- function(my.data.slake, my.qst.opendata) {
  ## data for testing the function
  ## my.data.slake <- data.slake.full.example
  ## my.qst.opendata <- qst.opendata
  ## max.time = my.qst.opendata[as.character(id.slakes), 'tmax']
  id.slakes <- my.data.slake$slake[1]
  max.time = my.qst.opendata[as.character(id.slakes), 'tmax']

  graph <- (my.data.slake %>%
            ggplot2::ggplot(aes(x =  difftime.num , y = pdiffmass))
    + geom_ribbon(aes(ymax = pdiffmass),ymin=0, xmin=0,
                  fill="gray60",colour=NA)
    + geom_ribbon(data = my.data.slake[8:11,],aes(ymax = pdiffmass),ymin=0, xmin=my.data.slake$difftime.num[8],
                  fill="gray80",colour=NA)
    + geom_ribbon(data = my.data.slake[1:8,],aes(ymax = pdiffmass),ymin=0,
                  fill="gray90",colour=NA)
    + geom_hline(yintercept = 1, col = "gray15", linewidth = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.data.slake$difftime.num[1:100], col = "gray5", linewidth = 0.3)
    + geom_vline(xintercept = my.data.slake$difftime.num[c(1,5,10,15,20,25,30,35)],
                 col = "gray5",
                 linewidth = 0.6)
    + geom_ribbon(aes(ymin=pdiffmass), ymax = 10, fill="white", colour=NA, alpha=0.8)
    + geom_line(colour = 'gray20', linewidth = 3)
    + geom_point(colour = 'gray80', size = 0.5)
    ## t0
    + geom_vline(xintercept = 0, col = "gray15", linewidth = 1)
    + geom_point(
        x =  0 , y = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray89", size = 3)
    + annotate(geom = "label", x = 0,
               y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
               label = "W[t0]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               hjust = 1.05,
               vjust = 1.2
               )
    ## Archimedes buoyancy
    + geom_vline(xintercept = my.data.slake$difftime.num[1],
                 col = "gray15", linewidth = 1)
    + geom_vline(xintercept = my.data.slake$difftime.num[8],
                 col = "gray15", linewidth = 1)
    + geom_segment(
        x =  my.data.slake$difftime.num[3],
        xend = my.data.slake$difftime.num[3],
        y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']),
        yend = my.qst.opendata[as.character(id.slakes), 'Wt0'],
        col = "gray15",
        linewidth = 1,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + geom_point(
        x =  my.data.slake$difftime.num[3],
        y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+(1/my.qst.opendata[as.character(id.slakes), 'Warchi']))/2,
        col = "gray15",
        size = 2
      )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[3],
               y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+(1/my.qst.opendata[as.character(id.slakes), 'Warchi']))/2,
               label = "Archimedes' \n buoyancy",
               ## parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               hjust = -0.05
               )
    ## tmax
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", linewidth = 1)
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax'], y = 1, col = "gray89", size = 3)
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 'tmax'],
               y = 1.1,
               ## label = paste("t[max]==", round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),"~sec"),
               label = "t[max]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               hjust = -0.05
               )
    ## data acquisition*
    ## no need as asked by reviewer
    ## + annotate(geom = "label", x = my.data.slake$difftime.num[1],
    ##            y = 0.2,
    ##            ## y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 1.3),
    ##            label = "Acquired data (*) :",
    ##            ## size=5,
    ##            label.size=0,
    ##            hjust = 0.1
    ##            )
    + annotate(geom = "text", x = my.data.slake$difftime.num[c(1,5,10,15,20,25,30,35)],
               y = 0,
               ## y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 1.2),
               label = c(1,5,10,15,20,25,30,35),
               ## size=5,
               angle = 90, hjust = 0, vjust = 0
               )
    ## Under water info
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[1],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.6),
               label = "1: Preparation phase, \nsoil sample is in the air",
               ## size=3,
               hjust = 0, vjust = 0
               )
    + geom_segment(
        x =  my.data.slake$difftime.num[1],
        xend = my.data.slake$difftime.num[8],
        y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 0.5),
        yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 0.5),
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[8],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.2),
               label = "2: Immersion phase",
               ## size=3,
               hjust = 0, vjust = 0
               )
    + geom_segment(
        x =  my.data.slake$difftime.num[8],
        xend = my.data.slake$difftime.num[11],
        y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.1),
        yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.1),
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[11],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.3),
               label = "3: QuantiSlakeTest, \nsoil sample is under the water",
               ## size=3,
               hjust = 0, vjust = 0
               )
    + geom_segment(
        x =  my.data.slake$difftime.num[11],
        xend = my.data.slake$difftime.num[50],
        y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.4),
        yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.4),
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[12],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.7),
               label = "Early increase in soil mass",
               ## size=3,
               hjust = 0
               )
    + geom_segment(
        x =  my.data.slake$difftime.num[11],
        xend = my.qst.opendata[as.character(id.slakes), 'tmax'],
        y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.8),
        yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.8),
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    ## Slope ?
    + geom_segment(
        x =  0,
        xend = my.qst.opendata[as.character(id.slakes), 'tmax'],
        y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
        yend = 1,
        col = "gray89",
        linewidth = 1
      )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax']/2,
        y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
        col = "gray89",
        size = 2
      )
    + annotate(geom = "label",
               x = my.qst.opendata[as.character(id.slakes), 'tmax']/2,
               y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
               label = "Slope[0-max]",
               parse = TRUE,
               col = "gray15", fill = "gray88",  label.r = unit(0, "lines"),
               label.size = 1, hjust = 0.5, vjust = 1.2
               )
    ## Legend
    + annotate("rect",
               xmin=240, xmax =1000,
               ymin=0.80,ymax =1.15,
               fill="white", colour="gray16"
               )
    + annotate(geom = "text", x = 250, y = 1.12,
               label = "bold('QuantiSlakeTest curve')",
               parse = TRUE,
               size = 5, hjust = 0
               )
    + annotate(geom = "text", x = 250, y = 1.07,
               label = "bold('and main indicators')",
               size=5, parse = TRUE,
               hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 1.01,  
               label = "'(i) Indicators related to the early increase in soil mass'",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), size = 1,
               label.size = 3, hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 0.95,
               label = "'(ii) Slopes in the decreasing part of the curve'",
               parse = TRUE,
               col = "gray16", size = 1,
               label.size = 3, hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 0.89,
               label = "(iii) Threshold values of mass loss", size = 3,
               fill = "white",
               hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 0.83,
               label = "bold('(iv) Global indicators')",
               parse = TRUE,
               col = "gray90",
               fill = "gray10",
               label.size = 3, hjust = 0
               )
    + labs(y = "Relative soil mass [-]",
           x = "Time [sec]"# , 
           )
    + coord_cartesian(xlim=c(NA,max.time+2),ylim=c(0,3))
    + scale_x_continuous(breaks=c(0,
                                  round(my.data.slake$difftime.num[1],1),
                                  round(my.data.slake$difftime.num[8],1),
                                  round(my.data.slake$difftime.num[20],1),
                                  ## round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 'tmax'],1)
                                  ))
    + scale_y_continuous(breaks=c(0,
                                  my.qst.opendata[as.character(id.slakes), 'Wt0'],
                                  round(1/my.qst.opendata[as.character(id.slakes), 'Warchi'],2),
                                  1))
    + theme_classic(base_size=16)
  )
  print(graph)
}
## IDEA with DERIV'
## graph <- (my.data.slake %>%
##           dplyr::select(difftime.num,pdiffmass,derivmass,derivsecmass) %>%
##           tidyr::pivot_longer(cols = -1) %>%          
##           ggplot2::ggplot(aes(
##   x =  difftime.num, y = value))
##   + geom_line(colour = 'gray20', linewidth = 3)
##   + geom_line(colour = 'gray20', linewidth = 3)
##   + geom_vline(xintercept = 0, col = "gray15", linewidth = 0.3)
##   + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", linewidth = 0.3)
##   ## Legend
##   + annotate("rect",
##              xmin=240, xmax =1000,
##              ymin=0.80,ymax =1.15,
##              fill="white", colour="gray16"
##              )
##   + annotate(geom = "text", x = 250, y = 1.12,
##              label = "bold('QuantiSlakeTest curve')",
##              parse = TRUE,
##              size = 5, hjust = 0
##              )
##   + annotate(geom = "text", x = 250, y = 1.07,
##              label = "bold('and main indicators')",
##              size=5, parse = TRUE,
##              hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 1.01,  
##              label = "'(i) Indicators related to the early increase in soil mass'",
##              parse = TRUE,
##              col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), size = 1,
##              label.size = 3, hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 0.95,
##              label = "'(ii) Slopes in the decreasing part of the curve'",
##              parse = TRUE,
##              col = "gray16", size = 1,
##              label.size = 3, hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 0.89,
##              label = "(iii) Threshold values of mass loss", size = 3,
##              fill = "white",
##              hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 0.83,
##              label = "bold('(iv) Global indicators')",
##              parse = TRUE,
##              col = "gray90",
##              fill = "gray10",
##              label.size = 3, hjust = 0
##              )
##   + labs(y = "Relative soil mass [-]",
##          x = "Time [sec]"# , 
##          )
##   + coord_cartesian(xlim=c(NA,max.time+2))
##   + scale_x_continuous(breaks=c(0,
##                                 ## round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),
##                                 round(my.qst.opendata[as.character(id.slakes), 'tmax'],0)
##                                 ))
##   + scale_y_continuous(breaks=c(0,
##                                 my.qst.opendata[as.character(id.slakes), 'Wt0'],
##                                 my.qst.opendata[as.character(id.slakes), 'Wend'],
##                                 1)
##   + facet_grid(name~., scales = "free_y")
##   + theme_classic(base_size=16)

Geek zone
illustr.qst.curve.indic.nb.nd <- function(my.data.slake, my.qst.opendata,  max.time = 1000) {

  id.slakes <- my.data.slake$slake[1]

  graph <- (ggplot2::ggplot(my.data.slake, aes(
                                             x =  difftime.num, y = pdiffmass))
    + geom_ribbon(aes(ymax =pdiffmass),ymin=0, fill="gray60",colour=NA)
    + geom_vline(xintercept = 0, col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't25'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't50'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't75'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't90'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't95'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't99'], col = "gray15", linewidth = 0.3)
    + geom_hline(yintercept = 1, col = "gray15", linewidth = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wend'], col = "black", linewidth = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray15", linewidth = 0.3)
    + geom_ribbon(aes(ymin=pdiffmass),ymax =2, fill="white",colour=NA,alpha=0.8)
    + annotate("rect", xmin=1000,xmax =1200,ymin=0.5,ymax =2, fill="white",colour=NA,alpha=0.8)
    + geom_line(colour = 'gray20', linewidth = 3)
    + geom_point(
        x =  0 , y = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray89", size = 3)
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax'], y = 1, col = "gray89", size = 3)
    + geom_segment(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax'],
        xend = my.data.slake$difftime.num[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 60)][1],
        y = 1,
        yend = my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 60)][1],
        col = "gray89",
        linewidth = 0.6
      )
    + geom_segment(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax'],
        xend = my.data.slake$difftime.num[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 300)][1],
        y = 1,
        yend = my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 300)][1],
        col = "gray66",
        linewidth = 0.6
      )
    + geom_segment(
        x =  my.data.slake$difftime.num[
                             my.data.slake$difftime.num > (
                               my.qst.opendata[as.character(id.slakes), 'tmax'
                                               ] + 60)][1],
        xend = my.data.slake$difftime.num[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 300)][1],
        y = my.data.slake$pdiffmass[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes), 'tmax'
                                              ] + 60)][1],
        yend = my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 300)][1],
        col = "gray66",
        linewidth = 0.6
      )
    + geom_segment(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax'],
        xend = my.data.slake$difftime.num[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 600)][1],
        y = 1,
        yend = my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 600)][1],
        col = "gray66",
        linewidth = 0.6
      )
    + geom_segment(
        x =  my.data.slake$difftime.num[
                             my.data.slake$difftime.num > (
                               my.qst.opendata[as.character(id.slakes), 'tmax'
                                               ] + 300)][1],
        xend = my.data.slake$difftime.num[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 600)][1],
        y = my.data.slake$pdiffmass[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes), 'tmax'
                                              ] + 300)][1],
        yend = my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'
                                                 ] + 600)][1],
        col = "gray66",
        linewidth = 0.6
      )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax']/2,
        y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
        col = "gray89",
        size = 1
      )
    + geom_segment(
        x =  0,
        xend = my.qst.opendata[as.character(id.slakes), 'tmax'],
        y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
        yend = 1,
        col = "gray89",
        linewidth = 0.3
      )
    + geom_segment(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax']/2,
        y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
        xend = -10,
        yend = 1.1,
        col = "gray15",
        linewidth = 0.3
      )
    + annotate(geom = "label", -10,
               y = 1.1,
               label = "Slope[0-max]",
               parse = TRUE,
               col = "gray15", fill = "gray88",  label.r = unit(0, "lines"),## , size = 1,
               label.size = 1, hjust = 1
               )
    + geom_segment(
        x =  -40,
        xend = -40,
        y = 1,
        yend = my.qst.opendata[as.character(id.slakes), 'Wt0'],
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + geom_segment(
        x =  my.qst.opendata[as.character(id.slakes), 't50'],
        xend = my.qst.opendata[as.character(id.slakes), 't75'],
        y = 0.1,
        yend = 0.1,
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + geom_segment(
        x =  my.qst.opendata[as.character(id.slakes), 't75'],
        xend = my.qst.opendata[as.character(id.slakes), 't90'],
        y = 0.05,
        yend = 0.05,
        col = "gray15",
        linewidth = 0.3,
        arrow = arrow(ends="both", length = unit(0.008, "npc"))
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't50'],
               y = 0.1,
               label = "dt[50-75]", 
               parse=TRUE,
               label.size = 1,
               fill = "white",
               hjust = 1.1
               )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't90'],
               y = 0.05,
               label = "dt[75-90]",
               parse=TRUE,
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(
        x =  -40,
        y = mean(c(1,my.qst.opendata[as.character(id.slakes), 'Wt0'])),
        col = "gray15",
        size = 1
      )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 'tmax'],
        y = 1.1,
        col = "gray15",
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 'tmax'],
               y = 1.1,
               label = "t[max]",
               label.size = 1,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), ## size = 1,
               hjust = -0.05
               )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 't25'],
        y = 0.30,
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't25'],
               y = 0.30,
               label = "t25", 
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 't50'],
        y = 0.25,
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't50'],
               y = 0.25,
               label = "t50",
               label.size = 1,
               fill = "white",             
               hjust = -0.05
               )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 't75'],
        y = 0.2,
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't75'],
               y = 0.2,
               label = "t75", 
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 't90'],
        y = 0.19,
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't90'],
               y = 0.19,
               label = "t90",
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 't95'],
        y = 0.19,
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't95'],
               y = 0.19,
               label = "t95", 
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(
        x =  my.qst.opendata[as.character(id.slakes), 't99'],
        y = 0.19,
        size = 1
      )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't99'],
               y = 0.19,
               label = "t99",
               label.size = 1,
               fill = "white",
               hjust = 1.05
               )
    + geom_point(
        x =  -100,
        y = 1,
        col = "gray5",
        size = 1
      )
    + annotate(geom = "label", x = -100,
               y = 1,
               label = "W[max]",
               label.size = 1,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), # size = 1,
               hjust = 1.05,
               vjust = -0.2
               )
    + geom_point(
        x =  -100,
        y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
        col = "gray5",
        size = 1
      )
    + annotate(geom = "label", x = -100,
               y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
               label = "W[t0]",
               label.size = 1,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), # size = 1,
               hjust = 1.05,
               vjust = 1.2
               )
    + geom_point(
        x =  my.data.slake$difftime.num[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1],
        y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1],
        col = "gray89", size = 3)
    + geom_point(
        x =  my.data.slake$difftime.num[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1],
        y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1],
        col = "gray89", size = 3)
    + geom_point(
        x =  my.data.slake$difftime.num[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1],
        y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1],
        col = "gray89", size = 3)
    + annotate(geom = "label", x = -50,
               y = mean(c(1,my.qst.opendata[as.character(id.slakes), 'Wt0'])),
               label = "W[max]-W[t0]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), ## size = 1,
               label.size = 1, hjust = 1.05,
               )
    + geom_point(
        x =  -50,
        y = my.qst.opendata[as.character(id.slakes), 'Wend'],
        col = "black",
        size = 1
      )
    + annotate(geom = "label", x = -50,
               y = my.qst.opendata[as.character(id.slakes), 'Wend'],
               label = "bold('W')[bold('end')]",
               label.size = 1,
               parse = TRUE,
               col = "gray90",
               ## fill = "lightgrey",
               fill = "gray10",
               hjust = 1.05,
               vjust = -0.05
               )
    + geom_point(
        x = mean(c(as.numeric(
          my.data.slake$difftime.num[
                          my.data.slake$difftime.num > (
                            my.qst.opendata[as.character(id.slakes),
                                            'tmax'] + 60)
                        ][1]),
          as.numeric(
            my.data.slake$difftime.num[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes),
                                              'tmax'] + 300)
                          ][1]))
          ),
        y = mean(
          c(my.data.slake$pdiffmass[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)
                          ][1],
            my.data.slake$pdiffmass[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                          ][1])
        ),
        col = "gray66",
        size = 1
      )
    + annotate(geom = "label",
               x = mean(c(as.numeric(
                 my.data.slake$difftime.num[
                                 my.data.slake$difftime.num > (
                                   my.qst.opendata[as.character(id.slakes),
                                                   'tmax'] + 60)
                               ][1]),
                 as.numeric(
                   my.data.slake$difftime.num[
                                   my.data.slake$difftime.num > (
                                     my.qst.opendata[as.character(id.slakes),
                                                     'tmax'] + 300)
                                 ][1]))
                 ),
               y = mean(
                 c(my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (
                                     my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)
                                 ][1],
                   my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (
                                     my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                                 ][1])
               ),
               label = "Slope[60-300]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = 1.1, vjust = 1.1
               )
    + geom_point(
        x = mean(c(as.numeric(
          my.data.slake$difftime.num[
                          my.data.slake$difftime.num > (
                            my.qst.opendata[as.character(id.slakes),
                                            'tmax'] + 600)
                        ][1]),
          as.numeric(
            my.data.slake$difftime.num[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes),
                                              'tmax'] + 300)
                          ][1]))
          ),
        y = mean(
          c(my.data.slake$pdiffmass[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)
                          ][1],
            my.data.slake$pdiffmass[
                            my.data.slake$difftime.num > (
                              my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                          ][1])
        ),
        col = "gray66",
        size = 1
      )
    + annotate(geom = "label",
               x = mean(c(as.numeric(
                 my.data.slake$difftime.num[
                                 my.data.slake$difftime.num > (
                                   my.qst.opendata[as.character(id.slakes),
                                                   'tmax'] + 600)
                               ][1]),
                 as.numeric(
                   my.data.slake$difftime.num[
                                   my.data.slake$difftime.num > (
                                     my.qst.opendata[as.character(id.slakes),
                                                     'tmax'] + 300)
                                 ][1]))
                 ),
               y = mean(
                 c(my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (
                                     my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)
                                 ][1],
                   my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (
                                     my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                                 ][1])
               ),
               label = "Slope[300-600]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = 0.5, vjust = 1.1
               )
    + geom_point(
        x = mean(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                   as.numeric(
                     my.data.slake$difftime.num[
                                     my.data.slake$difftime.num > (
                                       my.qst.opendata[as.character(id.slakes),
                                                       'tmax'] + 60)
                                   ][1]))
                 ),
        y = mean(
          c(1, my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (
                                 my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)
                             ][1])
        ),
        col = "gray89",
        size = 1
      )
    + annotate(geom = "label",
               x = mean(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                          as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (
                            my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1]))),
               y = mean(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (
                 my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1])),
               label = "Slope[max-60]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = -0.05, vjust = -0.05
               )
    + geom_point(
        x =  quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                        as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (
                          my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1])),0.5),
        y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1]),0.5),
        col = "gray66",
        size = 1
      )
    + annotate(geom = "label",
               x = quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                              as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (
                                my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1])),0.5),
               y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (
                 my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1]),0.5),
               label = "Slope[max-300]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = -0.05, vjust = -0.05
               )
    + geom_point(
        x =  quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                        as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (
                          my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1])),0.75),
        y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (
          my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1]),0.25),
        col = "gray66",
        size = 1
      )
    + annotate(geom = "label",
               x = quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                              as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (
                                my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1])), 0.75),
               y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (
                 my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1]),0.25),
               label = "Slope[max-600]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = -0.05, vjust = -0.05
               )
    + annotate(geom = "label", x = 500, y = 0.25,
               label = "bold('Area under curve (AUC)')",
               parse = TRUE,
               col = "gray90",
               fill = "gray10",
               label.size = 1, hjust = 0
               )
    ## Legend
    + annotate("rect",
               xmin=260, xmax =1000,
               ymin=0.80,ymax =1.15,
               fill="white", colour="gray16"
               )
    + annotate(geom = "text", x = 270, y = 1.12,
               label = "bold('QuantiSlakeTest curve')",
               parse = TRUE,
               ## label.size = 5,
               hjust = 0
               )
    + annotate(geom = "text", x = 270, y = 1.08,
               label = "bold('and main indicators')",
               parse = TRUE,
               ## label.size = 5,
               hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 1.01,  
               label = "'(i) Indicators of the early increase in soil mass'",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), # size = 1,
               label.size = 1, hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 0.95,
               label = "'(ii) Slopes in the decreasing part of the curve'",
               parse = TRUE,
               col = "gray16", # size = 1,
               label.size = 0.3, hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 0.89,
               label = "(iii) Threshold values of mass loss", label.size = 1,
               fill = "white",
               hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 0.83,
               label = "bold('(iv) Global indicators')",
               parse = TRUE,
               col = "gray90",
               fill = "gray10",
               label.size = 1, hjust = 0
               )
    + labs(y = "Relative soil mass [-]",
           x = "Time [sec]"# , 
           )
    + coord_cartesian(xlim=c(-200,max.time), ylim = c(0, 1.1))
    + scale_x_continuous(breaks=c(0,
                                  ## round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't50'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't90'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't95'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't99'],0)
                                  ))
    + scale_y_continuous(breaks=c(0,
                                  my.qst.opendata[as.character(id.slakes), 'Wt0'],
                                  my.qst.opendata[as.character(id.slakes), 'Wend'],
                                  1))
    + theme_classic(base_size=16)
  )

  print(graph)
}

Geek zone
f.df.corr <- function(df, var.l, var.c,
                      my.round = 3
                      ){

  ## keep only numeric data
  df.num <- df[, sapply(df, is.numeric)] 

  ## na.omit(df.num)
  mcor <- cor(df.num)

  df.corr <- as.data.frame(round(mcor,
                                 my.round))[var.l, var.c]
  return(df.corr)
}

f.gt.corr.color <- function(df,
                            my.pal = RColorBrewer::brewer.pal(n=8, name="RdYlBu")){
  tab <- df %>%
    dplyr::mutate(" " = "    ") %>%
    tibble::rownames_to_column("rownames") %>%
    gt(rowname_col = "rownames") %>%
    data_color(
      columns = names(df),
      fn = scales::col_numeric(
                     palette = my.pal,
                     domain = c(-1,1)
                   )) 
  return(tab)
}

Geek zone
plot.compare.trt <- function(trial = "P-K mineral fertiliser",
                             my.color="black"){ 
  df <- qst.opendata %>%
    dplyr::filter(essai_uk == trial) %>%
    dplyr::select(modalite_uk_abb, parcelle, bloc,
                  Wmax.Wt0,
                  tmax,
                  Sl.max.30,
                  Sl.max.60,
                  Sl.max.300,
                  ## Sl.max.600,              
                  t75,
                  t95,
                  Wend,
                  AUC) %>%
    dplyr::rename(
             "(a) Wmax-Wt0 [-]" = Wmax.Wt0,
             "(b) tmax [sec]" = tmax,
             "(c) Slope max-30 [sec-1]" = Sl.max.30,
             "(d) Slope max-60 [sec-1]" = Sl.max.60,
             "(e) Slope max-300 [sec-1]" = Sl.max.300,
             ## "(e) Slope max-600 [sec-1]" = Sl.max.600,
             "(f) t75 [sec]" = t75,
             "(g) t95 [sec]" = t95,
             "(h) Wend [-]" = Wend,
             "(i) AUC [sec]" = AUC) %>% 
    tidyr::pivot_longer(!c(modalite_uk_abb,parcelle,bloc),
                        names_to = "indicator",
                        values_to = "value"
                        )

  df$indicator <- factor(df$indicator,
                         levels=c(
                           "(a) Wmax-Wt0 [-]",
                           "(b) tmax [sec]",
                           "(c) Slope max-30 [sec-1]",
                           "(d) Slope max-60 [sec-1]",
                           "(e) Slope max-300 [sec-1]",
                           ## "(e) Slope max-600 [sec-1]" ,
                           "(f) t75 [sec]",
                           "(g) t95 [sec]",
                           "(h) Wend [-]",
                           "(i) AUC [sec]"
                         )
                         )

  levels(df$indicator) <- c(
    expression("(a)"~"W"["max"]~"-W"["t0"]~"[-]"),
    expression("(b)"~"t"["max"]~"[sec]"),
    expression("(c)"~"Slope"["max-30"]~"[sec"^{-1}~"]"),
    expression("(d)"~"Slope"["max-60"]~"[sec"^{-1}~"]"),
    expression("(e)"~"Slope"["max-300"]~"[sec"^{-1}~"]"),
    expression("(f)"~"t75"~"[sec]"),
    expression("(g)"~"t95"~"[sec]"),
    expression("(h)"~"W"["end"]~"[-]"),
    expression("(i)"~"AUC"~"[sec]")
  )

  g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                              y=value))
    + geom_boxplot(colour="black",fill="gray80")
    + facet_wrap(.~indicator,
                 nrow=3, scales="free_y",
                 labeller = label_parsed
                 )
    + labs(title = paste0("Trial - ",trial),
           x = "Treatment",
           y = "Value of the QuantiSlakeTest indicator")
    + theme_bw()
  )
}

Geek zone
plot.compare.trt.v2 <- function(trial = "P-K mineral fertiliser",
                                my.color="black"){ 
  df <- qst.opendata %>%
    dplyr::filter(essai_uk == trial) %>%
    dplyr::select(modalite_uk_abb, parcelle, bloc,
                  tmax,
                  local.slope.30.60,
                  delta.t.50.75,
                  Wend) %>%
    dplyr::rename(
             "(a) tmax [sec]" = tmax,
             "(b) Slope 30-60 [sec-1]" = local.slope.30.60,
             "(c) dt 50-75 [sec]" = delta.t.50.75,
             "(d) Wend [-]" = Wend) %>% 
    tidyr::pivot_longer(!c(modalite_uk_abb,parcelle,bloc),
                        names_to = "indicator",
                        values_to = "value"
                        )

  df$indicator <- factor(df$indicator,
                         levels=c(
                           "(a) tmax [sec]",
                           "(b) Slope 30-60 [sec-1]",
                           "(c) dt 50-75 [sec]",
                           "(d) Wend [-]"
                         )
                         )

  levels(df$indicator) <- c(
    expression("(a)"~"t"["max"]~"[sec]"),
    expression("(b)"~"Slope"["30-60"]~"[sec"^{-1}~"]"),
    expression("(c)"~"dt"["50-75"]~"[sec]"),
    expression("(d)"~"W"["end"]~"[-]")
  )

  g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                              y=value))
    + geom_boxplot(colour="black",fill="gray80")
    + facet_wrap(.~indicator,
                 nrow=3, scales="free_y",
                 labeller = label_parsed
                 )
    + labs(title = paste0("Trial - ",trial),
           x = "Treatment",
           y = "Value of the QuantiSlakeTest indicator")
    + theme_bw()
  )
}

Geek zone
plot.compare.trt.v3 <- function(trial = "P-K mineral fertiliser",
                                my.color="black"){ 
  df <- qst.opendata %>%
    dplyr::filter(essai_uk == trial) %>%
    dplyr::select(modalite_uk_abb, parcelle, bloc,
                  tmax,
                  local.slope.60.300,
                  delta.t.50.75,
                  Wend) %>%
    dplyr::rename(
             "(a) tmax [sec]" = tmax,
             "(b) Slope 60-300 [sec-1]" = local.slope.60.300,
             "(c) dt 50-75 [sec]" = delta.t.50.75,
             "(d) Wend [-]" = Wend) %>% 
    tidyr::pivot_longer(!c(modalite_uk_abb,parcelle,bloc),
                        names_to = "indicator",
                        values_to = "value"
                        )

  df$indicator <- factor(df$indicator,
                         levels=c(
                           "(a) tmax [sec]",
                           "(b) Slope 60-300 [sec-1]",
                           "(c) dt 50-75 [sec]",
                           "(d) Wend [-]"
                         )
                         )

  levels(df$indicator) <- c(
    expression("(a)"~"t"["max"]~"[sec]"),
    expression("(b)"~"Slope"["60-300"]~"[sec"^{-1}~"]"),
    expression("(c)"~"dt"["50-75"]~"[sec]"),
    expression("(d)"~"W"["end"]~"[-]")
  )

  g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                              y=value))
    + geom_boxplot(colour="black",fill="gray80")
    + facet_wrap(.~indicator,
                 nrow=3, scales="free_y",
                 labeller = label_parsed
                 )
    + labs(title = paste0("Trial - ",trial),
           x = "Treatment",
           y = "Value of the QuantiSlakeTest indicator")
    + theme_bw()
  )
}

Geek zone
plot.compare.trt.biss <- function(
                                  trial = "P-K mineral fertiliser",
                                  my.color="black"){ 
  df <- qst.mean.sd.opendata %>%
    dplyr::filter(essai_uk == trial) %>%
    dplyr::select(modalite_uk_abb, parcelle.c, bloc,
                  MWD_T1,
                  MWD_T3, 
                  MWD_T2
                  ) %>%
    dplyr::rename(
             "(a) MWD 1" = MWD_T1,
             "(b) MWD 2" = MWD_T3,
             "(c) MWD 3" = MWD_T2
           ) %>% 
    tidyr::pivot_longer(!c(modalite_uk_abb,parcelle.c,bloc),
                        names_to = "indicator",
                        values_to = "value"
                        )

  df$indicator <- factor(df$indicator,
                         levels=c(
                           "(a) MWD 1",
                           "(b) MWD 2",
                           "(c) MWD 3"
                         )
                         )

  levels(df$indicator) <- c(
    expression("(a)"~"MWD"~"1"),
    expression("(b)"~"MWD"~"2"),
    expression("(c)"~"MWD"~"3")
  )


  g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                              y=value))
    + geom_boxplot(colour="black",fill="gray80")
    + facet_wrap(.~indicator,
                 nrow=1, scales="free_y",
                 labeller = label_parsed
                 )
    + labs(title = paste0("Trial - ",trial),
           x = "Treatment",
           y = "Value of the Le Bissonnais indicator")
    + theme_bw()
  )
}

## g <- plot.compare.trt.biss("Organic matter")
## print(g)

Geek zone
f.qst.anova  <- function(my.campagne, indic="Wend"){
  res <- list()

  df <- qst.opendata %>%
    dplyr::filter(campagne == my.campagne) %>%
    dplyr::mutate(trmt = as.factor(as.character(modalite_uk_abb))) %>%
    dplyr::mutate(bloc = as.factor(bloc)) %>%
    dplyr::mutate(parcelle = as.factor(parcelle)) %>%
    dplyr::select(indic,bloc,trmt,parcelle)

  ## Test normality
  test.norm <- shapiro.test(df[[indic]])
  res[["Shapiro"]] <-  test.norm$p.value

  test.var <- bartlett.test(df[[indic]] ~ trmt,
                            data=df)

  ## Test homoscedacity
  res[["Bartlett"]] <-  test.var$p.value

  ## ANOVA
  rmodel <- lmer(
    get(indic) ~ trmt + (1|bloc) + (1|parcelle),
    data=df
  )
  mod <- car::Anova(rmodel, test = "F")
  res[["Anova"]] <- mod

  return(res)
}

f.qst.means.comp <- function(my.campagne,indic="Wend"){
  df.lf <- qst.opendata %>%
    dplyr::filter(campagne == my.campagne) %>%
    dplyr::mutate(trmt = as.factor(as.character(modalite_uk_abb))) %>%
    dplyr::mutate(bloc = as.factor(bloc)) %>%
    dplyr::mutate(parcelle = as.factor(parcelle))

  rmodel <- lmer(
    get(indic) ~ trmt + (1|bloc) + (1|parcelle),
    data=df.lf
  )

  t <- emmeans::contrast(
                  emmeans::emmeans(rmodel, ~ trmt),
                  alpha = 0.05,
                  method = "pairwise"
                )  %>%
    as_tibble() %>%
    mutate(across(where(is.double), signif, digits = 3))
  return(t)
}

s.f.qst.means.comp <- function(my.campagne,my.indic="Wend"){
  t <- f.qst.means.comp(my.campagne,my.indic) %>%
    tidyr::separate(contrast,c("group1","group2"),sep = " - ") %>%
    dplyr::mutate(p.value = signif(p.value,2)) %>%
    dplyr::mutate(campagne = my.campagne, indic = my.indic) %>%
    dplyr::select(campagne,indic,group1,group2,p.value) %>%
    as.data.frame()
}

f.qst.anova.biss  <- function(my.campagne="Les Fonds", indic="MWD_T1"){
  res <- list()

  df <- qst.mean.sd.opendata %>%
    dplyr::filter(campagne == my.campagne) %>%
    dplyr::mutate(trmt = as.factor(as.character(modalite_uk_abb))) %>%
    dplyr::mutate(bloc = as.factor(bloc)) %>%
    dplyr::mutate(parcelle = as.factor(parcelle.c)) %>%
    dplyr::select(indic,bloc,trmt,parcelle.c)

  ## ANOVA
  rmodel <- lmer(
    get(indic) ~ trmt + (1|bloc),
    data=df
  )
  mod <- car::Anova(rmodel, test = "F")
  res[["Anova"]] <- mod

  return(res)
}

Geek zone
f.scatter.fred <- function(var_x,var_y,database,databasemeansd) {
  g <- (ggplot(data=database,
               aes(x=.data[[var_x]],
                   y=.data[[var_y]])
               )
    + geom_smooth(data=databasemeansd,
                  aes(x=.data[[paste0(var_x,"_Mean")]],
                      y=.data[[paste0(var_y,"_Mean")]],
                      color = essai_uk),
                  method=lm,
                  formula = y~x,
                  ## color="black",
                  se=F)
    + geom_smooth(data=databasemeansd,
                  aes(x=.data[[paste0(var_x,"_Mean")]],
                      y=.data[[paste0(var_y,"_Mean")]],
                      ),
                  method=lm,
                  formula = y~x,
                  color="black",
                  se=F)
    + stat_summary(
        fun = 'mean',
        fun.min = function(x) 0,
        geom = 'point', size = 6,
        aes(color=essai_uk))
    + stat_summary(
        fun = mean,
        fun.min = function(y) mean(y) - sd(y), 
        fun.max = function(y) mean(y) + sd(y), 
        geom ="pointrange", aes(color=essai_uk),
        show.legend = FALSE,linewidth=1)
    + scale_color_manual(
        values = c(viridis::viridis(4)[1:3]),
        name = "Trials",
        drop = FALSE
      )
    ## + geom_jitter(aes(group=essai_uk, color=essai_uk),
    ##               alpha=0.5,size=1, width=0.0005)
    ## + annotate("text",label = "r=0.925",x=0.05,y=0.3,linewidth=6,hjust = 0)
    ## + labs(
    ##     y=expression("W"["max"] - "W"["t0"] ~ "[-]"),
    ##     x="SOC:Clay ratio [-]")
    + theme_classic(base_size = 16)
  )
  return(g)
}