Bayesian influence diagnostics using normalizing functional Bregman divergence

Bayesian influence diagnostics using normalizing functional Bregman divergence

Ian M Danilevicz Department of Statistics, Federal University of Minas Gerais, Brazil Ricardo S Ehlers Department of Applied Mathematics and Statistics, University of São Paulo, Brazil

Ideally, any statistical inference should be robust to local influences. Although there are simple ways to check about leverage points in independent and linear problems, more complex models require more sophisticated methods. Kullback-Leiber and Bregman divergences were already applied in Bayesian inference to measure the isolated impact of each observation in a model. We extend these ideas to models for dependent data and with non-normal probability distributions such as time series, spatial models and generalized linear models. We also propose a strategy to rescale the functional Bregman divergence to lie in the (0,1) interval thus facilitating interpretation and comparison. This is accomplished with a minimal computational effort and maintaining all theoretical properties. For computational efficiency, we take advantage of Hamiltonian Monte Carlo methods to draw samples from the posterior distribution of model parameters. The resulting Markov chains are then directly connected with Bregman calculus, which results in fast computation. We check the propositions in both simulated and empirical studies.

Key-words: Bayesian inference, functional Bregman divergence, influential observations, Hamiltonian Monte Carlo.

1 Introduction

After fitting a statistical model we need to investigate whether the model assumptions are supported. In particular, inference about parameters would be weak if it is influenced by a few individual results. In this paper, we make use of a new diagnostic analysis tool for detecting influential points. The idea is to adapt the functional Bregman divergence to compare two or more likelihoods (Goh and Dey (2014)) to the context of measuring how influent is each observation in a given model. An influential point consists of an observation which strongly changes the estimation of parameters. The classical example is a point which drastically alters the slope parameter in a linear regression. In Bayesian inference, our focus lies into the whole posterior distribution instead of a single parameter. Indeed, seeking for leverage effect in many parameters from a complex model seems unfeasible.

Bayesian inference should produce a posterior distribution based on Bayes theorem. So, if there is a function which measures the distance between two probability densities we can measure distance between two posterior distributions or between a Bayesian model and its perturbed version. The perturbed case could consist of the same sample without an element if we have identical and independent observations (Goh and Dey (2014)) or it may be a sample with an imputed element if we work with dependent models (Hong-Xia et al. (2016)). We can use a well know function such as Kullback-Leiber to measure the divergence between two posterior distributions as well as the functional Bregman divergence, which is a generalization of the previous one.

In the applications we have in mind, the posterior distributions are not available in closed form and we resort to Markov chain Monte Carlo (MCMC) methods to obtain approximations for parameter estimates and detection of influential observations. All the necessary computations in this paper were implemented using the open-source statistical software R (R Core Team (2015)). In particular, the rstan package which is an interface to the open-source Bayesian software Stan (Stan Development Team (2016)) was used to draw samples from the joint posterior distributions. Stan is a computing environment for implementing Hamiltonian Monte Carlo methods (HMC, Neal (2011)) coupled with the no-U-turn sampler (NUTS) which are designed to improve speed, stability and scalability compared to standard MCMC schemes. Typically, HMC methods result in high acceptance rates and low serial correlations thus leading to efficient posterior sampling.

The remainder of this paper is structured as follows. In Section 2, the models used to illustrate the application of our propositions are briefly reviewed and the associated prior distributions are described. The Hamiltonian Monte Carlo sampling scheme is also described here. In Section 3 we introduce the functional Bregman divergence and describe its use to detect influential observations in models for both independent and dependent data. Section 4 consists of simulation studies where we perform sensitivity analysis and investigate how accurately we can detect influential observations. Section 5 summarizes empirical studies in which we illustrate our proposed methodology applied in real data. A discussion in Section 6 concludes the paper.

2 Models

This section is dedicated to describe the models which we used around the paper and the HMC sampling scheme adopted.

2.1 Generalized Linear Models

Generalized linear models (GLM, Nelder and Wedderburn (1972)) are used here to illustrate applications of our methods in models for dependent data with non-normal distributions. Let conditionally independent, where the distribution of each belongs to the exponential family of distributions, i.e


The density in equation (1) is parameterized by the canonical parameter and and are known functions. Also, is related to regression coefficients by a monotone differentiable link function such that . The linear predictor is , where is the design matrix and is a vector of regression coefficients.

The likelihood function based on model (1) is given by,

This class of models includes several well kown distributions such as Poisson, Binomial, Gamma, Normal and inverse Normal.

2.2 Spatial Regression Models

We chose to illustrated our methods using spatial regression models (SRM) as a kind of geostatistical data model (Gaetan and Guyon (2010)). The model can be represented as,


where is the response of the observation , is the value of at the coordinate x, is the value of at the coordinate y and is an error, usually assumed . Most commonly, x and y are latitude and longitude however they could express as angles. If we assume normality of the errors the likelihood function can be expressed as follows,

where is a design matrix with the following columns: ones, coordinate x, coordinate y, interaction of x and y, squared x and squared y. The matrix describes the covariance structure between the observations and in the response vector. A suitable set of priors consists of assuming that and the variance-covariance matrix follows a inverse Wishart .

