1. Introduction

User complaints regarding public transit fare must be addressed as a primary issue in transport policy. When these problems affect the daily commute of specific groups that depend mostly on public transit, response is needed from the competent authorities to prevent inequity. Fare inequities can be the result of flaws in the operational characteristics of the transit network, which affect the fare charged for the service.

Equity can be referred as the distribution of impacts (benefits and costs) and whether that distribution is considered fair and appropriate [1]. Equity in bus public transportation fares is difficult to accomplish when control and information are insufficient, particularly when the area of service is extensive and presents high operational route diversity. The bus public transportation network in Costa Rica is regulated nationwide and consists of nearly 750 routes not separated in regions or urban areas; therefore, due to the regulation structure, there is a considerable lack of information and control from the local authorities.

Previous studies in Costa Rica [2-4] showed that some routes present operational characteristics that affect the price charged to the users, increasing inequity for vulnerable population.

The public transit system in Costa Rica is private-operated and the fare is stablished by the Public Service Regulator Authority. The fare model is complex and each private operator request revisions of the fare individually, which has resulted in significant differences on the fare by kilometer between routes that are otherwise similar. To address this issue, a methodology to identify the routes that have fares much higher than expected, after considering operational parameters, is proposed. Solving these differences can contribute to increase equity and efficiency in the bus public transportation system.

Clearly, to compare the fare between different transit routes, these should have similar characteristics. Route grouping by similar characteristics, using existing information, allows the identification of outliers within groups. Multilevel or hierarchical modeling was chosen as the method for the group evaluation, since it allows to control for the natural correlation existing between routes of the same group [5] while also taking advantage of the nested data structure to improve model estimation [6]. More importantly, non-hierarchical models are usually inappropriate for hierarchical data because, with few parameters, they generally cannot fit accurately large datasets as the one use in this study, whereas with many parameters, they tend to overfit the data [5].

To identify inequities in the transit routes, the difference between the observed fare per kilometer and what is expected in similar routes (after controlling for several operational characteristics) is proposed. The routes that have higher differences between observed and expected fare per kilometer are the ones that charge higher operational costs into the fare, even after controlling for different operational characteristics. Each category or group sets its own value for comparison, according to the estimates. The larger outliers within each group are the ones to be addressed first as they impose the most expensive fares after controlling for operational characteristics. This procedure is analogous to the one use to identify sites with excess crash frequency in highway safety [7]. As in the case of highway safety, the routes with the largest fare excess are the ones expected to have the highest probability of reduction.

A Bayesian approach is chosen because it provides the flexibility to model complex correlation structures as the ones included in multilevel models. Further, the Bayesian approach allows to easily compare different modeling approaches within the same framework. Bayesian models have been gaining popularity in several fields as the approach of choice when modeling multiple levels and incorporating random effects or complicated dependence structures [8,9].

The purpose of this study is to propose a method to identify transit routes with significantly higher fares compared to similar routes by applying Bayesian modeling and creating different route groups or clusters. These differences could represent inequities or inefficiencies in the transit system that should be studied further.

This paper is organized as follows. First, a literature review regarding equity, multilevel modeling, and cluster analysis is presented. Then, the methodology used for the research is described, followed by the presentation of the data, which includes information about the background of the bus transit network in Costa Rica, the group creation and the parameters used for the models. Finally, the results, which conduct a comparison of goodness-of-fit and estimates between generalized linear models and random coefficients models, are presented and discussed; followed by the conclusions and recommendations for future research.

2. Literature review

Equity has been widely explored regarding transit policies. It has important political implications, and it needs to be considered so that the whole transportation system [10] and the individual bus system revenue can be improved [11]. The proposed methodology doesn’t involve the formulation of policies but seeks to identify the routes in which equity can be achieved while improving efficiency [12].

The proposed model estimates the fare-per-kilometer variable with covariates such as ridership, vehicle occupancy and route length. Other studies have also considered the effects of operational characteristics of the bus public network on fare structure. Ling [13] estimated, given the fares and ridership in the flat fare system, how the total ridership, operator revenue, passenger-km, and consumer surplus would change if the fare structure changes to fare differential system. Fairness and equity, increased revenue and increased ridership were identified as some reasons for promoting differentiated fares.

