Analytic posteriors for Pearson's correlation coefficient.

Pearson's correlation is one of the most common measures of linear dependence. Recently, Bernardo (11th International Workshop on Objective Bayes Methodology, 2015) introduced a flexible class of priors to study this measure in a Bayesian setting. For this large class of priors, we show that the (marginal) posterior for Pearson's correlation coefficient and all of the posterior moments are analytic. Our results are available in the open-source software package JASP.

Introduction

Pearson's product–moment correlation coefficient ρ is a measure of the linear dependency between two random variables. Its sampled version, commonly denoted by r, has been well studied by the founders of modern statistics such as Galton, Pearson, and Fisher. Based on geometrical insights, fisher (1915, 1921) was able to derive the exact sampling distribution of r and established that this sampling distribution converges to a normal distribution as the sample size increases. Fisher's study of the correlation has led to the discovery of variance‐stabilizing transformations, sufficiency (FISHER, 1920), and, arguably, the maximum likelihood estimator (FISHER, 1922; STIGLER, 2007). Similar efforts were made in Bayesian statistics, which focus on inferring the unknown ρ from the data that were actually observed. This type of analysis requires the statistician to (i) choose a prior on the parameters, thus, also on ρ, and to (ii) calculate the posterior. Here we derive analytic posteriors for ρ given a large class of priors that include the recommendations of JEFFREYS (1961), LINDLEY (1965), BAYARRI (1981), and, more recently, BERGER and SUN (2008) and BERGER et al. (2015). Jeffreys's work on the correlation coefficient can also be found in the second edition of his book (JEFFREYS, 1961), originally published in 1948; see ROBERT et al. (2009) for a modern re‐read of Jeffreys's work. An earlier attempt at a Bayesian analysis of the correlation coefficient can be found in JEFFREYS (1935). Before presenting the results, we first discuss some notations and recall the likelihood for the problem at hand.

Notation and result

Let (X 1,X 2) have a bivariate normal distribution with mean and covariance matrix where and are the population variances of X 1 and X 2, and where ρ is Pearson's correlation coefficient ρ measures the linear association between X 1 and X 2. In brief, the model is parametrized by the five unknowns θ=(μ 1,μ 2,σ 1,σ 2,ρ).

Bivariate normal data consisting of n pairs of observations can be sufficiently summarized as , where is the sample correlation coefficient, the sample mean, and the average sums of squares. The bivariate normal model implies that the observations y are functionally related to the parameters by the following likelihood function: For inference we use the following class of priors: where η denotes the hyperparameters, that is, η=(α,β,γ,δ). This class of priors is inspired by the one that José Bernardo (2015) used in his talk on reference priors for the bivariate normal distribution at the ‘11th International Workshop on Objective Bayes Methodology in honor of Susie Bayarri’. This class of priors contains certain recommended priors as special cases.

If we set α=1,β=γ=δ=0 in Equation (3), we retrieve the prior that Jeffreys recommended for both estimation and testing (JEFFREYS, 1961, pp. 174–179 and 289–292). This recommendation is not the prior derived from Jeffreys's rule based on the Fisher information (e.g., ly et al., 2017), as discussed in BERGER and SUN (2008). With α=1,β=γ=δ=0, thus, a uniform prior on ρ, Jeffreys showed that the marginal posterior for ρ is approximately proportional to h (n,r|ρ), where represents the ρ‐dependent part of the likelihood Equation (2) with θ 0=(μ 1,μ 2,σ 1,σ 2) integrated out. For n large enough, the function h a is a good approximation to the true reduced likelihood h given below.

We thank an anonymous reviewer for clarifying how Jeffreys derived this approximation.

If we set α=β=γ=δ=0 in Equation (3), we retrieve Lindley's reference prior for ρ. lindley (1965, pp. 214–221) established that the posterior of is asymptotically normal with mean and variance n −1, which relates the Bayesian method of inference for ρ to that of Fisher. In Lindley's (1965, p. 216) derivation, it is explicitly stated that the likelihood with θ 0 integrated out cannot be expressed in terms of elementary functions. In his analysis, Lindley approximates the true reduced likelihood h with the same h a that Jeffreys used before. BAYARRI (1981) furthermore showed that with the choice γ=δ=0, the marginalization paradox (DAWID et al., 1973) is avoided.

In their overview, BERGER and SUN (2008) showed that for certain a,b with α=b/2−1,β=0,γ=a−2, and δ=b−1, the priors in Equation (3) correspond to a subclass of the generalized Wishart distribution. Furthermore, a right‐Haar prior (e.g., sun and berger, 2007) is retrieved when we set α=β=0,γ=−1,δ=1 in Equation (3). This right‐Haar prior then has a posterior that can be constructed through simulations, that is, by simulating from a standard normal distribution and two chi‐squared distributions (berger and sun, 2008, Table 1). This constructive posterior also corresponds to the fiducial distribution for ρ (e.g., fraser, 1961; hannig et al., 2006). Another interesting case is given by α=0,β=1,γ=δ=0, which corresponds to the one‐at‐a‐time reference prior for σ 1 and σ 2; see also Jeffreys (1961 p. 187).

The analytic posteriors for ρ follow directly from exact knowledge of the reduced likelihood h (n,r|ρ), rather than its approximation used in previous work. We give full details, because we did not encounter this derivation in earlier work.

