A Bayesian analysis of the time through the order penalty in baseball

Ryan S. Brill; Sameer K. Deshpande; Abraham J. Wyner

doi:10.1515/jqas-2022-0116

Published online by De Gruyter June 27, 2023

A Bayesian analysis of the time through the order penalty in baseball

Ryan S. Brill , Sameer K. Deshpande and Abraham J. Wyner

From the journal Journal of Quantitative Analysis in Sports

https://doi.org/10.1515/jqas-2022-0116

Showing a limited preview of this publication:

Abstract

As a baseball game progresses, batters appear to perform better the more times they face a particular pitcher. The apparent drop-off in pitcher performance from one time through the order to the next, known as the Time Through the Order Penalty (TTOP), is often attributed to within-game batter learning. Although the TTOP has largely been accepted within baseball and influences many managers’ in-game decision making, we argue that existing approaches of estimating the size of the TTOP cannot disentangle continuous evolution in pitcher performance over the course of the game from discontinuities between successive times through the order. Using a Bayesian multinomial regression model, we find that, after adjusting for confounders like batter and pitcher quality, handedness, and home field advantage, there is little evidence of strong discontinuity in pitcher performance between times through the order. Our analysis suggests that the start of the third time through the order should not be viewed as a special cutoff point in deciding whether to pull a starting pitcher.

Keywords: baseball; Bayesian statistics; mathematical modeling; pitching; time through the order penalty

Corresponding author: Ryan S. Brill, Graduate Group in Applied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, PA, USA, E-mail: ryguy123@sas.upenn.edu

Funding source: Wisconsin Alumni Research Foundation

Acknowledgments

The authors thank Tom Tango for his comments on an early draft of this paper. The authors acknowledge the High Performance Computing Center (HPCC) at The Wharton School, University of Pennsylvania for providing computational resources that have contributed to the research results reported within this paper.

Author contributions: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.
Research funding: Support for S.K.D. was provided by the University of Wisconsin–Madison, Office of the Vice Chancellor for Research and Graduate Education with funding from the Wisconsin Alumni Research Foundation.
Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

Appendix A: Our code and data

Our code is available on Github.^[3] The data_wrangling folder of the Github repository contains our dataset processing, including the Retrosheet data scraper. The data folder further processes the full dataset into a smaller dataset relevant for this paper. Finally, the model_positive_slope_prior folder contains our data analysis, including our Stan model.

The final datasets used in this paper are available for download.^[4] The cleaned dataset of all MLB plate appearances from 1990 to 2020 is retro_final_PA_1990-2020d.csv. The datasets T T O _ d a t a s e t _ 4 1 0 . c s v and T T O _ d a t a s e t _ 5 1 0 . c s v are processed subsets of the large dataset which we use to fit our models.

Appendix B: Model simulation study

We conduct a simulation study to assess the capacity of our model (Equation (2)) to estimate time through the order penalties of various sizes. Specifically, we simulate data consistent with different TTOPs and verify that our posterior estimates are close to the data generating parameters.

B.1 Simulation setup

For our first simulation, we generate data consistent with continuous pitcher fatigue and no TTOP for any of the plate appearance outcomes by setting β _2k = β _3k = 0 for each k ≠ 1. In our second simulation, for each k ≠ 1, we set the β _2k and β _3k so that the resulting xwOBA curves display TTOPs consistent with Tango, Lichtman, and Dolphin (2007)’s findings of about 10 expected wOBA points between successive times through the order. Finally, for our third simulation, we set β _2k and β _3k so that there is no 2TTOP (in terms of xwOBA) but a large 3TTOP of about 50 wOBA points. For each simulation, we set the values of the α _0k’s, α _1k’s, and η _k’s in a way that is consistent with observed data. Additional details about the simulation setup, including the data generating parameter values, are available in Appendix C.

For each simulation, we generate 225 full seasons worth of data. We fit our model to 80 % of the data from each simulated season and evaluate our fitted model’s predictive performance on the remaining 20 %. We further assess how well our fitted model recovers the function xwOBA(t, x ) for a set of average confounder values.

