We wish to determine which alpha level is equivalent to a Bayes
factor of 1. I.e. only reject the null if the data is at least at likely
under the null and under the alternative. To do this, we need a way to
connect the \(p\)-value to the Bayes
factor. The alphaN package does this for tests of
coefficients in regression models.
Basic functionality
This vignette provides an introduction to the basic functionality of
alphaN. For full details on methodology, please refer
to Wulff & Taylor (2023).
Setting the alpha level
Using the alphaN function, we can get the alpha level we
need to use to obtain a desired level of evidence when testing a
regression coefficient in regression model.
Here is an example: We are planning to run a linear regression model
with 1000 observations. We thus set n = 1000. The default
BF is 1 meaning that we want to avoid Lindley’s paradox,
i.e. we just want the null and the alternative to be at least equally
likely when we reject the null.
alpha <- alphaN(n = 1000, BF = 1)
alpha
#> [1] 0.008582267
Therefore, to obtain evidence of at least 1, we should set our alpha
to 0.0086.
Plotting the relationship between the Bayes factor and
p-value
The alphaN function works by mapping the
p-value to the Bayes factor. This relationship can be shown
using the JAB_plot. For instance:
JAB_plot(n = 1000, BF = 1)
The alpha level needed to achieve a Bayes factor of 1 is shown with a
red triangle in the plot. Lines for achieving Bayes factors of 3
(moderate evidence) and 10 (strong evidence) are also shown by default.
As it is evident a lower alpha level is needed to achieve higher
evidence.
Alpha as a decreasing function of N
An important point of the procedure is that alpha will be set as a
function of sample size. The larger the sample size, the lower the alpha
needed such that a significant result can be interpreted as evidence for
the alternative.
The graph below illustrates this relationship for previous
example:
seqN <- seq(50, 1000, 1)
plot(seqN, alphaN(seqN), type = "l",
xlab = "n", ylab = "Alpha")
Setting the prior
To set the alpha level as a function of sample size, we need to
choose the prior carefully. alphaN allows the user to
choose from four sensible prior options based on suggestions from the
previous literature: Jeffreys’ approximate BF
(method = "JAB"), the minimal training sample
(method = "min"), the robust minimal training sample
(method = "robust"), and balanced Type-I and Type-II errors
(method = "balanced"). method = "JAB" is a
good choice for users who want to be conservative against small effects,
method = "min" is for when the MLE is misspecified,
method = "robust" is for when the MLE is misspecified and
the sample size is small, and method = "balanced" is for
when Type-II errors are costly.
For instance, to achieve evidence of 3 for 1,000 observations while
we ensure balanced error rates, we run
alphaN(1000, BF = 3, method = "balanced")
#> [1] 0.024221
The package contains the convenience function
alphaN_plot that allows a quick comparison of alpha as a
function of sample size for the four different methods:
Example
In this section, we illustrate how the package may be used on a
dataset. In this case, we use a dataset on getting into graduate school
from UCLA.
df <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
The dataset contains four variables with 400 observations. The
variables are Graduate Record Exam scores (gre), grade
point average (gpa) and the rank the undergraduate
institution (`rank``).
str(df)
#> 'data.frame': 400 obs. of 4 variables:
#> $ admit: int 0 1 1 1 0 1 1 0 1 0 ...
#> $ gre : int 380 660 800 640 520 760 560 400 540 700 ...
#> $ gpa : num 3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ...
#> $ rank : int 3 3 1 4 4 2 1 2 3 2 ...
We imagine we are interested in testing the coefficient on
gre and we want to estimate the model
admit ~ gre + gpa + rank where we are interested in testing
the coefficient on hp. We set n = 400 because
we have 400 observations.
Let us also say that we would like it to be just as likely that the
alternative is true compared to the null if we reject the null. This
means that we want to know which alpha corresponds to a Bayes factor of
1 (If we instead would want it to be 3 times more likely that the
alternative is true than the null if we reject the null, we would find
the alpha corresponding to a Bayes factor of 3). Thus, we set
BF = 1. Because we wish to remain skeptical of trivial
effects, we use the default method = "JAB:
alpha_gre <- alphaN(n = 400, BF = 1, method = "JAB")
alpha_gre
#> [1] 0.01437526
The p-value that corresponds to a Bayes factor of 1 for this
particular model and sample size is 0.0144. We therefore set alpha to
0.0144 and estimate our model.
glm1 <- glm(admit ~ gpa + factor(rank) + gre, data = df, family = "binomial")
summary(glm1)
#>
#> Call:
#> glm(formula = admit ~ gpa + factor(rank) + gre, family = "binomial",
#> data = df)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -1.6268 -0.8662 -0.6388 1.1490 2.0790
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -3.989979 1.139951 -3.500 0.000465 ***
#> gpa 0.804038 0.331819 2.423 0.015388 *
#> factor(rank)2 -0.675443 0.316490 -2.134 0.032829 *
#> factor(rank)3 -1.340204 0.345306 -3.881 0.000104 ***
#> factor(rank)4 -1.551464 0.417832 -3.713 0.000205 ***
#> gre 0.002264 0.001094 2.070 0.038465 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 499.98 on 399 degrees of freedom
#> Residual deviance: 458.52 on 394 degrees of freedom
#> AIC: 470.52
#>
#> Number of Fisher Scoring iterations: 4
The p-value for the coefficient on gre is about
0.0385. Because this is larger than 0.0144, we cannot reject the null of
no relationship and conclude that the null is more likely to be true
conditional on this data.
Quantifying the evidence based on the data
Next, we can compute the actual Bayes factor for gre. We
can do this using the JAB function. It takes as an argument
the glm object. We specify that we are interest in the
gre variable and set method = JAB adj to
adjust for the number of parameters. The function automatically counts
the number of parameters based on the glm object:
JAB_gre <- JAB(glm1, covariate = "gre", method = "JAB")
JAB_gre
#> [1] 0.425894
We can see that the Bayes factor is 0.4259, which indeed does
indicate that it is more likely that the null is true. The Bayes factor
directly quantifies the evidence and suggests that it is 2.347969 times
more likely that the null is true compared to the compared, which is
just anecdotal evidence.
We could also have computed the Bayes factor manually using the
JABt function by plugging in the sample size and the
z-statistic from the regression:
JABt(400, 2.070, method = "JAB")
#> [1] 0.4260143
or by plugging in the \(p\)-value in
the JABp function while making sure to tell the function
that the \(p\)-value is based on a
z-test:
JABp(400, 0.038465, z = TRUE, method = "JAB")
#> [1] 0.4258952
The difference to the result from JAB is solely due to
rounding errors because JAB uses the exact values from the
glm object instead of the rounded values that we supplied
to the functions. The functions JABt and JABp
are useful in situations where the dataset may not avaliable, for
instance when for results printed in a journal article.