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

Table of Contents

1. Introduction

This is the public repository linked to the paper "The QuantiSlakeTest, dynamic weighting of soil under water to measure soil structural stability" submitted to the SOIL journal. You will find open data, meta data and commented R-scripts.

[to be amended]

This README file is done in a reproducible research spirit based on the powerful combo of emacs, org and babel. All the code blocks of this file have been extracted in one file, init.R in the root of the repository. The raw data are stored in the data-raw folder of the repository.

2. Initialisation

2.1. Cleaning session and workspace

Geek zonerm(list=ls())

2.2. Loading packages

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

2.3. Read data

2.3.1. QST data (1 line = 1 qst experiment)

Geek zoneqst.opendata <- read.csv(
  file = c(
    "./data-raw/qst-opendata.csv")
)
row.names(qst.opendata) <- qst.opendata$slake

qst.opendata$modalite_uk_abb <- factor(
  qst.opendata$modalite_uk_abb,
  levels = c(
    "P1 K0","P1 K1","P1 K2","P2 K0","P2 K1","P2 K2",
    "P","RT",
    "RE","FYM","RR"
  )
)

head(qst.opendata, n = 10) %>%
  mutate(across(where(is.double), signif, digits = 3))
  slake parcelle.x sample.id series.id duration Tmax Wt0 Wend Wmax.Wt0 Warchi Sl.0.max Sl.max.30 Sl.max.60 Sl.max.300 Sl.max.600 T25 T50 T75 T90 T95 T99 AUC campagne essai.x culture subzone X modalite essaib parcellea parcelleb plateforme traitement.x bloc ligne colonne traitment2 traitement3 traitement4 newbloc tmax t25 t50 t75 t90 t95 t99 essai.b parcelle.b parcelle.c sampleid corg humus pHeau sampleidgranulo clay siltf siltc silt sandf sandc sand ib nlisting essai.y parcelle.y traitement.y oc.clay tare tarepfrais tarepsec tarep105 pfrais psec p105 watercontent pctms watercontentcorr dapp rhob Wroot sample parcelle mode profondeur dryingtemp dryingj eau essaiuk modaliteuk modaliteukabb
382 382 39 ErnageA-39-PK-39-P2K0-2019-D-C Ernage A-39 971 4.69 0.94 0.061 0.06 0.42 0.013 -0.017 -0.011 -0.0029 -0.0015 21.4 32.5 67.2 146 360 785 126 Ernage A Ernage A Froment 39 5 P2 K0 PK PK39 E-39-2019 PK P2K0 2     P2 K0     4.69 21.4 32.5 67.2 146 360 785 PK 39 PK39 T193414 11.5 2.3 6.95 G190839 15 28 50.3 78.3 6.14 0.583 6.73 2.1 31 PK E-39-2019 P2 K0 0.077 8.96 164 136 137 155 127 128 18.2 0.986 18.5 1.27 1.25   PK-39-P2K0-2019-D-C 39 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K0 P2 K0
383 383 39 ErnageA-39-PK-39-P2K0-2019-D-B Ernage A-39 719 2.11 0.93 0.21 0.07 0.44 0.033 -0.014 -0.0091 -0.0026 -0.0013 19.9 29.4 96.5 265 298 331 211 Ernage A Ernage A Froment 39 5 P2 K0 PK PK39 E-39-2019 PK P2K0 2     P2 K0     2.11 19.9 29.4 96.5 265 298 331 PK 39 PK39 T193414 11.5 2.3 6.95 G190839 15 28 50.3 78.3 6.14 0.583 6.73 2.1 31 PK E-39-2019 P2 K0 0.077 8.96 164 136 137 155 127 128 18.2 0.986 18.5 1.27 1.25   PK-39-P2K0-2019-D-B 39 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K0 P2 K0
384 384 40 ErnageA-40-PK-40-P2K2-2019-D-B Ernage A-40 972 25.2 0.9 0.25 0.1 0.44 0.004 -0.0005 -0.0008 -0.0015 -0.0012 162 296 383 558 725 911 449 Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 PK P2K2 2     P2 K2     25.2 162 296 383 558 725 911 PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-B 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
386 386 16 ErnageA-16-PK-16-P1K0-2019-D-C Ernage A-16 971 12.4 0.93 0.64 0.07 0.43 0.0057 -0.0013 -0.0013 -0.0009 -0.00054 85.5 155 328 600 755 940 648 Ernage A Ernage A Froment 16 2 P1 K0 PK PK16 E-16-2019 PK P1K0 1     P1 K0     12.4 85.5 155 328 600 755 940 PK 16 PK16 T193411 11.9 2.38 6.86 G190836 19.9 27.5 47.4 74.9 5.06 0.214 5.27 1.76 28 PK E-16-2019 P1 K0 0.06 8.96 169 140 133 160 131 124 18.3 0.982 18.7 1.31 1.28   PK-16-P1K0-2019-D-C 16 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P1 K0 P1 K0
387 387 16 ErnageA-16-PK-16-P1K0-2019-D-B Ernage A-16 1130 32.5 0.93 0.62 0.07 0.44 0.0022 -0.002 -0.0017 -0.00093 -0.00059 86 195 349 569 726 881 642 Ernage A Ernage A Froment 16 2 P1 K0 PK PK16 E-16-2019 PK P1K0 1     P1 K0     32.5 86 195 349 569 726 881 PK 16 PK16 T193411 11.9 2.38 6.86 G190836 19.9 27.5 47.4 74.9 5.06 0.214 5.27 1.76 28 PK E-16-2019 P1 K0 0.06 8.96 169 140 133 160 131 124 18.3 0.982 18.7 1.31 1.28   PK-16-P1K0-2019-D-B 16 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P1 K0 P1 K0
388 388 41 ErnageA-41-PK-41-P1K2-2019-D-C Ernage A-41 726 11.8 0.85 0.25 0.15 0.43 0.013 -0.0022 -0.0045 -0.0024 -0.0012 65.8 77.7 177 286 319 471 275 Ernage A Ernage A Froment 41 7 P1 K2 PK PK41 E-41-2019 PK P1K2 2     P1 K2     11.8 65.8 77.7 177 286 319 471 PK 41 PK41 T193416 12.5 2.5 6.85 G190841 16.9 31.9 46.2 78.1 4.38 0.661 5.04 1.97 33 PK E-41-2019 P1 K2 0.074 8.03 176 143 143 168 135 135 19.5 0.986 19.8 1.35 1.33   PK-41-P1K2-2019-D-C 41 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P1 K2 P1 K2
391 391 40 ErnageA-40-PK-40-P2K2-2019-D-D Ernage A-40 1560 23.5 0.9 0.46 0.1 0.44 0.0043 -0.0024 -0.0033 -0.0015 -0.00082 72.6 91.5 252 535 941 1530 512 Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 PK P2K2 2     P2 K2     23.5 72.6 91.5 252 535 941 1530 PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-D 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
392 392 40 ErnageA-40-PK-40-P2K2-2019-D-E Ernage A-40 626 50.3 0.83 0.74 0.17 0.41 0.0034 -0.00091 -0.0012 -0.00082   110 150 275 297 440 473 508 Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 PK P2K2 2     P2 K2     50.3 110 150 275 297 440 473 PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-E 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
401 401 40 ErnageA-40-PK-40-P2K2-2019-D-F Ernage A-40 1500 25.3 0.87 0.49 0.13 0.43 0.0051 -0.0006 -0.00096 -0.0011 -0.00076 185 284 447 661 816 1190 597 Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 PK P2K2 2     P2 K2     25.3 185 284 447 661 816 1190 PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-F 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
403 403 39 ErnageA-39-PK-39-P2K0-2019-D-D Ernage A-39 1180 14.7 0.89 0.47 0.11 0.43 0.0075 -0.0024 -0.0038 -0.0017 -0.00088 57.7 84.9 133 198 308 429 468 Ernage A Ernage A Froment 39 5 P2 K0 PK PK39 E-39-2019 PK P2K0 2     P2 K0     14.7 57.7 84.9 133 198 308 429 PK 39 PK39 T193414 11.5 2.3 6.95 G190839 15 28 50.3 78.3 6.14 0.583 6.73 2.1 31 PK E-39-2019 P2 K0 0.077 8.96 164 136 137 155 127 128 18.2 0.986 18.5 1.27 1.25   PK-39-P2K0-2019-D-D 39 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K0 P2 K0