B.2 Simulation results

In all three simulation studies, we reliably recover the data generating parameters: averaged across all parameters, the estimated frequentist coverage of the marginal 95 % posterior credible intervals exceeds 92 % in each study. Importantly, the coverage of the 95 % posterior credible intervals for the discontinuity parameters β _2k and β _3k exceeds 91 % in each study. That is, for each simulated dataset, the 95 % credible intervals for the β _2k’s and β _3k’s usually contain the true data generating parameters. Furthermore, our model demonstrates good predictive capabilities (see Appendix C for details).

B.3 Simulation visualization

In each simulation, we visualize the trajectory of posterior expected wOBA over the course of the game for an average batter on the road facing an average pitcher with the same handedness. That is, we plot the sequence { xwOBA ( t , x ̃ ) } t = 1 27 where

(18) x ̃ ⊤ = ( x ( b ) ̄ , x ( p ) ̄ , 1,0 ) .

Figure 6 shows the sequence of posterior means, 50 %, and 95 % credible intervals of xwOBA ( t , x ̃ ) based on a single simulated dataset from each simulation setting. We overlay the true values of xwOBA ( t , x ̃ ) , computed from the data generating parameters, to each plot. We see that in each of the three simulation studies, we recover the true underlying expected wOBA trajectory.

Figure 6:

Trend in xwOBA over the course of a game from our first, second, and third simulation studies. The red dots indicate the true underlying expected wOBA values, the white dots indicate the posterior means of the xwOBA values, the thick black error bars denote the 50 % posterior credible intervals, and the thin black error bars denote the 95 % posterior credible intervals.

Appendix C: Simulation details

C.1 Data generating parameters

The exact data generating parameter values of β _2k and β _3k for our three simulation studies are shown in Table 4.

Table 4:

The data generating parameter values of β _2k and β _3k in each of our three simulations.

	k = BB	k = HBP	k = 1B	k = 2B	k = 3B	k = HR
β _2k for sim 1	0	0	0	0	0	0
β _3k for sim 1	0	0	0	0	0	0
β _2k for sim 2	2/65	0	4/65	2/65	0	2/65
β _3k for sim 2	1/15	0	2/15	1/15	0	1/15
β _2k for sim 3	0	0	0	0	0	0
β _3k for sim 3	1/10	1/10	3/10	1/10	1/10	3/20

Furthermore, in each of our simulation studies, we assume that pitchers fatigue linearly over the course of a game. The particular true parameter values of α _0k and α _1k used in each of our simulation studies are shown in Table 5.

Table 5:

The data generating parameter values of α _0k and α _1k in each of our three simulations.

	k = BB	k = HBP	k = 1B	k = 2B	k = 3B	k = HR
α _0k	−0.601	−1.804	−0.475	−0.943	−1.510	−0.565
α _1k	0.00271	0.0122	0.00354	0.00635	0.0223	0.00926

Finally, in each of our simulation studies, we set the value of η to mimic fitted values from observed data. The particular true parameter values of η used in each of our simulation studies are shown in Table 6.

Table 6:

The data generating parameter values of η in each of our three simulations.

	k = BB	k = HBP	k = 1B	k = 2B	k = 3B	k = HR
η _{bat_quality}	0.865	1.408	0.371	0.856	1.399	1.525
η _{pit_quality}	1.128	1.987	1.050	1.472	3.286	1.850
η _hand	−0.201	0.166	−0.0164	−0.0420	−0.462	−0.0958
η _home	0.0792	−0.0776	0.0245	−0.00103	0.107	0.0230

C.2 Predictive performance on simulated data

Our model demonstrates good predictive capabilities. To get a general sense of our model’s performance, we use out-of-sample cross entropy loss, given by

(19) − 1 n ∑ i = 1 n ∑ k = 1 7 1 { y i = k } ⋅ log P ( y i = k ) .

For each of our three simulations, the average cross entropy loss over each of our 25 datasets is 1.05, 1.06, and 1.07, respectively. Using the empirical outcome probabilities yields an average out-of-sample cross-entropy loss of 1.06, 1.08, and 1.08, respectively, for each of our three simulations. It is reassuring that our model (barely) outperforms the observed base rates.

