Extending the Behrens-Fisher Problem to Testing Equality of Slopes in Linear Regression: The Bayesian Approach

Mohamed Shoukri^*, Futwan Al-Mohanna

Department of Cell Biology, National Biotechnology Center, Research Centre, King Faisal Specialist Hospital and Research Center, Riyadh, Saudi Arabia.

DOI: 10.4236/ojs.2018.82018                

Abstract

Testing the equality of means of two normally distributed random variables when their variances are unequal is known in the statistical literature as the “Behrens-Fisher problem”. It is well-known that the posterior distributions of the parameters of interest are the primitive of Bayesian statistical inference. For routine implementation of statistical procedures based on posterior distributions, simple and efficient approaches are required. Since the computation of the exact posterior distribution of the Behrens-Fisher problem is obtained using numerical integration, several approximations are discussed and compared. Tests and Bayesian Highest-Posterior Density (H.P.D) intervals based upon these approximations are discussed. We extend the proposed approximations to test of parallelism in simple linear regression models.

Keywords

Bayesian Inference, Posterior Distributions, Behrens-Fisher Problem, Posterior Moments, Edgeworth Expansion, Monte-Carlo Integration

1. Introduction

Suppose that x₁ and x₂ are two independent normal random variables with means μ₁ and μ₂, and variances ${\sigma }_1^2$ and ${\sigma }_2^2$ , respectively. Samples of sizes n₁ and n₂ drawn from the corresponding populations are denoted by $x_{ij}$ ( $i=1,2$ and $j=1,2,\cdots ,n_i$ ). It is desired to test the hypothesis H₀: μ₁ = μ₂ when it cannot be taken as known that the variance ratio $\lambda ={\sigma }_1^2/{\sigma }_2^2$ is one. This problem, known as the Behrens-Fisher (BF) problem, has been investigated by many authors and each has proposed some particular solution. It is not the purpose of this paper to list or survey these solutions. However, the Bayesian approach to this problem, viewed as one of the most fascinating approaches of statistical inference on the means of heterogeneous normal populations, will be the main focus of this paper, and we shall also give special attention to the problem of testing parallelism of two linear regression lines when the variances of the error terms are not equal.

The exact Bayesian solution to the BF problem was given by Jeffreys [1] which was identical to the solution of Behrens [2] and Fisher [3] . The resulting posterior distribution of $U={\mu }_1-{\mu }_2$ , which is the primitive of a valid Bayesian procedure, possesses a form requiring an integration that cannot be performed analytically, and hence, direct and routine implementation of this test is impossible. Although evaluation of such an integral can be achieved numerically, using the approaches described by Reilly [4] and Naylor and Smith [5] , and by Monte-Carlo integration as given in Robert and Casella [6] it is of interest to provide expressions for routine applications under the Bayesian approach.

The paper has two chief objectives. First we present a comparison among several approximations to the tails of the posterior distribution of the variate $U={\mu }_1-{\mu }_2$ , where μ_j is the mean of the j^th population. Second we extend the methodologies to address the question of parallelism or equality of slopes when the variances of the error terms of the two regression lines are not equal. Recommendations will follow the examples in the last sections.

2. Testing Equality of Normal Means: Posterior Analysis

Based on the samples outcome, the usual sufficient statistics are defined as follows:

${\bar{x}}_i={\displaystyle {\sum }_{j=1}^{n_i}{x}_{ij}}/n_i$ and ${\Sigma }_i=\left(n_i-1\right)s_i^2={\displaystyle {\sum }_{j=1}^{n_i}{\left(x_{ij}-{\bar{x}}_i\right)}^2}$ , $i=1,2$ .

The likelihood function of the combined data from two samples drawn independently from two normal populations with means μ₁ and μ₂, and variances ${\sigma }_1^2$ and ${\sigma }_2^2$ , respectively, is proportional to:

${\left(1/{\sigma }_1\right)}^{n_1}{\left(1/{\sigma }_2\right)}^{n_2}\exp \left\{-\frac{1}{2}\left[\frac{{\Sigma }_1+n_1{\left({\bar{x}}_1-{\mu }_1\right)}^2}{{\sigma }_1^2}+\frac{{\Sigma }_2+n_2{\left({\bar{x}}_2-{\mu }_2\right)}^2}{{\sigma }_2^2}\right]\right\}$ . (2.1)

Following Lee [7] , we use the non-informative prior for the means and the standard deviations