2.3.2. QST aggregated data (mean and sd by plot)

Geek zoneqst.mean.sd.opendata <- read.csv(
  file = c(
    "./data-raw/qst-mean-sd-opendata.csv"
  )
)

head(qst.mean.sd.opendata, n = 10) %>%
  mutate(across(where(is.double), signif, digits = 3))
  series.id parcelle.c bloc campagne modaliteukabb slakeMean slakeSD durationMean durationSD TmaxMean TmaxSD Wt0Mean Wt0SD WendMean WendSD Wmax.Wt0Mean Wmax.Wt0SD WarchiMean WarchiSD Sl.0.maxMean Sl.0.maxSD Sl.max.30Mean Sl.max.30SD Sl.max.60Mean Sl.max.60SD Sl.max.300Mean Sl.max.300SD Sl.max.600Mean Sl.max.600SD T25Mean T25SD T50Mean T50SD T75Mean T75SD T90Mean T90SD T95Mean T95SD T99Mean T99SD AUCMean AUCSD tmaxMean tmaxSD t25Mean t25SD t50Mean t50SD t75Mean t75SD t90Mean t90SD t95Mean t95SD t99Mean t99SD corgMean corgSD humusMean humusSD pHeauMean pHeauSD clayMean claySD siltfMean siltfSD siltcMean siltcSD siltMean siltSD sandfMean sandfSD sandcMean sandcSD sandMean sandSD ibMean ibSD nlistingMean nlistingSD oc.clayMean oc.claySD tareMean tareSD tarepfraisMean tarepfraisSD tarepsecMean tarepsecSD tarep105Mean tarep105SD pfraisMean pfraisSD psecMean psecSD p105Mean p105SD watercontentMean watercontentSD pctmsMean pctmsSD watercontentcorrMean watercontentcorrSD dappMean dappSD rhobMean rhobSD WrootMean WrootSD dryingjMean dryingjSD X essai traitement MWDT1 MWDT2 MWDT3 MAT1 MAT2 MAT3 MWDmean MWD2.MWD1 MWD3.MWD1 MWD3.MWD2