The reduced likelihood hγ,δ(n,r|ρ) If |r|<1,n>γ+1 , and n>δ+1 , then the likelihood f(y|θ) times the prior Equation  (3) with θ 0=(μ 1,μ 2,σ 1,σ 2) integrated out is a function f that factors as The first factor is the marginal likelihood with ρ fixed at zero, which does not depend on r nor on ρ , that is, where . We refer to the second factor as the reduced likelihood, a function of ρ which is given by a sum of an even function and an odd function, that is, h =A +B , where where and where 2 F 1 denotes Gauss' hypergeometric function.

□

To derive f (y|ρ), we have to perform three integrals: (i) with respect to , (ii) , and (iii) .

The integral with respect to yields where we abbreviated . The factor p (y 0) follows directly by setting ρ to zero in Equation (8) and two independent gamma integrals with respect to σ 1 and σ 2 resulting in Equation (5). These gamma integrals cannot be used when ρ is not zero. For f (y|ρ), which is a function of ρ, we use results from special functions theory.

For the second integral, we collect only that part of Equation (8) that involves σ 1 into a function g, that is, The assumption n>γ+1 and the substitution allow us to solve this integral using Lemma A.1, which we distilled from the Bateman manuscript project (ERDéLYI et al., 1954), with and c=n−γ−1. This yields where and where 1 F 1 denotes the confluent hypergeometric function. The functions Å and B¨ are the even and odd solutions of Weber's differential equation in the variable , respectively.

With , we see that f (y|ρ) follows from integrating σ 2 out of the following expression: where Hence, the last integral with respect to σ 2 only involves the functions k and l in Equation (12). The assumption n>δ+1 and the substitution , thus, allow us to solve this integral using Equation (7.621.4) from gradshteyn and ryzhik (2007, p. 822) with . This yields After we combine the results, we see that , where Hence, f (y|ρ) is of the asserted form. Note that is even, while is an odd function of ρ.

This main theorem confirms Lindley's insights; h (n,r|ρ) is indeed not expressible in terms of elementary functions, and the prior on ρ is updated by the data only through its sampled version r and the sample size n. As a result, the marginal likelihood for data y then factors into p (y)=p (y 0)p (n,r;γ,δ), where is the normalizing constant of the marginal posterior of ρ. More importantly, the fact that the reduced likelihood is the sum of an even function and an odd function allows us to fully characterize the posterior distribution of ρ for the priors Equation (3) in terms of its moments. These moments are easily computed, as the prior π (ρ) itself is symmetric around zero. Furthermore, the prior π (ρ) can be normalized as where denotes the beta function. The case with β=0 is also known as the (symmetric) stretched beta distribution on (−1,1) and leads to Lindley's reference prior when we ignore the normalization constant, that is, , and, subsequently, let .

Characterization of the marginal posteriors of ρ If n>γ+δ−2α+1 , then the main theorem implies that the marginal likelihood with all the parameters integrated out factors as p (y)=p (y 0)p (n,r;γ,δ) where defines the normalizing constant of the marginal posterior for ρ . Observe that the integral involving B is zero, because B is odd on (−1,1) . More generally, the k th posterior moment of ρ is These posterior moments define the series where is the normalization constant of the prior Equation  (13) , W (n) is the ratios of gamma functions as defined under Equation  (7) , and refers to the Pochhammer symbol for rising factorials. The terms a and b are The series defined in Equation  (16) are hypergeometric when β is a non‐negative integer.

□

The series E(ρ |n,r) result from term‐wise integration of the hypergeometric functions in A and B . The assumption n>γ+δ−2α+1 and the substitution x=ρ 2 allow us to solve these integrals using Equation (3.197.8) in gradshteyn and ryzhik (2007, p. 317) with their and when k is even, while we use when k is odd. A direct application of the ratio test shows that the series converge when |r|<1.

Analytic posteriors for the case β=0

For most of the priors discussed earlier, we have β=0, which leads to the following simplification of the posterior.

(Characterization of the marginal posteriors of ρ, when β=0) If n>γ+δ−2α+1 and |r|<1 , then the marginal posterior for ρ is where p (n,r;γ,δ) refers to the normalizing constant of the (marginal) posterior of ρ , which is given by More generally, when β=0 , the k th posterior moment is when k is even, and when k is odd.

□

The assumption n>γ+δ−2α+1 and the substitution x=ρ 2 allow us to use Equation (7.513.12) in gradshteyn and ryzhik (2007, p. 814) with and when k is even, while we use when k is odd. The normalizing constant of the posterior p (n,r;γ,δ) is a special case with k=0.

The marginal posterior for ρ updated from the generalized Wishart prior, the right‐Haar prior, and Jeffreys's recommendation then follow from a direct substitution of the values for α,γ, and δ as discussed under Equation (3). Lindley's reference posterior for ρ is given by which follows from Equation (17) by setting γ=δ=0 and, subsequently, letting .

Lastly, for those who wish to sample from the posterior distribution, we suggest the use of an independence‐chain Metropolis algorithm (tierney, 1994) using Lindley's normal approximation of the posterior of as the proposal. This method could be used when Pearson's correlation is embedded within a hierarchical model, as the posterior for ρ will then be a full conditional distribution within a Gibbs sampler. For α=1,β=γ=δ=0,n=10 observations and r=0.6, the acceptance rate of the independence‐chain Metropolis algorithm was already well above 75%, suggesting a fast convergence of the Markov chain. For n larger, the acceptance rate further increases. The R code for the independence‐chain Metropolis algorithm can be found on the first author's home page. In addition, this analysis is also implemented in the open‐source software package jasp (https://jasp-stats.org/).