$\Pi \left({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2\right)\text{d}{\mu }_1\text{d}{\mu }_2\text{d}{\sigma }_1\text{d}{\sigma }_2\propto \text{d}{\mu }_1\text{d}{\mu }_2\frac{\text{d}{\sigma }_1\text{d}{\sigma }_2}{{\sigma }_1{\sigma }_2}$ (2.2)

From Box & Tiao [8] , the joint posterior density of $\left({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2\right)$ is the product of (2.1) and (2.2) and is then given by:

$\begin{array}{l}\Pi \left({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2|x\right)\\ \propto {\sigma }_1^{-\left(n_1+1\right)}{\sigma }_2^{-\left(n_2+1\right)}\exp \left\{-\frac{1}{2}\left[\frac{{\Sigma }_1+n_1{\left({\bar{x}}_1-{\mu }_1\right)}^2}{{\sigma }_1^2}+\frac{{\Sigma }_2+n_2{\left({\bar{x}}_2-{\mu }_2\right)}^2}{{\sigma }_2^2}\right]\right\}\end{array}$ (2.3)

Integrating ${\sigma }_1,{\sigma }_2$ out of (2.3), the marginal posterior density of $U={\mu }_1-{\mu }_2$ is shown to be

$\begin{array}{c}\Pi \left(u|x\right)\propto {\displaystyle {\int }_{-\infty }^{\infty }{\left[1+\frac{n_1{\left(y+u-{\bar{x}}_1\right)}^2}{{\Sigma }_1}\right]}^{-\frac{n_1}{2}}{\left[1+\frac{n_2{\left(y-{\bar{x}}_2\right)}^2}{{\Sigma }_2}\right]}^{-\frac{n_2}{2}}\text{d}y}\\ ={\displaystyle {\int }_{-\infty }^{\infty }\Pi \left(u|y,x\right)\cdot \Pi \left(y|x\right)\text{d}y}\end{array}$ (2.4)

The conditioning in Equation (2.4) is on the represents the data vector x.

Equation (2.4) was obtained by Jeffreys [1] . In the Bayesian approach, inferences about U are completely determined by (2.4), which is not amenable to simple manipulation in order to have the tests on U conducted in a routine fashion. Before applying some of the suggested approximations, we should understand the nature of the problem, and in order to draw safe conclusions, we consider not only the marginal posterior density of u, but also we need to examine the other components of (2.3). For this purpose, we state the following results without proof, since they are easily obtained by application of the calculus of probability:

1) The conditional posterior pdf for U, given μ₂, is the univariate student’s $t_{\left(n_1-1\right)}$ , with mean ${\bar{x}}_1-{\mu }_2$ and variance ${\text{Σ}}_1/\left(n_1\left(n_1-3\right)\right)$ , i.e.,

$\Pi \left(u|{\mu }_2,x\right)\propto {\left[1+\frac{n_1{\left(u-\left({\bar{x}}_1-{\mu }_2\right)\right)}^2}{{\Sigma }_1}\right]}^{-\frac{n_1}{2}}$ . (2.5)

2) The marginal posterior pdf for the variance ratio λ is such that $\left(s_2^2/s_1^2\right)\lambda$ has the well-known Snedecor’s F-distribution with $\left(n_2-1\right)$ , and $\left(n_1-1\right)$ degrees of freedom, i.e.,

$\Pi \left(\lambda |x\right)\propto {\lambda }^{\frac{n_2-3}{2}}{\left(1+\frac{{\Sigma }_2}{{\Sigma }_1}\lambda \right)}^{-\left(\frac{n_1+n_2-1}{2}\right)}$ (2.6)

3) The conditional posterior pdf for U, given λ, is such that

$t=\Pi \left(u|\lambda ,x\right)\propto {\left[1+a\left(\lambda \right){\left(u-{\bar{x}}_1+{\bar{x}}_2\right)}^2\right]}^{-\left(\frac{n_1+n_2-2}{2}\right)}$ (2.7)

where $\left(\lambda \right)=\frac{n_1{n}_2\lambda }{\left(n_1+n_2\lambda \right)\left({\text{Σ}}_1+\lambda {\text{Σ}}_2\right)}$ , i.e., it has student’s t density with

$\left(n_1+n_2-2\right)$ degrees of freedom. Hence as was indicated by Barnard [9] , if a range of λ values is considered plausible, λ can be varied over this range to see what effect this variation has on the value of t. He also postulated that when n₁ and n₂ are nearly equal, and moderately large, and $s_1^2$ and $s_2^2$ do not differ much, the conditional distribution in (2.7) will remain nearly constant over a wide range of λ values. In the next section we derive different approximations to deal with the situation when the range of λ is wide enough to have a considerable effect on t.

3. Approximate Posterior Inference

3.1. Monte Carlo Integration (Exact Approximation)

We shall write the posterior density of U as follows:

$\Pi \left(u|x\right)={\displaystyle {\int }_0^{\infty }\Pi \left(u|\lambda ,x\right)\cdot \Pi \left(\lambda |x\right)\text{d}\lambda }$ (3.1a)

In (3.1a) we face the same problem as we have with (2.4), however, this equation will be the tool of discussion in the remainder of this section. As can be seen evaluation of the posterior density of U is just evaluating the integral in (3.1a) which can be written as:

$E\left[\Pi \left(u|\lambda ,x\right)\right]={\displaystyle {\int }_0^{\infty }\Pi \left(u|\lambda ,x\right)\cdot \Pi \left(\lambda |x\right)\text{d}\lambda }$ (3.1b)

Referring to Robert and Casella [6] , the principle of the Monte Carlo method for approximating (3.1b) is to generate sample $\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\right)$ from their density $\text{Π}\left(\lambda |x\right)$ and suggest as an approximation the empirical average

$\Pi {\left(u|\lambda ,x\right)}^{*}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\Pi \left(u|{\lambda }_i,x\right)}$ (3.2)

3.2. Moments (Scale) Approximation

The second approximation to be considered is by the method of moments. While Patil [10] derived the posterior moments of U using (2.4), they are identical to the posterior moments

$E\left(u^r|x\right)={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_0^{\infty }{u}^r\cdot \Pi \left(u|\lambda ,x\right)\cdot \Pi \left(\lambda |x\right)\text{d}\lambda \text{d}u}}$ . (3.4)

Denoting the r^th central moment of U by $m_r$ , and its 4^th cumulant by $l_4$ , we can show from Equation (3.4) that

$m_2=\frac{{\Sigma }_1}{n_1\left(n_1-3\right)}+\frac{{\Sigma }_2}{n_2\left(n_2-3\right)}$ ,

$m_4=3\left[\frac{{\Sigma }_1^2}{n_1^2\left(n_1-3\right)\left(n_1-5\right)}+\frac{2{\Sigma }_1{\Sigma }_2}{n_1{n}_2\left(n_1-3\right)\left(n_2-3\right)}+\frac{{\Sigma }_2^2}{n_2^2\left(n_2-3\right)\left(n_2-5\right)}\right]$

The fourth cumulant $l_4$ is:

$l_4=m_4-3m_2^2$ (3.5)

Following Patil [10] , we suggest the following approximation to the posterior distribution of δ:

$\delta ~a{t}_{\left(b\right)}\text{where}\left(\delta =U-{\bar{x}}_1+{\bar{x}}_2\right)$ . (3.6)

Equating m₂ and m₄ to the second and fourth central moments of the R.H.S. of ~ in (3.6) we have:

$a=m_2{m}_4{\left(m_4+l_4\right)}^{-\frac{1}{2}}$ and $b=\left[2\left(1+\frac{m_4}{l_4}\right)\right]$ .

In the above equation [s] means the smallest integer larger than s. Thus one may use tables of student’s t-distribution with [b] degrees of freedom to make tests and construct approximate H.P.D. intervals for U. In other words:

$1-\alpha =P_r\left(t_{\alpha /2}<T_b<t_{1-\alpha /2}\right)=P_r\left(a{t}_{b,\frac{\alpha }{2}}<\delta <a{t}_{b,1-\frac{\alpha }{2}}\right)$

3.3. Modal Approximation

Barnard [9] suggested that if the intention is to make inference on U alone by integrating out λ and it is found that $\Pi \left(u|\lambda ,x\right)$ is very sensitive to changes in λ, then $\Pi \left(\lambda |x\right)$ should be carefully examined. If $\Pi \left(\lambda |x\right)$ were sharp with most of its probability mass concentrated over a small region about its mode denoted by λ^*, performing the numerical integration in (3.1a) will be nearly equivalent to assigning the model value to λ in the conditional density, $\Pi \left(u|\lambda ,x\right)$ so that:

$\Pi \left(u|\lambda ,x\right)\doteq \Pi \left(u|{\lambda }^{*},x\right)$ . (3.7)

Now, since ${\lambda }^{*}=\frac{\left(n_2-3\right){\Sigma }_1}{\left(n_1+1\right){\Sigma }_2}$ , then substituting in (3.7) we have:

$\Pi \left(u|x\right)\propto {\left[1+a\left({\lambda }^{*}\right){\left(u-{\bar{x}}_1+{\bar{x}}_2\right)}^2\right]}^{-\left(\frac{n_1+n_2-1}{2}\right)}$ , (3.8)

where

$a\left({\lambda }^{*}\right)={\left[\left(n_1+n_2-2\right)\left(\frac{{\Sigma }_1}{\left(n_1+1\right)n_1}+\frac{{\Sigma }_2}{n_2\left(n_2-3\right)}\right)\right]}^{-1}$ .

Therefore, as a modal approximation to the posterior distribution of we take:

$t^{*}=\frac{u-\left({\bar{x}}_1-{\bar{x}}_2\right)}{\sqrt{\frac{{\Sigma }_1}{n_1\left(n_1+1\right)}+\frac{{\Sigma }_2}{n_2\left(n_2-3\right)}}}~t_{\left(n_1+n_2-2\right)}$ . (3.9)

3.4. Approximation Based on Averaging

The following lemma due to Feller is quite appealing and may be used to approximate the integral (3.1a):

Lemma: Feller ( [11] , p. 219). For $n=1,2,\cdots$ , consider a family of distributions ${\varphi }_n\left(y\right)$ with expectation Θ and variance ${\sigma }_n^2\left(\Theta \right)$ (that is the variance is a function of the mean). If $g(\cdot )$ is bounded and continuous, and ${\sigma }_n^2\left(\Theta \right)\to 0$ for each Θ, then

$E\left(g(\cdot )\right)={\displaystyle {\int }_{-\infty }^{\infty }g\left(y\right){\varphi }_n\left(y\right)\text{d}y}\to g\left(\Theta \right)$ . (3.10)

Since $\bar{\lambda }=E\left(\lambda |x\right)\equiv \Theta =\frac{\left(n_1-1\right)s_1^2}{\left(n_1-3\right)s_2^2}$ and $\mathrm{var}\left(\lambda |x\right)={\sigma }_n^2\left(\theta \right)=2{\theta }^2\left[\frac{n_1^{-1}+n_2^{-1}-4{\left(n_1{n}_2\right)}^{-1}}{\left(n_1^{-1}-{\left(n_1{n}_2\right)}^{-1}\right)\left(n_2^{-1}-5{\left(n_1{n}_2\right)}^{-1}\right)}\right]$ with ${\sigma }_n^2\left(\Theta \right)\to 0$ as $\left(n_1,n_2\right)\to \infty$ , then by taking $\Pi \left(u|\lambda ,x\right)$ as $g(\cdot )$ and $\Pi \left(\lambda |x\right)$ as ${\varphi }_n\left(y\right)$ the right hand integral of (3.1a) can be approximated by:

$\Pi \left(u|x\right)\propto {\left[1+a\left(\bar{\lambda }\right){\left(u-{\bar{x}}_1+{\bar{x}}_2\right)}^2\right]}^{-\left(\frac{n_1+n_2-1}{2}\right)}$ ,

where

$a\left(\bar{\lambda }\right)={\left\{\left(n_1+n_2-4\right)\left[\frac{s_2^2}{n_2}+\frac{\left(n_1-1\right)s_1^2}{n_1\left(n_1-3\right)}\right]\right\}}^{-1}$ .

Accordingly, one may take

$E=\frac{U-\left({\bar{x}}_1-{\bar{x}}_2\right)}{\sqrt{\left(\frac{n_1+n_2-4}{n_1+n_2-2}\right)\left[\frac{{\Sigma }_1}{n_1\left(n_1-3\right)}+\frac{{\Sigma }_2}{n_2\left(n_2-1\right)}\right]}}~t_{\left(n_1+n_2-2\right)}$ . (3.11)

as an approximating distribution to the posterior distribution of U.

3.5. Edgeworth Expansion

The Edgeworth expansion has been used extensively by many authors in order to approximate the density function of any statistics ${\nu }_n$ . We refer to the paper by Barndorff-Nielson and Cox [12] for an overview and comparison between these techniques. The Edgeworth expansion, due to Edgeworth [13] [14] , for the

density function of any statistic ${\nu }_n=\frac{v_n-E\left(v_n\right)}{\sqrt{\text{Var}\left(v_n\right)}}$ at a point c is given by the general formula:

$\begin{array}{c}\psi \left(y\right)\doteq [1+\frac{\sqrt{{\vartheta }_1}}{6}\left(c^3-3c\right)+\frac{{\theta }_2-3}{24}\left(c^4-6c^2+3\right)\\ +\frac{{\theta }_1}{72}\left(c^6-15c^4+45c^2-15\right)]\cdot \frac{1}{\sqrt{2}\Pi }\exp \left(-\frac{1}{2}{c}^2\right),\end{array}$

where ${\theta }_1$ and $\left({\theta }_2-3\right)$ are respectively the coefficients of skewness and Kurtosis for the density of ${\nu }_n$ _. As can be seen when ${\vartheta }_1=0$ , the terms of order $n^{-1/2}$ disappear and the modified Edgeworth expansion of order $n^{-1}$ for the density function of the statistic ${\nu }_n=\frac{U-{\bar{x}}_1+{\bar{x}}_2}{\sqrt{{\tau }_2}}$ at c is given by:

$\psi \left(y\right)\doteq \left[1+\frac{{\kappa }_4}{24{\tau }_2^2}\left(c^4-6c^2+3\right)\right]\cdot \frac{1}{\sqrt{2\Pi }}\exp -\left(\frac{1}{2}{c}^2\right)$ . (3.12)

The Edgeworth expansion for the density of U can easily be obtained, from Equation (3.12), by using linear transformation, $U={\bar{x}}_1-{\bar{x}}_2+\sqrt{m_2}c$ .

Although we are not approximating the distribution function directly, in practice these approximations may be used for calculating tail areas. Thus it is of interest to see how these approximations for the posterior distributions of $U={\mu }_1-{\mu }_2$ perform at the tail.

Example 1: “Mean tumor recurrence scores in breast cancer patients stratified by tumor grades.”

Oncotype DX is a commercial assay used for making decisions regarding the treatment of breast cancer. The results are reported as a tumor recurrence score ranging from 0 to 100, Klein et al. [15] showed that the recurrence score correlated (among other genetic factors) with the tumor grade. That is there was a significant difference in the mean recurrence scores in patients with tumor grade 1 as compared to the mean recurrence score of patients with tumor grade 3. In this example we present the summary data for 40 and 37 breast cancer women samples with mean recurrence scores in tumor rage 1 and tumor grade 3. This data is a good example that we can use to apply the proposed methodology. The following results were obtained:

$n_1=40$ , ${\bar{x}}_1=11.55$ , $s_1^2=18.3$

$n_2=37$ , ${\bar{x}}_2=34.57$ , $s_2^2=171.25$

Approximations, discussed in section 3, to the posterior density of $U={\mu }_1-{\mu }_2$ are applied to the above data. In addition we also consider the usual asymptotic normal approximation which for large n₁ and n₂, is given as suggested by Welch [16] as:

$z=\frac{U-\left({\bar{x}}_1-{\bar{x}}_2\right)}{\sqrt{{\tau }_2}}\stackrel{p}{\to }N\left(0,1\right)$ ,

The results of the data analyses based on the proposed approximations are presents in Table 1.

As we can see, all approximations have confidence limits close to the exact limits, probably because the sample sizes are moderately large. We provide the R-codes for the calculations of these limits in the Appendix within the applications of linear regression.

4. Extension of the Behrens-Fisher Problem to Testing Equality of Slopes of Two Independent Regression Lines

Linear and nonlinear regression models are ubiquitous in medical and biological research. Testing equality of slopes of two linear regression lines is of special interest. This is illustrated in the following application which is a continuation of example 1.

Table 1. Approximate 95% HPD for the difference in two means under heterogeneity using the breaking strength data.

Example 2: continuation of example 1.

Ki67 is a commonly used marker of cancer cell proliferation, and has significant prognostic value in tumor recurrence of breast cancer. In this illustration which is a continuation to example 1, we use a sub-sample of women who had Oncotype DX testing performed and available Ki67 indices which correlated with tumor grades (1 versus 3). Literature documented that Ki67 scores contribute significantly to models that predict risk of recurrence in breast cancer, for example see; Cuzick et al. [17] , Klein et al. [14] , and more recently Thakur et al. [18] . In this example we examine the relationship between Ki67 as predictor of tumor recurrence scores of breast for tumor grade 3 and grade 1 separately. Of interest is to test the equality of the slopes of the linear regressions of tumor recurrence score on the Ki67 for both grades. We use Bayesian methodology to answer this question.

We shall use the following notations. First we denote the independent variable by x to predict values of the dependent variable denoted by y.

Suppose that we have two linear regression lines, then conditional on $x_{ij}$ we assume that

$y_{ij}~N\left({\alpha }_j+{\beta }_j\left(x_{ij}-{\bar{x}}_j\right),{\varphi }_j\right)$ , $i=1,2,\cdots ,n_j$ & $j=1,2$ .

In general we assume that we have two conditions, and we would like to estimate $\delta ={\beta }_1-{\beta }_2$ .

In our Bayesian analysis we shall take a reference prior that is independently uniform in ${\alpha }_j,{\beta }_j$ and $\log {\varphi }_j$ , such that:

$\pi \left({\alpha }_j,{\beta }_j,{\varphi }_j\right)\propto 1/{\varphi }_j,j=1,2$

Let us define the following quantities:

$S_{jee}=S_{jyy}-S_{jxy}^2/S_{jxx}$ , $a_j={\bar{y}}_j$ , $b_j=S_{jxy}/S_{jxx}$ , $e_{j0}={\bar{y}}_j-b_j{\bar{x}}_j$

where

$S_{jxy}={\displaystyle {\sum }_{i=1}^{n_j}\left(x_{ij}-{\bar{x}}_j\right)\left(y_{ij}-{\bar{y}}_j\right)}$ , $S_{jyy}={\displaystyle {\sum }_{i=1}^{n_j}{\left(y_{ij}-{\bar{y}}_j\right)}^2}$ , $S_{jxx}={\displaystyle {\sum }_{i=1}^{n_j}{\left(x_{ij}-{\bar{x}}_j\right)}^2}$ ,