1 Ernage A-16 PK16 1 Ernage A P1 K0 413 23.8 1080 95.4 10.3 13.2 0.926 0.0472 0.342 0.279 0.074 0.0472 0.45 0.0158 0.0194 0.0158 -0.00666 0.00588 -0.00536 0.00402 -0.00197 0.00102 -0.00107 0.000474 49 35.3 103 68.5 194 134 378 199 539 216 763 177 386 250 10.3 13.2 49 35.3 103 68.5 194 134 378 199 539 216 763 177 11.9 0 2.38 0 6.86 0 19.9 0 27.5 0 47.4 0 74.9 0 5.06 0 0.214 0 5.27 0 1.76 0 28 0 0.06 0 8.96 0 169 0 140 0 133 0 160 0 131 0 124 0 18.3 0 0.982 0 18.7 0 1.31 0 1.28 0     30 0 2 PK P1 K0 0.264 1.22 0.768 0.443 0.802 0.816 0.75 4.62 2.91 0.631
2 Ernage A-2 PK02 1 Ernage A P1 K1 424 19.9 2000 2210 2.49 0.541 0.942 0.0386 0.214 0.227 0.0575 0.0386 0.435 0.0252 0.0236 0.0154 -0.0136 0.00427 -0.00872 0.00293 -0.0025 0.000707 -0.00143 0.000208 19.4 2.18 33.9 4.21 81.7 19.6 234 128 339 218 625 298 209 115 2.49 0.541 19.4 2.18 33.9 4.21 81.7 19.6 234 128 339 218 625 298 9.75 0 1.95 0 6.91 0 20 0 28.3 0 46 0 74.2 0 5.52 0 0.286 0 5.8 0 1.95 0 27 0 0.0488 0 8.94 0 171 0 142 0 151 0 162 0 133 0 142 0 17.7 0 0.981 0 18 0 1.33 0 1.31 0     30 0 1 PK P1 K1 0.225 0.997 0.537 0.358 0.762 0.751 0.586 4.44 2.39 0.538
3 Ernage A-20 PK20 1 Ernage A P1 K2 423 18.9 1200 239 5.95 6.94 0.93 0.0339 0.248 0.108 0.07 0.0339 0.448 0.0148 0.0199 0.0115 -0.0066 0.00297 -0.00566 0.00168 -0.00208 0.000327 -0.00113 0.000207 37.9 17 77.8 25.2 258 225 526 238 724 263 1070 161 352 81.5 5.95 6.94 37.9 17 77.8 25.2 258 225 526 238 724 263 1070 161 11.4 0 2.28 0 7 0 19.9 0 26.8 0 47.4 0 74.2 0 5.57 0 0.398 0 5.96 0 1.77 0 29 0 0.0574 0 8.96 0 171 0 141 0 140 0 162 0 132 0 131 0 18.1 0 0.983 0 18.5 0 1.32 0 1.3 0     30 0 3 PK P1 K2 0.279 1.17 0.762 0.407 0.803 0.828 0.737 4.2 2.73 0.651
4 Ernage A-35 PK35 2 Ernage A P2 K1 425 13.7 3170 2200 17.2 11.6 0.896 0.027 0.278 0.122 0.104 0.027 0.436 0.0152 0.0122 0.0147 -0.00266 0.000699 -0.00358 0.0015 -0.00224 0.000404 -0.00117 0.000189 72.4 24.9 115 29.6 210 72.6 264 75.8 352 58.8 1620 1630 351 101 17.2 11.6 72.4 24.9 115 29.6 210 72.6 264 75.8 352 58.8 1620 1630 10.8 0 2.17 0 6.61 0 16.2 0 27.9 0 49.8 0 77.7 0 5.61 0 0.513 0 6.12 0 2.09 0 30 0 0.0668 0 9 0 178 0 146 0 146 0 169 0 137 0 137 0 19 0 0.986 0 19.2 0 1.37 0 1.35 0     30 0 4 PK P2 K1 0.295 0.868 0.53 0.492 0.734 0.758 0.565 2.94 1.79 0.611
5 Ernage A-39 PK39 2 Ernage A P2 K0 395 11.8 1090 266 8.49 5.5 0.898 0.0432 0.276 0.222 0.102 0.0432 0.43 0.01 0.0168 0.0111 -0.00936 0.00627 -0.0075 0.00328 -0.00226 0.000666 -0.00116 0.000347 35.5 15.6 52.1 23.6 101 26.5 253 102 481 223 777 396 308 183 8.49 5.5 35.5 15.6 52.1 23.6 101 26.5 253 102 481 223 777 396 11.5 0 2.3 0 6.95 0 15 0 28 0 50.3 0 78.3 0 6.14 0 0.583 0 6.73 0 2.1 0 31 0 0.077 0 8.96 0 164 0 136 0 137 0 155 0 127 0 128 0 18.2 0 0.986 0 18.5 0 1.27 0 1.25 0     30 0 5 PK P2 K0 0.319 1.24 0.628 0.508 0.763 0.833 0.728 3.88 1.97 0.508
6 Ernage A-40 PK40 2 Ernage A P2 K2 392 6.98 1160 444 31.1 12.8 0.875 0.0332 0.485 0.201 0.125 0.0332 0.43 0.0141 0.0042 0.000707 -0.0011 0.000882 -0.00156 0.00117 -0.00123 0.000332 -0.000927 0.000239 133 50.9 205 101 339 91.6 513 154 730 213 1020 446 517 60.9 31.1 12.8 133 50.9 205 101 339 91.6 513 154 730 213 1020 446 11.7 0 2.35 0 6.8 0 16.6 0 30.3 0 47.3 0 77.5 0 5.16 0 0.651 0 5.81 0 2.02 0 32 0 0.0705 0 8.57 0 184 0 150 0 149 0 175 0 141 0 140 0 19.5 0 0.986 0 19.8 0 1.41 0 1.39 0     30 0 6 PK P2 K2 0.312 1.33 0.779 0.482 0.801 0.79 0.807 4.26 2.5 0.587
7 Ernage A-41 PK41 2 Ernage A P1 K2 414 20.7 1290 418 10.1 1.3 0.864 0.0472 0.302 0.1 0.136 0.0472 0.458 0.0192 0.0139 0.00565 -0.00678 0.00268 -0.00646 0.00184 -0.00208 0.000327 -0.00108 0.000187 40 14.8 61.7 17.2 141 55.9 418 239 651 228 1060 384 354 84.1 10.1 1.3 40 14.8 61.7 17.2 141 55.9 418 239 651 228 1060 384 12.5 0 2.5 0 6.85 0 16.9 0 31.9 0 46.2 0 78.1 0 4.38 0 0.661 0 5.04 0 1.97 0 33 0 0.074 0 8.03 0 176 0 143 0 143 0 168 0 135 0 135 0 19.5 0 0.986 0 19.8 0 1.35 0 1.33 0     30 0 7 PK P1 K2 0.297 1.64 0.642 0.526 0.834 0.844 0.859 5.51 2.16 0.392
8 Ernage A-47 PK47 2 Ernage A P1 K0 427 1.41 1140 157 31.1 8.34 0.82 0.0707 0.63 0.141 0.18 0.0707 0.44 0 0.0057 0.000707 -0.00175 0.000354 -0.00155 0.000636 -0.000985 0.000163 -0.000595 0.000205 80.6 76.9 127 143 277 193 444 304 510 256 655 182 634 46.1 31.1 8.34 80.6 76.9 127 143 277 193 444 304 510 256 655 182 11 0 2.21 0 6.69 0 12.8 0 28.7 0 52.4 0 81 0 5.76 0 0.413 0 6.18 0 2.36 0 34 0 0.0862 0 8.06 0 172 0 140 0 142 0 164 0 132 0 134 0 19.6 0 0.973 0 20.1 0 1.32 0 1.28 0     30 0 8 PK P1 K0 0.272 1.09 0.574 0.472 0.761 0.778 0.647 4.02 2.11 0.525
9 Ernage A-51 PK51 2 Ernage A P1 K1 434 14.5 967 165 15.5 5.03 0.86 0.0216 0.393 0.247 0.14 0.0216 0.442 0.025 0.00948 0.00256 -0.006 0.00428 -0.00595 0.00297 -0.0019 0.000821 -0.000995 0.000438 46.6 13.4 69.5 24.2 116 42.6 266 141 374 123 733 239 405 202 15.5 5.03 46.6 13.4 69.5 24.2 116 42.6 266 141 374 123 733 239 9.69 0 1.94 0 6.69 0 13 0 31.3 0 50.3 0 81.6 0 5.07 0 0.313 0 5.38 0 2.61 0 35 0 0.0743 0 8.09 0 166 0 137 0 139 0 158 0 128 0 131 0 18.5 0 0.987 0 18.7 0 1.28 0 1.27 0     30 0 9 PK P1 K1 0.271 1.15 0.519 0.509 0.781 0.772 0.647 4.24 1.91 0.451
10 Les Fonds-13NT LF13 5 Les Fonds RT 579 7.54 3150 2930 46.8 35.3 0.8 0.0548 0.76 0.0898 0.2 0.0548 0.47 0.00816 0.00628 0.00407 -0.00128 0.000483 -0.00121 0.000416 -0.000535 0.000222 -0.000332 0.000135 97.7 55.7 206 112 589 159 846 154 1160 450 1650 949 740 51.1 46.8 35.3 97.7 55.7 206 112 589 159 846 154 1160 450 1650 949 13 0 2.6 0 7.17 0 13.1 0 30.4 0 49.2 0 79.6 0 5.86 0 1.49 0 7.35 0 2.08 0 3 0 0.0994 0 8.08 0 159 0 131 0 129 0 150 0 123 0 121 0 18.4 0 0.988 0 18.6 0 1.23 0 1.21 0 55.8 10.2 30 0 30 Tvl du sol TCS 0.28 1.21 1 0.473 0.729 0.794 0.832 4.34 3.57 0.824