Appendix D: Observed model fit details

D.1 The impact of pitcher decline on the outcome of a plate appearance

In this Section, we quantify the effect size of pitcher decline over the course of a game, again using the 2017 season as our primary example.

In particular, we examine how the probability of each outcome of a plate appearance changes over the course of a game. Specifically, we use the posterior distribution of P ( y = k | t , x ) , defined in Formula (6), to characterize the amount by which pitchers decline within a game. In particular, we compute the posterior distribution of the change in the probability of outcome k ≠ 1 from 1TTO to 2TTO, over average,

(20) D 12 ( k , x ) = 1 9 ∑ t = 10 18 P ( y = k | t , x ) − 1 9 ∑ t = 1 9 P ( y = k | t , x ) ,

and the similarly defined D 23 ( k , x ) , which captures the change in the probability of outcome k ≠ 1 from 2TTO to 3TTO, over average.

In Figure 7 we plot the posterior distribution of D 12 ( k , x ̃ ) , using plate-appearance-state vector x ̃ from Formula (12). From 1TTO to 2TTO, the probability of a single increases by about 0.005, the probability of a home run increases by about 0.003, and the probability of the other non-out categories change negligibly. With this, the probability of an out decreases by about 0.01. So, there is a small decrease in pitcher performance on average from 1TTO to 2TTO.

$Figure 7: The difference in probability of each plate appearance outcome between 2TTO and 1TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 12 ( k , x ̃ ) ${\mathcal{D}}_{12}(k,\tilde{\boldsymbol{x}})$ . The red line denotes the mean, and the blue line denotes 0.$

Figure 7:

The difference in probability of each plate appearance outcome between 2TTO and 1TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 12 ( k , x ̃ ) . The red line denotes the mean, and the blue line denotes 0.

Similarly, in Figure 8 we plot the posterior distribution of D 23 ( k , x ̃ ) . From 2TTO to 3TTO, the probability of a single increases by about 0.005, the probability of a double increases by about 0.004, and the probability of the other non-out categories change negligibly. With this, the probability of an out decreases by about 0.01. So, there is a small decrease in pitcher performance on average from 2TTO to 3TTO.

$Figure 8: The difference in probability of each plate appearance outcome between 3TTO and 2TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 23 ( k , x ̃ ) ${\mathcal{D}}_{23}(k,\tilde{\boldsymbol{x}})$ . The red line denotes the mean, and the blue line denotes 0.$

Figure 8:

The difference in probability of each plate appearance outcome between 3TTO and 2TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 23 ( k , x ̃ ) . The red line denotes the mean, and the blue line denotes 0.

Additionally, we examine how the expected wOBA of each outcome of a plate appearance changes over the course of a game. In particular, we compute the posterior distribution of the change in the expected wOBA of outcome k ≠ 1 from 1TTO to 2TTO, over average,

(21) D 12 ′ ( k , x ) = 1 9 ∑ t = 10 18 1000 ⋅ w k ⋅ P ( y = k | t , x ) − 1 9 ∑ t = 1 9 1000 ⋅ w k ⋅ P ( y = k | t , x ) ,

where w _k is the wOBA weight for outcome k as discussed in Section 2.4. Similarly, we define D 23 ′ ( k , x ) , which captures the change in the expected wOBA of outcome k ≠ 1 from 2TTO to 3TTO, over average.

In Figure 9 we plot the posterior distribution of D 12 ′ ( k , x ̃ ) , using plate-appearance-state vector x ̃ from Formula (12). From 1TTO to 2TTO, the expected wOBA points of a home run increases by about six, the expected wOBA points of a single increases by about four, and the other non-out categories change negligibly. Note that the expected wOBA of an out doesn’t change because an out is worth zero wOBA.

$Figure 9: The difference in xwOBA of each plate appearance outcome between 2TTO and 1TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 23 ′ ( k , x ̃ ) ${\mathcal{D}}_{23}^{\prime }(k,\tilde{\boldsymbol{x}})$ . The red line denotes the mean, and the blue line denotes 0.$