and

${\Sigma }_j=\left(n_j-2\right)S_{jee}/S_{jxx},j=1,2$ .

The joint posterior distribution of the model parameters $\left({\alpha }_j,{\beta }_j,{\varphi }_j\right)$ are proportional to

$\begin{array}{l}f\left({\alpha }_j,{\beta }_j,{\varphi }_j\text{|}{x}_{ij},y_{ij}\right)\\ \propto {\pi }_j={\varphi }_j^{-\left(n_j+2\right)/2}\exp \left[-\frac{1}{2}\left\{S_{jee}+n_j{\left({\alpha }_j-a_j\right)}^2+S_{jxx}{\left({\beta }_j-b_j\right)}^2\right\}/{\varphi }_j\right]\end{array}$ (4.1)

For the two regression lines, the joint posterior of $\left({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\varphi }_1,{\varphi }_2\right)$ is ${\pi }_1{\pi }_2$ , or

$\begin{array}{l}\pi \left({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\varphi }_1,{\varphi }_2|x,y\right)\\ \propto {\varphi }_1^{-\left(\frac{n_1+2}{2}\right)}{\varphi }_2^{-\left(\frac{n_2+2}{2}\right)}\underset{j=1}{\overset{2}{{\displaystyle \prod }}}\exp \left[-\frac{1}{2}\left\{S_{jee}+n_j{\left({\alpha }_j-a_j\right)}^2+S_{jxx}{\left({\beta }_j-b_j\right)}^2\right\}/{\varphi }_j\right]\end{array}$ (4.2)

Integrating out ${\alpha }_1$ and ${\alpha }_2$ , the joint posterior of $\left({\beta }_1,{\beta }_2,{\varphi }_1,{\varphi }_2\right)$ is thus given as:

$\propto {\varphi }_1^{-\left(\frac{n_1+1}{2}\right)}{\varphi }_2^{-\left(\frac{n_2+1}{2}\right)}\underset{j=1}{\overset{2}{{\displaystyle \prod }}}\exp \left[-\frac{1}{2}\left\{\frac{S_{jee}+S_{jxx}{\left({\beta }_j-b_j\right)}^2}{{\varphi }_j}\right\}\right]$

Integrating out ${\varphi }_1$ and ${\varphi }_2$ we get:

$\pi \left({\beta }_1,{\beta }_2|x,y\right)\propto {\left[1+\frac{n_1-2}{{\Sigma }_1}{\left({\beta }_1-b_1\right)}^2\right]}^{-\left(\frac{n_1-1}{2}\right)}{\left[1+\frac{n_2-2}{{\Sigma }_2}{\left({\beta }_2-b_2\right)}^2\right]}^{-\left(\frac{n_2-1}{2}\right)}$ (4.3)

Under the transformation $\delta ={\beta }_1-{\beta }_2$ , one can show that the posterior density of δ is given by:

$\begin{array}{c}\pi \left(\delta |x,y\right)\propto \underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}{\left[1+\frac{r_1}{{\Sigma }_1}{\left(\delta +{\beta }_2-b_1\right)}^2\right]}^{-\left(\frac{r_1+1}{2}\right)}{\left[1+\frac{r_2}{{\Sigma }_2}{\left({\beta }_2-b_2\right)}^2\right]}^{-\left(\frac{r_2+1}{2}\right)}\text{d}{\beta }_2\\ =\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\pi \left(\delta |{\beta }_2,x,y\right),\pi \left({\beta }_2|x,y\right)\text{d}{\beta }_2\end{array}$ (4.4)

where $r_j=n_j-2$ .

The Bayesian inferences on δ are completely determined by Equation (4.4) which we cannot easily manipulate in order to have statistical inferences test on δ conducted in a routine fashion. Moreover, it is clear from (4.4) that the exact

marginal posterior distribution of $\frac{\left({\beta }_j-b_j\right)\sqrt{r_j}}{{\Sigma }_j}$ is that of a Student t-distribution with r_j degrees of freedom.

Similar to testing the equality of means of two normally distributed distributions and as shown in the first part of the paper, we use the suggested approximations to the integral given in (4.4) to find the marginal posterior distribution of δ. However, it is quite helpful to not only examine the posterior density of δ, but also examine the components of the joint posterior density given in (4.4). For this purpose, we state the following results without proof, since they can be easily obtained by applications of the calculus of probability.

1) The conditional posterior p∙d∙f of δ, given ${\beta }_2$ is the univariate Student-t with ( $n_1-2$ ), with conditional mean $b_1-{\beta }_2$ and conditional variance $\frac{{\Sigma }_1}{\left(n_1-2\right)\left(n_1-4\right)}$ .

The unconditional posterior mean and variance of δ are:

$E\left(\delta |x,y\right)=b_1-b_2$

${\mu }_2\left(\delta \right)=\mathrm{var}\left(\delta |x,y\right)=\frac{{\Sigma }_1}{\left(n_1-2\right)\left(n_1-4\right)}+\frac{{\Sigma }_2}{\left(n_2-2\right)\left(n_2-4\right)}$

$\begin{array}{c}{\mu }_4\left(\delta \right)=\frac{3{\Sigma }_1^2}{\left(n_1-2\right)\left(n_1-2\right)\left(n_1-4\right)\left(n_1-6\right)}+\frac{3{\Sigma }_2^2}{\left(n_2-2\right)\left(n_2-2\right)\left(n_2-4\right)\left(n_2-6\right)}\\ +6\frac{{\Sigma }_1{\Sigma }_2}{\left(n_1-2\right)\left(n_2-2\right)\left(n_1-4\right)\left(n_2-4\right)}\end{array}$

2) The marginal posterior p∙d∙f of the variance ratio λ is a scale multiplicative of the F-distribution. That is:

$\lambda ={\varphi }_1/{\varphi }_2=\frac{S_{1ee}/\left(n_1-2\right)}{S_{2ee}/\left(n_2-2\right)}{F}_{n_2-2,n_1-2}$

Therefore, the posterior moments of λ are:

$E\left[\lambda |x,y\right]=\frac{S_{1ee}}{S_{2ee}}\cdot \frac{\left(n_2-2\right)}{\left(n_1-4\right)}$

$\text{Mode}\left[\lambda |x,y\right]=\frac{S_{1ee}}{S_{2ee}}\cdot \frac{\left(n_2-4\right)}{n_1}$

$\mathrm{var}\left[\lambda |x,y\right]={\left(S_{1ee}/S_{2ee}\right)}^2\left[\frac{2\left(n_2-2\right)\left(n_1+n_2-6\right)}{{\left(n_1-4\right)}^2\left(n_1-6\right)}\right]$

Using the inverse moments of the F-distribution we can show that:

$E\left(\frac{1}{\lambda }|x,y\right)=\left(\frac{S_{2ee}}{S_{1ee}}\right)\left(\frac{n_1-2}{n_2-4}\right)$ (4.5)

and

$E\left(\frac{1}{{\lambda }^2}|x,y\right)={\left(\frac{S_{2ee}}{S_{1ee}}\right)}^2\left[\frac{n_1\left(n_1-2\right)}{\left(n_2-4\right)\left(n_2-6\right)}\right]$

3) The conditional posterior p∙d∙f of δ, given λ is such that

$\pi \left(\delta |\lambda ,x,y\right)\propto {\left[1+\frac{AB}{A+B}{\left(\delta -\left(b_1-b_2\right)\right)}^2\right]}^{-\left(\nu -1/2\right)}$ (4.6)

$E\left(\delta |\lambda ,x,y\right)=b_1-b_2$

$\mathrm{var}\left(\delta |\lambda ,x,y\right)=\frac{1}{n_1+n_2-6}\left[\frac{1}{A}+\frac{1}{B}\right]=\frac{1}{n_1+n_2-6}\left[\frac{S_{1ee}+\lambda S_{2ee}}{S_{1xx}}+\frac{S_{1ee}+\lambda S_{2ee}}{\lambda S_{2xx}}\right]$

These results are derived based on the fact that conditional on λ, the posterior distribution of

$t=\frac{\left(\delta -D\right)\sqrt{n_1+n_2-4}}{\sqrt{A^{-1}+B^{-1}}}$

Is that of a t-distribution with $\left(n_1+n_2-4\right)$ degrees of freedom.

5. Approximating the Posterior Distribution of δ

5.1. Moments Approximation (Scale Approximation)

As a first approximation to the posterior distribution of $\Delta =\delta -\left(b_1-b_2\right)$ is to assume that

$\Delta \stackrel{d}{=}at\left(\nu \right)$ (5.1)

Equating ${\mu }_2\left(\Delta \right)$ and ${\mu }_4\left(\Delta \right)$ to the second and fourth control moments of the R∙H∙S of $\stackrel{d}{=}$ in (5.1) we get:

$a={\left[\frac{{\mu }_2\left(\Delta \right){\mu }_4\left(\Delta \right)}{{\mu }_4\left(\Delta \right)+{\kappa }_4\left(\Delta \right)}\right]}^{1/2}$

$\nu =\left[2\left(1+\frac{{\mu }_4\left(\Delta \right)}{{\kappa }_4\left(\Delta \right)}\right)\right]$

Here, ${\kappa }_4\left(\Delta \right)={\mu }_4\left(\Delta \right)-3{\mu }_2^2\left(\Delta \right)$ and $\left[\xi \right]$ mean the smallest integer larger than ξ. Thus one may use tables of student’s t-distribution with $\left[\nu \right]$ degrees of freedom to construct H.P.D. intervals on δ. As can be seen this result is identical to the moment or the scale approximation of the posterior distribution of the difference between means.

5.2. Modal Approximation

As suggested by Box and Tiao [8] , we approximate the posterior density on replacing λ by its modal value so that:

$\pi \left(\delta |x,y\right)\circeq \pi \left(\delta |{\lambda }^{*},x,y\right)$ (5.2)

Hence we take

$T^{*}=\frac{\left(\delta -\left(b_1+b_2\right)\right)\sqrt{\left(n_1+n_2-4\right)}}{\sqrt{\left(A+B\right)/AB}}$ (5.3)

As a t-statistic with $\left(n_1+n_2-4\right)$ degrees of freedom, where $A=S_{1xx}/\left(S_{1ee}+{\lambda }^{*}{S}_{2ee}\right)$ and $B={\lambda }^{*}{S}_{2xx}/\left(S_{1ee}+{\lambda }^{*}{S}_{2ee}\right)$ , and ${\lambda }^{*}=\frac{S_{1ee}}{S_{2ee}}\cdot \frac{\left(n_2-4\right)}{n_1}$ is the modal value of λ.

5.3. Approximation Based on Averaging

Here we find that the conditions of Feller’s [8] lemma are satisfied by the two central moments of $\bar{\lambda }$ . This is because: $\bar{\lambda }=E\left(\lambda |x,y\right)=\phi =\frac{n_2-2}{n_1-4}\frac{S_{1ee}}{S_{2ee}}$ and $\mathrm{var}\left(\lambda |x,y\right)={\sigma }_n^2\left(\phi \right)=2{\phi }^2\left[\frac{n_1+n_2-6}{\left(n_1-6\right)\left(n_2-2\right)}\right]$ with ${\sigma }_n^2\left(\phi \right)\to 0$ as $\left(n_1,n_2\right)\to \infty$ .

By taking $\pi \left(\delta |x,y\right)$ as $g(\cdot )$ and $\pi \left(\lambda |x,y\right)$ as ${\psi }_n\left(\xi \right)$ , the R.H.S. of (5.2) can be approximated by