A particular case of SRM consists of assuming an independence variance structure, i.e. , where is an identity matrix. In this case default prior distributions for the could be Inverse Gamma, or Gamma distributions if we are not restricted to conjugate priors.

2.3 GARCH Model

The generalized autoregressive conditional heteroscedasticity (GARCH) model (Bollerslev (1986)) is the most used class of models to study the volatility in financial markets. The GARCH() model is typically presented as the following sequence of equations,

where is the observed return at time and and are unknown parameters. The are independent and identically distributed error terms with mean zero and variance one. Also, , and define the positivity constraints and , ensures covariance stationarity of .

Given an observed time series of returns the conditional likelihood function is given by,


where and represents the set of all model parameters. In practice, to get this recursive definition of the volatility off the round in Stan we need to impute non-negative initial values for .

Prior distributions for the GARCH parameters were proposed by Deschamps (2006) and also used in Ardia and Hoogerheide (2010), who suggest a multivariate Normal distribution for and truncated to satisfy the associated constraints. However, to avoid truncation we propose a simpler approach and specify the following priors, , and for and respectively.

2.4 Hamiltonian Monte Carlo

Our approach to detect influential observations relies on MCMC methods that should produce Markov chains which efficiently explore the parameter space. This motivates seeking for sampling strategies that aim at reducing correlation within the chains thus improving convergence to the posterior distribution. Hamiltonian Monte Carlo (HMC) comes as a recent and powerful simulation technique when all the parameters of interest are continuous. HMC uses the gradient of the log posterior density to guide the proposed jumps in the parameter space and reduces the random walk effect in the traditional Metropolis-Hastings algorithm (Duane et al. (1987) and Neal (2011)).

For a -dimensional vector of parameters and denoting the posterior density of , the idea is to augment the parameter space whereas the invariant distribution is now a Hamiltonian density given by,

for a normalizing constant . The Hamiltonian function is decomposed as, , where is the potential energy, is the position vector, is the kinetic energy and is the momentum vector in the physics literature. In a Bayesian setup we set .

Trajectories between points are defined theoretically by some differential equations which in practice cannot be solved analytically. So, in terms of simulation a method is required to approximately integrate the Hamiltonian dynamics. The leapfrog operator (Leimkuhler and Reich (2004)) is typically used to discretize the Hamiltonian dynamics and it updates at time as the following steps,

where is a user specified small step-size and is the gradient of with respect to . Then, after a given number of time steps this results in a proposal and a Metropolis acceptance probability is employed to correct the bias introduced by the discretization and ensure convergence to the invariant posterior distribution.

So, using Hamiltonian Monte Carlo involves specifying the number of leapfrogs by iteration, the step-size length and the initial distribution of the auxiliary variable . The choice of an appropriated which associated with will not produce a constant periodicity may be done using the No-U-Turn sampler (NUTS, Hoffman and Gelman (2014)), which aims at avoiding the need to hand-tune and in practice. During the warmup the algorithm will test different values of leapfrogs and step-size and automatically judges the best range to sample. The basic strategy is to double until increasing the leapfrog will no longer enlarge the distance between an initial value of and a proposed value . The criterion is the derivative with respect to time of half the squared distance between the and . To define an efficient value of , NUTS constantly checks if the acceptance rate is sufficiently high during the warmup. If it is not, the algorithm just shortens the step-size at next iteration (see Nesterov (2009) and Hoffman and Gelman (2014)).

It is also worth noting that the Stan programming language provides a numerical gradient using reverse-mode algorithmic differentiation so that obtaining the gradient analytically is not necessary. Finally, the distribution of is a multivariate normal with either a diagonal or a full variance-covariance matrix. The former is usually selected because the precision increase is almost irrelevant compared to the computational memory costs (Stan Development Team (2016)).

3 Functional Bregman divergence

The functional Bregman divergence aims at measuring dissimilarities between functions, and in particular we are interested in comparing posterior distributions. The method is briefly described here and adapted to our models for detection of influential observations. We define as a finite measure space and and as two non-negative functions.

Definition 1.

Let us consider being a strictly convex and differentiable function on . Then the functional Bregman divergence is defined under the marginal density as


where represents the derivative of .

This divergence has some well-known properties (see for example Goh and Dey (2014)), the proofs of which appear in Frigyik et al. (2008a, b).

Clearly, if then and the functional Bregman divergence is zero for any and . However, if we choose a strictly convex then the Bregman divergence will be always greater than zero, except for the trivial case . Furthermore, works as a tuning parameter and increasing its distance from the identity we would have as large as desired no matter the functions and . In this paper, we follow the suggestion in Goh and Dey (2014) and restrict to the class of convex functions defined by Eguchi and Kano (2001),


Three popular choices of are: (Itakura-Saito distance), (Kullback-Leibler divergence), and (squared Euclidean distance or ).

3.1 Perturbation in dependent models