Figure 9:

The difference in xwOBA of each plate appearance outcome between 2TTO and 1TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 23 ′ ( k , x ̃ ) . The red line denotes the mean, and the blue line denotes 0.

Similarly, in Figure 10 we plot the posterior distribution of D 23 ′ ( k , x ̃ ) . From 2TTO to 3TTO, the expected wOBA of a double and single increases by about five, the xwOBA of a home run increases by about three, and the other categories change negligibly.

$Figure 10: The difference in xwOBA of each plate appearance outcome between 3TTO and 2TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 23 ′ ( k , x ̃ ) ${\mathcal{D}}_{23}^{\prime }(k,\tilde{\boldsymbol{x}})$ . The red line denotes the mean, and the blue line denotes 0.$

Figure 10:

The difference in xwOBA of each plate appearance outcome between 3TTO and 2TTO on average (assuming a batter of average quality on the road faces a pitcher of average quality with a handedness match during each plate appearance). Equivalently, the posterior distribution of D 23 ′ ( k , x ̃ ) . The red line denotes the mean, and the blue line denotes 0.

Furthermore, we aggregate the increase in the probability of each non-out plate appearance outcome k from one TTO to the next via expected wOBA, defined in Equation (8). In particular, recall from Section 3.2 that a pitcher declines by about 13 wOBA points from one TTO to the next, over average, which is consistent with the effect sizes from Figures 7 and 8. Figure 11 illustrates this via a histogram of the posterior samples of D 12 ( x ̃ ) and D 23 ( x ̃ ) . We see that virtually all of these samples are positive, suggesting that average pitcher performance declines from one TTO to the next, and that the means of these distributions are around 13 wOBA points, which are consistent with Tango, Lichtman, and Dolphin (2007)’s findings. Specifically, our model suggests that the expected wOBA points of an average plate appearance increases by 13.4 (with a 95 % credible interval of [7.78, 19.0]) from the first TTO to the second, and by 12.5 (with a 95 % credible interval of [5.98, 18.7]) from the second TTO to the third.

$Figure 11: The posterior distribution of the mean batter improvement, or mean pitcher decline in xwOBA, from 1TTO to 2TTO (left) and from 2TTO to 3TTO (right). Equivalently, the posterior distributions of D 12 ( x ̃ ) ${\mathcal{D}}_{12}(\tilde{\boldsymbol{x}})$ (left) and D 23 ( x ̃ ) ${\mathcal{D}}_{23}(\tilde{\boldsymbol{x}})$ (right) (see Formula (13)). The red line denotes the mean, and the blue line denotes 0. We see that batters improve relative to the pitcher by about 13 wOBA points on average from one TTO to the next.$

Figure 11:

The posterior distribution of the mean batter improvement, or mean pitcher decline in xwOBA, from 1TTO to 2TTO (left) and from 2TTO to 3TTO (right). Equivalently, the posterior distributions of D 12 ( x ̃ ) (left) and D 23 ( x ̃ ) (right) (see Formula (13)). The red line denotes the mean, and the blue line denotes 0. We see that batters improve relative to the pitcher by about 13 wOBA points on average from one TTO to the next.

D.2 Predictive performance on observed data

To get a general sense of our model’s performance on observed data, we run a five-fold cross validation to predict the probability of each plate appearance outcome for each plate appearance in 2017. The out-of-sample cross entropy loss, given by Formula (19), is 1.035. We compare our model’s cross entropy loss to that of other prediction strategies to better understand its performance. Consider a five-fold cross validation using the base rates of each plate appearance outcome. So, for each fold, find the proportion of plate appearances in which each outcome occurs, and compute the cross entropy loss using these base rates on the remaining out-of-sample plate appearances. For reference, in 2017, an out occurs in 67.6 % of plate appearances, an uBB 7.8 %, an HBP 0.9 %, a 1B 14.9 %, a 2B 4.8 %, a 3B 0.45 %, and an HR in 3.5 % of plate appearances. The out-of-sample cross entropy loss of the base rates of each outcome is 1.042. So, our model very slightly outperforms the base rates. Finally, note that our model using raw batter and pitcher quality covariates, rather than logit-transformed batter and pitcher quality covariates, has a cross-validated out-of-sample cross entropy loss of 1.040. That the logit-transformed player quality covariates have better out-of-sample predictive performance helps justify using the logit transform.