More recently, Liu et al [14] analyzed the effects of differential fare strategies on social welfare. The authors found that all differential fare strategies produced higher social welfare than the flat fare. Similarly, Tang et al [15] proposed an optimization of bus line fares using elastic demand. The authors maximize social welfare by proposing several operational strategies. The proposed methodology seeks to improve these and other aspects by identifying operational parameters that contribute to inefficiency within groups.

In order to avoid the biases due to analyzing systems with different sizes and operating environments, Karlaftis and McCarthy [16] used cluster analysis to classify 256 transit systems into homogenous groups. The 742 routes used in this analysis required a similar grouping in order to improve model performance. For regression of clustered data, random coefficients models are usually considered, as they control for between-cluster variation [17].

As mentioned above, the route classification allows a multilevel approach by analyzing two levels of data: route-level and category-level. Multilevel modeling has been applied to transportation in different areas of study, but none have explored the effect of operational parameters in fare estimation. Cervero and Kang [18] analyzed the impact in land-use changes and land values due to the upgrade of BRT services in Seoul Korea. The analysis was performed with multilevel modeling, in which the first level corresponded to the parcels and the second level corresponded to the neighborhood groups. Similar individual and neighborhood levels were used by Yavuz et al [19] and Paez and Mercado [20]. Yavuz et al [19] performed a multilevel approach to analyze how perceptions of bus and train safety in Chicago, Illinois vary as a function of person-level characteristics and neighborhood-level characteristics. Paez and Mercado [21] determined individual and neighborhood characteristics that affected the distance traveled and the variability of these factors on each mode type using multilevel analysis. Wang et al. [21] analyzed a Demand Responsive Transport (DRT) System in Manchester; they explored the relationship between a range of socioeconomic variables and service area factors and the demand for DRT using a linear multilevel model. In this model the first level was the census track and the second level the service area.

The studies mentioned above applied multilevel modeling to population or land distribution parameters as parcels, neighborhoods, census tracks or service areas. This research applied multilevel modelling to the transportation system itself, in which the reference level is the whole network, while subnetworks represent the subsequent levels. Ma and Lebacque [22] applied multilevel modeling using similar reference and subsequent levels to study system optimal routing for public transit systems.

Generalized linear models and random coefficients models have been applied to clustered observations. Similar approaches were used by Laird and Ware [23] for longitudinal data and by Aguero-Valverde and Jovanis [24] for crash frequency models. Laird and Ware [23] explored a general family of two-stage random-effects models in two examples from an epidemiological study of the health effects of air pollution. Aguero-Valverde and Jovanis [24] implemented multilevel models in road segments of different functional types by using a full Bayes hierarchical approach to analyze spatial correlation in road crash models in Pennsylvania and Washington. The research concluded that random effects significantly improved the precision, particularly for small sample sizes and low sample means.

Methods

The models are estimated using a full Bayes hierarchical approach in order to estimate not only the parameters but also the estimated excess fare of each transit route and its statistical significance. In Bayesian inference, the posterior distribution of the parameters of interest is estimated as the product of the prior distribution times the likelihood of the model up to a constant, based on the Bayes Theorem. For more details on Bayesian inference, interested readers can refer to Gelman et al [5], Congdon [6] or Koch [25].

At the first level of hierarchy, the logarithm of the fare per kilometer is assumed normally distributed:

where FKM _ik is the fare per kilometer for route i of group k, is the expected value (i.e. the mean) for route i of group k and T is the precision for the normal distribution (i.e. the inverse of the variance).

For this model, the second level of the hierarchy defines the mean:

This formulation is equivalent to a linear model:

where the error term is normally distributed:

The hyper-prior for 3fb7" T is supposed gamma:

The selection of the prior distribution for the betas differentiates between and independent and correlated model. The prior distributions in the case of the independent coefficient model are:

The coefficients for each group of routes are completely independent and the model is almost equivalent to estimate a single model for each group. The only “shared” information among groups is the random variability between transit routes, since they share the same T as shown in (5).

For the random coefficients model (i.e. correlated model), the coefficients change for each route group, but belong to the same normal distribution:

This correlated or random coefficients model allows for “shared” information among route groups through the precision parameter, which in turn, improves parameter estimation. Finally, to identify inequities and inefficiencies in the transit routes the difference between the observed fare per kilometer and what is expected based on similar routes (i.e. belonging to the same group and controlling for the covariates) is estimated:

Clearly , the highest deltas express the highest excess or inequities between transit routes thatshould, in theory, have very similar fares per kilometer. Here, another of the advantages of FullyBayesian models is evident since the full posterior distribution of the deltas is known, which allowsto explore the statistical significance of the excess fare.

4. DATA

The data for the study was provided by ARESEP (Regulatory Authority of Public Services). Several details about public transport background in Costa Rica are needed to have a better understanding of the methods applied to the dataset of the bus routes. Urban development in Costa Rica is concentrated in the Central Valley, the natural barriers serve as borders for the Great Metropolitan Area (GAM), which main city is San José, the capital.

In terms of transport network, cities in Costa Rica can be hierarchized in four categories:

1. 1. San José
2. 2. Main districts of the provinces within the GAM: Heredia, Cartago and Alajuela
3. 3. Main districts of the main cities outside the GAM
4. 4. Main districts of the secondary cities outside the GAM

The transport network revolves around the capital; San José is the distribution center for all the country’s regions, with about 33 % of the total routes connecting San José directly with other cities, within or outside the GAM. Each of the cities mentioned above present their own transport system, this represents about 65 % of the routes of the country. The last 2 % is composed of the routes that communicate main districts between themselves, without stopping in San José. As noted, the transport system in Costa Rica presents substantial size and operational differences according to its location.

To perform the analysis adequately, 742 bus public transportation routes were classified in 30 groups according to four main operation characteristics: location, route structure, terrain and length. Location was the main parameter for the classification. The other three parameters depend on the route configuration. For route structure, according to Molinero and Sánchez [27] and Vuchic [28], routes were classified in radial or circular configuration. Almost 95 % of the routes in Costa Rica are radial. The terrain was classified according to the Central American Manual for Geometric Design of Highways and Streets [29] in Level – Rolling (0 % to 10 % ) and Rolling – Mountainous (> 10% ). Last, the routes were classified by length, in order to differentiate the ones that operate in the inner city, from the ones that reach the outer city or travel between cities. They were classified in short (0 km to 20 km), medium (20 km to 40 km), long (40 km to 125 km) and very long (more than 125 km).

Figura 1 shows the summary of the classification method described above. Location is the firstparameter for grouping the routes, then each unit is assorted to a different category according toits characteristics from the subsequent three levels. For example, one of the final groups consistsof the units of San José, with radial structure, level – rolling terrain and short length. Not everypossible combination was defined as a group, each category was evaluated by size and relevanceto decrease the total number of groups.

The classification parameters returned 30 categories from the 742 routes analyzed, some groups consist of less than 10 routes, in contrast with groups with more than 50 routes per category. Therefore, models such as random coefficients models, that controls for between-cluster variation, must be considered [17]. For more details about the classification of bus routes interested readers can check [4].Once the categories were established, the variable fare per kilometer was selected to model inequities within each group. As a result of exploratory analysis and linear regression models, three independent variables were included in the model:

1. 1. Route length: route length measured in kilometers.
2. 2. Monthly ridership: total average passengers per month for each route.
3. 3. Vehicle occupancy: monthly ridership divided by the total average number of trips per month.

Several other variables can be included in the models to explain the differences in fare per kilometer. In the case of this study, these three operation variables were used as they were available for almost all the bus transit routes in Costa Rica. In addition to the variables used in the categorization of the bus routes, these variables should explain most of the variability in the data. Additional variability in the data should be represented by the excess fare per kilometer or delta, previously introduced in (9). Further analysis of deltas will reflect on the routes with the highest operational costs per user.

From the summary statistics in Tabla 1, it is observed that monthly ridership has high dispersion, since the standard deviations are higher than the mean frequencies for many of the groups and the total routes. The bus route length also presented high dispersion on the aggregated statistics but the dispersion inside each group is significantly reduced. This was expected since the route length was one of the criteria used to create the groups. The fare per kilometer and the vehicle occupancy presented small variance both at the aggregated level and at the group level.

Results

Models were estimated using the open source software OpenBUGS [30]. OpenBUGS is a popular software used for Full Bayesian estimation in fields ranging from transportation to medicine. For more details on Bayesian estimation using OpenBUGS interested readers can check Lawson et al [31], Congdon [32] or Ntzoufras [33].

For the models, 1.000 iterations were discarded as burn-in. The following 175.000 iterations were used to obtain summary statistics of the posterior distribution of parameters. Convergence was assessed by visual inspection of the Markov chains for the parameters. Furthermore, the number of iterations was selected so that the Monte Carlo error for each parameter in the model would be less than 5 % of the value of the standard deviation of that parameter.