$\pi \left(\delta |x,y\right)\alpha {\left[1+a\left(\bar{\lambda }\right){\left(\delta -\left(b_1-b_2\right)\right)}^2\right]}^{-\left(\nu -1/2\right)}$

where

$a\left(\bar{\lambda }\right)={\left(\frac{1}{\bar{A}}+\frac{1}{\bar{B}}\right)}^{-1}={\left[\frac{S_{1ee}}{S_{1xx}}+\bar{\lambda }\frac{S_{2ee}}{S_{1xx}}+\frac{S_{2ee}}{S_{2xx}}+\frac{1}{\bar{\lambda }}\frac{S_{1ee}}{S_{2xx}}\right]}^{-1}$

Hence, we take

$\bar{T}=\frac{\delta -\left(b_1-b_2\right)}{\sqrt{\left({\bar{\Sigma }}_1+{\bar{\Sigma }}_2\right)}}~t_{\left(n_1+n_2-4\right)}$

where

${\bar{\Sigma }}_1=\frac{\left(n_1+n_2-6\right)}{\left(n_1-2\right)\left(n_1-4\right)}{\Sigma }_1$

${\bar{\Sigma }}_2=\frac{n_1+n_2-6}{{\left(n_2-2\right)}^2}{\Sigma }_2$

We now analyze the data in this example. We take the Ki67 to be the explanatory variable (x), while the recurrence score is the dependent variable (y). The tumor grades 1 and 3 form the two groups whose slopes are to be compared. The summary data are:

Tumor grade 1        Tumor grade 3

$n_1=40$                       

                  

                  

                    

The scatter plots are given for tumor grade 1 (Figure 1) and for tumor grade 3 (Figure 2)

Figure 1. Scatter plot of recurrence score against Ki67 for tumor grade 1.

Figure 2. Scatter plot of recurrence score against Ki67 for tumor grade 3.

As we can see, for tumor grade 1, the association between Ki67 and tumor recurrence is quite weak, but the association is stronger for tumor grade 3.

Remarks:

For all the proposed methods, our data analysis approach is simulation-based. The number of replications used is sufficient. For the Monte-Carlo integration, which we consider to be the exact we monitored the simulation. As we can see in Figure 3 the sequence of simulations tends to stabilize as we approach the specified number of simulations.

We should note that the red line bands in Figure 3 are not 95% confidence bands in classical sense, but they correspond to the confidence assessment that is produced for every number of iteration, if we decide to stop at this number of iterations. The 2.5% and 97.5% quantiles for the methods are shown in Table 2.

6. Discussion

In the data analytic part, one may be interested in the shape of the density of the approximating distributions and how they deviate from the exact density. We did not discuss this issue since most of the time we are interested in the tail area

Figure 3. Monitoring the approximation of the function with mean ± two standard deviations against the iterations for a single sequence of simulations.

Table 2. Approximate 95% HPD for the difference in two slopes of linear regression lines heterogeneity using the Tumor recurrence scores data.

of the distribution in order to construct confidence intervals on the mean difference or the difference of slopes. From Table 1 and Table 2, we can see that all the approximations perform well when compared to the exact limits. However, the scale approximation, which uses the first 4 central moments, provides more accurate confidence limits. The modal approximation seems to err in the upper limits of difference between the target parameters. Our general conclusion is that when the sample sizes are reasonably large as in the breast cancer example all the approximations (except the modal) may be used. For sample sizes and routine implementation of the proposed test, the scale (the four moments) approximation is recommended.

Conflict of Interest

None declared by all authors.

Appendix: R Codes for the Regression Model. The Codes for the Comparison between Means Is Quite Similar and Not Included

R-CODES.

#Summary data

n1=40, n2=37, s1xx=1019.6, s2xx=11167.57, s1ee=700.13, s2ee=3931.172

b1=.124, b2=.447.

#Monte Carlo integration

sig1=(n1-2)*s1ee/s1xx

sig2=(n2-2)*s2ee/s2xx

T2=(sig1/((n1-2)*(n1-4)))+(sig2/((n2-2)*(n2-4)))

ll=rf(10000,n2-2,n1-2) #simulating from the F-distribution

l=ll*(s1ee/s2ee)*((n2-2)/(n1-2))

A=s1xx/(s1ee+l*s2ee)

B=l*s2xx/(s1ee+l*s2ee)

nn=n1+n2-4

beta1=b1+rt(10000,n1-2)*sqrt(sig1/((n1-2)*(n1-4)))

beta2=b2+rt(10000,n2-2)*sqrt(sig2/((n2-2)*(n2-4)))

delta=(beta1-beta2)

q1=quantile(delta,prob=c(0.025,0.975))

q1

#Monitoring the simulation

h=function(x){(1/(1+((A*B)/(A+B))*(x-b1+b2)^2)^(((nn+1)/2)))}

x=h(l)

estint=cumsum(x)/(1:10^4)

mean(estint)

esterr=sqrt(cumsum((x-estint)^2))/(1:10^4)

mean(esterr)

plot(estint,xlab="Mean and error range",lwd=2,

ylim=mean(x)+20*c(-esterr[10^4],esterr[10^4]),ylab="")

lines(estint+2*esterr,col="red",lwd=2)

lines(estint-2*esterr,col="red",lwd=2)

#Moments (Scale) Approximation

sig1=(n1-2)*s1ee/s1xx

sig2=(n2-2)*s2ee/s2xx

T2=(sig1/((n1-2)*(n1-4)))+(sig2/((n2-2)*(n2-4)))

T41=3*sig1^2/((n1-2)^2*(n1-4)*(n1-6))

T42=3*sig2^2/((n2-2)^2*(n2-4)*(n2-6))

T412=6*sig1*sig2/((n1-2)*(n1-4)*(n2-2)*(n2-4))

T4=T41+T42+T412

K4=T4-3*T2^2

a=sqrt(T2*T4/(T4+K4))

b=ceiling(2*(1+(T4/K4))) #integer degrees of freedom