D.3 The trend is persistent across years

In Figure 12 we show boxplots of the posterior distributions of the discontinuity parameters β _2k and β _3k − β _2k from our model (Equation (2)) fit separately on data from each season from 2012 to 2019. For some outcomes (e.g., walks), the posterior distributions are tightly concentrated around 0, and for other outcomes (e.g., triples and hit-by-pitches, which are rare events), the posterior distributions are quite wide, which is compatible with a large effect in either direction. Overall, the posterior distributions of the discontinuity parameters cover both positive and negative values, and most of them are centered around 0. In particular, we don’t see what we would expect to see if there were strong evidence for a TTOP (i.e., we don’t see the posterior distributions tighly concentrated around a positive number). Ultimately, we do not find the posterior distributions in Figure 12 to be consistent with large, systematic time through the order penalties.

Figure 12:

Posterior boxplots of the TTOP discontinuity parameters from Model (2), fit separately on data from each year from 2012 to 2019. The blue line denotes 0. We see that each posterior distribution covers both positive and negative values.

In Figure 13 we plot the posterior distribution of xwOBA over the course of a game according to our model fit separately on data from each year from 2012 to 2019. We see that expected wOBA increases steadily over the course of a game, without significant discontinuity (in particular, significant upward discontinuity) between times through the order. The 2018 season is the only season in which we see an upward discontinuity in the posterior means, which occurs between 2TTO and 3TTO. This discontinuity, however, lies inside of the credible intervals and so is not significant.

Figure 13:

Trend in expected wOBA over the course of a game for an average batter facing an average pitcher of the same handedness on the road, according to the model from Equation (2) fit on separately on data from each year from 2012 to 2019. The white dots indicate the posterior means of the expected wOBA values, the thick black error bars denote the 50 % credible intervals, and the thin black error bars denote the 95 % credible intervals.

Appendix E: Alternative models

E.1 A more flexible model: the indicator model

In Equation (2) we model pitcher decline over the course of a game as the combination of discontinuous decline from each TTO to the next and continuous linear pitcher decline across all the batters. A more flexible model would not enforce a particular functional form on within-game pitcher decline. In particular, the most flexible model has a separate coefficient for each batter t ∈ {1, …, 27},

(22) log P ( y i = k ) P ( y i = 1 ) = ∑ t = 1 27 α t k I t i = t + x i ⊤ η k .

With this more flexible model, the qualitative results of our study don’t change. For instance, as in Figure 4, in Figure 14 we plot the posterior distribution of the trajectory of expected wOBA over the course of a game, according to the indicator model from Equation (22) fit on data from 2017. We do not see a significant discontinuity in pitcher performance from one TTO to the next. In other words, we don’t find evidence of a strong batter discontinuity between times through the order. This trend is persistent across each year from 2012 to 2019.

Figure 14:

Trend in expected wOBA over the course of a game in 2017 for an average batter facing an average pitcher of the same handedness on the road, according to the indicator model from Equation (22). The white dots indicate the posterior means of the expected wOBA values, the thick black error bars denote the 50 % credible intervals, and the thin black error bars denote the 95 % credible intervals.

E.2 A more elaborate model: pitcher-specific and batter-specific effects

In our model from Equation (2), we make the simplifying assumption that the trajectory of within-game pitcher deterioration is the same across all pitchers and batters. Nonetheless, it is likely that pitcher performance declines at different rates for different players. To account for such heterogeneity, we extend our model by introducing player-specific rates of decline. Specifically, we model

(23) log P ( y i = k ) P ( y i = 1 ) = α 0 k p ( i ) + α 1 k p ( i ) t i + β 2 k b ( i ) I t i ∈ 2 TTO + β 3 k b ( i ) I t i ∈ 3 TTO + x i ⊤ η k ,