Independent and correlated or random parameters models

Table 2 presents the estimates and the standard deviation for each coefficient for the independent and correlated random coefficients model. The main interest is to determine which model performs better in terms of goodness-of-fit and coefficient significance.

The goodness-of-fit measures commonly used in full Bayesian statistics are presented in thetable: the posterior mean of the deviance and the deviance information criterion (DIC)[26]. The deviance is estimated in the same way for frequentist and Bayesian statistics, while the DIC is the Bayesian equivalent of the Akaike information criterion. As in the case of their frequentist counterparts, the deviance and the DIC quantify the relative goodness-of-fit of the models; therefore, they are useful for comparing models.

Regarding the results in Table 2, the estimates, in both the independent and the correlated models, corresponding to the independent variables have a negative effect over the fare per kilometer. As expected, when the route length (\beta_1), average occupancy (\beta_2) and monthly ridership (\beta_3) increase, the fare per kilometer decreases. According to (3) , the coefficients of the variables correspond to the constant elasticities since they relate the natural logarithm of the fare per kilometer with the natural logarithm of the covariates. Therefore, it is observed that route length has higher influence over the fare per kilometer than any other independent variable, since in general it has the highest coefficents. Concerning the estimates significance, both models present a similar behavior. The random coefficients model performs slightly better with the route length estimates, which, as mentioned previously, have more influence in the general model. On the other hand, \beta_2 and \beta_3, are significant just in one more group in the independent model compared to the correlated or random parameters model.

The main difference between models can be observed at the estimates for the mean and standard deviation. The random coefficients model return values that are closer to the mean, as they consider the whole sample for modelling and not just each group sample. For the same reason, the standard deviation is also reduced, resulting in a better estimation of the fare per kilometer. In terms of the goodness-of-fit measures, the DIC for the random parameters model is 22 points lower, which means that the model is s significantly better [26].

From Table 2it is observed that the estimates for groups 13, 23 and 30 are not significant. Even though, these categories have a small sample size (17, 14 and 3 respectively) the groups 1, 10 and 15 have also small sample sizes (5, 9, and 22 respectively) but they return one or more significant variables. The random coefficients model, as mentioned before, has a stronger influence in reducing the small sample groups mean and standard deviation.

The estimated kernel density of the posterior distribution of the parameters further shows the benefits of random coefficients models over univariate models. The estimated density of \beta_0 for the groups with the largest and the smallest sample size are presented in Figure 2. The density in group 2, with a sample size of 90, shows almost no changes between the univariate and the multivariate model. On the other hand, the density in Group 30, with a sample size of 3, shows how by considering the population and not just the within-group variation sample, the multivariate model improves the model prediction compared to the univariate model, particularly for groups with small samples.

Significant Deltas Obtained from the Independent and Random Coefficients Models

Bayesian modeling allows to estimate the significance of the deltas obtained after running both models. Table 3 shows the 5 % of the higher deltas obtained for all the population in the random coefficients model, and as seen, the independent model has seven not significant values, while the correlated model has three. A positive delta indicates a value of fare per kilometer higher than expected, these are the routes that have to be revised first in order to set an adequate fare in terms of operation, equity and the category in which the route is included.

The deltas are ranked from high to low and compared to the ranking obtained in the univariate model. Coincidentally, the highest 5 % deltas in the random parameters model are the same as in the independent model, though they are not in the same order of hierarchy. The delta ranked in the 4th position in the random parameters model serves as an example as how not considering the significance can mislead the interpretation of the results.

Once the outliers are identified, it is necessary to perform individual analysis to each route to identify which are the operational characteristics that have price implications for the served population. This method is intended to be used by the competent authorities as a first approach to perform a diagnostic of the whole network, so that, the resources for specific audits can be focused to solve the most problematic routes in terms of equity and efficiency.

Ridership and origin-destination studies can be performed to determine the necessity of improving the route’s design in terms of structure, configuration, frequency or schedule. The change of this parameter will imply a variation in the fare charged for the service and therefore modify the ridership and average occupancy. The route improvement must be accompanied by model optimization, seeking to conform categories with more homogeneous behavior.

Figure 3 makes a comparison between both model rankings around the unitary diagonal, showing that the models present slight differences in the deltas ranking, as a considerable number of points are close to the diagonal; nevertheless, several values, specifically those of lower ranking for the univariate model are far from the unitary reference.

