Phosphorus fertilisation rules particularly rely on a diagnosis step that consists of comparing the level of plant available P in the soil by analysis with a "critical P threshold value (P crit ). P crit are generally calculated from data obtained on Long Term Fertilisation Experiments (LTFEs), providing yield vs available-P relationships. According to the statistical functions used to fit such relationships, the P crit may vary up to two-fold, therefore introducing a tremendous variability in recommendations between users of different methods. The purpose of this study is to provide objective arguments to help choose the "right" model that can handle data sets of varying quality. Our analysis is based on six LTFEs providing complete environmental description, regular Olsen P measurements of available P, and annual yields for wheat, maize, barley and durum wheat. The dataset is composed of 4500 point pairs "yield-Olsen P". Four models were tested: Linear Plateau (LP), Quadratic-Plateau (QP), Mitscherlich 95% (Mi95) and Cate-Nelson clustering algorithm (CN). The model assessment was mainly based on the three following criteria: quality of adjustment fitting (RMSE), risk of high interannual variability (Coefficient of Variation, CV), risk of providing out of range P crit , standing for accuracy, certainty and robustness, respectively. The P crit calculation of a given crop on a given trial was tested using either pooled data relying on relative yields, or by averaging on annual trials. Our results clearly show that the two first indicators (RMSE, CV) could not discriminate between models. However, the ability of models to cope with scarce data strongly favoured the LP function. We also proved that using pooled relative yields or averaging annual trails not only led to the same conclusions, but provided also the same orders of magnitude of P crit . As mentioned in the literature, we confirmed the strong impact of the method on the final P crit , with CN and QP functions exhibiting the lowest and highest values, respectively. Finally, we recommend the use of the LP method, which combines many advantages: simple algorithm, strong ability to provide sensible thresholds even on trials relying on few data, equivalent model fitting and interannual variations than more complex models, no need to rely on an arbitrary yield gap, and intermediate values of thresholds. The relative yield method should be preferred wherever annual data are too sparse for model fitting; alternatively, averaging annual data seemed simpler and not subject to the arbitrary choice of a "reference yield".