Here we extend the ideas in Goh and Dey (2014) where perturbation was defined in models for independent and identically distributed observations to dependent models. A general perturbation is defined as the ratio of unnormalized posterior densities,


where indicates that likelihood and/or prior suffers some perturbation. Particularly, to assess potential influence of any observation the perturbation is restricted to the likelihood function while keeping the prior unaltered. The associated perturbation is then given by,

where denotes the vector without the th case. In models for dependent data however we can not exclude an observation without modifying the likelihood structure. In any case, the general rule to measure the local influence of the th point is to compute the divergence between and ,

The integral in Equation (5) however is analytically intractable for most practical situations and an approximation is needed. It is convenient to define the normalizing constant for as,

where is the marginal density and is any probability density function. So, given a sample from the posterior distribution (which could be generated by HMC) we can estimate the normalizing constant as,

This is the so called Importance-Weighted Marginal Density Estimate (IWMDE, Chen (1994)). Denoting the resulting posterior distribution as , the approximate perturbed posterior is given by,


where . Consequently, we can approximate the functional Bregman divergence between and by,

which for the convex functions in (6) is simplified as,

where . In particular, for which corresponds to the Kullback-Leibler divergence, we can simplify many terms of the above expression and obtain,


3.2 Normalizing Bregman divergence

When using a functional Bregman divergence to evaluate influential points each and in this scale we might have doubts about one or more values being substantially higher than the others. To facilitate comparison, McCulloch (1989) proposed a calibration which compresses the scale between 0.5 and 1 by making an analogy with the comparison between two Bernoulli distributions one of which with success probabilities equal to 1/2. However, extending this idea to any functional Bregman divergence and comparing any probability distribution with a Bernoulli sounds difficult to justify theoretically. Therefore, we propose a different route to compare the Bregman divergence between two densities, which we call a normalizing Bregman divergence.

Proposition 1.

Given n+1 probability functions we have n divergences between and , which we write as . Then, there is a function for which the sum of the divergences returns one,

and is called a normalizing Bregman divergence.


There is a sequence of functions , which tunes the divergence intensity between any two density functions and . Suppose we have a full probability density and we wish to compare it with each likelihood without th element as to check local influence. We already know that each divergence is positive, so the sum of divergences belongs to positive real domain,

where if and only if is the identity, but also could be arbitrarily high as becomes more and more convex. In particular may be one. ∎

Proposition 2.

Given n+1 probability functions we have n divergences between and , which we write as . There is a mean operator which transforms any Bregman divergence in a normalizing Bregman divergence.


where is called a normalizing Bregman operator.


By the generalized Pythagorean inequality (Frigyik et al. (2008b)), it is natural to supposed order maintenance as for any and under the restriction of strictly convexity. If the above order relation is maintained then we can guarantee that all Bregman divergence with any consists of the same divergence just with a different location scale. ∎

We gather these two arguments together as follows. A finite sum of divergences is finite and just tunes the scale but not the order of Bregman. So there is a special case of , let us call it , for which the sum with respect to a set of densities results in one, and we call this divergence a normalizing Bregman divergence. This is so because all Bregman divergence preserves the same order.

So, the attractiveness of our proposal is that and it is quite intuitive to work in this scale to compare divergences in the context of identifying influential observations. Also, one possible caveat is that a result being high or low would depend on the sample size so that any cut-off point should take into account. In this paper we argue that uUnder the null hypothesis that there is no influential observation in the sample, a reasonable expected normalizing Bregman would be , i.e.

so that we expected each observation would present the same divergence. This bound becomes our starting point to identify influential observations. If any observation returns , then it is a natural candidate to be an influential point, which we must investigate. This should be seen as a useful practical device to seek for influence rather than a definitive theoretical constant which can separate influential from not influential cases. Finding a better cut-off point other than is still as a open problem for future research. Finally, we note that using the Kullback-Leibler divergence which is approximated using Equation (9) leads to faster computations.

4 Simulation Study

In this section, we assess the performance of the algorithms and methods proposed by conducting a simulation study. In particular, we verify whether reliable results are produced and which parameters are the most difficult to estimate. We also check sensitivity to prior specification and the performance for detecting influential observations. We concentrate on the performance of posterior expectations as parameter estimators using Hamiltonian Monte Carlo methods and the Stan package. For all combinations of models and prior distributions we generate replications of data and the performances were evaluated considering the bias and the square root of the mean square error (SMSE), which are defined as,

where denotes the point estimate of a parameter in the th replication, . Finally, for each data set we generated two chains of 4000 iterations using Stan and discarded the first 2000 iterations as burn-in.

4.1 Performance for Estimation and Sensitivity Analysis

We begin with a logistic regression with an intercept and two covariates and simulate data using two parametric sets. The values of the two covariates and were generated independently from a standard normal distribution. The first model (Model 1) has true parameters given by , , while for the second one (Model 2) the true parameters were set to , , and each model was tested for two different sample sizes, and . Finally, inspired by Gelman et al. (2008), we adopted three different prior distributions for the coefficients , as follows. Prior 1: , Prior 2: , Prior 3: . The main results of this exercise are summarized in Table 1. Model 1 with and Cauchy prior presents less bias for most estimations, except for . This prior also leads to lowest SMSE for all parameters. When we observed the same set but with , all estimations get better and again the Cauchy prior provides the least biased estimation. The mean SMSE falls from around 0.100 to approximately 0.025 when the sample size incrases. For Model 2, the results are quite similar.