2.3.3. Raw data from one QST experiment

Geek zonedata.slake.example <- read.csv(
   file = c(
     "./data-raw/data-slake-example.csv"
   )
)

head(data.slake.example, n = 10) %>%
   mutate(across(where(is.double), signif, digits = 3))
  X.1 sample.id slake i campaign serie sample.x time mass mass.2 time.2 derivmass derivsecmass diffmass.full pdiffmass.full difftime.full difftime difftime.num diffmass pdiffmass sample.y campagne parcelle mode profondeur dryingtemp dryingj eau operator starttime device device2 essai culture series.id subzone X modalite essaib parcellea parcelleb corg a cn hum xxx1 xxx2 xxx3 xxx4 xxx5 xxx6 xxx7 plateforme traitement bloc ligne colonne traitment2 traitement3 traitement4 newbloc slake.sample
1 1 ErnageA-20-PK-20-P1K2-2019-D-F 410 11 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:22.739357 54.3 54.3 2019-07-01 10:46:22 -3.88 49.9 -78.7 0.408 8.62 0 0 -7.55 0.878 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
2 2 ErnageA-20-PK-20-P1K2-2019-D-F 410 12 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:23.496483 58.7 58.7 2019-07-01 10:46:23 4.42 8.3 -74.3 0.442 9.37 0.757 0.757 -3.13 0.949 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
3 3 ErnageA-20-PK-20-P1K2-2019-D-F 410 13 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:24.333913 60 60 2019-07-01 10:46:24 1.31 -3.11 -73 0.451 10.2 1.59 1.59 -1.82 0.971 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
4 4 ErnageA-20-PK-20-P1K2-2019-D-F 410 14 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:25.300481 59.9 59.9 2019-07-01 10:46:25 -0.09 -1.4 -73 0.451 11.2 2.56 2.56 -1.91 0.969 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
5 5 ErnageA-20-PK-20-P1K2-2019-D-F 410 15 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:26.260727 60.3 60.3 2019-07-01 10:46:26 0.33 0.42 -72.7 0.453 12.1 3.52 3.52 -1.58 0.974 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
6 6 ErnageA-20-PK-20-P1K2-2019-D-F 410 16 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:27.217738 60.5 60.5 2019-07-01 10:46:27 0.25 -0.08 -72.5 0.455 13.1 4.48 4.48 -1.33 0.978 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
7 7 ErnageA-20-PK-20-P1K2-2019-D-F 410 17 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:28.187114 60.8 60.8 2019-07-01 10:46:28 0.3 0.05 -72.2 0.457 14.1 5.45 5.45 -1.03 0.983 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
8 8 ErnageA-20-PK-20-P1K2-2019-D-F 410 18 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:29.140590 60.9 60.9 2019-07-01 10:46:29 0.08 -0.22 -72.1 0.458 15 6.4 6.4 -0.95 0.985 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
9 9 ErnageA-20-PK-20-P1K2-2019-D-F 410 19 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:30.096646 61.1 61.1 2019-07-01 10:46:30 0.21 0.13 -71.9 0.46 16 7.36 7.36 -0.74 0.988 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F
10 10 ErnageA-20-PK-20-P1K2-2019-D-F 410 20 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:30.901920 61.2 61.2 2019-07-01 10:46:30 0.15 -0.06 -71.7 0.461 16.8 8.16 8.16 -0.59 0.99 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       PK P1K2 1     P1 K2     410 - PK-20-P1K2-2019-D-F

3. Personal functions

3.1. Illustrative plot of one QST curve (bw,no data)

