1. Introduction

2. Literature Review

3. Materials and Methods

3.4. Models

Both models are generalized linear models. We have defined linear predictor as $\eta =X\beta$ where X denotes vector of features described in Table 1 and B is vector of coefficients. Each coefficient’s distribution is described in appropriate Model section. Both vectors are size Nx1 [18].

We then used logarithmic link function to transform linear predictor’s domain to positive real numbers. It was one of the possible options, but nevertheless necessary, as both models are defined by inverse gamma function. This way we obtained explanatory variable ${\mu }_i$, representing mean of the outcome variable [18]. We defined prior distribution for standard deviation of our models, denoted as $\sigma$, to be exponential distribution with rate parameter equal 0.5.

Lastly, we defined likelihood as inverse gamma distribution with parameters shape ($\alpha$) and scale ($\beta$), computed from $\mu$ and $\sigma$ in such way, that the resulting distribution had mean and standard deviation of $\mu$ and $\sigma$. Reasoning behind this particular distribution was to model skewness of the data effectively. Also, time has to be strictly positive continuous, which inverse gamma also provides.

While models itself were defined in Stan, the experiments were done via CmdStanPy, which is one of Python’s interface for it [15].

3.4.1. Model 1

First model is defined as follows:

$$delivery\_tim{e}_i\sim InverseGamma({\alpha }_i,{\beta }_i)$$

(5)

$${\alpha }_i=\frac{{\mu }_i^2}{{\sigma }_i^2}+2$$

(6)

$${\beta }_i=\frac{{\mu }_i^3}{{\sigma }_i^2}+{\mu }_i$$

(7)

$${\sigma }_i\sim Exponential\left(0.5\right)$$

(8)

$$\begin{array}{cc}\hfill {\mu }_i=exp(& \mathit{distance}\_{\mathit{coeff}}_{\mathit{i}}·distanc{e}_i+\mathit{traffic}\_\mathit{level}\_\mathit{coeff}[\mathit{traffic}\_{\mathit{level}}_{\mathit{i}}]+\hfill \\ \hfill & +\mathit{meal}\_\mathit{prep}\_{\mathit{coeff}}_{\mathit{i}}·meal\_preparation\_tim{e}_i+mea{n}_i)\hfill \end{array}$$

(9)

$$mea{n}_i\sim N(3,0.1)$$

(10)

$$\mathit{distance}\_{\mathit{coeff}}_{\mathit{i}}\sim Normal(0,0.3)$$

(11)

$$\mathit{meal}\_\mathit{prep}\_{\mathit{coeff}}_{\mathit{i}}\sim Normal(0,0.3)$$

(12)

$$\mathit{traffic}\_\mathit{level}\_\mathit{coeff}\left[\mathit{1}\right]\sim Normal(0,0.3)$$

(13)

$$\mathit{traffic}\_\mathit{level}\_\mathit{coeff}\left[\mathit{2}\right]\sim Normal(0,0.3)$$

(14)

$$\mathit{traffic}\_\mathit{level}\_\mathit{coeff}\left[\mathit{3}\right]\sim Normal(0,0.3)$$

(15)

$$\mathit{traffic}\_\mathit{level}\_\mathit{coeff}\left[\mathit{4}\right]\sim Normal(0,0.3)$$

(16)

3.4.2. Model 2

Second model is extension of the first model by two predictors: number of deliveries and standardized delivery person rating. As such we only present changes necessary for creating model 2 out of model 1:

$$\begin{array}{cc}\hfill {\mu }_i=exp(& \mathit{distance}\_{\mathit{coeff}}_{\mathit{i}}·distanc{e}_i+\mathit{traffic}\_\mathit{level}\_\mathit{coeff}[\mathit{traffic}\_{\mathit{level}}_{\mathit{i}}]+\hfill \\ \hfill & +\mathit{meal}\_\mathit{prep}\_{\mathit{coeff}}_{\mathit{i}}·meal\_preparation\_tim{e}_i+\hfill \\ \hfill & +\mathit{deliveries}\_\mathit{number}\_\mathit{coeff}[\mathit{number}\_\mathit{of}\_{\mathit{deliveries}}_{\mathit{i}}]+\hfill \\ \hfill & +\mathit{person}\_\mathit{rating}\_{\mathit{coeff}}_{\mathit{i}}·delivery\_person\_ratin{g}_i+mea{n}_i)\hfill \end{array}$$