[ Table 1 around here ]

We now turn to the analysis of the spatial regression model given in (2). Data from two models with an intercept and four covariates were generated where the true coefficients are given by = (Model 1) and = (Model 2), both with the same variance . Each one was tested for two different sample sizes, and . Both latitude and longitude were generated by independent standard normal distributions without truncation, as an hypothetical surface without borders.

Each model was estimated under three different prior specifications for the coefficients , and as follows. Prior 1: , , Prior 2: , , Prior 3: , , where Prior 1 is more flat than the others and follows the suggestion in Chung et al. (2013), Prior 3 is more informative following Stan Development Team (2016) warnings to avoid eventual computational errors and Prior 2 is an intermediate one. The main results of this exercise are summarized in Table 2. In this table, 0.0000 means smaller than 0.0001.

In Model 1 with the priors present bias of the same order for most parameters except for , where Prior 3 is the best, and , where Prior 3 is the worst. The SMSEs are quite similar across all prior specifications. When increases there is some changes in bias order: has the best performance with Prior 3, however , and show a one order decrease with Prior 2, then this is the best prior. For Model 2 and , we see Prior 2 again with better bias results for and , but Prior 3 is better for estimating . With the bias results show an advantage of Prior 1 to estimate and , as well as Prior 2 is better to estimate and , and likewise Prior 3 for . However, the SMSEs were very similar, so that the differences between priors were not so relevant in this last set. Overall, Prior 2 presents the best results.

[ Table 2 around here ]

Our last exercise concerns to GARCH(1,1) models where we generate artificial time series with Normal errors and two different sets of parameters: , , and , , . However we propose to estimate both cases with Normal and Student error terms, even though all series were built using Normal errors. Replacing a Normal by a Student is a commonly used strategy to control overdispersed data. The prior distributions were assigned as follows. Prior 1: , , , Prior 2: , , and Prior 3: , , ,

The results for the GARCH(1,1) with Normal errors are presented in Table 3. We notice that Prior 3 attained the best results for the parameter set 1, but Prior 1 was better in set 2. This outcome happens because Prior 1 is perhaps too informative about and even a value of as large as 900 was not enough for the model to learn from data. However, the different priors assigned to and do not imply in any drastic output change. Table 4 summarizes the output from a GARCH(1,1) model estimated with Student errors. From this table we notice that Prior 1 returned the best results for the parameter set 2, but in set 1 the three priors share similar performances. We now look at both tables in tandem to compare Student and Normal errors. The Normal GARCH presented better results than Student for the parameter set 2, but they were similar in set 1, so that there is no need of a robust model in this case (the series were generated with normal errors).

[ Table 3 around here ]

[ Table 4 around here ]

4.2 Influence Identification

To evaluate the normalizing Bregman divergence as a useful tool to identify point influence we proceed with three simulation sets, each refering to a different model. We use the same models presented in the previous subsection. Within each model we created four scenarios: I without any kind of perturbation, II where there is one perturbed observation, III where there are two influential points and IV with three influential points. The point contamination in time series and spatial models follows the scheme proposed by Hong-Xia et al. (2016) and Cho et al. (2009), i.e., , where is the standard deviation of the observed sample . For the logistic regression however we need a different approach to contaminate data. In this case we simply exchange the output, i.e. if is to be contaminated and then we set , otherwise if we set .

The results for the case influence diagnostic using normalizing Bregman divergence in logistic regression are shown in Table 5. For this table, the true parameter values are , , and the prior distributions are ,  . Also, the perturbation schemes are: I no perturbation, II observation 64 has an additional noise, III observations 44 and 64 present perturbation and IV observations 19, 44 and 64 have an extra noise. The table then shows the estimated (mean and standard deviation) divergences for the three observations, 19, 44 and 64.

We first notice that in the no perturbation scenario the estimated divergences are mostly as expected on average, i.e. and for and respectively. On the other hand, when the output is perturbed the average divergence was between 0.028 and 0.031 for and between 0.008 and 0.009 for . Finally, there is a correlation between mean and standard deviation in the sense that a small value of one corresponds to low estimation of the other.

[ Table 5 around here ]

The results are even more emphatic in spatial regression models as shown in Table 6. In this table, the true parameter values are , , , and we chose Prior 2. Also, the influence scenarios are: I without perturbation, II observation 19 has an additional noise, III observations 15 and 19 present perturbation and IV where 3, 15 and 19 have an extra noise.