Geek zoneillustr.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", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't25'], col = "gray15", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't50'], col = "gray15", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't75'], col = "gray15", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't90'], col = "gray15", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't95'], col = "gray15", size = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't99'], col = "gray15", size = 0.3)
    + geom_hline(yintercept = 1, col = "gray15", size = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wend'], col = "black", size = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray15", size = 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', size = 3)
    + geom_point(x = 0 , y = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray89", size = 3)
    ## + annotate("label",x = 0 ,
    ##              y = (my.qst.opendata[as.character(id.slakes), 'Wt0']-0.047),
    ##              label.size = 0, fill = "white",
    ##              label = "0",
    ##              col = "gray15", size = 3, vjust = 0, hjust = 1.5 )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 'tmax'], y = 1, col = "gray89", size = 3)
    ## + annotate("label",x = my.qst.opendata[as.character(id.slakes), 'tmax'],
    ##             y = 1.01,
    ##             label.size = 0, fill = "white",
    ##             label = "max", col = "gray15", size = 3, vjust = 0, hjust = -0.1
    ##             )
    ## + 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'] + 30)][1],
    ##                y = 1,
    ##                yend = my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 30)][1],
    ##                col = "gray89",
    ##                size = 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'] + 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",
                   size = 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",
                   size = 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",
                   size = 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",
                   size = 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",
                   size = 0.3
                   )
    + annotate(geom = "label", -10,
               y = 1.1,
               ## label = paste("Slope[0-max]==", my.qst.opendata[as.character(id.slakes), 'Sl.0.max'],"~sec^{-1}"),
               label = "Slope[0-max]",
               parse = TRUE,
               col = "gray15", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               size = 3, hjust = 1
               )
    + geom_segment(x = -90,
                   xend = -90,
                   y = 1,
                   yend = my.qst.opendata[as.character(id.slakes), 'Wt0'],
                   col = "gray15",
                   size = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + geom_point(x = -90,
                 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 = paste("t[max]==", round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),"~sec"),
               label = "t[max]",
               size = 3,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.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", 
               size = 3,
               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",
               size = 3,
               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", 
               size = 3,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't90'],
                 y = 0.15,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't90'],
               y = 0.15,
               label = "t90",
               size = 3,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't95'],
                 y = 0.12,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't95'],
               y = 0.12,
               label = "t95", 
               size = 3,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't99'],
                 y = 0.10,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't99'],
               y = 0.10,
               label = "t99",
               size = 3,
               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]",
               size = 3,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.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]",
               size = 3,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.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'] + 30)][1],
    ##              y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 30)][1],
    ##              col = "gray89", size = 3)
    ## + annotate("label",
    ##            x = my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 30)][1],
    ##            y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 30)][1],
    ##            label.size = 0, fill = "white",
    ##            label = "t[max]+30", hjust = -0.1, vjust = -0.2, size=3,
    ##            parse=TRUE,
    ##            col = "gray16")
    ## + annotate("label",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],
    ##              label.size = 0, fill = "white",
    ##              label = "t[max]+60", hjust = -0.1, vjust = -0.2, size=3,
    ##              parse=TRUE,
    ##              col = "gray16")
    + 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)
    ## + annotate("label",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],
    ##             label.size = 0, fill = "white",
    ##             label = "t[max]+300", hjust = -0.1, vjust = -0.2, size=3,
    ##             parse=TRUE,
    ##             col = "gray16")
    + 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)
    ## + annotate("label",
    ##            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],
    ##            label.size = 0, fill = "white",
    ##            label = "t[max]+600", hjust = -0.1, vjust = -0.2, size=3,
    ##            parse=TRUE,
    ##            col = "gray16")
    + 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 = -100,
               y = mean(c(1,my.qst.opendata[as.character(id.slakes), 'Wt0'])),
               label = "W[max-Wt0]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               size = 3, 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')]",
               size = 3,
               parse = TRUE,
               col = "gray90",
               ## fill = "lightgrey",
               fill = "gray10",
               hjust = 1.05,
               vjust = -0.05
               )
    + 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", label.size = 1,
               size = 3, hjust = -0.05, vjust = -0.05
               )
    ## + annotate(geom = "text",
    ##            x = 500,
    ##            y = 0.75,
    ##             ## 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'] + 30)][1]))),
    ##            ## y = mean(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 30)][1])),
    ##            label = paste("Slope[max-30]==", my.qst.opendata[as.character(id.slakes), "Sl.max.30"],"~sec^{-1}"),
    ##            parse = TRUE,
    ##            col = "gray16",
    ##            size = 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", label.size = 1,
               size = 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", label.size = 1,
               size = 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",
               size = 3, hjust = 0
               )
    ## 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')",
               parse = TRUE,
               size = 5, 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"), label.size = 1,
               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", label.size = 1,
               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",
               size = 3, 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)
}

3.2. Matrix of boxplots

Geek zoneplot.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()
    )
}

3.3. Statistical analyses, model, tests, ANOVA and comparisons of means

Geek zonef.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()
}

4. Figures and Tables of the article

4.1. Material & Methods

4.1.1. Figure 1 - Material & Methods - Illustration of the material (see montage photo)

illustration_device.png

4.1.2. Figure 2 - Material & Methods - Illustration of a QST curve and main indicators

Geek zonefig2 <- illustr.qst.curve.indic.nb.nd(
  data.slake.example,
  qst.opendata
)
ggsave("./fig/fig02.pdf",fig2, width=2000, height=1750, units="px")

illustration_indicators.png

4.2. Results

4.2.1. Table 1 - Results - Soil properties

  1. Data cooking
    Geek zonedf.soil.prop <- qst.opendata %>%
  dplyr::select(
           campagne,
           essai_uk,
           modalite_uk,
           parcelle,
           clay,
           silt_f,
           silt_c,
           silt,
           sand_f,
           sand_c,
           sand,
           corg,
           oc.clay,
           pH_eau,
           rho_b,
           ) %>%
    dplyr::group_by(essai_uk,parcelle) %>%
  dplyr::slice_head() %>%
  dplyr::ungroup() %>%
  dplyr::mutate(plot = 1:length(parcelle)) %>%
  dplyr::rename(
           Treatment = modalite_uk,
           Plot = plot,
           Clay = clay,
           "Silt (fine)" = silt_f,
           "Silt (coarse)" = silt_c,
           "Silt (total)" = silt,
           "Sand (fine)" = sand_f,
           "Sand (coarse)" = sand_c,
           "Sand (total)" = sand,    
           "Organic C (g.kg^-1)" = corg,
           "SOC:Clay" = oc.clay,
           "pH [-]" = pH_eau,
           "Bulk density [-]" = rho_b
         )

write.csv(df.soil.prop,file="./data-output/check_up_table_soil_prop.csv")

df.soil.prop <- df.soil.prop %>%
  dplyr::select(-c(campagne,essai_uk,parcelle))
  2. Table as an image with gt
    Geek zone    tab1 <- df.soil.prop %>%
      mutate(empty = "") %>%
      gt() %>%
      tab_spanner(
        label = md("Texture (%)"),
        columns = names(df.soil.prop)[2:8]
      ) %>%
      cols_move(
        columns = "empty",
        after = "Treatment"
      ) %>%
      cols_move(
        columns = "Plot",
        after = "empty"
      ) %>%
      cols_move(
        columns = "Treatment",
        after = "Plot"
      ) %>%
      tab_row_group(
        label = "P-K mineral fertiliser",
        rows = 27:35
      ) %>%
      tab_row_group(
        label = "Tillage",
        rows = 19:26
      ) %>%
      tab_row_group(
        label = "Organic matter",
        rows = 1:18
      ) %>%
      fmt_number(
        columns =  c("Clay","Silt (fine)","Silt (coarse)","Silt (total)",
                     "Sand (fine)",
                     "Sand (coarse)",
                     "Sand (total)"),
             decimals = 1
           ) %>%
      fmt_number(
             columns = "SOC:Clay",
       decimals = 3
      ) %>%
          fmt_number(
             columns = c("pH [-]","Bulk density [-]"),
             decimals = 2
           ) %>%
      tab_options(table.font.size = 10)
gtsave(tab1,"./fig/tab01.pdf", zoom = 1)

gtsave(tab1,"./fig/dataframe_soil_properties.png",
             expand=30)

    dataframe_soil_properties.png

