# Creating a prior for a logistic model in brms

**URL:** <https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307>\
**Category:** Modeling\
**Tags:** brms\
**Created:** [May 30, 2024, 6:00pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307 "2024-05-30T18:00:09Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![JulianChristov](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianchristov/32/18657_2.png) [@JulianChristov](https://discourse.mc-stan.org/u/JulianChristov)\
**Post date:** [May 30, 2024, 6:00pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/1 "2024-05-30T18:00:09Z")

</div>

I am a recent graduated from George Mason in statistics. My thesis paper was to create a model that could predict match outcomes by using player performance indicators, specifically transfer value. I have only 80 observations from World Cup knockout-stage matches. The first time I ran my model with loo\_predict I received a lot of warnings, numerous high pareto k values. This error gave me unreliable estimates for predictions, and predictions were skewed in the true outcomes favor. I presume this may have been because I specified the prior incorrectly. I can send you my data if you need and would really appreciate to join you on a zoom call to talk about it. Thank you so much for your help in advance.

1. Would you please help me specifying a prior

a. Do power priors exist in BRMS to weight variables based on my beliefs?  
b. What should I specify. Just the fixed effects. The fixed and random effects or something more?  
c. How should I evaluate the distribution of best fit.  
My data: Random effects seem to represent something in between a bell curve and uniform distribution. Most fixed effects represent a compressed bell curve with extremely heavy tail but one variable is binomial -1 and 1? Would using sensitivity analysis and WAIC be the best way to go to fit data to distrubtutions?

Sincerely,

Julian Christov 312-945-9463 [julianchri@gmail.com](mailto:julianchri@gmail.com)

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 31, 2024, 11:31am UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/2 "2024-05-31T11:31:02Z")

</div>

What is the frequency breakdown of the 80 observations’ outcomes? How many predictors do you have? You will get more benefit from limiting the predictor space to what the sample size will support rather than tinkering with smooth priors. For example do a good deal of unsupervised learning e.g. variable clustering, redundancy analysis, possibly sparse principal component analysis. If your outcome is binary you need 96 observations just to estimate the **intercept** in a logistic model.

---

<div class="post-metadata">

**Author:** ![JulianChristov](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianchristov/32/18657_2.png) [@JulianChristov](https://discourse.mc-stan.org/u/JulianChristov)\
**Post date:** [May 31, 2024, 1:17pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/3 "2024-05-31T13:17:05Z")

</div>

Thank you for the timely response. Below, I provided answers to your questions and included the formula and model, as well as histograms for the predictors. I have minimal experience with unsupervised learning, most of which is PCA. Is it possible you could provide code and further insight into applying unsupervised learning techniques to my Bayesian logistic regression model. Thank you so much.

What is the frequency breakdown of the 80 observations’ outcomes?

54 instances of 1 and 26 instances of 0 (54 wins for the home team and 26 losses for the home team).

How many predictors do you have?

11 predictors, 6 fixed and 5 random. Just to clarify I used “0” as the intercept for my model because my variables show the difference in market value and age from one team and the other team.

My model looks at the differences between the two competing teams and I needed to make adjustments to transfer value to better improve the model. In terms of the variables used in my model, ‘TV’ represents the proportional difference in transfer value across player positions. Due to the closed market structure for international soccer, positional value is relevant. ‘AGE’ indicates the annual difference in ages across positions. Age is critical since transfer value only looks to measure projected long-term performance of the player and value to a club and does not evaluate immediate player performance, meaning older players near or slightly above their playing peak are undervalued. When looking at values, ‘CB\_TV = -1.2’ and ‘CB\_AGE = 2’ would mean the team’s starting Centre Backs (CBs) are, on average, 2 years older, and the opposing team’s starting CBs have a transfer value 120% higher, on average.

Code:

formula ← bf(Outcome ~ 0 + ELO + 0 + GK\_TV + (0 + GK\_TV | GK\_AGE) + 0 + CB\_TV + (0 + CB\_TV | CB\_AGE) + 0 + CM\_TV + (0 + CM\_TV | CM\_AGE) + 0 + CF\_TV + (0 + CF\_TV | CF\_AGE) + 0 + OUT\_TV + (0 + OUT\_TV | OUT\_AGE))

model ← brm(  
formula = formula,  
data = soccer,  
family = bernoulli(),  
prior = prior,  
cores = num.cores,  
control = list(max\_treedepth = 15, adapt\_delta = 0.99))

 ![image.png](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/e/9/e97ad4377f48d4c775c8cc6fe630de43eedd8ebf.png)

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 31, 2024, 3:01pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/4 "2024-05-31T15:01:55Z")

