4.2.2. Goodness of Fit and Prediction Including "dataset2"

To assess the goodness of fit of the selected models, we explored the uncertainty of parameters and the prediction quality by comparing the predicted and observed values of *Vpt*. We evaluated the prediction quality using two methods: cross-validation and training validation. In both methods, model performance was evaluated on data not included in the sample used to estimate model parameters.

For increasing the domain where the uncertainty of the parameters and the prediction quality of the models could be explored, a larger dataset ("dataset2") was added to the first dataset ("dataset1") with which the model parameters were estimated.

In the cross-validation method, observations of the hives of a given apiary were removed from the database, the model was fitted to the remaining data, and estimated parameters were plugged in to predict *Vpt* for the hives of the apiary whose observations were removed (this case corresponds to predicting *Vpt* for a new apiary based on observations collected from other apiaries). This procedure was repeated for each apiary of the dataset (i.e., 54 times for model A and 40 times for model B) and allowed us to provide averaged cross-validation assessments of the prediction performance.

The prediction performance was assessed with respect to the two following criteria:


These criteria (CI and quantiles) were empirically calculated from 1000 simulations of the zero-inflated beta distribution in which the estimated values of μ, σ, and ν were inserted (random factors incorporated into μ, σ, and ν were randomly drawn at each simulation from centered normal distributions with standard deviations equal to their estimated values). Note that estimation uncertainty was neglected in this simulation procedure; this choice may lead to un-calibrated confidence intervals and quantiles. The comparisons of quantiles Q97.5%, Q85%, Q75%, and Q50% with the threshold of 3 Varroa mites for 100 bees can be used as indicators to assess whether *Vpt* will exceed this problematic threshold. The efficiency of these indicators was assessed with the error rate *τerror*(*α*) calculated as the ratio between (i) the number of hives for which the predicted quantile Qα(*Vpt*) is less than or equal to 3 at time *t*, whereas the actual observation *Vpt* is greater than 3, and (ii) the number of observations *Vpt* greater than 3:

$$\tau^{error}(a) = \frac{\sum\_{R=1}^{K} \sum\_{r=1}^{N\_R} 1(\lfloor Vp\_{t,R,r} > 3 \rfloor \ Q\_a^{-R}(\lfloor Vp\_{t,R,r} \rfloor \le 3) }{\sum\_{R=1}^{K} \sum\_{r=1}^{N\_R} 1(\lfloor Vp\_{t,R,r} > 3 \rfloor \tag{4}$$

where *K* is the number of apiaries (54 for model A and 40 for model B), *NR* is the number of hives in the apiary *R* (which ranges between 7 and 51), *Vpt*,*R*,*<sup>r</sup>* is the observed value of phoretic Varroa at time *t* for the hive *r* of the apiary *R*, and *Q*−*<sup>R</sup> <sup>α</sup>* ( *V pt*,*R*,*r*) is the predicted quantile at α% of *Vpt,R,r* for the hives of the apiary *R* whose observations were removed (−*R*). The indicator function E→**1I** takes the value of 1 if event E is true, or otherwise 0.

Equation (4) presents the case in which observations are greater than 3 and predictions are less than or equal to 3, and *τerror*(*α*) was also calculated when observations are less than or equal to 3 and predictions are greater than 3. The error rate *τerror*(*α*) can be computed in other specific conditions, for example, conditions related to the number of phoretic Varroa at *t*−x (*Vpt–x*,*R*,*r*).

In the training validation method, observations of the hives of all apiaries after a specific date, *tA* = 31 (15 April) for model A, and *tB* = 118 (11 July) for model B, were removed from the database, the model was fitted to remaining data, and estimated parameters were plugged in to predict *Vpt* for the hives of all apiaries after *tA* or *tB* (this case corresponds to predicting *Vpt* for an apiary already installed, based on observations collected beforehand from this apiary). This procedure was repeated for each year of the dataset (i.e., from 2014 to 2016 and 2018 for both models A and B) and allowed us to provide averaged training validation assessments of the prediction performance already introduced in the cross-validation approach.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/pathogens10060678/s1, Table S1: Comparisons of the tested models investigating the influence of phoretic Varroa mite (per 100 bees) numbers, capped brood cell numbers, varbrood, and date on next phoretic Varroa numbers as a function of the estimation length, using the AICc criterion. N = 867 for data adjustment at one month (x = 1) and N = 93 for data adjustment at three months (x = 3), Table S2: R graphic output of model A summary, Table S3: R graphic output of model B summary, Table S4: Performance comparisons between different quantiles (Q97.5, Q85, Q75, Q50) for models A and B, depending on the observed phoretic Varroa numbers at t−x and the observed phoretic Varroa numbers at *t*. Error rates represent the colony percentage for which the Varroa load at the horizon of prediction was badly predicted with the cross-validation method. For each quantile, the number of hives to treat depends on the error rate for which their percentages were reported, Table S5: Performance comparisons between different quantiles (Q97.5, Q85, Q75, Q50) for models A and B, depending on the observed phoretic Varroa numbers at t-x and the observed phoretic Varroa numbers at t. Error rates represent the colony percentage for which the Varroa load at the horizon of prediction was badly predicted with the training validation method. For each quantile, the number of hives to treat depends on the error rate for which their percentages were reported.

**Author Contributions:** Conceptualization, F.M. and A.K.; methodology, F.M., A.K., L.M., H.D. and S.S.; investigation, P.M., A.M., C.V., B.B., F.M. and A.K.; data curation, H.D., L.M., S.S., Y.P., M.P. and A.K.; writing—original draft preparation, H.D., F.M. and A.K.; supervision, F.M. and A.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data are available upon request to the corresponding author.

**Acknowledgments:** We wish to thank all beekeepers that gave us access to their apiaries to record data.

**Conflicts of Interest:** The authors declare no conflict of interest.