4.2.2. Table 2 - Results - Correlation Matrix QST v.s. Le Bissonnais & Soil Properties

  1. Data cooking
    Geek zonedf.comp <- qst.mean.sd.opendata %>%
  dplyr::select(
           MWD_T1,
           MWD_T2,
           MWD_T3,
           MA_T1,
           MA_T2,
           MA_T3,
           ## MWD2.MWD1,
           ## MWD3.MWD1,
           ## MWD3.MWD2,
           corg_Mean,
           clay_Mean,
           oc.clay_Mean,
           pH_eau_Mean,
           rho_b_Mean,
           Wmax.Wt0_Mean,
           tmax_Mean,
           Sl.0.max_Mean,
           Sl.max.30_Mean,
           Sl.max.60_Mean,
           Sl.max.300_Mean,
           Sl.max.600_Mean,
           t25_Mean,
           t50_Mean,
           t75_Mean,
           t90_Mean,
           t95_Mean,
           t99_Mean,
           Wend_Mean,
           AUC_Mean
         ) %>%
  dplyr::rename(
           "MWD 1" = MWD_T1,
           "MWD 3" = MWD_T2, ## NB INVERSION 2 et 3 !!!
           "MWD 2" = MWD_T3, ## NB INVERSION 2 et 3 !!!
           "MA 1" = MA_T1,
           "MA 3" = MA_T2, ## NB INVERSION 2 et 3 !!!
           "MA 2" = MA_T3, ## NB INVERSION 2 et 3 !!!
           ## "MWD3:MWD1" = MWD2.MWD1,
           ## "MWD2:MWD1" = MWD3.MWD1,
           ## "MWD2:MWD3" = MWD3.MWD2,
           SOC = corg_Mean,
           Clay = clay_Mean,
           "SOC:Clay" = oc.clay_Mean,
           pH = pH_eau_Mean,
           "Bulk density" = rho_b_Mean,
           "Slope 0-max" = Sl.0.max_Mean,
           tmax = tmax_Mean,
           "Wmax-Wt0" = Wmax.Wt0_Mean,
           "Slope max-30" = Sl.max.30_Mean,
           "Slope max-60" = Sl.max.60_Mean,
           "Slope max-300" = Sl.max.300_Mean,
           "Slope max-600" = Sl.max.600_Mean,
           t25 = t25_Mean,
           t50 = t50_Mean,
           t75 = t75_Mean,
           t90 = t90_Mean,
           t95 = t95_Mean,
           t99 = t99_Mean,
           "Wend" = Wend_Mean,
           AUC = AUC_Mean
         )

df.comp.num <- df.comp[,sapply(df.comp, is.numeric)] 

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

df.corr <- as.data.frame(round(mcor,3))[
  c(
    "Slope 0-max",
    "tmax",
    "Wmax-Wt0",
    "Slope max-30",
    "Slope max-60",
    "Slope max-300",
    "Slope max-600",
    "t25",
    "t50",
    "t75",
    "t90",
    "t95",
    "t99",
    "Wend",
    "AUC"
    ),
  c(
    "MWD 1",
    "MWD 2",
    "MWD 3",
    "MA 1",
    "MA 2",
    "MA 3",
    "SOC",
    "Clay",
    "SOC:Clay",
    "pH",
    "Bulk density"
  )
]
  2. Table as an image with gt
    Geek zonetab2 <- df.corr %>%
  dplyr::mutate(" " = "    ") %>%
  dplyr::add_rownames("qst") %>%
  gt(rowname_col = "qst") %>%
  data_color(
    columns = names(df.corr),
    colors = scales::col_numeric(
                       palette = RColorBrewer::brewer.pal(n=8, name="RdYlBu"),
                       domain = c(-1,1)
                     )) %>%
  tab_spanner(
    label = md("Le Bissonnais *et al.* (1996)"),
    columns = names(df.corr)[1:6]
  ) %>%
  tab_spanner(
    label = md("Soil properties"),
    columns = names(df.corr)[7:11]
  ) %>%
  cols_move(
    columns = " ",
    after = "MA 3"
  ) %>%
  tab_row_group(
    label = "QST indicators",
    rows = row.names(df.corr)
  )

gtsave(tab2,"./fig/tab02.pdf", zoom = 1)

gtsave(tab2,"./fig/dataframe_cor_qst_biss_soil.png",
         expand=30)

    dataframe_cor_qst_biss_soil.png

4.2.3. Figure 3 - Results - Scatterplot correlation SOC:CLay X Wmax.Wt0

Geek zoneqst.opendata$essai_uk <- factor(
  qst.opendata$essai_uk,
  levels = c("Organic matter","Tillage","P-K mineral fertiliser")
)