</div>

Since your effective sample size is very low and you don’t have enough observations to even estimate an intercept-only model (the frequentist Wilson 0.95 compatibility interval for the probability of Y=1 is (0.57, 0.77), which is relevant event if you force the covariate-adjusted intercept to be zero), analysis may be futile. For small datasets you have to have an extremely high-resolution Y.

What is a random predictor?

It is often safe to do data reduction outside of the main model fitting stage so that need not necessarily involve Bayesian methods. In the future I will have a way to insert PCA restrictions into logistic regression so all modeling will be on the original predictor scales.

It is often dangerous to omit the intercept in the best of conditions.

---

<div class="post-metadata">

**Author:** ![JulianChristov](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianchristov/32/18657_2.png) [@JulianChristov](https://discourse.mc-stan.org/u/JulianChristov)\
**Post date:** [May 31, 2024, 7:29pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/5 "2024-05-31T19:29:14Z")

</div>

My apologies there are no random predictors. Instead the age variables are just random effects on transfer value to adjust to immediate performance. Please let me know if my data would be helpful to you for the PCA.

Based on my soccer knowledge I strongly believe the intercept is close to if not zero because this is when all things are equal therefore the teams should have the same shot at winning, at 50/50 odds. If one were to say the model would have a intercept it would imply a home team or away team have an advantage regardless of player composition (transfer value) or past results (Elo) which does not seem logical to me, especially where games are played at a neutral location with neutral referees. I understand that there may  
additional complexities to evaluate, such as different coach, but for the most part it seems close to if not zero.

I would really appreciate to get on a zoom call with you. I can send a link whenever you are available. Thank you so much for taking so much time to help me on this.

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [June 1, 2024, 1:19pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/6 "2024-06-01T13:19:12Z")

</div>

Still not understanding ‘random’ here. But to the main issue pick one binary fixed predictor and fit a simple model with it. Get the 0.95 HPD uncertainty interval for the odds ratio for that factor. You’ll find it’s so wide that you haven’t learned anything from the data. It just gets worse with more predictors.

---

<div class="post-metadata">

**Author:** ![JulianChristov](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianchristov/32/18657_2.png) [@JulianChristov](https://discourse.mc-stan.org/u/JulianChristov)\
**Post date:** [June 2, 2024, 2:52am UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/7 "2024-06-02T02:52:54Z")

</div>

In the formula, age indicates that the influence of Transfer Value (TV) on the outcome varies across different age groups. I do not mean age is ‘random’ just that the effect of age on transfer value is not necessarily constant.

I can omit ELO in the modeling process and see if it explains the incorrect predictions. However, it becomes tough to omit or combine other variables since I am looking to measure disparities in players from the two competing teams. Please let me know if you would like my Capstone paper as it explains much of the context and my thought process behind for the topic.

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [June 2, 2024, 12:30pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/9 "2024-06-02T12:30:04Z")

</div>

I don’t mean to collapse all predictors to one measure. Pick one good measure and do a complete analysis. You’ll see that the sample size is inadequate as judged by the extreme width of an uncertainty interval for an odds ratio. Therefore the sample is inadequate for everything you are doing.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [June 3, 2024, 7:02pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/10 "2024-06-03T19:02:18Z")

</div>

> [@harrelfe](#):
>
> You’ll see that the sample size is inadequate as judged by the extreme width of an uncertainty interval for an odds ratio. Therefore the sample is inadequate for everything you are doing.

I think you are too pessimistic here. I know your justifications for binary target sample size, but you are leaving out here in this post the background that you are assuming no prior and requiring quite narrow interval. You claim the posterior interval would have extreme width, but then it’s much narrower than without any data.

> [@harrelfe](#):
>
> Since your effective sample size is very low and you don’t have enough observations to even estimate an intercept-only model (the frequentist Wilson 0.95 compatibility interval for the probability of Y=1 is (0.57, 0.77),

Even with this computation the interval width is only 21% from the uniform interval.

> [@JulianChristov](#):
>
> My thesis paper was to create a model that could predict match outcomes

Also, the original goal was to make predictions and not to infer any posterior intervals. There are plenty of papers showing that we can get useful predictions with even fewer observation (and binary target).

> [@harrelfe](#):
>
> Get the 0.95 HPD uncertainty interval for the odds ratio for that factor. You’ll find it’s so wide that you haven’t learned anything from the data. It just gets worse with more predictors.

Of course, you learn from the data, and usually you learn more information you have. There are plenty of papers showing you can get useful predictions with fewer observations and thousands of predictors. In this case the predictors are likely to be correlating and thus providing similar information, and thus adding more of them is likely to improve predictive performance. Of course, this is assuming you are using Bayesian inference and not using priors that force overfitting.

Couple references

- Juho Piironen, Markus Paasiniemi, and **Aki Vehtari** (2020). Projective inference in high-dimensional problems: prediction and feature selection. _Electronic Journal of Statistics_ , 14(1):2155-2197. [Online](https://doi.org/10.1214/20-EJS1711).
- Juho Piironen and **Aki Vehtari** (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors. _Electronic Journal of Statistics_ , 11(2):5018-5051. [Online](https://projecteuclid.org/euclid.ejs/1513306866)

> [@JulianChristov](#):
>
> I have only 80 observations from World Cup knockout-stage matches.

But, these are not just binary data, as there is also score for both teams, and instead of Bernoulli model, it would be better to use bivariate Poisson model for both scores or Poisson-difference model for the score difference. You can’t do these with brms, but would need to write them in Stan language.

> [@JulianChristov](#):
>
> formula ← bf(Outcome ~ 0 + ELO + 0 + GK\_TV + (0 + GK\_TV | GK\_AGE) + 0 + CB\_TV + (0 + CB\_TV | CB\_AGE) + 0 + CM\_TV + (0 + CM\_TV | CM\_AGE) + 0 + CF\_TV + (0 + CF\_TV | CF\_AGE) + 0 + OUT\_TV + (0 + OUT\_TV | OUT\_AGE))

You have many redundant `0`s here and it would be good to read the documentation on what `0` means in the formula.

While I did say above that for binary target, fewer than 80 can be plenty for getting a useful predictive model (for example, in repeated betting beating a model with no data at all), and I think something like R2D2 type prior would be sensible (instead of improper uniform or indenpendent normals that favor overfitting), I don’t think there would be enough information in these data to improve the predictive performance with varying slope components.

Sorry, I don’t have more time for this, but wanted to counter Franks statements which seemed to state too universal claims about what can be done with binary outcomes.

---

<div class="post-metadata">

**Author:** ![JulianChristov](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianchristov/32/18657_2.png) [@JulianChristov](https://discourse.mc-stan.org/u/JulianChristov)\
**Post date:** [June 5, 2024, 2:10am UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/11 "2024-06-05T02:10:58Z")

</div>

Thank you so much for the advice. First off, I took away the 0’s that were not needed in the formula:

formula ← bf(Outcome ~ 0 + ELO + GK\_TV + (0 + GK\_TV | GK\_AGE) + CB\_TV + (0 + CB\_TV | CB\_AGE) + CM\_TV + (0 + CM\_TV | CM\_AGE) + CF\_TV + (0 + CF\_TV | CF\_AGE) + OUT\_TV + (0 + OUT\_TV | OUT\_AGE))

I know a lot of successful soccer model’s have used a bivariate Poisson model to predict outcomes for soccer matches, like Dixon and Coles.  
However, personally I do not think score line is best in predicting knock-out stage soccer matches and the information it provides may not be as accurate as the true outcome itself.

- There are only 2 outcomes that ultimately determine success in knock-out matches: win or go home.
- Certain teams consistently win by small margins while other teams win less consistently but by large margins.
- Game scenarios can greatly affect score lines even when teams are closely matched.
- The GK position may be more important in knockout-stage matches since games can go penalty kicks.

These sources seemed very useful for what I am working on but a lot of what was talked about is beyond my current level of expertise.

Predictive projection in order to reduce the number of features may not be as useful as regulation techniques to decrease the magnitude of those features, especially since the variables seem key to understanding the outcome in there own ways. I could be wrong though and this may help me mitigate overfitting the training data.

The paper having to do with the regularized horseshoe prior seems very useful in my case.

I will follow your advice using the R2D2 prior though. It seems to fit best with the model I am trying to build since:

- Use R2D2 prior when you are interested in controlling the proportion of variance explained by the predictors and have a moderate number of predictors.
- Use regularized horseshoe prior when you are dealing with high-dimensional data with many predictors, where you expect only a few of them to be truly relevant.

Please let me know if this code is okay for computing LOOCV to determine accuracy of this kind of model since it is my first time using rstan and code came from ChatGPT (Imay end up taking ELO out and also seeing how it performs). Also would PCA be helpful at all to reduce the dimensionality (I am guessing the shrinkage prior should be fine without).

Here is my code:

stan\_model\_code\_r2d2 ← "  
data {  
int\<lower=0\> n; // number of observations  
int\<lower=0\> d; // number of predictors  
int\<lower=0\> k; // number of age groups  
int\<lower=0, upper=1\> y[n]; // binary outcome  
matrix[n, d] x; // predictors  
matrix[n, k] age\_group\_indices; // age group indices  
real\<lower=0\> scale\_icept; // prior std for the intercept  
}

parameters {  
real beta0; // intercept  
real\<lower=0, upper=1\> r2; // R-squared parameter  
vector[d] beta\_raw; // raw coefficients  
matrix[d, k] u; // random effects for each predictor by age group  
}

transformed parameters {  
vector[d] beta; // regression coefficients  
vector[n] eta; // linear predictor

beta = beta\_raw \* sqrt(r2 / (1 - r2));

// Compute linear predictor  
eta = beta0 + x \* beta;  
for (i in 1:n) {  
eta[i] += u[1,] \* age\_group\_indices[i, 1] \* x[i, 2]; // GK\_TV | GK\_AGE  
eta[i] += u[2,] \* age\_group\_indices[i, 2] \* x[i, 3]; // CB\_TV | CB\_AGE  
eta[i] += u[3,] \* age\_group\_indices[i, 3] \* x[i, 4]; // CM\_TV | CM\_AGE  
eta[i] += u[4,] \* age\_group\_indices[i, 4] \* x[i, 5]; // CF\_TV | CF\_AGE  
eta[i] += u[5,] \* age\_group\_indices[i, 5] \* x[i, 6]; // OUT\_TV | OUT\_AGE  
}  
}

model {  
// Prior for R-squared  
r2 ~ beta(2, 2);

// Prior for coefficients  
beta\_raw ~ normal(0, 1);

// Prior for random effects  
for (j in 1:k) {  
u[, j] ~ normal(0, 1);  
}

// Prior for intercept  
beta0 ~ normal(0, scale\_icept);

// Logistic regression model  
y ~ bernoulli\_logit(eta);  
}  
"

library(rstan)  
library(dplyr)

# Load the data

soccer ← read.table(“/mnt/data/luluboyo.txt”, header = TRUE, sep = “\t”)

# Define the number of observations

n\_obs ← nrow(soccer)

# Initialize a vector to store the predicted probabilities

predicted\_probs ← numeric(n\_obs)

# Perform LOOCV

for (i in 1:n\_obs) {

# Prepare training data leaving out the ith observation

soccer\_train ← soccer[-i,]  
soccer\_test ← soccer[i,]

# Extract predictors and outcome for training data

x\_train ← as.matrix(soccer\_train[, c(“ELO”, “GK\_TV”, “CB\_TV”, “CM\_TV”, “CF\_TV”, “OUT\_TV”)])  
age\_group\_indices\_train ← as.matrix(soccer\_train[, c(“GK\_AGE”, “CB\_AGE”, “CM\_AGE”, “CF\_AGE”, “OUT\_AGE”)])  
y\_train ← soccer\_train$Outcome

# Prepare data for Stan

stan\_data\_train ← list(  
n = nrow(x\_train),  
d = ncol(x\_train),  
k = ncol(age\_group\_indices\_train),  
y = as.integer(y\_train),  
x = x\_train,  
age\_group\_indices = age\_group\_indices\_train,  
scale\_icept = 10  
)

# Compile and fit the model

stan\_model ← stan(model\_code = stan\_model\_code\_r2d2, data = stan\_data\_train,  
iter = 2000, chains = 4, cores = parallel::detectCores(),  
control = list(max\_treedepth = 15, adapt\_delta = 0.99))

# Extract posterior samples

posterior\_samples ← extract(stan\_model)

# Prepare the new data (left-out observation)

new\_x ← as.matrix(soccer\_test[, c(“ELO”, “GK\_TV”, “CB\_TV”, “CM\_TV”, “CF\_TV”, “OUT\_TV”)])  
new\_age\_group\_indices ← as.matrix(soccer\_test[, c(“GK\_AGE”, “CB\_AGE”, “CM\_AGE”, “CF\_AGE”, “OUT\_AGE”)])

# Number of posterior samples

n\_samples ← length(posterior\_samples$beta0)

# Initialize a vector to store the predictions

predictions ← numeric(n\_samples)

# Compute the linear predictor (eta) for each posterior sample

for (j in 1:n\_samples) {  
eta ← posterior\_samples$beta0[j] + sum(new\_x \* posterior\_samples$beta[j,])  
eta ← eta + posterior\_samples$u[1, j] \* new\_age\_group\_indices[1, 1] \* new\_x[1, 2]; // GK\_TV | GK\_AGE  
eta ← eta + posterior\_samples$u[2, j] \* new\_age\_group\_indices[1, 2] \* new\_x[1, 3]; // CB\_TV | CB\_AGE  
eta ← eta + posterior\_samples$u[3, j] \* new\_age\_group\_indices[1, 3] \* new\_x[1, 4]; // CM\_TV | CM\_AGE  
eta ← eta + posterior\_samples$u[4, j] \* new\_age\_group\_indices[1, 4] \* new\_x[1, 5]; // CF\_TV | CF\_AGE  
eta ← eta + posterior\_samples$u[5, j] \* new\_age\_group\_indices[1, 5] \* new\_x[1, 6]; // OUT\_TV | OUT\_AGE  
predictions[j] ← 1 / (1 + exp(-eta)) # Apply logistic function to get probabilities  
}

# Store the mean prediction for the ith observation

predicted\_probs[i] ← mean(predictions)

cat(“Finished LOOCV for observation”, i, “with predicted probability:”, predicted\_probs[i], “\n”)  
}

# Print all predicted probabilities

predicted\_probs

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [June 5, 2024, 11:41am UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/12 "2024-06-05T11:41:32Z")

</div>

I’m not surprised you can get a model that is better than no model. But I would be surprised if the posterior distribution of a typical risk estimate is usefully narrow unless the prior information is so extreme that you almost don’t need data.

---

<div class="post-metadata">

**Author:** ![JulianChristov](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianchristov/32/18657_2.png) [@JulianChristov](https://discourse.mc-stan.org/u/JulianChristov)\
**Post date:** [June 6, 2024, 8:33pm UTC](https://discourse.mc-stan.org/t/creating-a-prior-for-a-logistic-model-in-brms/35307/13 "2024-06-06T20:33:21Z")

</div>

Last question. I would like to determine which team is more likely to win based on the transfer value of their starters, adjusted for age. Could I perform a sensitivity analysis on the age groups across different positions to assess the impact on this prediction?