For and scenario I the normalizing Bregman divergence has mean around 0.020 in the three observed points, which corresponds to the expected . In scenario II, the estimates for observations 3 and 15 fall to 0.011 and 0.013, because observation 19 was perturbed and its divergence estimate rises to 0.423. In scenario III the estimate for observation 3 falls even more, because both 15 and 19 were perturbed and both have a 0.283 estimate for the divergence. Finally, when the three observations were perturbed they share the impact between 0.210 and 0.214 estimates. For and scenario I, we have again results precisely as expected, i.e. 0.005 compared to . Furthermore, scenarios II, III and IV are quite similar, the mean values are slightly smaller, but this is expected for a larger sample.

[ Table 6 around here ]

The effect is still clear in time series with moderate sample sizes, as shown in Table 7. This table shows results for a GARCH(1,1) model with normal errors and true parameter values , and . Influence scenarios for : I without perturbation, II observation 64 has an additional noise, III observations 44 and 64 present perturbation and IV where 19, 44 and 64 have an extra noise. Influence scenarios to : I without perturbation, II observation 464 has an additional noise, III observations 344 and 464 present perturbation and IV where 119, 344 and 464 have an extra noise.

For and scenario I the normalizing Bregman divergence has mean equal to 0.009 in the three observed points, which is slightly bellow the expected . In scenario II the estimated divergence for observations 19 and 44 fall to 0.007, because observation 64 was perturbed and its estimated divergence rises to 0.247, which might not seem a large value but it is more than 20 times the expected value . In III the estimate for observation 19 falls even more, because both 44 and 64 were perturbed and have estimated divergences 0.188 and 0.192. Finally, when all three observations were perturbed they share the impact with similar estimated divergences. For and scenario I, we have again results around 0.002, i.e. . Furthermore, scenarios II, III and IV show quite similar results, the estimated values being slightly smaller, but this is expected for a larger sample.

[ Table 7 around here ]

5 Empirical Analysis

In this section, we investigate influential points in real data sets using the normalizing Bregman divergence. In all examples, convergence assessment of the Markov chains were based on visual inspection of trace plots, autocorrelation plots and the statistic since we ran two chains for each case. All results indicated that the chains reached stationarity relatively fast.

5.1 Binary Regression for Alpine Birds

A study about an endemic coastal alpine bird was conducted in Vancouver Island (Southwest coast of British Columbia, Canada) for more than a decade and the results were published in Jackson et al. (2015). The presence or absence of birds in a grid of space was registered over the years together with other environmental characteristics as covariates. The authors proposed an interpretation of data by a Random Forest model. Here we extend their model to a Bayesian framework and consider a binary logistic regression with other covariates.

For illustration, we selected the following covariates: elevation (1000 meters) and average temperature in summer months (in Celcius degrees) and the model also includes an intercept. We then ran a HMC with two chains, each one with 4,000 iterations where the first half was used as burn-in. This setup was used to fit models with probit and logit link functions. The normalizing Bregman divergences estimated for each observation are displayed in Figure 1. From this figure, it is hard to judge what is a high value of divergence, because there are more than one thousand observations in the sample. However, we can easily conclude that the model with logit link performs better because the highest values of logit are lower than the highest values of probit. This is to say that the most influential points in the logit model are not so influent as in the probit model.

[ Figure 1 around here ]

5.2 Spatial models for rainfall in South of Brazil

Here we illustrate a spatial regression approach to analyze the data on precipitation levels in Paraná State, Brazil. This data is freely available in the geoR package and was previously analyzed by for example Diggle and Ribeiro Jr (2002) and Gaetan and Guyon (2010). The data refers to average precipitation levels over 33 years of observation during the period May-June (dry-season) in 143 recording stations throughout the state. The original variable was summarized in 100 millimeters of precipitation per station. We changed it to 10,000 millimeters of rain per station, which seems more intuitive once the average local precipitation is around 1,000 millimeter per year and the observation time was larger than 10 years.

We the fitted three models for the average rainfall. These are the full SRM presented in Equation (2), the same model but without the squared components and , and the smallest one without squared components and neither the interaction term , which we refer to as full, middle and small models respectively. All the models were fitted from two HMC chains, each one with 20,000 iterations where the first half was used as burn-in.

We estimated the normalizing Bregman divergence for each recording station and compared the results in the same way as in the previous example. We conclude that the small model was the best one in the sense that it shows the smallest peaks. For example, the maximum values for each model were 0.072, 0.068 and 0.039 respectively, which are already pretty high relative to the expected .

We chose to display only the results for this one best model in Figure 2 from which we can see that the largest values of normalizing Bregman (largest circles) are scattered around the map, notwithstanding the rainy region is concentrated in the southwest.

[ Figure 2 around here ]

5.3 GARCH for Bitcoin exchange to US Dollar

The cryptocurrencies were born in the new millennium dawn as an alternative as governments and banks. As such, they changed the rules of financial market and they appreciated very fast, although high fluctuation and sharp falls are common. In particular, the Bitcoin is likely the most famous cryptocurrency and shows the largest volume of crypt transactions.

