`library(alphaN)`

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.

You can install the development version of alphaN from GitHub with:

```
# install.packages("devtools")
::install_github("jespernwulff/alphaN") devtools
```

This vignette provides an introduction to the basic functionality of
**alphaN**. For full details on methodology, please refer
to Wulff & Taylor (2023).

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.

```
<- alphaN(n = 1000, BF = 1)
alpha
alpha#> [1] 0.008582267
```

Therefore, to obtain evidence of at least 1, we should set our alpha to 0.0086.

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.

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:

```
<- seq(50, 1000, 1)
seqN plot(seqN, alphaN(seqN), type = "l",
xlab = "n", ylab = "Alpha")
```

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:

`alphaN_plot(BF = 3)`

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.

`<- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv") df `

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`

:

```
<- alphaN(n = 400, BF = 1, method = "JAB")
alpha_gre
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.

```
<- glm(admit ~ gpa + factor(rank) + gre, data = df, family = "binomial")
glm1 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.

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(glm1, covariate = "gre", method = "JAB")
JAB_gre
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.