As a general rule we do four analyses fairly commonly in the Bridges Lab, summarized here based on the nature of the independent and dependent variables:
Dependent Variable
Independent Variable
Analysis
Continuous
Continuous
Linear Regression
Continuous
Counts (Yes/No or Groups)
Binomial Regression
Counts (Yes/No or Groups)
Continuous
Pairwise Test (t-test/Mann-Whitney/ANOVA)
Counts (Yes/No or Groups)
Counts (Yes/No or Groups)
\(\chi^2\) or Fisher’s Test
Lets take a Bayesian approach to each of these using the Mtcars dataset, which looks like this after a bit of fiddling.
Lets start by testing if there is a relationship between one quantitative variables, the mpg (miles per gallon) and a categorical variable - transmission (manual or automatic transmission) using the default priors
library(broom.mixed)tidy(pairwise.fit) %>%kable(caption="Summary of model fit for mpg versus transmission")
Summary of model fit for mpg versus transmission
effect
component
group
term
estimate
std.error
conf.low
conf.high
fixed
cond
NA
(Intercept)
17.089490
1.1440067
14.874103
19.262258
fixed
cond
NA
transmissionmanual
7.277792
1.8424683
3.713655
10.912369
ran_pars
cond
Residual
sd__Observation
5.010808
0.6489197
3.932789
6.461536
ran_pars
cond
Residual
prior_sigma__NA.NA.prior_sigma
5.893894
6.9432256
0.161323
23.833297
plot(pairwise.fit)
hypothesis(pairwise.fit, "transmissionmanual>0") # testing for whether manual tramsmission has higher mpg
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (transmissionmanual) > 0 7.28 1.84 4.23 10.32 Inf
Post.Prob Star
1 1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
As you can see, this analysis estimates that the manual transmission has a 7.2777915 higher mpg (\(\pm\) 1.8424683) with a very high Bayes Factor () and posterior probability (1).
Looking at the model diagnostics three things are relevant:
Was there convergence in the model. This is measured by Rhat (also known as potential scale reduction factor (PSRF), or the Gelman-Rubin statistic). It should be between 0.99 and 1.01.
What was the ESS (effective sample size, not a great name, but a measure of the number of samples of the posterior distribution). Generally, an \(ESS > 400\) is considered good for reliable estimates. This is most important for the random effects rather than the fixed effects.
Lets take a look at these:
kable(data.frame(Parameter =names(rhat(pairwise.fit)),Rhat =format(rhat(pairwise.fit), nsmall =5)),caption="Rhat values for model testing association between mpg and transmission (should b between 0.99 and 1.01 for convergence",row.names = F)
Rhat values for model testing association between mpg and transmission (should b between 0.99 and 1.01 for convergence
Parameter
Rhat
b_Intercept
1.000238
b_transmissionmanual
1.000444
sigma
1.000516
Intercept
1.000369
prior_Intercept
1.000925
prior_sigma
1.000403
lprior
1.000087
lp__
1.001772
tidy(pairwise.fit, effects ="ran_pars", ess =TRUE) %>%kable(caption="Model random effects for mpg vs transmission")
Model random effects for mpg vs transmission
effect
component
group
term
estimate
std.error
conf.low
conf.high
ess
ran_pars
cond
Residual
sd__Observation
5.010808
0.6489197
3.932789
6.461536
3072.006
ran_pars
cond
Residual
prior_sigma__NA.NA.prior_sigma
5.893894
6.9432256
0.161323
23.833297
3859.087
The posterior distribution for the effect of a manual transmission is here:
We didnt specify the model type (gausian is the default). We also used the default priors here, which were
prior_summary(pairwise.fit) %>%kable(caption="Default priors for a pairwise analysis of mpg vs transmission in the mtcars data")
Default priors for a pairwise analysis of mpg vs transmission in the mtcars data
prior
class
coef
group
resp
dpar
nlpar
lb
ub
source
b
default
b
transmissionmanual
default
student_t(3, 19.2, 5.4)
Intercept
default
student_t(3, 0, 5.4)
sigma
0
default
This includes flat priors for the beta coefficients, student’s t distributions for intercept (centered around the mean mpg) and residual error (centered around zero).
What if I Want an ANOVA
Perhaps I only am interested in whether there is a difference between groups, but not pairwise differences (i.e an ANOVA). This can be done as well (again using default priors).
anova.fit <-brm(mpg~0+cylindars,data=mtcars.data, sample_prior =TRUE) #zero sets there to be no intercept, shows the mean for each groupanova.fit.null <-brm(mpg~0,data=mtcars.data, sample_prior =TRUE) #null model#estimate the bayes factorbayes_factor(anova.fit,anova.fit.null) -> anova.bfpost_prob(anova.fit,anova.fit.null) -> anova.pp
Now lets look at the results of those analyses comparasons. The bayes factor for the ANOVA is 3.4898088^{28} which is very extreme evidence for a difference between groups. The posterior probabilities are 1 (very hight) for the fitted model and 2.8829788^{-29} (very low) for the null model.
prior_summary(anova.fit) %>%kable(caption="Summary of priors for comparason between cylindar numbers on mpg")
Summary of priors for comparason between cylindar numbers on mpg
prior
class
coef
group
resp
dpar
nlpar
lb
ub
source
b
default
b
cylindars4
default
b
cylindars6
default
b
cylindars8
default
student_t(3, 0, 5.4)
sigma
0
default
#model fitting parameterskable(data.frame(Parameter =names(rhat(anova.fit)),Rhat =format(rhat(anova.fit), nsmall =5)),caption="Rhat values for model testing association between mpg and cylindars (should b between 0.99 and 1.01 for convergence",row.names = F)
Rhat values for model testing association between mpg and cylindars (should b between 0.99 and 1.01 for convergence
Parameter
Rhat
b_cylindars4
0.9998450
b_cylindars6
1.0016103
b_cylindars8
1.0012933
sigma
0.9997617
prior_sigma
0.9997447
lprior
0.9997617
lp__
1.0007648
tidy(anova.fit, effects ="ran_pars", ess =TRUE) %>%kable(caption="Model random effects for mpg vs cylindars")
Model random effects for mpg vs cylindars
effect
component
group
term
estimate
std.error
conf.low
conf.high
ess
ran_pars
cond
Residual
sd__Observation
3.332159
0.460378
2.5931131
4.377553
3015.966
ran_pars
cond
Residual
prior_sigma__NA.NA.prior_sigma
5.905686
7.034321
0.1814052
22.768172
3510.893
#create a teable of the pairwise hypothesesrbind(hypothesis(anova.fit, "cylindars4>cylindars6")$hypothesis,hypothesis(anova.fit, "cylindars4>cylindars8")$hypothesis,hypothesis(anova.fit, "cylindars6>cylindars8")$hypothesis) %>%kable(caption="Pairwise hypothesis tests for cylindars on mpg")
Pairwise hypothesis tests for cylindars on mpg
Hypothesis
Estimate
Est.Error
CI.Lower
CI.Upper
Evid.Ratio
Post.Prob
Star
(cylindars4)-(cylindars6) > 0
6.925078
1.592501
4.362640
9.505415
Inf
1.00000
*
(cylindars4)-(cylindars8) > 0
11.540329
1.367975
9.291548
13.788081
Inf
1.00000
*
(cylindars6)-(cylindars8) > 0
4.615250
1.517746
2.095993
7.108223
443.4444
0.99775
*
tidy(anova.fit) %>%kable(caption="Summary of model fit for mpg versus cylindars")
Lets start by testing if there is a relationship between two quantitative variables, the mpg (miles per gallon) and the weight (wt) using the default priors
prior_summary(linear.fit) %>%kable(caption="Prior summary for effects of weight on mpg")
Prior summary for effects of weight on mpg
prior
class
coef
group
resp
dpar
nlpar
lb
ub
source
b
default
b
wt
default
student_t(3, 19.2, 5.4)
Intercept
default
student_t(3, 0, 5.4)
sigma
0
default
tidy(linear.fit) %>%kable(caption="Summary of model fit for mpg versus weight")
Summary of model fit for mpg versus weight
effect
component
group
term
estimate
std.error
conf.low
conf.high
fixed
cond
NA
(Intercept)
37.281958
1.9944452
33.329685
41.141723
fixed
cond
NA
wt
-5.347938
0.5933914
-6.509334
-4.197569
ran_pars
cond
Residual
sd__Observation
3.146442
0.4190495
2.454004
4.105158
ran_pars
cond
Residual
prior_sigma__NA.NA.prior_sigma
5.951883
6.5251426
0.193961
22.365550
plot(linear.fit)
hypothesis(linear.fit, "wt<0") # testing for whether weight has a negative effect on mpg
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (wt) < 0 -5.35 0.59 -6.32 -4.36 Inf 1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
Notice that the priors were the same here (they are based on the dependent variable), and again there is a highly probable inverse relationship between mpg and weight (heavier cars get lower fuel economy).
What if I am more interested in the R-squared?
To do this we can use the bayes_R2 function.
kable(bayes_R2(linear.fit),caption="Estimates for R2 between weight and mpg")
Estimates for R2 between weight and mpg
Estimate
Est.Error
Q2.5
Q97.5
R2
0.7416091
0.0472108
0.6223189
0.7978628
r2.probs <-bayes_R2(linear.fit, summary=F) #summary false is to get the posterior probabilitiesggplot(data=r2.probs, aes(x=R2)) +geom_density(fill="blue") +geom_vline(xintercept=0,color="red",lty=2) +labs(y="",x="R2",title="Posterior probabilities") +lims(x=c(0,1)) +theme_classic(base_size=16)
What is a Chi-Squared Equivalent
First lets do an example with a standard \(\chi^2\) test (and a Fisher’s test since the counts are quite low) on whether there is a relationship bewteen engine type and transmission type.
library(tidyr) #for pivot widerlibrary(tibble) #for column to rownameengine.trans.counts <- mtcars.data %>%group_by(engine,transmission) %>%count() %>%pivot_wider(names_from=transmission,values_from=n) %>%column_to_rownames('engine')chisq.test(engine.trans.counts) %>% tidy %>%kable(caption="Chi-squared test for engine/transmission")
Chi-squared test for engine/transmission
statistic
p.value
parameter
method
0.3475355
0.5555115
1
Pearson’s Chi-squared test with Yates’ continuity correction
fisher.test(engine.trans.counts) %>% tidy %>%kable(caption="Fisher's test for engine/transmission")
Fisher’s test for engine/transmission
estimate
p.value
conf.low
conf.high
method
alternative
0.511233
0.4726974
0.0944144
2.614145
Fisher’s Exact Test for Count Data
two.sided
Both agree, not much of a relationship here. For the brms modelling we need to make a few changes. First, we will use a bernoulli distribution since our data is only zeros and ones, again we will use all the default priors.
# Fit the modelcounts.model <-brm(engine ~ transmission, data = mtcars.data, family =bernoulli())
Lets look at these results
prior_summary(counts.model) %>%kable(caption="Prior summary for effects of transmission on engine type")
Prior summary for effects of transmission on engine type
prior
class
coef
group
resp
dpar
nlpar
lb
ub
source
b
default
b
transmissionmanual
default
student_t(3, 0, 2.5)
Intercept
default
tidy(counts.model) %>%kable(caption="Summary of model fit for transmission versus engine type")
Summary of model fit for transmission versus engine type
effect
component
group
term
estimate
std.error
conf.low
conf.high
fixed
cond
NA
(Intercept)
0.5564968
0.4787228
-0.3322844
1.5439365
fixed
cond
NA
transmissionmanual
-0.7260799
0.7635121
-2.2286394
0.7678828
plot(counts.model)
hypothesis(counts.model, "transmissionmanual<0") # testing for whether weight has a negative effect on mpg
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (transmissionmanual) < 0 -0.73 0.76 -1.97 0.52 4.87
Post.Prob Star
1 0.83
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
Now in this case there is moderate evidence for a negative relationship between transmnission and engine type (\(\beta\)=-0.7260799 \(\pm\) 0.7635121) with a Bayes Factor of 4.8651026 and a Posterior Probability of 0.8295.
What is a Binomial Regression Equivalent?
Lets now look at the relationsbip between transmission type (binomial variable) and weight (continuous variable). Again we will use a bernouli distribution (which will use, by default a logit link function)
# Fit the modelbinomial.fit <-brm(transmission ~ wt, data = mtcars.data, family =bernoulli())
prior_summary(binomial.fit) %>%kable(caption="Prior summary for effects of transmission on engine type")
Prior summary for effects of transmission on engine type
prior
class
coef
group
resp
dpar
nlpar
lb
ub
source
b
default
b
wt
default
student_t(3, 0, 2.5)
Intercept
default
tidy(binomial.fit) %>%kable(caption="Summary of model fit for transmission versus engine type")
Summary of model fit for transmission versus engine type
effect
component
group
term
estimate
std.error
conf.low
conf.high
fixed
cond
NA
(Intercept)
14.516094
5.252784
6.173255
26.369970
fixed
cond
NA
wt
-4.800125
1.668465
-8.550785
-2.122508
plot(binomial.fit)
hypothesis(binomial.fit, "wt<0") # testing for whether weight has a negative effect on mpg
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (wt) < 0 -4.8 1.67 -7.78 -2.47 Inf 1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
as_draws_df(binomial.fit) %>%ggplot(aes(x=b_wt)) +geom_density(fill="blue") +geom_vline(xintercept=0,color="red",lty=2) +labs(y="",x="Beta coefficient for weight on transmission",title="Posterior probabilities") +theme_classic(base_size=16)
Again there is strong evidence that a higher weight makes an automatic transmission more likely.
Hopefully these examples helped, but a couple things since we used defaults.
Think closely about your priors, if you have a good reason to set them to something other than the default you should. In this case the defaults worked pretty well and as you can see are set based on the data that is input. These are generally very weakly informative priors and the modelling could be improved on by setting your own.
Also remember, nowhere in here was there a p-value. We could consider a posterior probability of >0.95 our cutoff for significance if we prefer, but make sure to state this in your methods (along with your choice of priors and the package used.)
Note this script used some examples generated by perplexity.ai and then modified further