We illustrate the statistical analysis with one year of daily data on the log-returns of Bitcoin (BTC) exchange to U.S. Dollar (USD) from August 5, 2017 to August 5, 2018. This data was produced from the CoinDesk price page (see We then fitted GARCH(1,1) models with normal and Student errors for the log-return of BTC to USD exchange. We ran the HMC with two chains, each one with 4,000 iterations where the first half was used as burn-in. The estimation of main parameter of Normal model are: has zero mean and SD, is 0.15 (0.06) and is 0.07 (0.08), on the other hand the Student is with zero mean and SD too, is 0.11 (0.05) and is 0.06 (0.07).

We estimated the normalizing Bregman divergence for each day, the result could be see in Figure 3. Here is not so trivial to choose between the Normal or the Student model. Because there is no clear dominance of one or another, even though the Student presents the highest value of divergence, both form a mixed cloud of values very closed. However the highest points are quite sure very influent observations, because they represent more than 20 times the expected mean of 1/364. Consequently, it is not a surprise that the observed high divergences in January correspond to what the Consumer News and Business Channel (CNBC) called a Bitcoin nightmare. A time of new regulations in South Korea as well as a Facebook currency policy change, which implied a devaluation.

[ Figure 3 around here ]

6 Discussion

In this article we explored the possibilities of using the functional Bregman divergence as a useful generalization of the Kullback-Leiber divergence to identify influential observations in Bayesian models, for both dependent and independent data. Kullback-Leiber is easier to estimate, but overall it is difficult to infer if a point represents an influential point or not. So we propose to normalizing the Bregman divergence based on the order maintenance of the functional. It has two intuitive advantages: firstly, it belongs to range between zero and one which is easier to interpret, secondly we can evaluate its intensity according to the sample size. In particular, the normalizing Bregman divergence for the Kullback-Leiber case avoids the need for heavy computations.

As we saw in the simulation study, the expected average of a normalizing Bregman divergence to any observation without perturbation is approximately . Of course that number of influential points is a relevant issue to evaluate the value of a normalizing Bregman divergence. The simulation study embraced three different fields of statistic: GLM, spatial models and time series, with similar conclusions in all of them. Besides, the empirical analysis explored three scientific fields: Ecology, Climatology and Finance. Finally, in all cases the Hamiltonian Monte Carlo was an efficient and fast way to obtain samples from the posterior distribution of parameters


Ricardo Ehlers received support from São Paulo Research Foundation (FAPESP) - Brazil, under grant number 2016/21137-2.


  • Ardia and Hoogerheide (2010) Ardia, D. and Hoogerheide, L. F. (2010). Bayesian estimation of the GARCH(1,1) model with Student-t innovations. The R Journal, 2(2):41–47.
  • Bollerslev (1986) Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3):307–327.
  • Chen (1994) Chen, M.-H. (1994). Importance-Weighted Marginal Bayesian Posterior Density Estimation. Journal of the American Statistical Association, 89:818–824.
  • Cho et al. (2009) Cho, H., Ibrahim, J. G., Sinha, D., and Zhu, H. (2009). Bayesian case influence diagnostic for survival models. Biometrics, 65(1):116–124.
  • Chung et al. (2013) Chung, Y., Rabe-Hesketh, S., Dorie, V., Gelman, A., and Liu, J. (2013). A nondegenerate penalized likelihood estimator for variance parameters in multilevel models. Psychometrika, 78(4):685–709.
  • Deschamps (2006) Deschamps, P. J. (2006). A flexible prior distribution for Markov switching autoregressions with Student-t errors. Journal of Econometrics, 133(1):153–190.
  • Diggle and Ribeiro Jr (2002) Diggle, P. and Ribeiro Jr, P. (2002). Bayesian inference in Gaussian model-based geostatistics. Geographical and Environmental Modelling, 6(2):129–146.
  • Duane et al. (1987) Duane, A., Kennedy, A., Pendleton, B., and Roweth, D. (1987). Hybrid Monte Carlo. Physics Letters B, 195(2):216–222.
  • Eguchi and Kano (2001) Eguchi, S. and Kano, Y. (2001). Robustifying maximum likelihood estimation. Institute of Statistical Mathematics. Technical Report, Tokyo Japan.
  • Frigyik et al. (2008a) Frigyik, B., Srivastava, S., and Gupta, M. (2008a). Functional Bregman divergence and Bayesian estimation of distributions. IEEE Transactions on Information Theory, 54:5130–5139.
  • Frigyik et al. (2008b) Frigyik, B., Srivastava, S., and Gupta, M. (2008b). An introduction to functional derivatives. Technical report [online], Department of Electrical Engineering, University of Washington, Seattle, WA.
  • Gaetan and Guyon (2010) Gaetan, C. and Guyon, X. (2010). Spatial Statistics and Modeling. Springer-Verlag.
  • Gelman et al. (2008) Gelman, A., Jakulin, A., Pittau, M. G., and Su, Y. S. (2008). A weakly informative default prior distribution for logistic and other regression models. The Annals of Applied Statistics, 2(4):1360–1383.
  • Goh and Dey (2014) Goh, G. and Dey, D. K. (2014). Bayesian model diagnostics using functional Bregman divergence. Journal of Multivariate Analysis, 124:371–383.
  • Hoffman and Gelman (2014) Hoffman, M. D. and Gelman, A. (2014). The no-u-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15:1351–1381.
  • Hong-Xia et al. (2016) Hong-Xia, H., Jin-Guan, L., Hong-Xia, W., and Xing-Fang, H. (2016). Bayesian case influence analysis for GARCH models based on Kullback-Leibler divergence. Journal of the Korean Statistical Society, 45(4):595–609.
  • Jackson et al. (2015) Jackson, M. M., Gergel, S. E., and Martin, K. (2015). Effects of climate change on habitat availability and configuration for an endemic coastal alpine bird. Plos One, 10(11):e0142110.
  • Leimkuhler and Reich (2004) Leimkuhler, B. and Reich, S. (2004). Simulating Hamiltonian dynamics, volume 14 of Cambridge monographs on applied and computational mathematics. Cambridge University Press, Cambridge, 2:18.
  • McCulloch (1989) McCulloch, R. E. (1989). Local model influence. Journal of American Statistical Association, 84(406):473–478.
  • Neal (2011) Neal, R. (2011). MCMC using Hamiltonian dynamics. In Brooks, S., Gelman, A., Jones, G. L., and Meng, X.-L., editors, Handbook of Markov Chain Monte Carlo, pages 116–162. Chapman and Hall/CRC.
  • Nelder and Wedderburn (1972) Nelder, J. A. and Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society A, 135(3):370–394.
  • Nesterov (2009) Nesterov, Y. (2009). Primal-dual subgradient methods for convex problems. Mathematical Programming, 120(1):221–259.
  • R Core Team (2015) R Core Team (2015). R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Stan Development Team (2016) Stan Development Team (2016). Stan modeling language users guide and reference manual. Version 2.14.0.
