Model that shrinks a set of coefficients towards their common mean I am interested in estimating the odds of a certain disease based on a medium sized group of correlated biological markers (roughly 20 markers). The model will also include several confounding variables as covariates (i.e. sex, age, BMI, education).
Typically this is done with logistic regression - either in a large model including all predictors, or in several separate models (one for each biological predictor including the confounders). The problem with the latter approach is that it neglects to account for the correlation amongst those markers.
An approach I have seen in the literature is to 'shrink' the regression coefficients for the biological markers towards a common mean. This makes sense in my case - these biological markers are related and I would assume their effect sizes would be roughly similar. A large effect size for a single marker is unlikely to be valid.
The problem is that I have no idea how to implement this! My exposure to Bayesian methods is limited (although I have used R package brms), so if someone would be able to point me in the direction or provide any insight, that would be excellent!
 A: I'll offer a partial answer since there has been little activity since you asked this.
A keyword here is hierarchical modeling. It is possible to specify a hierarchical model for your biological markers to shrink their coefficients toward one another while letting the other predictors vary independently. A Bayesian hierarchical model would be something like
$$
\beta_i \sim \mathrm{Normal}(\beta, \sigma_\beta),\\
\beta \sim \mathrm{Normal}(0, 1),\\
\sigma_\beta \sim \mathrm{exponential}(1),\\
\gamma_i \sim \mathrm{Normal}(0, 1).
$$
In this example the $\beta_i$ are coefficients for the correlated biomarkers, and the $\gamma_i$ are coefficients for the other predictors. What we are assuming here is that the several $\beta_i$ vary around a common average $\beta$ (no subscript) with a standard deviation of $\sigma_\beta$. Both of these to be estimated from the data, but with some prior assumption about roughly where they will lie.
Depending on how comfortable you are with Bayesian data analysis this will seem more or less natural and straight-forward. Were it me I would have written something like this directly in Stan (other samplers are available), which lets you specify the whole model explicitly. It might also be possible to do this in brms, which builds a Stan program in the background, but you would have to wrestle the model formulation syntax (borrowed from lme4), which may or may not be suited to your needs.
As you say I don't think shrinkage on only some coefficients is possible in glmnet. Such an approach could sidestep the need for priors at the cost of estimating a shrinkage parameter and sacrificing interpretability. I think it is possible in rms (not to be confused with brms) somehow, but I am no expert.