Conclusions and recomendations for future research

• Univariate and multivariate models were estimated using a full Bayes hierarchical approach. The models estimated the effect of route length, average occupancy and monthly ridership in fare per kilometer, for 742 clustered routes of the bus network in Costa Rica. The estimates in both models have a negative effect over the fare per kilometer. An increase of the independent variables leads to a decrease in the fare per kilometer, being the route length the parameter with the higher influence overall.

• The random coefficients (correlated) model performed better than the univariate model in terms of goodness-of-fit statistics; the DIC for the multivariate model was 22 points lower, which means that the multivariate model is significantly better. The random coefficients model estimates usually returned smaller standard deviations, because it considers the whole sample for modelling and not just each group sample, resulting in a better estimation of the fare per kilometer.

• The classification of the data in groups enables a better prediction of the within-group mean, standard deviation and outliers. Outliers allow to identify routes that behave significantly different from those in which they are grouped; therefore, individual analyses should be performed to determine the operation parameters that need to be adjusted to improve performance. Positive deltas, that show the difference between the observed values and the estimated ones, indicate a value of fare per kilometer higher than expected, hence, these are the routes that should be revised first in order to set an adequate fare in terms of efficiency and equity.

• Even though fare is just one of many parameters in a transit system, it has a direct effect on the users in terms of equity, and as stated by Ling [13], improving the fare structure can result in increasing ridership and revenue. With better data availability, similar analyses can be conducted in different bus transit networks as a tool for evaluating the efficiency of each system. Other parameters can be introduced as independent variables for determining specific efficiency aspects.

Acknowledgements

The Authors thank ARESEP (Regulatory Authority of Public Services) for providing data on bus public transportation routes.

References

T. Litman, Evaluating transportation equity - guidance for incorporating distributional impacts in transportation planning. Victoria Transport Policy Institute, Victoria, BC, 2016.

Programa de Investigación en Desarrollo Urbano Sostenible. Elaboración de auditorías de demanda y cálculo de indicadores y parámetros operativos del servicio de transporte remunerado de personas, modalidad autobús. Universidad de Costa Rica, San José, Costa Rica, 2014.

Programa de Investigación en Desarrollo Urbano Sostenible. Elaboración de auditorías de demanda y cálculo de indicadores y parámetros operativos del servicio de transporte remunerado de personas, modalidad autobús que brindan el servicio en los corredores San José – Heredia y San José-Moravia. Universidad de Costa Rica, San José, Costa Rica, 2015.

C.S. Quirós-Calderón and J. Aguero-Valverde, “Clasificación de las rutas de la red de transporte público modalidad autobús de Costa Rica”, Revista Ingeniería, Vol. 28, no. 2, pp.: 112-136, 2018

A. Gelman, J. B. Carlin, H. S. Stern and D. B. Rubin., Bayesian Data Analysis, Chapman y Hall/CRC Texts in Statistical Science, 2003.

P. Congdon, Bayesian statistical modelling. Chichester, UK: John Wiley y Sons, 2007.

AASHTO, Highway Safety Manual. Washington DC, USA. 2010.

E. Hauer, “Identification of sites with promise”, Transportation Research Record, 1542, pp. 54-60. ,1996

S. Banerjee, B.P. Carlin, A.E. Gelfand, Hierarchical modeling and analysis for spatial data. Chapman y Hall/CRC, Florida, 2004.

M. Garrett, and B. Taylor, “Reconsidering social equity in public transit”, Berkeley Planning Journal, 13.1, pp. 6-27, 1999.

C. Nuworsoo, A. Golub and E. Deakin, “Analyzing equity impacts of transit fare changes: Case study of Alameda–Contra Costa Transit, California”, Evaluation and Program Planning, vol. 32, no 4, pp. 360-368, 2009. https://doi.org/10.1016/j.evalprogplan.2009.06.009.

R. Cervero, “Flat versus differentiated transit pricing: what’s a fair fare?”, Transportation, vol. 10, no 3, pp. 211-232, 1981..

J. H. Ling, “Transit fare differentials: A theoretical analysis”, Journal of advanced transportation, vol. 32, no 3, pp. 297-314, 1998. http://dx.doi.org/10.1002/atr.5670320304.