d=rt(10000,b)

delta.m=(b1-b2)+a*d

q2=quantile(delta.m,prob=c(.025,.975))

q2

#Welch Normal Approximation

zz=rnorm(10000)

u=(b1-b2)+zz*sqrt(T2)

mean.u=mean(u)

sd.u=sd(u)

qqnorm(u)

q3=quantile(u,prob=c(0.025,0.975))

q3

#Edgeworth expansion

library(distr)

aa=K4/(24*T2^2)

p <- function(x) (1/sqrt(2*pi) *(1/(exp(x^2/2)*

(1+aa*(x^4-6*x^2+3)))))

# probability density function

dist <-AbscontDistribution(d=p) # signature for a distribution with pdf ~ p

rdist <- r(dist) # function to create random variates from p

set.seed(1)

XX <- rdist(10000) # sample from X ~ p

x <- seq(-10,10, .01)

hist(XX, freq=F, breaks=50, xlim=c(-5,5))

lines(x,p(x),lty=2, col="red")

mean(XX)

sd(XX)

edgeworth=b1-b2+sqrt(T2)*XX

hist(edgeworth)

q4=quantile(edgeworth,prob=c(0.025,.975))

q4

#Averaging approximation based on Feller’s lemma

nu=n1+n2-4

f=rf(10000,n2-2,n1-2)

t=rt(10000,nu)

l=(s1ee/s2ee)*f

mean.l=mean(l)

denom=s1ee+mean.l*s2ee

A=s1xx/denom

B=mean.l*s2xx/denom

d.mean=(b1-b2)+ (sqrt((A+B)/(A*B))*t)/sqrt(nu)

q5=quantile(d.mean,prob=c(.025,.975))

q5

#Modal approximation

ff=rf(10000,n2-2,n1-2)

f=ff*s1ee*(n2-2)/s2ee*(n1-2)

m=(s1ee/s2ee)*(n2-4)/n1

nu=n1+n2-4

t=rt(10000,nu)

denom=s1ee+m*s2ee

AA=s1xx/denom

BB=m*s2xx/denom

d.mod=(b1-b2)+ (sqrt((AA+BB)/(AA*BB))*t)/sqrt(nu)

q6=quantile(d.mod,prob=c(.025,.975))

q6

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Jeffreys, B. (1940) Note on the Behrens-Fisher Formula. Annals of Eugenics, 10, 48-51. https://doi.org/10.1111/j.1469-1809.1940.tb02236.x
[2]	Behrens, W.V. (1929) Ein Beitragzur Fehlerberechunung bei Weniger Beoachtungen. Landw. Jb., 68, 807-837.
[3]	Fisher, R.A. (1939) The Comparison of Samples with Possibly Unequal Variances. Annals of Eugenics, 9, 174-180. https://doi.org/10.1111/j.1469-1809.1939.tb02205.x
[4]	Reilly, P.M. (1976) The Numerical Computation of Posterior Distributions. Applied Statistics, 25, 201-209. https://doi.org/10.2307/2347227
[5]	Naylor, J.C. and Smith, A.F.M. (1982) Applications of a Method for the Efficient Computation of Posterior Distributions. Applied Statistics, 31, 214-225. https://doi.org/10.2307/2347995
[6]	Robert, C. and Casella, G. (2009) Introducing Monte Carlo Methods with R. Springer, New York.
[7]	Lee, P.M. (1989) Bayesian Statistics: An Introduction. Oxford University Press.
[8]	Box, G.E.P. and Taio, G.C. (1973) Bayesian Inference in Statistical Analysis. Addison-Wesley Publishing Company, Inc., Reading, Mass.
[9]	Barnard, G.A. (1984) Comparing the Means of Two Independent Samples. Applied Statistics, 33, 266-271. https://doi.org/10.2307/2347702
[10]	Patil, V.H. (1964) The Behrens-Fisher Problem and Its Bayesian Solution. Journal of the Indian Society for Probability and Statistics, 2, 21-31.
[11]	Feller, W. (1971) An Introduction to Probability Theory and Its Applications. 2nd Edition, John Wiley & Sons, Inc., New York.
[12]	Barndorff-Nielson, O. and Cox, D.R. (1979) Edgeworth and Saddle-Point Approximations with Statistical Applications. Journal of the Royal Statistical Society: Series B, 41, 279, 312.
[13]	Edgeworth, F.Y. (1896) The Asymmetrical Probability Curve. Philosophical Magazine, 5th Series, 41, 90-99.
[14]	Edgeworth, F.Y. (1907) On the Representation of Statistical Frequency by Series. Journal of the Royal Statistical Society: Series A, 70, 102-106. https://doi.org/10.2307/2339504
[15]	Klein, M.E., Dabbs, D.J., Shuai, Y., Brufsky, A.M., Jankowitz, R., Puhalla, S.L., et al. (2013) Prediction of the Oncotype DX Recurrence Score: Use of Pathology-Generated Equations Derived by Linear Regression Analysis. Modern Pathology, 26, 658-664.
[16]	Welch, B.L. (1937) The Significance of the Difference between Two Means when the Population Variances Are Unequal. Biometrika, 34, 28-35.
[17]	Cuzick, J., Dowsett, M., Pineda, S., Wale, C., Salter, J., Quinn, E., et al. (2011) Prognostic Value of a Combined Estrogen Receptor, Progesterone Receptor, Ki-67, and Human Epidermal Growth Factor Receptor 2 Immuno-Histochemical Score and Comparison with the Genomic Health Recurrence Score in Early Breast Cancer. Journal of Clinical Oncology, 29, 4273-4278.
[18]	Thakur, S.S., Li, H.C., Angela, M.Y., Chan, R.T., Bigras, G., Morris, D., Enwere, E.K. and Yang, H. (2018) The Use of Automated Ki67 Analysis to Predict Oncotype DX Risk-of-Recurrence Categories in Early-Stage Breast Cancer: January 5, 2018. PLOS ONE.