Simple Bayesian testing of scientific expectations in linear regression models
Joris Mulder Anton Olsson-Collentine



1,2]J. Mulder
3]A. Olsson-Collentine
[1]Department of Methodology and Statistics, Tilburg University, Warrandelaan 1, Tilburg, the Netherlands. Email: j.mulder3@tilburguniversity.edu
[2]Jheronimus Academy of Data Science, Sint Janssingel 92, 5211 DA 's-Hertogenbosch, The Netherlands
[3]Department of Methodology and Statistics, Tilburg University, Warrandelaan 1, Tilburg, the Netherlands. Email: j.a.e.olssoncollentine@tilburguniversity.edu
                     





Simple Bayesian Testing of Scientific Expectations in Linear Regression Models
    [
    
==============================================================================




Scientific theories can often be formulated using equality and order constraints on the relative effects in a linear regression model. For example, it may be expected that the effect of the first predictor is larger than the effect of the second predictor, and the second predictor is expected to be larger than the third predictor. The goal is then to test such expectations against competing scientific expectations or theories. In this paper a simple default Bayes factor test is proposed for testing multiple hypotheses with equality and order constraints on the effects of interest.  
The proposed testing criterion can be computed without requiring external prior information about the expected effects before observing the data. The method is implemented in R-package called `lmhyp' which is freely downloadable and ready to use. The usability of the method and software is illustrated using empirical applications from the social and behavioral sciences.


Key words: Bayes factors, Bayesian hypothesis testing, equality and order constraints, regression modeling



§ INTRODUCTION

The linear regression model is the most widely used statistical method for assessing the relative effects of a given set of predictors on a continuous outcome variable. This assessment of the relative effects is an essential part when testing, fine-graining, and building scientific theories. For example, in work and organizational psychology the regression model has been used to better understand the effects of discrimination by coworkers and managers on workers' well-being <cit.>; in sociology to assess the effects of the different dimensions of socioeconomic status on one's attitude towards immigrants <cit.>; and in experimental psychology to make inferences regarding the effects of gender when hiring employees (Carlson and Sinclair, 2017). Despite the extensive literature on statistical tools for linear regression analysis, methods for evaluating multiple hypotheses with equality and order constraints on the relative effects in a direct manner are still limited. This paper presents a Bayes factor testing procedure with accompanying software for testing such hypotheses with the goal to aid researchers in the development and evaluation of scientific theories.

As an example, let us consider the following linear regression model where a dependent variable is regressed on three predictor variables, say, X_1, X_2, and X_3:

    y_i   =   β_0 +β_1 X_i,1 + β_2 X_i,2 + β_3 X_i,3 + ϵ_i,

where y_i is the dependent variable of the i-th observation, X_i,k denotes the k predictor variable of the i-th observation, β_k is the regression coefficient of the k-th predictor, for k=1,…,3, β_0 is the intercept, and ϵ_i are independent normally distributed errors with unknown variance σ^2, for i=1,…,n.

In exploratory studies the interest is typically whether each predictor has an effect on the dependent variable, and if there is evidence of a nonzero effect, we would be interested in whether the effect is positive or negative. In the proposed methodology such an exploratory analysis can be executed by simultaneously testing whether an effect is zero, positive, or negative. For the first predictor, the exploratory multiple hypothesis test would be formulated as

    H_0   :   β_1=0
    H_1   :   β_1>0
     H_2   :   β_1<0.

The proposed Bayes factor test will then provide a default quantification of the relative evidence in the data between these hypotheses.

In confirmatory studies, the interest is typically in testing specific hypotheses with equality and order constraints on the relative effects based on scientific expectations or psychological theories <cit.>. Contrasting regression effects against each other using equality or order constraints can be more informative than interpreting them at certain benchmark values (e.g., standardized effects of .2, .5, and 1, are sometimes interpreted as `small', `medium', and `large' effects, respectively) because effects are not absolute but relative quantifications; relative to each other and relative to the scientific field and context <cit.>. For example, a standardized effect of .4 may be important for an organizational psychologist who is interested in the effect of discrimination on well-being on the work floor but less so for a medical psychologist who wishes to predict the growth of a tumor of a patient through a cognitive test. As such, interpreting regression effects relative to each other using equality and order constraints would be more insightful than interpreting the effects using fixed benchmarks.

In the above regression model for instance, let us assume that β_1, β_2, and β_3 denote the effects of a strong, medium and mild treatment, respectively. It may then be hypothesized that the effect of the strong treatment is larger than the effect of the medium treatment, the effect of the medium treatment is expected to be larger than the effect of the mild treatment, and all effects are expected to be positive. Alternatively it may be expected that all treatments have an equal positive effect. These hypotheses can then be tested against a third hypothesis which complements the other hypotheses. This comes down to the following multiple hypothesis test:

    H_1   :   β_1 > β_2>β_3>0
    H_2   :   β_1 = β_2=β_3>0
     H_3   :   neither H_1, nor H_2.

Here the complement hypothesis H_3 covers the remaining possible values of β_1, β_2, and β_3 that do not satisfy the constraints under H_1 and H_2. Subsequently the interest is in quantifying the relative evidence in the data for these hypotheses.

A general advantage of Bayes factors for testing statistical hypotheses is that we obtain a direct quantification of the evidence in the data in favor of one hypothesis against another hypothesis. Furthermore, Bayes factors can be translated to the posterior probabilities of the hypotheses given the observed the data and the hypotheses of interest. These probabilities give a direct answer to the research question which hypothesis is most likely to be true and to what degree given the data. These posterior probabilities can be used to obtain conditional error probabilities of drawing an incorrect conclusion when `selecting' a hypothesis in light of the observed data. These and other properties have greatly contributed to the increasing popularity of Bayes factors for testing hypotheses in psychological research <cit.>.

The proposed testing criterion is based on the prior adjusted default Bayes factor <cit.>. The method has various attractive properties. First, the proposed Bayes factor has an analytic expression when testing hypotheses with equality and order constraints in a regression model. Thus computationally demanding numerical approximations can be avoided resulting in a fast and simple test. Furthermore, by allowing users to formulate hypotheses with equality as well as ordinal constraints a broad class of hypotheses can be tested in an easy and direct manner. Another useful property is that no proper (subjective) prior distribution needs to be formulated based on external prior knowledge, and therefore the method can be applied in an automatic fashion. This is achieved by adopting a fractional Bayes methodology <cit.> where a default prior is implicitly constructed using a minimal fraction of the information in the observed data and the remaining (maximal) fraction is used for hypothesis testing <cit.>. This default prior is then relocated to the boundary of the constrained space of the hypotheses. In the confirmatory test about the unconstrained default prior for (β_1,β_3,β_3) would be centered around 0. Because this Bayes factor can be computed without requiring external prior knowledge it is called a `default Bayes factor'. Thereby, these default Bayes factors differ from regular Bayes factors where a proper prior is specified reflecting the anticipated effects based on external prior knowledge <cit.>. Other default Bayes factors that have been proposed in the literature are the fractional Bayes factor <cit.>, the intrinsic Bayes factor <cit.>, and the Bayes factor based on expected-posterior priors <cit.>.

Although various alternative testing procedures are available for hypothesis testing for linear regression analysis, these methods are limited to some degree. First, classical significance tests are only suitable for testing a null hypothesis against a single alternative, and unsuitable for testing multiple hypotheses with equality as well as order constraints <cit.>. Second, traditional model comparison tools (e.g., the AIC, BIC, or CFI) are generally not suitable for evaluating models (or hypotheses) with order constraints on certain parameters <cit.>. Third, currently available Bayes factor tests cannot be used for testing order hypotheses <cit.>, are not computationally efficient <cit.>, or are based on large sample approximations <cit.>. The proposed Bayes factor, on the other hand, can be used for testing hypotheses with equality and/or order constraints, is very fast to compute due to its analytic expression, and is an accurate default quantification of the evidence in the data in the case of small to moderate samples because it does not rely on large sample approximations. Other important properties of the proposed methodology are its large sample consistent behavior and its information consistent behavior <cit.>.

The Bayesian test is implemented in the R-package `lmhyp', which is freely downloadable and ready for use in R. The main function `test_hyp' needs a fitted modeling object using the `lm' function together with a string that formulates a set of hypotheses with equality and order constraints on the regression coefficients of interest. The function computes the Bayes factors of interest as well as the posterior probabilities that each hypothesis is true after observing the data.

The paper is organized as follows. Section 2 presents the derivation of the default Bayes factor between hypotheses with equality and order hypotheses on the relative effects in a linear regression model. Section 3 presents the `lmhyp' package and explains how it can be used for testing scientific expectations in psychological research. Section 4 shows how to apply the new procedure and software for testing scientific expectations in work and organizational psychology and social psychology. The paper ends with a short discussion.




§ A DEFAULT BAYES FACTOR FOR EQUALITY AND ORDER HYPOTHESES IN A LINEAR REGRESSION MODEL



 §.§ Model and hypothesis formulation

For a linear regression model,

    y=Xβ+N(0,σ^2I_n),

where y is a vector of length n of outcome variables, X is a n× k matrix with the predictor variables, and β is a vector of length k containing the regression coefficients, consider a hypothesis with equality and inequality constraints on certain regression coefficients of the form

    H_t:R_Eβ=r_E    &   R_Iβ>r_I,

where [R_E | r_E] and [R_I|r_I] are the augmented matrices with q_E and q_I rows that contain the coefficients of the equality and inequality constraints, respectively, and k+1 columns. For example, for the regression model from the introduction, with β=(β_0,β_1,β_2,β_3)', and the hypothesis H_1:β_1 > β_2>β_3>0 in (<ref>), the augmented matrix of the inequalities is given by

    [R_I|r_I] = [
    [  0  1 -1  0  0;  0  0  1 -1  0;  0  0  0  1  0 ]]

and for the hypothesis H_2:β_1 = β_2=β_3>0, the augmented matrices are given by

    [R_E|r_E]    =   [
    [  0  1 -1  0  0;  0  0  1 -1  0;    ]]
    [ R_I|r_I]     =   [
    [ 0 0 0 1 0 ]]


The prior adjusted default Bayes factor will be derived for a constrained hypothesis in (<ref>) against an unconstrained alternative hypothesis, denoted by H_u:β∈ℝ^k, with no constraints on the regression coefficients. First we transform the regression coefficients as follows

    ξ=[[ ξ_E; ξ_I ]]=[[ R_E;   D ]]β=Tβ,

where D is a (k-q_E)× k matrix consisting of the unique rows of I_k-R_E'(R_ER_E)^-1R_E. 
Thus, ξ_E is a vector of length q_E and ξ_I is a vector of length k-q_E. Consequently, model (<ref>) can be written as

    y = XR_E^-1ξ_E + XD^-1ξ_I + N(0,σ^2I_n),

because

    Xβ=XT^-1ξ=X[R_E^-1 D^-1][[ ξ_E; ξ_I ]]=XR_E^-1ξ_E+XD^-1ξ_I,

where R_E^-1 and D^-1 are the (Moore-Penrose) generalized inverse matrices of R_E and D,
and the hypothesis in (<ref>) can be written as

    H_t:ξ_E=r_E    &   R̃_Iξ_I>r̃_I,

because

    R_Eβ=R_IT^-1ξ=
    R_E[R_E^-1 D^-1]ξ=[I_q_E 0]ξ=ξ_E=r_E

and

    R_Iβ=R_IT^-1ξ=R_I[R_E^-1 D^-1]
    [[ r_E; ξ_I ]]=R_IR_E^-1r_E+R_ID^-1ξ_I>r_I
          ⇔R̃_Iξ_I>r̃_I, with r̃_I=r_I-R_IR_E^-1r_E and R̃_I=R_ID^-1.




 §.§ A default Bayes factor for testing hypotheses

The Bayes factor for hypothesis H_1 against H_2 is defined as the ratio of their respective marginal likelihoods,

    B_12=p_1(y)/p_2(y).

The marginal likelihood quantifies the probability of the observed data under a hypothesis <cit.>. For example, if B_12 = 10 this implies that the data were 10 times more likely to have been observed under H_1 than under H_2. Therefore, the Bayes factor can be seen as a relative measure of evidence in the data between two hypotheses. The marginal likelihood under a constrained hypothesis H_t in (<ref>) is obtained by integrating the likelihood over the order constrained subspace of the free parameters weighted with the prior distribution,

    p_t(y) = ∬_R_Iβ>r_I p_t(y|β,σ^2)π_t(β,σ^2)dβ dσ^2,

where p_t(y|β,σ^2) denotes the likelihood of the data under hypothesis H_t given the unknown model parameters, and π_t denotes the prior distribution of the free parameters under H_t. The prior quantifies the plausibility of possible values that the model parameters can attain before observing the data.

Unlike in Bayesian estimation, the choice of the prior can have a large influence on the outcome of the Bayes factor. For this reason ad hoc or arbitrary prior specification should be avoided when testing hypotheses using the Bayes factor. However, specifying a prior that accurately reflects one's uncertainty about the model parameters before observing the data can be a time-consuming and difficult task <cit.>. A complicating factor in the case of testing multiple, say, 3 or more, hypotheses, is that priors need to be carefully formulated for the free parameters under all hypotheses separately. Because noninformative improper priors also cannot be used when computing marginal likelihoods, there has been an increasing interest in the development of default Bayes factors where ad hoc or subjective prior specification is avoided. In these default Bayes factors a proper default prior is often (implicitly) constructed using a small part of the data while the remaining part is used for hypothesis testing. An example is the fractional Bayes factor <cit.> where the marginal likelihood is defined by

    p_t(y) = ∬_R_Iβ>r_I p_t(y|β,σ^2)^1-bπ_t(β,σ^2|y^b)dβ dσ^2,

where the (subjective) proper prior in (<ref>) is replaced by a proper default prior based on a (minimal) fraction “b” of the observed data[The proper default prior in (<ref>) is obtained by updating the noninformative improper (independence) Jeffreys' prior, π^N(β,σ^2)∝σ^-2, with a fraction b of the data: π_t(β,σ^2|y^b)∝π^N(β,σ^2)f_t(y|β,σ^2)^b <cit.>.], and the likelihood is raised to a power equal to the remaining fraction “1-b”, which is used for hypothesis testing.

In this paper an adjustment of fractional Bayes factor is considered where the default prior is centered on the boundary (or null value) of the constrained space. The motivation for this adjustment is two-fold. First when testing a precise hypothesis, say, H_0:β=0 versus H_a:β≠0, Jeffreys argued that a default prior for β under H_1 should be concentrated around the null value because if the null would be false, the true effect would likely to be close to the null, otherwise there would be no point in testing H_0. Second, when testing hypotheses with inequality or order constraints, the prior probability that the constraints hold serves as a measure of the relative complexity (or size) of the constrained space under a hypothesis <cit.>. This quantification of relative complexity of a hypothesis is important because the Bayes factor balances fit and complexity as an Occam's razor. This implies that simpler hypotheses (i.e., hypotheses having “smaller” parameter spaces) would be preferred over more complex hypotheses in the case of an approximately equal fit. Only when centering the prior at 0 when testing H_1:β<0 versus H_2:β>0, both hypotheses would be considered as equally complex with prior probabilities of .5 corresponding to half of the complete parameter space of β of all real values (ℝ).

Given the above considerations, the fractional Bayes factor is adjusted such that the default prior is (i) centered on the boundary of the constrained parameter space and (ii) contains minimal information by specifying a minimal fraction. Because the model consists of k+1 unknown parameters (k regression coefficients and an unknown error variance), a default prior is obtained using a minimal fraction[Updating the noninformative Jeffreys prior π^N(β,σ)∝σ^-2 with a sample of k+1 observations yields a proper marginal distribution for β having a multivariate Student t distribution with 1 degree of freedom, which is equivalent to a multivariate Cauchy distribution.] of b=k+1/n.

In order to satisfy the prior property (i) when testing a hypothesis (<ref>), the prior for β under the alternative should thus be centered at R^-1r, where R'=[R_E' R_I'] and r'=(r_E',r_I'), which is equivalent to centering the prior for ξ at μ^0=(μ_E^0',μ_I^0')'=TR^-1r=(r'_E,μ_I^0')', with R̃_Iμ_I^0=r̃_I. The following lemma gives the analytic expression of the default Bayes factor of a hypothesis with equality and order constraints on the regression coefficients versus an unconstrained alternative.


The prior adjusted default Bayes factors for an equality-constrained hypothesis, H_1:R_Eβ=r_E, an order-constrained hypothesis, H_2:R_Iβ>r_I, and a hypothesis with equality and order constraints, H_3:R_Eβ=r_E, R_Iβ>r_I, against an unconstrained hypothesis H_u:β∈ℝ^k are given by

    B_1u   =   f^E_1/c^E_1=
    t(r_E;R_Eβ̂,s^2(n-k)^-1R_E(X'X)^-1R_E',n-k)/t(r_E;r_E,s^2R_E(X'X)^-1R_E',1),
    
    B_2u   =   f^I_2/c^I_2=
    Pr(R_Iβ>r_I|y,H_u)/Pr(R_Iβ>r_I|y^b,H_u),
    
    B_3u   =   f^E_3/c^E_3×f^I|E_3/c^I|E_3=t(r_E;R_Eβ̂,s^2(n-k)^-1R_E(X'X)^-1R_E',n-k)/t(r_E;r_E,s^2R_E(X'X)^-1R_E',1)
          ×Pr(R̃_Iξ_I>r̃_I|ξ_E=r_E,y,H_u)/Pr(R̃_Iξ_I>r^*_I|ξ_E=ξ̂_E,y^b,H_u),

where r^*_I=R̃_Iξ̂_I, t(ξ;μ,S,ν) denotes a Student t density for ξ with location parameter μ, scale matrix S, and degrees of freedom ν, β̂=(X'X)^-1X'y is the maximum likelihood estimate (MLE) of β and s^2=(y-Xβ̂)'(y-Xβ̂) is the sums of squares, and the (conditional) distributions are given by

    π(β|y,H_u)   =   t(β;β̂,s^2(X'X)^-1/(n-k),n-k)
    π(β|y^b,H_u)   =   t(β;R^-1_Ir_I,s^2(X'X)^-1,1)
    π(ξ_I|ξ_E=r_E,y,H_u)   =   t(ξ_I;μ_I^N,S_I^N,n-k)
    π(ξ_I|ξ_E=r_E,y^b,H_u)   =   t(ξ_I;μ_I^0,S_I^0,1)

with

    μ_I^N   =   Dβ̂+D(X'X)^-1R_E'(R_E(X'X)^-1R_E')^-1(r_E-R_Eβ̂)
    S_I^N   =   (1+s^-2(r_E-R_Eβ̂)'(R_E(X'X)^-1R_E')^-1(r_E-R_Eβ̂))(n-k+q_E)^-1s^2
          (D(X'X)^-1D'-D(X'X)^-1R_E'(R_E(X'X)^-1R_E')^-1R_E(X'X)^-1D')
    S_I^0   =   s^21+q^E(D(X'X)^-1D'-D(X'X)^-1R_E'(R_E(X'X)^-1R_E')^-1R_E(X'X)^-1D'),


Proof: Appendix A.


Note that the factors in (<ref>) and (<ref>) are multivariate Savage-Dickey density ratio's <cit.>. These ratio's have an analytic expression because the marginal posterior and default prior have multivariate Student t distributions. In R these can be computed using the dmvt function in the mvtnorm-package <cit.>. 

The ratio's of (conditional) probabilities in (<ref>) and (<ref>) can also be computed in a straightforward manner. If R̃_I is of full row-rank then the transformed parameter vector, say, η_I=R̃_Iξ_I has a Student t distribution so that Pr(R̃_Iξ_I>r̃_I|ξ_E=r_E,y,H_u)=Pr(η_I>r̃_I|ξ_E=r_E,y,H_u) can be computed using the pmvt function from the mvtnorm-package <cit.>. If the rank of R̃_I is lower than q^I, then the probability can be computed as the proportion of draws from an unconstrained Student t distribution satisfying the order constraints.

The posterior quantities in the numerators reflect the relative fit of a constrained hypothesis, denoted by “f”, relative to the unconstrained hypothesis: a larger posterior probability implies a good fit of the order constraints and a large posterior density at the null value indicates a good fit of a precise hypothesis. The prior quantities in the denominators reflect the relative complexity of a constrained hypothesis, denoted by “c”, relative to the unconstrained hypothesis: a small prior probability implies a relatively small inequality constrained subspace, and thus a `simple' hypothesis, and a small prior density at the null value corresponds to a large spread (variance) of possible values under the unconstrained alternative implying the null hypothesis is relatively simple in comparison to the unconstrained hypothesis.

Figure <ref> gives more insight about the nature of the expressions in (<ref>) to (<ref>) in Lemma <ref> for an equality constrained hypothesis, H_1:β_1=β_2=0 (upper panels), an inequality constrained hypothesis, H_2:β>0 (middle panels), and hypothesis with an equality constraint and an inequality constraint, H_3:β_1>β_2=0 (lower panels). The Bayes factor for H_1 against the unconstrained hypothesis H_u in (<ref>) corresponds to the ratio of the unconstrained posterior density and the unconstrained default prior (which has a multivariate Cauchy distribution centered at the null value) evaluated at the null value. The Bayes factor for H_2 against H_u in (<ref>) corresponds to the ratio of posterior and default prior probabilities that the constraints hold under H_u. In the case of independent predictors, for example, the prior probability would be equal .25 as a result of centering the default prior at 0. The inequality constrained hypothesis would then be quantified as 4 times less complex than the unconstrained hypothesis. Finally for a hypothesis with equality and inequality constraints, H_3:β_1>β_2=0, the Bayes factor in (<ref>)-(<ref>) corresponds to the ratio of the surfaces of cross section of the posterior and prior density on the line β_1>0, β_2=0.



The default Bayes factors between these hypotheses are computed for a simulated data set with MLEs (β̂_1,β̂_2)=(.7,.03) (Appendix B) that results in  B_1u=f^E_1/c^E_1=0.061/0.159=0.383, B_2u=f^I_2/c^I_2=0.546/0.250=2.183, B_3u=f^E_3/c^E_3×f^I|E_3/c^I|E_3=1.608/0.318×0.996/0.500=10.061. As will be explained in the next section, it is recommendable to include the complement hypothesis in an analysis. The complement hypothesis covers the subspace of ℝ^2 that excludes the subspaces under H_1, H_2, and H_3. In this example the Bayes factor of the complement hypothesis against the unconstrained hypothesis equals B_cu=1-f^I_2/1-c^I_2=0.454/0.750=0.606.

After having obtained the default Bayes factor of each hypothesis against the unconstrained hypothesis, Bayes factors between the hypotheses of interest can be obtained through the transitivity property of the Bayes factor, e.g., B_31=B_3u/B_1u=10.061/.383=26.299. This implies there is strong evidence for H_3 relative to H_1, as the data were approximately 26 times more likely to have been produced under H_3 than under H_1.























Once the default Bayes factors of the hypotheses of interest against the unconstrained hypothesis are computed using Lemma 1, posterior probabilities can be computed for the hypotheses. In the case of, say, four hypotheses of interest against, the posterior probability that hypothesis H_t is true can be obtained via

    Pr(H_t|y)= B_tuPr(H_t)/B_1uPr(H_1)+B_2uPr(H_2)+B_3uPr(H_3)
    +B_cuPr(H_c),

for t=1, 2, or 3, where Pr(H_t) denotes the prior probability of hypothesis H_t, i.e., the probability that H_t is true before observing the data. As can be seen, the posterior probability is a weighted average of the Bayes factors weighted with the prior probabilities. Throughout this paper we will work with equal prior probabilities, but other choices may be preferred in specific applications <cit.>. For the example data from Appendix B and Figure <ref>, the posterior probabilities would be equal to P(H_1|y)=0.029, P(H_2|y)=0.165, P(H_3|y)=0.760, and P(H_c|y)=0.046. Based on these outcomes we would conclude that there is most evidence for H_3 that the effect of the first predictor is positive and the effect of the second predictor is zero with a posterior probability of .76. In order to draw a more decisive conclusion (e.g., when obtaining a posterior probability for a hypothesis larger than, say, .99) more data are needed.


























§ SOFTWARE

The Bayes factor testing criterion for evaluating inequality and order constrained hypotheses was implemented in a new R package called `lmhyp' to ensure general utilization of the methodology[Run `devtools::install_github("jaeoc/lmhyp")' in R to install the package.]. As input the main function `test_hyp' needs a fitted linear regression modeling object from the lm-function as well as a string that specifies the constrained hypotheses of interest.

As output the function provides the default Bayes factors between all pairs of hypotheses. By default a complement hypothesis is also included in the analysis. For example, when testing the hypotheses, say, H_1:β_1>β_2>β_3>0 versus H_2:β_1=β_2=β_3>0, a third complement hypothesis H_3 will be automatically added which covers the remaining parameter space, i.e., ℝ^3 excluding the subspaces under H_1 and H_2. The reason for including the complement hypothesis is that Bayes factors provide relative measures of evidence between the hypotheses. For example, it may be that H_2 receives, say, 30 times more evidence than H_1, i.e., B_21=30, which could be seen as strong evidence for H_2 relative to H_1, yet it may be that H_2 still badly fits to the data in an absolute sense. In this case the evidence for the complement hypothesis H_3 against H_2 could be very large, say, B_32=100.

Besides the default Bayes factor the function also provides the posterior probabilities of the hypotheses. Posterior probabilities may be easier for users to interpret than Bayes factors because the posterior probabilities sum up to 1. Note that when setting equal prior probabilities between two hypotheses, the posterior odds of the hypotheses will be equal to the Bayes factor. By default all hypotheses receive equal prior probabilities. Thus, in the case of T hypotheses, then P(H_t)=1/T, for t=1,…,T. Users can manually specify the prior probabilities by using the `priorprobs' argument. In the remaining part of the paper we will work with the default setting of equal prior probabilities. A step-by-step guide for using the software will be provided in the following section.




§ APPLICATION OF THE NEW TESTING PROCEDURE USING THE SOFTWARE PACKAGE `LMHYP'

In this section we illustrate how to use the `lmhyp' package to test hypotheses, by applying the procedure to two empirical examples from psychology. We begin by describing the published research. In the two following subsections we then formulate hypotheses for each example and test these using the function  from our R-package .[The R-script used to produce the results in this section is available at https://osf.io/g8c9p/]

For the first example we use data from a study of mental health workers in England <cit.>. The data of Johnson et al. measured health workers' well-being and its correlates, such as perceived discrimination from managers, coworkers, patients and visitors. Well-being was operationalized by scales measuring anxiety, depression and job dissatisfaction, the first two scales consisting of three items and the latter of five. The perceived discrimination variables are binary variables  that were meant to capture whether the worker believed they had been discriminated against from the four different sources in the last 12 months. This example demonstrates hypothesis testing in regards to single variables and the "exploratory" option of the test_hyp function.

Our second empirical example comes from research by <cit.>. Over four experiments, Carlson and Sinclair compare two theoretical explanations for perceptions of gender discrimination in hiring, although we use data from only the first experiment (available at https://osf.io/qcdgp/). In this study Carlson and Sinclair showed university students two fictive job applications from a man and a woman for a position as either a computer specialist or nurse. Participants were told that the fictive job applications had been sent to real companies as part of a previous study, but that only one of the two applicants had been invited to a job interview despite being equally qualified. A two-item scale was then used to measure participants' belief the outcome was due to gender discrimination. Several potential correlates were also measured using two-item scales, such as the individual's belief that (wo)men are generally discriminated against, their expectation that they are gender-stereotyped by others (`stigma consciousness') and the extent to which they identify as feminists. This example demonstrates testing hypotheses involving multiple variables.




 §.§ Hypothesis testing of single effects in organizational psychology

In our first example we illustrate how our approach might be used to explore competing hypotheses for single variables. It is common when testing the effect of an independent variable in regression to look at whether it is significantly different from zero, or to do a one-sided test of a positive versus a negative effect. When using a Bayes factor test we can test all these hypotheses directly against each other and compare the relative evidence for each hypothesis.

<cit.> theorized that work-place discrimination has a negative impact on workers' well-being. Here we are testing this expectation against a positive effect and a zero effect, while controlling for discrimination from different sources. For example, in the case of discrimination by managers we have

    H_1: β_manager < 0 
    
    H_2: β_manager = 0 
    
    H_3: β_manager > 0,

while controlling for discrimination by coworkers, patients, and visitors through the following regression model

    y_anxiety,i   =   β_0 + β_manager  X_manager,i + β_coworkers  X_coworker,i
           + β_patient  X_patient,i + β_visitor  X_visitor,i + error_i

where the β's are the regression effects of the various sources of discrimination on anxiety.

Evaluating these three hypotheses in R is straightforward with the  function from our R-package . This function takes as arguments `object', a fitted object using the lm function, `', a string vector specifying one or several hypotheses (separated by semicolons), `priorprob', specifying the prior probabilities of each hypotheses (by default equal, priorprob = 1),  and `mcrep', an integer that specifies the number of draws to compute the prior and posterior probabilities in the (unusual) case the matrix with the coefficients of the order constraints is not of full row rank (by default mcrep=1e6). In addition, the argument hyp also allows as input the string "exploratory", which will test the likelihood of the data for a zero, positive, or negative effect of all variables in the regression model, including the intercept. We will make use of this functionality below, after first discussing how to test the three hypotheses for a single variable. To test the hypotheses, we first fit a linear model on the variables as usual:



Next,  hypotheses are specified in R as character strings using the variable names from the fitted linear model. It is possible to test the traditional null hypothesis of β_manager=0 against the two-sided alternative example β_manager≠0 by writing



Note that the complement hypothesis, β_manager≠0, is automatically included. However, by testing whether the effect is zero, positive, or negative simultaneously, we obtain a more complete picture of the possible existence and direction of the population effect. This can be achieved by specifying all hypotheses as a single character vector in which the hypotheses are separated by semicolons:



Note that spacing does not matter. Once the hypotheses have been specified, they are tested by simply inputting them together with the fitted linear model object into the function :



This will compute the default Bayes factors from Lemma 1 between the hypotheses, as well as the posterior probabilities for the hypotheses. The posterior probabilities are printed as the primary output:



As can be seen the evidence is overwhelmingly in favor of a positive effect of discrimination from managers on anxiety amongst health workers.  In fact when concluding that H_3:β_manager>0 is true, we would have a conditional error probability of drawing the wrong conclusion of approximately zero. To perform this test for all regression effects one simply needs to set the second hyp argument equal to "exploratory": 



This option assumes that each hypothesis is equally likely a priori. In the current example we then get the following output:



The posterior probabilities for discrimination by managers are the same as when tested separately. In regards to the other variables, there seems to be positive evidence that there is no effect of discrimination by coworkers, patients, or visitors on anxiety. Note that the evidence for this is not as compelling as for the effect of discrimination by managers, as can be seen from the conditional error probabilities of .216, .369, and .082, respectively, which are quite large. Therefore more data are needed in order to draw more decisive conclusions. Note here that classical significance tests cannot be used for quantifying the evidence in the data in favor of the null; the classical test can only be used to falsify the null. When a null hypothesis cannot be rejected we are left in a state of ignorance because we cannot reject the null but also not claim there is evidence for the null <cit.>. 

Because the prior probabilities of the hypotheses are equal, the ratio of the posterior probabilities of two hypotheses corresponds with the Bayes factor, e.g., B_23=Pr(H_2|y)/Pr(H_3|y)=.780/.216=3.615, for the effect of discrimination by coworkers. By calling , we obtain the default Bayes factors between all pairs of hypotheses. For convenience the printed Bayes factors are rounded to three digits, though exact values can be calculated from the posterior probabilities (unrounded posterior probabilities are available by calling ). The Bayes factor matrix for discC (discrimination from coworkers) can be obtained by calling



Hence, the null hypothesis of no effect is 162 times more likely than hypothesis H_1 which assumes a negative effect (B_21 = 162.367), but only 3.6 times more likely than hypothesis H_3 which assumes a positive effect (B_23 = 3.615). Similar Bayes factor matrices can be printed for all variables when using the "exploratory" option.

To summarize the first application, regressing the effects of perceived discrimination from managers, coworkers, patients, and visitors on the anxiety levels of English health workers, we found very strong evidence for a positive effect of perceived discrimination from managers on anxiety, mild to moderate evidence for no effect of discrimination from coworkers, patients, and visitors on anxiety. More research is needed to draw clearer conclusions regarding the existence of a zero or positive effect of these latter three variables.




 §.§ Hypothesis testing of multiple effects in social psychology

In our second example we illustrate how our testing procedure can be used when testing multiple hypotheses with competing equal and order constraints on the effect of different predictor variables. Carlson and Sinclair (2017) compared two different theoretical explanations for perceptions of gender discrimination in hiring for the roles of computer specialist and nurse. To test individual differences they regressed perceptions of discrimination towards female victims on belief in discrimination against women, stigma consciousness and feminist identification, while controlling for gender and belief in discrimination against men. As a regression equation this can be expressed as


    y_discriminationW,i    =   β_0 + β_beliefW  X_beliefW,i + β_stigma  X_stigma,i + β_feminist  X_feminist
          + β_gender  X_gender,i + β_beliefM  X_beliefM,i + error_i.


where the β's are standardized regression effects of the variables on perceived discrimination. Since in this subsection we will compare the beta-coefficients of different variables against each other, it facilitates interpretation if they are on the same scale. As such we standardize all variables before entering them in the model.

The two theories that Carlson and Sinclair (2008) examined make different explanations for what individual characteristics are most important to perceptions of gender discrimination. The `prototype explanation' suggests that what matters are the individual's beliefs that the gender in question is discriminated against, whereas the `same-gender bias explanation' suggests that identification with the victim is most important. In our example, the victim of discrimination is female and Carlson and Sinclair operationalize identification with the victim as stigma consciousness and feminist identity. Note that neither theory makes any predictions regarding the control variables (gender and general belief that men are discriminated against). A first hypothesis, based on the prototype explanation, might thus be that belief in discrimination of women in general is positively associated with the belief that the female applicant has been discriminated against, whereas stigma consciousness and feminist identity have no effect on this belief. Formally, this can be expressed as


    H_1: β_beliefW > β_stigma= β_feminist = 0


which is equivalent to:


    H_1: β_beliefW > (β_stigma, β_feminist) = 0


Alternatively, we might expect all three variables to have a positive effect on the dependent variable (all β's > 0), but that, in accordance with the prototype explanation, a belief that women are generally discriminated against should have a larger effect on perceptions of discrimination than identifying with the job applicant. Formally this implies:


    H_2: β_beliefW > (β_stigma, β_feminist) > 0


A third hypothesis, based on the same-gender bias explanation, would be the reverse of the H_1, namely that stigma consciousness and feminist identity are positively associated with the outcome while a general belief in discrimination against women has no impact on the particular case. That is:


    H_2:  (β_stigma, β_feminist) > β_beliefW = 0


In this example we have thus specified three contradicting hypotheses regarding the relationships between three variables and wish to know which hypothesis receives most support from the data at hand. However, there is one additional implied hypothesis in this case: the complement. The complement, H_c, is the hypothesis that none of the specified hypotheses are true. The complement exists if the specified hypotheses are not exhaustive, that is, do not cover the entire parameter space. In other words, the complement exists if there are possible values for the regression coefficients which are not contained in the hypotheses, for example, (β_1, β_2, β_3) = (-1, -1, -1) is a combination of effects which do not satisfy the constraints of either H_1, H_2, or H_3. Thus, the interest is in testing the following hypotheses:


    H_1   :   β_beliefW > (β_stigma, β_feminist) = 0 
    
    H_2   :   β_beliefW > (β_stigma, β_feminist) > 0  
    
    H_3   :     (β_stigma, β_feminist) > β_beliefW = 0 
    
    H_c   :   not H_1, H_2, H_3


As before, we begin by fitting a linear regression on the (standardized) variables:



Next, we specify the hypotheses separated by semicolons as a character vector, here on separate lines for space reasons:



The complement does not need to be specified, as the function will include it automatically if necessary. For this example we get the following output:



From the output posterior probabilities we see that H_1 and H_2, both based on the prototype explanation, received the most support, whereas H_3, which was derived from the same-gender bias model, and the complementary hypothesis are both highly unlikely. These results can be succinctly reported as: "Using a default Bayes factor approach, we obtain overwhelming evidence that either hypothesis H_1 or H_2 is true with posterior probabilities of approximately .637, .359, .000,  and .004 for H_1, H_2, H_3, and H_4, respectively." Printing the Bayes Factor matrix yields:



We see that the evidence for both H_1 and H_2 is very strong compared to the complement and in particular compared to H_3, but that H_1 is only 1.8 times as likely as H_2 (B_13 = 1.777).

To summarize the second application, our data demonstrated strong evidence for the prototype explanation and a lack of support for the same-gender bias explanation in explaining perceptions of discrimination against female applicants in the hiring process of computer specialist and nurses. The relative evidence for the prototype explanation depended on its exact formulation, but was at least 919 times stronger than for the same-gender bias explanation, and 91 times stronger than for the complement. However, further research is required to determine whether identification with a female victim has zero or a positive effect on perceived discrimination.



 §.§ Supplementary output

When saving results from the  function to an object it is possible to print additional supplementary output. This output is provided to support a deeper understanding of the method and the primary output outlined in the above subsections. We illustrate these two additional commands using the example in section 4.2. Calling  prints the measures of relative fit “f” and complexity “c” in (<ref>)-(<ref>) of the Bayes factor of each hypothesis against the unconstrained hypothesis. Thus, for the data and hypotheses of section 4.2 we get



where  is the prior density at the null value,  the prior probability that the constraints hold,  the product of these two, and the columns labeled as , , and  have similar interpretations for the posterior quantities.  is the Bayes factor of hypothesis H_t against the unconstrained (H_u) and  is the posterior probability of hypothesis H_t. We rounded the output to three decimals for convenience. Cells with “NA” indicate that a column is “Not Available” to a particular hypothesis. For example, because H_2 contains only inequality comparisons  it has a prior (and posterior) probability but no prior density evaluated at a null value. Hypothesis H_1 and H_3 contain both equality and inequality comparisons and thus has both prior and posterior densities and probabilities. The Bayes factor for H_1 and H_3 against H_u can thus be calculated as B_1u=4.398/0.075=58.265 and B_3u=0.002/0.055=0.036 (see column ). The posterior hypothesis probabilities are calculated using (<ref>) by setting equal prior probabilities, i.e., Pr(H_t|y)=B_tu/∑_t' B_t'u, yielding, for example, Pr(H_1|y) = 58.265/58.265+32.525+0.036+0.357 = 0.639 (as indicated in column ).

If R_I is not of full row rank, the posterior and prior that the inequality constraints hold are computed as the proportion of draws from unconstrained Student t distributions. Under these circumstances there will be a, typically small, numerical Monte Carlo error. The 90% credibility intervals of the numerical estimate of the Bayes factors of the hypotheses against the unconstrained hypothesis can be obtained by calling



where  is the Bayes factor of hypothesis t against the unconstrained (u),  is the lower bound of the 90% credibility interval estimate of the Bayes factor and  is the upper bound. Credibility intervals are only printed when the computed Bayes factors have numerical errors. If the user finds the Monte Carlo error to be too large they can increase the number of draws from the Student t distributions by adjusting the input value for the  argument (default 10^6 draws).



§ DISCUSSION

The paper presented a new Bayes factor test for evaluating hypotheses on the relative effects in a linear regression model. The proposed testing procedure has several useful properties such as its flexibility to test multiple equality and/or order constrained hypotheses directly against each other, its intuitive interpretation as a measure of the relative evidence in the data between the hypotheses, and its fast computation. Moreover no prior information needs to be manually specified about the expected magnitude of the effects before observing the data. Instead, a default procedure is employed where a minimal fraction of the data is used for default prior specification and the remaining fraction is used for hypothesis testing. A consequence of this choice is that the statistical evidence cannot be updated using Bayes' theorem when observing new data. This is common in default Bayes factors <cit.>. Instead, the statistical evidence needs to be recomputed when new data are observed. This however is not a practical problem because of the fast computation of the default Bayes factor due to its analytic expression.

Furthermore the readily available lmhyp-package can easily be used in combination with the popular lm-package for linear regression analysis. The new method will allow researchers to perform default Bayesian exploratory analyses about the presence of a positive, negative or zero effect and to perform default Bayesian confirmatory analyses where specific relationships are expected between the regression effects which can be translated to equality and order constraints. The proposed test will therefore be a valuable contribution to the existing literature on Bayes factor tests <cit.>, which are gradually winning ground as alternatives to classical significance tests in social and behavioral research. Due to this increasing literature, a thorough study about the qualitative and quantitative differences between these Bayes factors is called for. Another useful direction for further research would be to derive Bayesian (interval) estimates under the hypothesis that receives convincing evidence from the data.





§ ACKNOWLEDGEMENTS

We would like to thank Stephen Wood for allowing us to use the empirical data in the first empirical application in this paper. These data are from the Service Delivery and Organisation (SDO) National Inpatient Staff Morale Study, which was funded by the National Institute for Health Research Service Delivery and Organisation (NIHR SDO) programme (Project Number 08/1604/142) and of which Stephen Wood was a co-investigator. The views and opinions expressed in the article are those of the authors and do not necessarily reflect those of the NIHR SDO programme or the Department of Health. The first author was supported by a NWO Vidi grant (452-17-006).






§ PROOF OF LEMMA 1

A derivation is given for the prior adjusted default Bayes factor for a hypothesis H_1:R_Eβ=r_E, R_Iβ>r_I against an unconstrained hypothesis H_u:β∈ℝ^k in (<ref>)-(<ref>). Based on the reparameterization ξ=[[ ξ_E; ξ_I ]]=[[ R_E;   D ]]β in (<ref>), the hypothesis is equivalent to H_1:ξ_E=r_E, R̃_Iξ_I>r̃_I against an unconstrained hypothesis H_u:ξ∈ℝ^k. The marginal likelihood under the constrained hypothesis H_1 is defined as in the fractional Bayes factor <cit.> with the exception that we integrate over an adjusted integration region <cit.>. This adjustment ensures that the implicit default prior is centered on the boundary of the constrained space. The marginal likelihood under H_1 is defined by

    p_1(y,b) = ∬_R̃_Iξ_I>r̃_I p(y|ξ_E=r_E,ξ_I,σ^2)π_u^N(ξ_I,σ^2)dξ_I dσ^2/∬_R̃_Iξ_I>r^*_I p(y|ξ_E=ξ̂_E,ξ_I,σ^2)^bπ_u^N(ξ_I,σ^2)dξ_I dσ^2.

As can be seen the adjustment implies that in the denominator the fraction of the likelihood is evaluated at ξ̂_E instead of r_E and the integration region equals

    R̃_I(ξ_I-ξ̂_I+μ_I^0)>r̃_I
    ⇔R̃_Iξ_I>R̃_Iξ̂_I=r_I^*,

because R̃_Iμ_I^0=r̃_I, instead of R̃_Iξ_I>r̃_I. Note that ξ̂_I=Dβ̂. This adjustment of the fractional Bayes factor ensures that the proposed default Bayes factor is computed using an implicit default prior that is centered on the boundary of the constrained space, following Jeffreys' heuristic argument <cit.>. Furthermore this ensures that the complexity of an order-constrained hypothesis is properly incorporated in the Bayes factor <cit.>.
The marginal likelihood under H_u is defined by

    p_u(y,b) = ∭ p(y|ξ_E,ξ_I,σ^2)π_u^N(ξ_E,ξ_I,σ^2)dξ_E dξ_I dσ^2/∭ p(y|ξ_E,ξ_I,σ^2)^bπ_u^N(ξ_E,ξ_I,σ^2)dξ_E dξ_I dσ^2.

The fraction b will be set to k+1/n because k+1 observations are needs to obtain a finite marginal likelihood when using a noninformative prior π^N_u(ξ_E,ξ_I,σ^2)=σ^-2 under H_u.

The default Bayes factor is then given by

    B_1u,b   =   p_1(y,b)/p_u(y,b)
    = ∬_R̃_Iξ_I>r̃_I p(y|ξ_E=r_E,ξ_I,σ^2)π_u^N(ξ_I,σ^2)dξ_I dσ^2/∬_R̃_Iξ_I>r^*_I p(y|ξ_E=ξ̂_E,ξ_I,σ^2)^bπ_u^N(ξ_I,σ^2)dξ_I dσ^2/
          ∭ p(y|ξ_E,ξ_I,σ^2)π_u^N(ξ_E,ξ_I,σ^2)dξ_E dξ_I dσ^2/∭ p(y|ξ_E,ξ_I,σ^2)^bπ_u^N(ξ_E,ξ_I,σ^2)dξ_E dξ_I dσ^2
       =   ∬_R̃_Iξ_I>r̃_I p(y|ξ_E=r_E,ξ_I,σ^2)π_u^N(ξ_I,σ^2)/∭ p(y|ξ_E,ξ_I,σ^2)π_u^N(ξ_E,ξ_I,σ^2)dξ_E dξ_I dσ^2dξ_I dσ^2
    /
          ∬_R̃_Iξ_I>r^*_I p(y|ξ_E=ξ̂_E,ξ_I,σ^2)^bπ_u^N(ξ_I,σ^2)/∭ p(y|ξ_E,ξ_I,σ^2)^bπ_u^N(ξ_E,ξ_I,σ^2)dξ_E dξ_I dσ^2
    dξ_I dσ^2
       =   ∬_R̃_Iξ_I>r̃_Iπ_u(ξ_E=r_E,ξ_I,σ^2|y)
    dξ_I dσ^2
    /
          ∬_R̃_Iξ_I>r^*_Iπ_u(ξ_E=ξ̂_E,ξ_I,σ^2|y^b)
    dξ_I dσ^2
       =   Pr(R̃_Iξ_I>r̃_I|y,ξ_E=r_E)/Pr(R̃_Iξ_I>r^*_I|y,ξ_E=ξ̂_E)×π_u(ξ_E=r_E|y)/π_u(ξ_E=ξ̂_E|y^b).

Furthermore, using standard calculus it can be shown that the marginal posterior for β for a fraction b of the data and a noninformative prior has a Student t distribution

    π_u(β|y^b)    ∝   ∫ p(y|β,σ^2)^b π_u^N(β,σ^2) dσ^2
       ∝    t(β;β̂,s^2(nb-k)^-1(X'X)^-1,nb-k),

and therefore, because ξ=Tβ|y^b, it holds that

    π_u(ξ|y^b) =
    t(ξ;ξ̂,s^2(nb-k)^-1T(X'X)^-1T',nb-k),

where T=[[ R_E;   D ]], ξ̂=(ξ̂_E',ξ̂_I')' with ξ̂_E=R_Eβ̂ and ξ̂_I=Dβ̂, and t(ξ;m,K,ν) denotes a Student t distribution for ξ with location parameters m, scale matrix K, and ν degrees of freedom.
Then it is well-known that the marginal distribution of ξ_E and the conditional distribution of ξ_I|ξ_E also have Student t distributions <cit.> given by

    ξ_E|y^b    ∼   
    t(ξ̂_E,s^2(nb-k)^-1R_E(X'X)^-1R'_E,nb-k)
    ξ_I|ξ_E,y^b    ∼   
    t(m_I|E,K_I|E,ν),

with

    m_I|E   =   Dβ̂+s^-2(nb-k)D(X'X)^-1R_E'(R_E(XX)^-1R_E)^-1(ξ_E-ξ̂_E)
    K_I|E   =   nb-k + s^-2(nb-k)(ξ_E-ξ̂_E)'(R_E'(X'X)^-1R_E)^-1(ξ_E-ξ̂_E)/nb-k+q^E
          (s^2(nb-k)^-1D(X'X)^-1D'-
           s^2(nb-k)^-1D(X'X)^-1R_E'(R_E(X'X)^-1R_E')^-1R_E(X'X)^-1D).


Thus, when plugging in b=1 and ξ_E=r_E in (<ref>) and (<ref>), and then in (<ref>), gives the numerators in (<ref>) and (<ref>), and plugging in b=k+1/n and ξ_E=ξ̂_E in (<ref>) and (<ref>), and then in (<ref>), gives the denominators in (<ref>) and (<ref>), which completes the proof.




§ EXAMPLE ANALYSIS FOR FIGURE 1




apacite



