Null hypothesis testing: does x1 exclude n?

Modeling
1

I have a dataset in which I’m interested in whether responses in condition b are more accurate than in condition a.

I’m also interested in whether they are higher than chance, which is not 50%/logit(.5) but rather n where n is a probability. So I want to find out my certainty about an across-condition difference and also whether a condition excludes n.

I’ll illustrate this working with the iris dataset in R and a linear model. Here, n can be any number and there is no transformation involved but the problem remains the same.

I have two questions:

Q1: Does sepal length differ significantly between setosa (the intercept) and versicolor?
Q2: Does the range of sepal length for species: versicolor exclude n=7?

For Q1, I can fit

lm1 = lm(Sepal.Length ~ 1 + Species, data = iris)

and then get confidence intervals for the predictor speciesversicolor: [0.73;1.13]. So yes, there is a significant difference here with p<0.05.

For Q2, I can fit

lm2 = lm(Sepal.Length ~ 0 + Species, data = iris)

and get confidence intervals for each level of species and see whether versicolor excludes 7. It does ([5.79;6.08]). I could also relevel the factor and refit the model but let’s say I want this for each level of species.

I’m new to Bayesian modelling. I feel like the answer to Q2 is available in the model:

stan_lm1 = stan_glm(Sepal.Length ~ 1 + Species, data = iris)

I mean there is no significance level but I can get a posterior interval for Species: versicolor and see if it includes 7. Not sure how to do that though. The posterior draws are about the difference from the intercept. brms::hypothesis can test ‘Speciesversicolor < 0’ but not, it seems, ‘Speciesversicolor < 7’.

I can of course fit

stan_lm2 = stan_glm(Sepal.Length ~ 0 + Species, data = d)

and use that one to test Q1 as well, so I can have 1 model instead of 2, but it feels like I’m missing something.

2

I can use tidybayes::add_epred_draws to get draws from the expectation of the posterior predictive distribution to get the average expectation for each level (see here).

I can then use these to visualise expected distributions:


draws = iris %>% 
  distinct(Species) %>% 
  add_epred_draws(stan_lm1, scale = 'response', re_formula = NA, ndraws = 500)

draws %>% 
  ggplot(aes(x = Species, y = .epred)) +
  geom_violin() +
  geom_boxplot(width = .1)
3

Both of your models contain the same information and can be used to test both Q1 and Q2. The way that the factors are coded changes the meaning of the intercept term, and you’ll need to express the hypotheses slightly differently to accomodate the different meanings between you two models. I also recommend fitting these with brms::brm instead of rstanarm::stan_glm because I’m not sure how to include the intercept term from an rstanarm fit in hypothesis.

library(brms)

stan_lm1 <- brm(Sepal.Length ~ 1 + Species, data = iris)
stan_lm2 <- brm(Sepal.Length ~ 0 + Species, data = iris)

# Q1
hypothesis(stan_lm1, "Speciesversicolor = 0")
hypothesis(stan_lm2, "Speciesversicolor - Speciessetosa = 0")

# Q2
hypothesis(stan_lm1, "Intercept + Speciesversicolor < 7")
hypothesis(stan_lm2, "Speciesversicolor < 7")

You should see that these return the exact same estimates and credible intervals for the hypothesis test. There’s also a plot() method for brmshypothesis objects if you want to visualize the posterior distribution of the value calculated by your hypothesis.

image
image1443×924 38.9 KB

4

A note to myself (and perhaps others): I was generally confused by how brms::hypothesis is printed.

See in the example above:

hyp1 = hypothesis(stan_lm1, “Intercept + Speciesversicolor < 7”)
hyp1

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (Intercept+Specie... < 0    -1.07      0.07    -1.18    -0.94        Inf         1    *

but

hyp1$hypothesis

                            Hypothesis  Estimate  Est.Error  CI.Lower   CI.Upper Evid.Ratio
1 (Intercept+Speciesversicolor)-(7) < 0 -1.065003 0.07419522 -1.184565 -0.9423582        Inf
  Post.Prob Star
1         1    *

so x < 7 is formatted as x - 7 < 0 and then some of it might get clipped.

5

Yeah the automatic algebraically re-arrangement can be a bit jarring. Intercept + Speciesversicolor < 7 becomes (Intercept + Speciesversicolor) - (7) < 0. At least the plot method displays the transformed hypothesis equations in the facet label.

You can avoid this behavior by expressing your hypothesis with 0 on the righthand side of the (in)equality of the hypothesis formula. This flexibility to re-arrange terms is also helpful if it’s more natural to interpret your hypothesis Estimate and posterior distribution with the sign flipped, e.g., (7 - (Intercept + Speciesversicolor)) > 0.