where p(i) is the index of the pitcher and b(i) is the index of the batter in at-bat i. The pitcher-specific continuous decline parameters and batter-specific discontinuity parameters have Gaussian priors,

(24) α 0 k p ( i ) ∼ N α 0 k , σ 0 k 2 , α 1 k p ( i ) ∼ N α 1 k , σ 1 k 2 , β 2 k b ( i ) ∼ N β 2 k , σ 2 k 2 , β 3 k b ( i ) ∼ N β 3 k , σ 3 k 2 ,

which themselves have priors,

(25) α 0 k , α 1 k , β 2 k , β 3 k ∼ N ( 0,25 ) , σ 0 k 2 , σ 1 k 2 , σ 2 k 2 , σ 3 k 2 ∼ half N ( 0,1 ) .

With this more flexible model, the qualitative results of our study don’t change. For instance, as in Figure 4, in Figure 15 we plot the posterior distribution of the trajectory of expected wOBA over the course of a game, according to the player-specific model from Equation (23) fit on data from 2017. In particular, we use the posterior distributions of the prior means α _0k, α _1k, β _2k, and β _3k to compute the xwOBA trajectory for an average pitcher facing an average batter. We do not see a significant upwards discontinuity in expected wOBA from one TTO to the next. In other words, we find little evidence for a strong batter discontinuity between times through the order. This trend is persistent across each year from 2012 to 2019.

Figure 15:

Trend in expected wOBA over the course of a game in 2017 for an average batter facing an average pitcher of the same handedness on the road, according to the model from Equation (23). The white dots indicate the posterior means of the expected wOBA values, the thick black error bars denote the 50 % credible intervals, and the thin black error bars denote the 95 % credible intervals.

References

Brown, L. D. 2008. “In-Season Prediction of Batting Averages: A Field Test of Empirical Bayes and Bayes Methodologies.” Annals of Applied Statistics 2 (1): 113–52. https://doi.org/10.1214/07-aoas138.Search in Google Scholar

Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell. 2017. “Stan: A Probabilistic Programming Language.” Journal of Statistical Software 76 (1): 1–32. https://doi.org/10.18637/jss.v076.i01.Search in Google Scholar PubMed PubMed Central

Fangraphs. 2021. wOBA and FIP Constants. https://www.fangraphs.com/guts.aspx?type=cn.Search in Google Scholar

Gelman, A., and D. B. Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science 7: 457–72. https://doi.org/10.1214/ss/1177011136.Search in Google Scholar

Greenhouse, J. 2011. Spitballing: Fourth Time’s the Harm. https://www.baseballprospectus.com/news/article/13117/spitballing-fourth-times-the-harm/.Search in Google Scholar

Laurila, D. 2015. Managers on the Third Time through the Order. https://blogs.fangraphs.com/managers-on-the-third-time-through-the-order/.Search in Google Scholar

Lichtman, M. 2013. Baseball ProGUESTus: Everything You Always Wanted to Know about the Times through the Order Penalty. https://www.baseballprospectus.com/news/article/22156/.Search in Google Scholar

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.Search in Google Scholar

Rivera, J. 2020. Rays’ Kevin Cash Explains Decision to Pull Blake Snell in World Series: ’I Regret it Because it Didn’t Work Out’. https://www.sportingnews.com/us/mlb/news/kevin-cash-blake-snell-world -series-explained/lfnyfc4nqwys1pcncc2lnyjho.Search in Google Scholar

Slowinski, P. 2010. wOBA. https://library.fangraphs.com/offense/woba/.Search in Google Scholar

Stan Development Team. 2022. RStan: The R Interaface for Stan.Search in Google Scholar

Tango, T., M. Lichtman, and A. Dolphin. 2007. The Book: Playing the Percentages in Baseball. Washington, D.C.: Potomac Books.Search in Google Scholar

Received: 2022-12-12

Accepted: 2023-05-17

Published Online: 2023-06-27