Hi,

I am trying to fit a binomial regression with a log link to some count data to estimate a relative risk (model posted last). The model seems to work OK, and retrieves the expected relative risk corrected for age (RR = exp(betaA)). However, I am sometimes running into the problem of bad initialisation. Namely that theta > 1. I’ve tried altering the init range from 0.0000000000001-5, with no luck.

I saw there was some previous discussion on binomial models with a log link (Log binomial random intercept model). From what I gathered, the conclusion was that the log-link is not good for probabilities. However, if not using the binomial-log model I am unsure how to model relative risks. I could, technically, fit any kind of model and then simulate the counts using the input variables and calculate the Cochran-Mantel-Haenzel adjusted RR (http://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/BS704-EP713_Confounding-EM/BS704-EP713_Confounding-EM7.html). This, however, poses two problems:

- Which alternative model to fit? A fully saturated poisson log-linear model would involve a lot of interactions and seems undesirable (especially if I want to examine a risk factor x age interaction at some point)
- I might end up with some zero division, if some values of c in the posterior are zero (notation according to the website)

Alternatively, I could provide some custom initialisation values and hope Stan figures out how to effectively explore the posterior – it seems to be running OK as long as the initialisation is valid. One option could be to pick a set of values which sum < 0. Would this still be considered a diffuse initialisation, so that the R-hat values are appropriate?

Any advice would be much appreciated.

Thanks in advance!

```
data {
int <lower=0> N; // Number of rows
int <lower=0> B[N]; // Count of patients in group
int <lower=0> Total[N]; // Count of patients in group
vector<lower=0, upper=1>[N] A;
vector<lower=0, upper=1>[N] Agebin6;
vector<lower=0, upper=1>[N] Agebin7;
vector<lower=0, upper=1>[N] Agebin8;
vector<lower=0, upper=1>[N] Agebin9;
vector<lower=0, upper=1>[N] Agebin10;
vector<lower=0, upper=1>[N] Agebin11;
vector<lower=0, upper=1>[N] Agebin12;
vector<lower=0, upper=1>[N] Agebin13;
vector<lower=0, upper=1>[N] Agebin14;
vector<lower=0, upper=1>[N] Agebin15;
vector<lower=0, upper=1>[N] Agebin16;
vector<lower=0, upper=1>[N] Agebin17;
}
parameters {
// Coefficients
vector[J] beta0;
vector[J] betaA;
vector[J] betaAge6;
vector[J] betaAge7;
vector[J] betaAge8;
vector[J] betaAge9;
vector[J] betaAge10;
vector[J] betaAge11;
//vector[J] betaAge12; // Reference
vector[J] betaAge13;
vector[J] betaAge14;
vector[J] betaAge15;
vector[J] betaAge16;
vector[J] betaAge17;
}
model {
beta0 ~ normal(0, 1);
betaA ~ normal(0, 1);
betaAge6 ~ normal(0, 1);
betaAge7 ~ normal(0, 1);
betaAge8 ~ normal(0, 1);
betaAge9 ~ normal(0, 1);
betaAge10 ~ normal(0, 1);
betaAge11 ~ normal(0, 1);
//betaAge12 ~ normal(0, 1);
betaAge13 ~ normal(0, 1);
betaAge14 ~ normal(0, 1);
betaAge15 ~ normal(0, 1);
betaAge16 ~ normal(0, 1);
betaAge17 ~ normal(0, 1);
/* Likelihood specification */
{
vector[N] theta;
theta = exp(beta0[j] +
betaA[j] .* A +
betaAge6[j] .* Agebin6 +
betaAge7[j] .* Agebin7 +
betaAge8[j] .* Agebin8 +
betaAge9[j] .* Agebin9 +
betaAge10[j] .* Agebin10 +
betaAge11[j] .* Agebin11 +
//betaAge12[j] .* Agebin12 +
betaAge13[j] .* Agebin13 +
betaAge14[j] .* Agebin14 +
betaAge15[j] .* Agebin15 +
betaAge16[j] .* Agebin16 +
betaAge17[j] .* Agebin17);
// Increment likelihood
B ~ binomial(Total, theta);
}
}```
```