Model 1
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
100 0.125 0.105 0.126 0.112 0.096 0.092
-0.084 0.113 -0.079 0.108 -0.060 0.095
0.040 0.088 0.045 0.089 0.095 0.083
300 0.031 0.024 0.031 0.026 0.030 0.025
-0.021 0.027 -0.026 0.027 -0.015 0.024
0.009 0.022 0.010 0.024 0.009 0.021
Model 2
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
100 -0.155 0.151 -0.171 0.172 -0.117 0.146
0.145 0.164 0.130 0.171 0.094 0.131
-0.055 0.106 -0.061 0.109 -0.037 -0.037
300 -0.056 0.038 -0.045 0.039 -0.026 0.036
0.051 0.040 0.033 0.041 0.022 0.036
-0.011 0.027 -0.014 0.030 -0.008 0.028
Table 1: Bias and square root of mean square error for parameter estimates in the logistic regression.
Model 1
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
50 -0.0026 0.0478 -0.0126 0.0465 -0.0024 0.0426
0.0024 0.0278 0.0013 0.0248 -0.0031 0.0254
0.0016 0.0272 0.0007 0.0242 -0.0063 0.0256
0.0035 0.0326 0.0014 0.0306 -0.0133 0.0319
0.0018 0.0178 0.0055 0.0160 0.0018 0.0157
-0.0038 0.0166 0.0020 0.0153 -0.0023 0.0159
0.0324 0.0137 0.0329 0.0135 0.0065 0.0117
200 -0.0022 0.0112 -0.0050 0.0105 -0.0006 0.0102
0.0021 0.0054 -0.0029 0.0056 0.0033 0.0054
-0.0056 0.0053 -0.0003 0.0054 0.0029 0.0051
0.0067 0.0055 -0.0034 0.0053 -0.0017 0.0058
-0.0008 0.0029 0.0000 0.0029 0.0001 0.0030
-0.0026 0.0031 0.0001 0.0029 0.0017 0.0029
0.0064 0.0026 0.0093 0.0027 0.0049 0.0025
Model 2
Prior 1 Prior 2 Prior 3
50 -0.0041 0.0474 0.0008 0.0459 0.0082 0.0488
-0.0031 0.0258 -0.0074 0.0265 -0.0088 0.0271
-0.0050 0.0283 -0.0028 0.0266 0.0056 0.0269
-0.0086 0.0313 -0.0092 0.0335 0.0010 0.0307
-0.0041 0.0169 0.0008 0.0167 -0.0070 0.0167
0.0018 0.0161 -0.0013 0.0173 -0.0012 0.0169
0.0344 0.0132 0.0266 0.0121 0.0099 0.0117
200 -0.0024 0.0101 0.0017 0.0108 0.0007 0.0105
0.0023 0.0053 -0.0019 0.0052 0.0055 0.0048
-0.0001 0.0054 0.0009 0.0054 0.0010 0.0053
0.0019 0.0058 0.0004 0.0057 -0.0012 0.0060
0.0000 0.0027 0.0018 0.0028 -0.0014 0.0030
0.0015 0.0027 -0.0001 0.0028 0.0007 0.0029
0.0053 0.0026 0.0084 0.0025 0.0037 0.0025
Table 2: Bias and square root of mean square error for parameter estimates in the spatial regression model.
Model 1
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
500 0.716 0.725 0.339 0.317 0.178 0.195
-0.111 0.013 -0.119 0.015 -0.121 0.015
-0.174 0.037 -0.095 0.015 -0.065 0.009
900 0.773 0.837 0.399 0.371 0.215 0.235
-0.131 0.018 -0.136 0.019 -0.136 0.019
-0.167 0.036 -0.095 0.016 -0.063 0.011
Model 2
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
500 0.031 0.021 -0.049 0.023 -0.056 0.022
-0.300 0.095 -0.319 0.106 -0.321 0.107
0.009 0.005 0.067 0.010 0.071 0.011
900 0.093 0.024 0.033 0.015 0.036 0.015
-0.297 0.091 -0.315 0.102 -0.315 0.102
-0.029 0.005 0.014 0.004 0.012 0.004
Table 3: Bias and square root of mean square error for parameter estimates in the GARCH(1,1) model with normal errors.
Model 1
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
500 0.133 0.136 -0.191 0.132 -0.262 0.147
-0.121 0.015 -0.128 0.017 -0.129 0.017
-0.174 0.036 -0.093 0.013 -0.078 0.010
900 0.192 0.152 -0.147 0.116 -0.224 0.126
-0.139 0.020 -0.143 0.021 -0.144 0.021
-0.171 0.035 -0.090 0.013 -0.072 0.009
Model 2
Prior 1 Prior 2 Prior 3
Parameter Bias SMSE Bias SMSE Bias SMSE
500 -0.218 0.060 -0.291 0.096 -0.276 0.087
-0.361 0.133 -0.383 0.149 -0.385 0.151
0.033 0.006 0.099 0.015 0.093 0.013
900 -0.154 0.032 -0.219 0.057 -0.214 0.054
-0.368 0.137 -0.378 0.145 -0.376 0.143
-0.009 0.003 0.040 0.006 0.036 0.005
Table 4: Bias and square root of mean square error for parameter estimates in the GARCH(1,1) model with Student errors.
Perturbation Obs Mean SD Perturbation Obs Mean SD
I 19 0.010 0.016 II 19 0.009 0.012
I 44 0.009 0.016 II 44 0.009 0.009
I 64 0.010 0.017 II 64 0.031 0.035
III 19 0.009 0.012 IV 19 0.028 0.025
III 44 0.028 0.027 IV 44 0.030 0.029
III 64 0.030 0.034 IV 64 0.031 0.030
Perturbation Obs Mean SD Perturbation Obs Mean SD
I 19 0.003 0.005 II 19 0.003 0.005
I 44 0.003 0.004 II 44 0.003 0.003
I 64 0.003 0.004 II 64 0.008 0.010
III 19 0.003 0.004 IV 19 0.009 0.010
III 44 0.008 0.010 IV 44 0.008 0.009
III 64 0.008 0.008 IV 64 0.008 0.010
Table 5: Case influence diagnostic in logistic regression by normalizing Bregman divergence.
Perturbation Obs Mean SD Perturbation Obs Mean SD
I 3 0.020 0.027 II 3 0.011 0.016
I 15 0.018 0.027 II 15 0.013 0.019
I 19 0.018 0.026 II 19 0.423 0.111
III 3 0.009 0.013 IV 3 0.211 0.074
III 15 0.283 0.089 IV 15 0.214 0.075
III 19 0.283 0.088 IV 19 0.210 0.078
Perturbation Obs Mean SD Perturbation Obs Mean SD
I 3 0.005 0.007 II 3 0.004 0.006
I 15 0.005 0.007 II 15 0.004 0.006
I 19 0.004 0.007 II 19 0.183 0.047
III 3 0.003 0.004 IV 3 0.133 0.035
III 15 0.154 0.040 IV 15 0.130 0.036
III 19 0.153 0.040 IV 19 0.133 0.038
Table 6: Case influence diagnostic for spatial regression model using normalizing Bregman divergence.
Perturbation Obs Mean SD Perturbation Obs Mean SD
I 19 0.009 0.014 II 19 0.007 0.010
I 44 0.009 0.014 II 44 0.007 0.011
I 64 0.009 0.014 II 64 0.247 0.072
III 19 0.005 0.006 IV 19 0.162 0.051
III 44 0.188 0.059 IV 44 0.156 0.048
III 64 0.192 0.060 IV 64 0.156 0.049
Perturbation Obs Mean SD Perturbation Obs Mean SD
I 119 0.001 0.002 II 119 0.001 0.002
I 344 0.002 0.004 II 344 0.001 0.003
I 464 0.001 0.003 II 464 0.063 0.029
III 119 0.001 0.002 IV 119 0.054 0.023
III 344 0.058 0.026 IV 344 0.054 0.025
III 464 0.059 0.028 IV 464 0.052 0.021
Table 7: Case influence diagnostic for GARCH using normalizing Bregman divergence.
Figure 1: Normalizing Bregman divergence for an endemic coastal alpine bird (1990-2013).
Figure 2: Normalizing Bregman divergence for average precipitation in Paraná, Brazil.
Figure 3: Normalizing Bregman divergence for Bitcoin exchange to US Dollar (August 5 2017 - August 5 2018).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description