(17)

$$\mathit{person}\_\mathit{rating}\_{\mathit{coeff}}_{\mathit{i}}\sim Normal(0,0.3)$$

(18)

$$\mathit{deliveries}\_\mathit{number}\_\mathit{coeff}\left[\mathit{1}\right]\sim Normal(0,0.3)$$

(19)

$$\mathit{deliveries}\_\mathit{number}\_\mathit{coeff}\left[\mathit{2}\right]\sim Normal(0,0.3)$$

(20)

$$\mathit{deliveries}\_\mathit{number}\_\mathit{coeff}\left[\mathit{3}\right]\sim Normal(0,0.3)$$

(21)

$$\mathit{deliveries}\_\mathit{number}\_\mathit{coeff}\left[\mathit{4}\right]\sim Normal(0,0.3)$$

(22)

Figure 2. Sampling check for prior distributions of model 1 link function parameters (parameters with _coeff suffix). X-axis represents coefficients values and Y-axis represents sample count. Each of the coefficients follows its distribution, which is necessary for prior check to be successful. (a) Prior distribution of distance coefficient, defined as $Normal(0,0.3)$. (b) Prior distribution of meal preparation time coefficient, defined as $Normal(0,0.3)$. (c) Prior distribution of $mean$ parameter, defined as $Normal(3,0.1)$. $mean$ parameter represents our belief of what mean delivery time should be in case all other parameters are 0. (d) Joint plot of prior distributions of traffic level coefficients, all defined as $Normal(0,0.3)$.

Figure 2. Sampling check for prior distributions of model 1 link function parameters (parameters with _coeff suffix). X-axis represents coefficients values and Y-axis represents sample count. Each of the coefficients follows its distribution, which is necessary for prior check to be successful. (a) Prior distribution of distance coefficient, defined as $Normal(0,0.3)$. (b) Prior distribution of meal preparation time coefficient, defined as $Normal(0,0.3)$. (c) Prior distribution of $mean$ parameter, defined as $Normal(3,0.1)$. $mean$ parameter represents our belief of what mean delivery time should be in case all other parameters are 0. (d) Joint plot of prior distributions of traffic level coefficients, all defined as $Normal(0,0.3)$.