B.Z. Liu, Y.E. Ge, K. Cao, X Jiang et al., “Optimizing a desirable fare structure for a bus-subway corridor”, PLoS ONE, vol 12, no 10, 2017. e0184815.

C. Tang, A.A. Ceder, and Y.E. Ge, “Integrated optimization of bus line fare and operational strategies using elastic demand”, Journal of Advanced Transportation, 2017. https://doi.org/10.1155/2017/7058789

M. Karlaftis, and P. McCarthy, “Cost structures of public transit systems: a panel data analysis”, Transportation Research Part E: Logistics and Transportation Review, vol. 38, no 1, pp. 1-18, 2002.. https://doi.org/10.1016/S1366-5545(01)00006-0.

N.T. Longford, “Logistic regression with random coefficients”, Computational Statistics y Data Analysis, vol. 17, no 1, pp. 1-15, 1994. https://doi.org/10.1016/0167-9473(92)00062-V.

R. Cervero, and C. D. Kang. “Bus rapid transit impacts on land uses and land values in Seoul, Korea” Transport Policy, vol. 18, no 1, pp. 102-116, 2011. https://doi.org/10.1016/j.tranpol.2010.06.005

N. Yavuz, E. Welch and P. Sriraj, “Individual and neighborhood determinants of perceptions of bus and train safety in Chicago, Illinois: Application of hierarchical linear modeling”, Transportation Research Record: Journal of the Transportation Research Board, no 2034, pp. 19-26, 2007. https://doi.org/10.3141/2034-03.

A. Paez, and R. G. Mercado, “Mobility of Canadian Elderly: Multilevel Analysis of Distance Traveled in the Hamilton Census Metropolitan Area, Ontario, Canada”. Presented at 86th Annual Meeting of the Transportation Research Board, Washington, D.C., 2007.

C. Wang, M. Quddus, M. Enoch, T. Ryley and L. Davison, “Multilevel modelling of Demand Responsive Transport (DRT) trips in Greater Manchester based on area-wide socio-economic data”, Transportation, vol. 41, no 3, pp. 589-610, 2014.. https://doi.org/10.1007/s11116-013-9506-1.

T.Y. Ma, and J. P. Lebacque. “Dynamic system optimal routing in multimodal transit network.” Transportation Research Record: Journal of the Transportation Research Board, no. 2351, 2013, pp. 76-84, 2013.. https://doi.org/10.3141/2351-09.

N.M. Laird, and J. H. Ware, “Random-effects models for longitudinal data”, Biometrics, pp. 963-974, 1982. https://doi.org/10.2307/2529876.

J. Aguero-Valverde, and P. Jovanis, “Spatial correlation in multilevel crash frequency models: Effects of different neighboring structures”, Transportation Research Record: Journal of the Transportation Research Board, no. 2165, pp. 21-32, 2010. https://doi.org/10.3141/2165-03.

K.R. Koch, Introduction to Bayesian statistics. Springer Science y Business Media, 2007.

D.J. Spiegelhalter, N.G. Best, B.P. Carlin and A. Van Der Linde, “Bayesian measures of model complexity and fit”, Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 64, no. 4, p. 583-639, 2002. https://doi.org/10.1111/1467-9868.00353.

A. Molinero, and L. Sánchez Arellano, Transporte público. Universidad Autónoma del Estado de México, Ciudad de México, 1997.

Vuchic, V., Urban transit systems and technology. John Wiley y Sons, Hoboken, NJ, 2005.

Secretaría de Integración Económica Centroamericana (SIECA), Manual Centroamericano de Normas para el Diseño Geométrico de las Carreteras Regionales con enfoque de Gestión de Riesgo y Seguridad Vial. 3ª Edición. Ciudad de Guatemala, Guatemala, 2011.

A. Thomas, B. O’Hara, U. Ligges and S. Sturtz. “Making BUGS open”, R news, vol. 6, pp. 12-17, 2006.

A.B. Lawson, W. J. Browne and C. L. V. Rodeiro. Disease mapping with WinBUGS and MLwiN. John Wiley y Sons, Chichester, UK, 2003.

P. Congdon, “The Basis for, and Advantages of, Bayesian Model Estimation via Repeated Sampling”, in Applied Bayesian Modelling. John Wiley y Sons, Chichester, UK, 2003, pp. 1-30.

I. Ntzoufras, Bayesian modeling using WinBUGS. John Wiley y Sons, Hoboken, NJ, 2009.