fig3 <- (ggplot(data=qst.opendata,
             aes(x=oc.clay,
                 y=Wmax.Wt0)
             )
  + geom_smooth(data=qst.mean.sd.opendata,
                aes(x=oc.clay_Mean,
                    y=Wmax.Wt0_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,size=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,size=6,hjust = 0)
  + labs(
      y=expression("W"["max"] - "W"["t0"] ~ "[-]"),
      x="SOC:Clay ratio [-]")
  + theme_classic(base_size = 16)
)
ggsave("./fig/fig03.pdf",fig3, width=3000, height=2000, units="px")
print(fig3)

xy_wmax_wt0_oc_clay_cropland.png

4.2.4. Figure 4 - Results - Comparison treatments : Organic Matter

Geek zoneg.bp.om <- plot.compare.trt("Organic matter", my.color = viridis::viridis(4)[2])
ggsave("./fig/fig04.pdf",g.bp.om, width=2300, height=1700, units="px")
print(g.bp.om)

boxplot_practices_OM.png

4.2.5. Figure 5 - Results - Curves and boxplot Tillage and Wend

Illustrative figure done using function from the "slaker" R-package, under devlopment : https://frdvnw.gitlab.io/slaker/dev/. Real data from the QST curves of the 3 long term field trials.

illustr_curves_wend_tillage.png

4.2.6. Figure 6 - Results - Comparison treatments : Tillage

Geek zoneg.bp.t <- plot.compare.trt("Tillage", my.color = viridis::viridis(4)[3])
ggsave("./fig/fig06.pdf",g.bp.t, width=2300, height=1700, units="px")
print(g.bp.t)

boxplot_practices_tillage.png

4.3. Discussion

4.3.1. Figure 7 - Discussion - Root Weights

  1. Data cooking
    Geek zoneqst.opendata.roots <- qst.opendata %>%
  dplyr::filter(essai_uk == "Tillage") %>%
  dplyr::select(modalite_uk, Wroot, parcelle,
                Wmax.Wt0,
                Sl.max.30,
                Sl.max.60,
                Sl.max.300,
                Sl.max.600,              
                Wend
                ) 

df.roots.long <- qst.opendata.roots %>%
  tidyr::pivot_longer(
           !c(
              modalite_uk, Wroot, parcelle
            ),
           names_to = "indicator",
           values_to = "value"
         )

df.roots.long$indicator <- factor(
  df.roots.long$indicator,
  levels=c(
    "Wmax.Wt0",
    "Sl.max.30",
    "Sl.max.60",
    "Sl.max.300",
    "Sl.max.600",
    "Wend"
  )
)

levels(df.roots.long$indicator) <- c(
  expression("(a) W"["max"]~"-W"["t0"]~"[-]"),
  expression("(b) Slope"["max-30"]~"[sec"^{-1}~"]"),
  expression("(c) Slope"["max-60"]~"[sec"^{-1}~"]"),
  expression("(d) Slope"["max-300"]~"[sec"^{-1}~"]"),
  expression("(e) Slope"["max-600"]~"[sec"^{-1}~"]"),
  expression("(f) W"["end"]~"[-]")
)
  2. Graph
    Geek zonefig7 <- (ggplot(data=df.roots.long,
             aes(x=Wroot,
                 y=value,
                 color=modalite_uk))
  + geom_point()
  + geom_smooth(method="lm")
  + ggpmisc::stat_poly_eq(
               aes(label = paste(after_stat(rr.label), sep = "**\", \"**")),
               label.x = "right",
               label.y = "bottom",
               vstep=0.07, size = 3
             )
  + facet_wrap(.~indicator,
               nrow=2, scales="free_y",
               ## labeller = labeller(indicator = as_labeller(indic_names,  label_parsed)),
               labeller = label_parsed
               )
  + scale_color_manual(values = c(viridis::turbo(8)[c(2,7)]), name="Treatment")
  + labs(title = paste0("Root biomass (mg)"),
         x = "Root biomass",
         y = "Value of the QuantiSlakeTest indicators")
  + theme_classic()
)
ggsave("./fig/fig07.pdf",fig7, width=3000, height=1500, units="px")
print(fig7)

    scatter_root_weight_indics.png

5. Statistical analyses

5.1. Calculation of some particular value

5.1.1. Correlation between mean value (oc:clay / Wmax-Wt0 indic)

Geek zoneround(cor(na.omit(qst.mean.sd.opendata[c("oc.clay_Mean","Wmax.Wt0_Mean")])),3)
  oc.clayMean Wmax.Wt0Mean
oc.clayMean 1 0.925
Wmax.Wt0Mean 0.925 1

5.1.2. Root biomass - mean calcultation by treatment

Geek zone## Exploration racines - brieuc
qst.opendata.roots %>%
  dplyr::group_by(modalite_uk) %>%
  dplyr::summarize(mean = round(mean(Wroot),2),
                   sd = round(sd(Wroot),2)
                   )
  modaliteuk mean sd
1 Ploughing 30.56 16.28
2 Reduced tillage 41.78 18.86

5.2. Linear Mixed Effect Model & ANOVA

5.2.1. Tillage trial (i.e. Les Fonds)

  1. Wmax.Wt0
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","Wmax.Wt0")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 10.4 1 2.09 0.0797
  2. Sl.max.30
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","Sl.max.30")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 1.04 1 2.27 0.404
  3. Sl.max.60
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","Sl.max.60")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 24.9 1 2.24 0.0301
  4. Sl.max.300
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","Sl.max.300")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 69.1 1 2.31 0.00905
  5. Sl.max.600
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","Sl.max.600")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 64.7 1 2.36 0.00914
  6. t75
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","t75")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 2.93 1 2.42 0.207
  7. t95
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","t95")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.00422 1 2.17 0.954
  8. tmax
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","tmax")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 4.33 1 2.85 0.134
  9. AUC
    Geek zonebroom::tidy(f.qst.anova(
         "Les Fonds","AUC"
       )[["Anova"]]
       ) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 66 1 2.27 0.0101
  10. Wend
    Geek zonebroom::tidy(f.qst.anova("Les Fonds","Wend")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 64.6 1 2.48 0.00781

5.2.2. Organic matter trial (i.e. Longs Tours)

  1. Wmax.Wt0
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","Wmax.Wt0")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 1.42 2 9.89 0.287
  2. Sl.max.30
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","Sl.max.30")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.923 2 9.93 0.429
  3. Sl.max.60
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","Sl.max.60")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 1.2 2 9.93 0.34
  4. Sl.max.300
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","Sl.max.300")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 2.74 2 9.88 0.113
  5. Sl.max.600
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","Sl.max.600")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 2.82 2 9.82 0.108
  6. t75
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","t75")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.957 2 9.84 0.417
  7. t95
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","t95")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.795 2 9.84 0.479
  8. tmax
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","tmax")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 4.24 2 9.83 0.047
  9. AUC
    Geek zonebroom::tidy(f.qst.anova(
         "Longs Tours","AUC"
       )[["Anova"]]
       ) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 1.95 2 9.92 0.193
  10. Wend
    Geek zonebroom::tidy(f.qst.anova("Longs Tours","Wend")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 2.96 2 9.91 0.0982

5.2.3. P-K mineral fertiliser trial (i.e. Ernage A)

  1. Wmax.Wt0
    Geek zonebroom::tidy(f.qst.anova("Ernage A","Wmax.Wt0")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.716 5 1.34 0.692
  2. Sl.max.30
    Geek zonebroom::tidy(f.qst.anova("Ernage A","Sl.max.30")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.253 5 1.34 0.9
  3. Sl.max.60
    Geek zonebroom::tidy(f.qst.anova("Ernage A","Sl.max.60")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.129 5 0.725 0.957
  4. Sl.max.300
    Geek zonebroom::tidy(f.qst.anova("Ernage A","Sl.max.300")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.743 5 1.28 0.685
  5. Sl.max.600
    Geek zonebroom::tidy(f.qst.anova("Ernage A","Sl.max.600")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.241 5 1.12 0.904
  6. t75
    Geek zonebroom::tidy(f.qst.anova("Ernage A","t75")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.0172 5 -0.109  
  7. t95
    Geek zonebroom::tidy(f.qst.anova("Ernage A","t95")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.0163 5 0.0363 0.994
  8. tmax
    Geek zonebroom::tidy(f.qst.anova("Ernage A","tmax")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 1.22 5 1.54 0.539
  9. AUC
    Geek zonebroom::tidy(f.qst.anova(
         "Ernage A","AUC"
       )[["Anova"]]
       ) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.443 5 1.41 0.802
  10. Wend
    Geek zonebroom::tidy(f.qst.anova("Ernage A","Wend")[["Anova"]]) %>%
  mutate(across(where(is.double), signif, digits = 3))
      term statistic df Df.res p.value
    1 trmt 0.582 5 1.29 0.743

5.2.4. Estimated marginal means & post-hoc comparisons

Organic matter trial only on indicators that show significant differences here before

  1. tmax
    Geek zonef.qst.means.comp("Longs Tours","tmax")
      contrast estimate SE df t.ratio p.value
    1 FYM - RE 2.95 1.2 9.44 2.46 0.0821
    2 FYM - RR -0.204 1.23 10 -0.166 0.985
    3 RE - RR -3.16 1.23 10 -2.57 0.0661
  2. Wend
    Geek zonef.qst.means.comp("Longs Tours","Wend")
      contrast estimate SE df t.ratio p.value
    1 FYM - RE -0.0524 0.0557 9.59 -0.941 0.629
    2 FYM - RR -0.137 0.0565 10.1 -2.42 0.0844
    3 RE - RR -0.0842 0.0565 10.1 -1.49 0.336

5.3. Summary tables

5.3.1. Define variables for the loops

Geek zoneall.campagne <- c("Les Fonds","Longs Tours")
all.indics <- c(
  "Wmax.Wt0",
  "tmax",
  "Sl.max.30",
  "Sl.max.60",
  "Sl.max.300",
  "Sl.max.600",
  "t75",
  "T95",
  "Wend",
  "AUC"
)
  x
1 Wmax.Wt0
2 tmax
3 Sl.max.30
4 Sl.max.60
5 Sl.max.300
6 Sl.max.600
7 t75
8 T95
9 Wend
10 AUC

5.3.2. Linear Mixed Effect Model & ANOVA

  1. Data cooking
    Geek zonet.pvalues.lmm <- cbind(data.frame(
  campagne="Les Fonds", indic="Sl.max.30"),
  f.qst.anova("Les Fonds","Sl.max.30")[["Anova"]])[0,]
for (my.indic in all.indics) {
  for (my.campagne in all.campagne){
    t.pvalues.lmm <- rbind(t.pvalues.lmm,
                       cbind(data.frame(campagne=my.campagne,indic=my.indic),f.qst.anova(my.campagne,my.indic)[["Anova"]])
                       )
  }
}

t.pvalues.lmm$indic.signif <- stats::symnum(
                                   t.pvalues.lmm[,"Pr(>F)"],
                                   cutpoints = c(0,0.001,0.01,0.05,0.1,1),
                                   symbols = c("***","**","*",".","n.s."),
                                   na = "") %>%
  as.character()

t.pvalues.lmm %>%
  dplyr::arrange(campagne,indic) %>%
  mutate(across(where(is.numeric), signif, digits = 3))  %>%
  dplyr::filter(indic.signif != "n.s.")
      campagne indic F Df Df.res Pr(>F) indic.signif
    trmt18 Les Fonds AUC 66 1 2.27 0.0101 *
    trmt8 Les Fonds Sl.max.300 69.1 1 2.31 0.00905 **
    trmt6 Les Fonds Sl.max.60 24.9 1 2.24 0.0301 *
    trmt10 Les Fonds Sl.max.600 64.7 1 2.36 0.00914 **
    trmt16 Les Fonds Wend 64.6 1 2.48 0.00781 **
    trmt Les Fonds Wmax.Wt0 10.4 1 2.09 0.0797 .
    trmt3 Longs Tours tmax 4.24 2 9.83 0.047 *
    trmt17 Longs Tours Wend 2.96 2 9.91 0.0982 .
  2. Table as an image with gt
    Geek zone
t.pvalues.lmm |> 
  dplyr::select(-indic.signif,-F,-Df,-Df.res) |>
  tidyr::pivot_wider(names_from=indic,
                     values_from=c("Pr(>F)")) |>
  mutate(across(where(is.double), signif, digits = 3)) |>
  gt() |> 
  data_color(
    columns = c(AUC,Sl.max.30,Sl.max.300,Sl.max.60,Sl.max.600,
                t75,T95,tmax,Wend,Wmax.Wt0
                ),
    colors = scales::col_numeric(
                       palette = rev(RColorBrewer::brewer.pal(n=6, name="Oranges")),
                       domain = c(0, 0.1)
                     )
  ) |>
  gtsave("./fig/all_p_values_lmm.png",
         expand=30)

    all_p_values_lmm.png

6. Supplementary figures not shown in the article

6.1. Results comparison treatments : PK full modality

Geek zoneg.bp.pk <- plot.compare.trt("P-K mineral fertiliser", my.color = viridis::viridis(4)[1])
print(g.bp.pk)

boxplot_practices_pk_full.png

6.2. Results comparison treatments : PK - K0-1-2

Geek zoneqst.opendata$modalite_uk <- as.character(qst.opendata$modalite_uk)
qst.opendata$modalite_uk[grepl("P[1,2] K[0,1,2]",qst.opendata$modalite_uk)] <-
  substr(qst.opendata$modalite_uk[grepl("P[1,2] K[0,1,2]",qst.opendata$modalite_uk)],
         start=4,stop=5)

qst.opendata$modalite_uk <- factor(
  qst.opendata$modalite_uk,
  levels = c(
    "K0","K1","K2",
    "Ploughing","Reduced tillage",
    "Farmyard manure", "Residue restitution", "Residue exportation"
  )
)


qst.opendata$modalite_uk_abb[grepl("P[1,2] K[0,1,2]",
                                   qst.opendata$modalite_uk_abb)] <-
  substr(qst.opendata$modalite_uk_abb[grepl("P[1,2] K[0,1,2]",
                                            qst.opendata$modalite_uk_abb)],
         start=4,stop=5)

g.bp.k <- plot.compare.trt("P-K mineral fertiliser", my.color = viridis::viridis(4)[1])
print(g.bp.k)

boxplot_practices_pk_K.png

7. Session info

The scripts have been run the Thu Oct 13 11:48:39 2022 in R version 4.2.1 (2022-06-23) on a platform x86_64-pc-linux-gnu (64-bit).