Figure 3. Computation and sample check of model 1 likelihood parameters. X-axis is time in minutes and Y-axis represents sample count. (a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 4.4 to 47 minutes. Mean of this distribution (at the top of the plot) is 23 minutes, which is reasonable value. (b) Prior distribution of standard deviation of the model, defined as $Exponential\left(0.5\right)$.

Figure 3. Computation and sample check of model 1 likelihood parameters. X-axis is time in minutes and Y-axis represents sample count. (a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 4.4 to 47 minutes. Mean of this distribution (at the top of the plot) is 23 minutes, which is reasonable value. (b) Prior distribution of standard deviation of the model, defined as $Exponential\left(0.5\right)$.

Figure 4. Prior predictive checks - model 1. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 3.2 to 47 minutes, which is broad range. Mean of this distribution (at the top of the plot) is 23 minutes, which is reasonable value. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data is included within generated data, which means that all observations are possible within prior model. This means that prior checks are successful.

Figure 4. Prior predictive checks - model 1. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 3.2 to 47 minutes, which is broad range. Mean of this distribution (at the top of the plot) is 23 minutes, which is reasonable value. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data is included within generated data, which means that all observations are possible within prior model. This means that prior checks are successful.

Figure 5. Sampling check for prior distributions of model 2 link function parameters (parameters with _coeff suffix). X-axis represents coefficients values and Y-axis represents sample count. Each of the coefficients follows its distribution, which is necessary for prior check to be successful. (a) Prior distribution of distance coefficient, defined as $Normal(0,0.3)$. (b) Prior distribution of meal preparation time coefficient, defined as $Normal(0,0.3)$. (c) Prior distribution of $mean$ parameter, defined as $Normal(3,0.1)$. $mean$ parameter represents our belief of what mean delivery time should be in case all other parameters are 0. (d) Joint plot of prior distributions of traffic level coefficients, all defined as $Normal(0,0.3)$. (e) Prior distribution of delivery person rating coefficient, defined as $Normal(0,0.3)$. (f) Joint plot of prior distributions of deliveries number coefficients, all defined as $Normal(0,0.3)$.

Figure 5. Sampling check for prior distributions of model 2 link function parameters (parameters with _coeff suffix). X-axis represents coefficients values and Y-axis represents sample count. Each of the coefficients follows its distribution, which is necessary for prior check to be successful. (a) Prior distribution of distance coefficient, defined as $Normal(0,0.3)$. (b) Prior distribution of meal preparation time coefficient, defined as $Normal(0,0.3)$. (c) Prior distribution of $mean$ parameter, defined as $Normal(3,0.1)$. $mean$ parameter represents our belief of what mean delivery time should be in case all other parameters are 0. (d) Joint plot of prior distributions of traffic level coefficients, all defined as $Normal(0,0.3)$. (e) Prior distribution of delivery person rating coefficient, defined as $Normal(0,0.3)$. (f) Joint plot of prior distributions of deliveries number coefficients, all defined as $Normal(0,0.3)$.

Figure 6. Computation and sample check of model 2 likelihood parameters. X-axis is time in minutes and Y-axis is sample count. (a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 2.3 to 57 minutes. Mean of this distribution (at the top of the plot) is 26 minutes, more than for model 1, but still within reasonable range. (b) Prior distribution of standard deviation of the model, defined as $Exponential\left(0.5\right)$.

Figure 6. Computation and sample check of model 2 likelihood parameters. X-axis is time in minutes and Y-axis is sample count. (a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 2.3 to 57 minutes. Mean of this distribution (at the top of the plot) is 26 minutes, more than for model 1, but still within reasonable range. (b) Prior distribution of standard deviation of the model, defined as $Exponential\left(0.5\right)$.

Figure 7. Prior predictive checks - model 2. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 1.3 to 58 minutes. It is very broad, improbable range, but for prior checks it is sufficient. Mean of this distribution (at the top of the plot) is 26 minutes, which is reasonable value. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data is included within generated data, which means that all observations are possible within prior model. This means that prior checks are successful.

Figure 7. Prior predictive checks - model 2. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 1.3 to 58 minutes. It is very broad, improbable range, but for prior checks it is sufficient. Mean of this distribution (at the top of the plot) is 26 minutes, which is reasonable value. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data is included within generated data, which means that all observations are possible within prior model. This means that prior checks are successful.

3.4.3. Priors and Prior Predictive Checks

We have decided to use unbounded weakly informative priors for all parameters for two main reasons. First, we lack expert knowledge on the influence of each feature. Second, the abundance of data reduces the influence of priors on the final distribution as more data points are added.

We chose a normal distribution with a mean of 0 and a standard deviation of 0.3 for our parameters. This distribution provides a value range of approximately -1 to 1, with most values likely clustering around 0. This choice reflects our initial assumption that each parameter is not highly influential, while still allowing for a small probability that they could be significantly influential. The standard deviation of 0.3 was selected to accommodate the exponential coming from inverse link function, as larger values from the linear combination could numerically destabilize the model.

The exception to this is the intercept parameter, $mea{n}_i$, for which we chose a strong prior: $N(3,0.1)$. It results in the base mean delivery time around 15-30 minutes. We opted for a strong prior here to ensure that our model accurately reflects the most average delivery time. Given that our features are standardized with a mean of 0, the linear combination for the most average case would be close to zero, leading to an unrealistic delivery time in the posterior distribution. The $mea{n}_i$ prior helps anchor the model, providing a trusted average time for an average case.

For both models prior predictive checks gave good results i.e. observed data was included within simulated data range and no outright impossible values was generated from either of the models.

4. Results

4.1. Posterior Predictive Checks for Model 1

Model 1 gave decent results. All of the observed data falls withing samples from posterior distribution, visual overlap is quite high. Posterior distribution exhibits long tail, which is drawback of using inverse gamma function. Data represented on Figure 10

Model coefficients for distance and meal preparation time ended with very narrow distributions, albeit positive ones, indicating that they impact output variable. Mean intercept parameter ended with mean closer to 3.1, which is also closer to the mean of the dataset ( $e^{3.1}\approx 22.18$ while mean of dataset is $\approx 27.05$ ). Traffic level coefficient represent trend, in which low traffic contributes to faster delivery times, and as traffic level increases the delivery times become longer. This interpretation is viable as it is not multiplicand, but sum component, so negative values will result in smaller mean and positive in larger mean. There is almost no distinction between influence of high traffic and jam. Data represented on Figure 8.

Linear model of mean delivery times ${\mu }_i$ results in probable values with respect to dataset. Standard deviation has completely changed it’s distribution and now follows normal distribution centered around 9.5. It is quite close to std of dataset, which is $\approx 8.99$. Data represented on Figure 9.

Figure 8. Sampling check for posterior distributions of model 1 link function parameters. X-axis represents coefficients values and Y-axis represents sample count. (a) Posterior distribution of distance coefficient. It is much narrower than prior distribution, but still follows normal distribution. Positive value indicates that it has impact on the output variable. (b) Posterior distribution of meal preparation time coefficient. Conclusions the same as for the distance coefficient. (c) Posterior distribution of $mean$ parameter. It has mean closer to 3.1 which more likely represents mean of the dataset. (d) Joint plot of posterior distributions of traffic level coefficients. The bigger the traffic level, the more impact it has on the outcome variable. Jam and High levels have the same impact.

Figure 8. Sampling check for posterior distributions of model 1 link function parameters. X-axis represents coefficients values and Y-axis represents sample count. (a) Posterior distribution of distance coefficient. It is much narrower than prior distribution, but still follows normal distribution. Positive value indicates that it has impact on the output variable. (b) Posterior distribution of meal preparation time coefficient. Conclusions the same as for the distance coefficient. (c) Posterior distribution of $mean$ parameter. It has mean closer to 3.1 which more likely represents mean of the dataset. (d) Joint plot of posterior distributions of traffic level coefficients. The bigger the traffic level, the more impact it has on the outcome variable. Jam and High levels have the same impact.

Figure 9. Computation and sample check of model 1 likelihood parameters. X-axis represents coefficients values and Y-axis represents sample count.(a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 21 to 32 minutes. Those are much more realistic values than the ones from prior distribution. Mean of this distribution (at the top of the plot) is 27 minutes, a reasonable value. (b) Posterior distribution of standard deviation of the model, it no longer resembles prior, now it follows normal distribution with mean ≈ 9.5.

Figure 9. Computation and sample check of model 1 likelihood parameters. X-axis represents coefficients values and Y-axis represents sample count.(a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 21 to 32 minutes. Those are much more realistic values than the ones from prior distribution. Mean of this distribution (at the top of the plot) is 27 minutes, a reasonable value. (b) Posterior distribution of standard deviation of the model, it no longer resembles prior, now it follows normal distribution with mean ≈ 9.5.

Figure 10. Posterior predictive checks - model 1. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 11 to 46 minutes. It is broad range, but realistic nevertheless. Mean of this distribution (at the top of the plot) is 27 minutes, which is reasonable value. It follows inverse gamma distribution as defined. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data has high overlap with sampled data from posterior distribution.

Figure 10. Posterior predictive checks - model 1. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 11 to 46 minutes. It is broad range, but realistic nevertheless. Mean of this distribution (at the top of the plot) is 27 minutes, which is reasonable value. It follows inverse gamma distribution as defined. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data has high overlap with sampled data from posterior distribution.

4.2. Posterior Predictive Checks for Model 2

Model 2 gave visually better results. All of the observed data falls withing samples from posterior distribution, as with model 1. Tail is shorter than in Model 1. Data represented on Figure 12. Since there is not much difference between posterior distributions for shared features of both models, we will only comment on new features, distinct to Model 2, as well as on likelihood related parameters.

Delivery person rating follows narrow normal distribution centered around -0.085. Since it is negative and has small std we can reason that delivery person rating is inversely related to delivery time. This is expected as couriers with higher scores are more likely to deliver food faster. Number of deliveries follows exactly the same trend as traffic level coefficients, but numerically is more important as the values range is greater. Data represented on Figure 11.

Linear model of mean delivery times ${\mu }_i$ has much longer tail with regards to Model 1, which results in different 94% HDI. Standard deviation has similar distribution to Model 1, although it’s mean value is smaller, around 7.85. Data represented on Figure 12.

Figure 11. Sampling check for posterior distributions of model 2 link function parameters. X-axis represents coefficients values and Y-axis represents sample count. (a) Posterior distribution of distance coefficient. It is much narrower than prior distribution, but still follows normal distribution. Positive value indicates that it has impact on the output variable. (b) Posterior distribution of meal preparation time coefficient. Conclusions the same as for the distance coefficient. (c) Posterior distribution of $mean$ parameter. It has mean closer to 3.1 which more likely represents mean of the dataset. (d) Joint plot of posterior distributions of traffic level coefficients. The bigger the traffic level, the more impact it has on the outcome variable. Jam and High levels have the same impact.(e) Posterior distribution of delivery person rating coefficient. It is much narrower than prior distribution, but still follows normal distribution. Negative values indicate inverse relationship between delivery time and rating; the bigger the courier rating the faster delivery will be made. (f) Joint plot of posterior distributions of deliveries number coefficients. The more deliveries, the more impact it has on the outcome variable. This is the same trend as for traffic level, but greater range translates to greater impact.

Figure 11. Sampling check for posterior distributions of model 2 link function parameters. X-axis represents coefficients values and Y-axis represents sample count. (a) Posterior distribution of distance coefficient. It is much narrower than prior distribution, but still follows normal distribution. Positive value indicates that it has impact on the output variable. (b) Posterior distribution of meal preparation time coefficient. Conclusions the same as for the distance coefficient. (c) Posterior distribution of $mean$ parameter. It has mean closer to 3.1 which more likely represents mean of the dataset. (d) Joint plot of posterior distributions of traffic level coefficients. The bigger the traffic level, the more impact it has on the outcome variable. Jam and High levels have the same impact.(e) Posterior distribution of delivery person rating coefficient. It is much narrower than prior distribution, but still follows normal distribution. Negative values indicate inverse relationship between delivery time and rating; the bigger the courier rating the faster delivery will be made. (f) Joint plot of posterior distributions of deliveries number coefficients. The more deliveries, the more impact it has on the outcome variable. This is the same trend as for traffic level, but greater range translates to greater impact.

Figure 12. Computation and sample check of model 2 likelihood parameters. X-axis represents coefficients values and Y-axis represents sample count.(a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 19 to 37 minutes. Those are much more realistic values than the ones from prior distribution. Mean of this distribution (at the top of the plot) is 27 minutes, a reasonable value. (b) Posterior distribution of standard deviation of the model, it no longer resembles prior, now it follows normal distribution with mean ≈ 7.85.

Figure 12. Computation and sample check of model 2 likelihood parameters. X-axis represents coefficients values and Y-axis represents sample count.(a) Computed $\mu$ represent mean delivery time for each sample. HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 19 to 37 minutes. Those are much more realistic values than the ones from prior distribution. Mean of this distribution (at the top of the plot) is 27 minutes, a reasonable value. (b) Posterior distribution of standard deviation of the model, it no longer resembles prior, now it follows normal distribution with mean ≈ 7.85.

Figure 13. Posterior predictive checks - model 2. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 12 to 45 minutes. It is slightly narrower than model 1 range, but realistic nevertheless. Mean of this distribution (at the top of the plot) is 27 minutes, which is reasonable value. It follows inverse gamma distribution as defined. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data has high overlap with sampled data from posterior distribution. Generated data has shorter tail than model 1, which is desirable.

Figure 13. Posterior predictive checks - model 2. (a) HDI 94% is represented as black bar at the bottom of the plot and tells us that 94% of shown mean times falls in range 12 to 45 minutes. It is slightly narrower than model 1 range, but realistic nevertheless. Mean of this distribution (at the top of the plot) is 27 minutes, which is reasonable value. It follows inverse gamma distribution as defined. (b) Real and simulated data overlay. Both are normalized, so that integral of the graph is 1. It was necessary for comparison. Measured data has high overlap with sampled data from posterior distribution. Generated data has shorter tail than model 1, which is desirable.

4.3. Models Comparison

Models were compared using WAIC and PSIS-LOO criteria using ArviZ library for exploratory analysis of Bayesian Models [19].

The WAIC (Widely Applicable or Watanabe-Akaike Information Criterion) is a statistical measure used to estimate the out-of-sample predictive accuracy of a model. It does this by evaluating the within-sample predictive accuracy and making necessary adjustments. WAIC calculates the log pointwise posterior predictive density (LPPD) and includes a correction for the effective number of parameters to account for overfitting. This correction is done by subtracting the sum of the posterior variances of the log predictive densities for each data point [20].

Table 2. Comparison result with WAIC criterion.

Table 2. Comparison result with WAIC criterion.

Model	Rank	waic	p_waic	d_waic	weight	SE	dSE
2	0	-117249.732164	14.307179	0.000000	0.997408	145.013769	0.000000
1	1	-122442.256308	9.314501	5192.524145	0.002592	137.412164	88.510163

Parameters reference [21].

Model 2 has higher ELPD score (denoted as waic) which indicates it’s better within-sample fit. WAIC also correctly states that it has higher number of effective parameters (p_waic). Weight parameter clearly states that Model 2 has near 1 probability within given data. It is slightly more uncertain than Model 1, which is indicated by SE (standard error) parameter, but when compared to differences in WAIC score and size of the dataset, it is not overly large. Overall WAIC clearly evaluated Model 2 as superior.

Figure 14. Comparison plot for WAIC criterion. Black dots indicate ELPD of each model with their standard error (black lines). Grey triangle represents standard error of difference in ELPD between Model 1 and top-ranked Model 2. Plot indicates that Model 2 performs better.

Figure 14. Comparison plot for WAIC criterion. Black dots indicate ELPD of each model with their standard error (black lines). Grey triangle represents standard error of difference in ELPD between Model 1 and top-ranked Model 2. Plot indicates that Model 2 performs better.

Figure 15. Comparison plot for PSIS-LOO criterion. Black dots indicate ELPD of each model with their standard error (black lines). Grey triangle represents standard error of difference in ELPD between Model 1 and top-ranked Model 2. Plot indicates that Model 2 performs better.

Figure 15. Comparison plot for PSIS-LOO criterion. Black dots indicate ELPD of each model with their standard error (black lines). Grey triangle represents standard error of difference in ELPD between Model 1 and top-ranked Model 2. Plot indicates that Model 2 performs better.

The PSIS-LOO (Pareto Smoothed Importance Sampling using Leave-One-Out validation) method is used to calculate out-of-sample predictive fit by summing the log leave-one-out predictive densities. These densities are evaluated using importance ratios (IS-LOO). However, the importance ratios can exhibit high or infinite variance, leading to instability in the estimates. To mitigate this issue, a generalized Pareto distribution is fitted to the largest 20% of the importance ratios [20].

Table 3. Comparison result with PSIS-LOO criterion.

Table 3. Comparison result with PSIS-LOO criterion.

Model	Rank	loo	p_loo	d_loo	weight	SE	dSE
2	0	-117249.732784	14.307799	0.000000	0.997408	145.013813	0.000000
1	1	-122442.256329	9.314522	5192.523545	0.002592	137.412164	88.510227

Parameters reference [21].

The PSIS-LOO evaluation provides near identical results as WAIC and the same as above conclusions can be drawn.

5. Conclusions

Abbreviations

References

1. Statista. Online Food Delivery - Worldwide. https://www.statista.com/outlook/emo/online-food-delivery/worldwide, 2024. accessed on 04.05.2024.
2. IMARC Group. India Online Food Delivery Market Report. https://www.imarcgroup.com/india-online-food-delivery-market, 2023. accessed on 04.05.2024.
3. Alalwan, A.A. Mobile food ordering apps: An empirical study of the factors affecting customer e-satisfaction and continued intention to reuse. International Journal of Information Management 2020, 36, 28–44. [Google Scholar] [CrossRef]
4. Abdi, A.; Amrit, C. A review of travel and arrival-time prediction methods on road networks: classification, challenges and opportunities. PeerJ Computer Science 2021, 37, e689. [Google Scholar] [CrossRef] [PubMed]
5. Lee, W.H.; Tseng, S.S.; Tsai, S.H. A knowledge based real-time travel time prediction system for urban network. Expert systems with Applications 2009, 36, 4239–4247. [Google Scholar] [CrossRef]
6. Li, X.; Cong, G.; Sun, A.; Cheng, Y. Learning travel time distributions with deep generative model. The World Wide Web Conference, 2019, pp. 1017–1027.
7. Asghari, M.; Emrich, T.; Demiryurek, U.; Shahabi, C. Probabilistic estimation of link travel times in dynamic road networks. Proceedings of the 23rd SIGSPATIAL International Conference on Advances in Geographic Information Systems, 2015, pp. 1–10.
8. Wang, Y.; Zheng, Y.; Xue, Y. Travel time estimation of a path using sparse trajectories. Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, 2014, pp. 25–34.
9. Wang, Z.; Fu, K.; Ye, J. Learning to estimate the travel time. Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, 2018, pp. 858–866.
10. Zhu, L.; Yu, W.; Zhou, K.; Wang, X.; Feng, W.; Wang, P.; Chen, N.; Lee, P. Order fulfillment cycle time estimation for on-demand food delivery. Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2020, pp. 2571–2580.
11. Li, Y.; Fu, K.; Wang, Z.; Shahabi, C.; Ye, J.; Liu, Y. Multi-task representation learning for travel time estimation. Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, 2018, pp. 1695–1704.
12. Wu, C.H.; Ho, J.M.; Lee, D.T. Travel-time prediction with support vector regression. IEEE transactions on intelligent transportation systems 2004, 5, 276–281. [Google Scholar] [CrossRef]
13. Wang, H.; Tang, X.; Kuo, Y.H.; Kifer, D.; Li, Z. A simple baseline for travel time estimation using large-scale trip data. ACM Transactions on Intelligent Systems and Technology (TIST) 2019, 10, 1–22. [Google Scholar] [CrossRef]
14. Gelman, A.; Carlin, J.B.; Stern, H.S.; Dunson, D.B.; Vehtari, A.; Rubin, D.B. Bayesian Data Analysis, 3rd ed.; Chapman and Hall/CRC: New York, USA, 2013; pp. 1–27. [Google Scholar]
15. Team, S.D. Stan Modeling Language Users Guide and Reference Manual, 2.34. https://mc-stan.org/docs/stan-users-guide/index.html. accessed on 13.05.2024.
16. Food Delivery Dataset. https://www.kaggle.com/datasets/gauravmalik26/food-delivery-dataset, 2023. accessed on 13.05.2024.
17. Gibas, J.; Pomykacz, J.; Baranowski, J. Bayesian modelling of travel times on the example of food delivery: Part 1 - Spatial data analysis and processing. Electronics, 2024. submitted.
18. Gelman, A.; Carlin, J.B.; Stern, H.S.; Dunson, D.B.; Vehtari, A.; Rubin, D.B. Bayesian Data Analysis, 3rd ed.; Chapman and Hall/CRC: New York, USA, 2013; pp. 405–411. [Google Scholar]
19. ArviZ. API Reference. https://python.arviz.org/en/stable/api/index.html, 2024. accessed on 20.05.2024.
20. Vehtari, A.; Gelman, A.; Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and computing 2017, 27, 1413–1432. [Google Scholar] [CrossRef]
21. ArviZ. arviz.compare. https://python.arviz.org/en/stable/api/generated/arviz.compare.html, 2024. accessed on 20.05.2024.