# Linear excess risk model

I am trying to fit a linear excess risk model as mentioned here:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3139969/

Briefly, this is the model:

``````   odds = exp(alpha)*(1+b1*x1+b2*x2+b3*x1*x2)
pr(y=1) = odds/(1+odds)

``````

Stan model:

``````data{
int N;               //sample size
int ncol;          //number of columns of data
matrix[N,ncol] X;   //matrix of data
int<lower=0, upper=1> y[N]; //outcome
}
parameters{
real a;
vector[ncol] b;
}
transformed parameters{
vector[N] xb1;
vector<lower=0>[N] odds;                       //odds cannot be negative
vector<lower=0,upper=1>[N] pr;

xb1 = X*b;
odds = exp(a)*(1+xb1);

for(n in 1:N){
pr[n] = odds[n]/(1+odds[n]);
}
}
model{
a ~ cauchy(0,2.5);
for(i in 1:ncol){
b[i] ~ cauchy(0,2.5);
}
y ~ bernoulli(pr);
}
generated quantities{
vector[N] yrep;
for(n in 1:N){
yrep[n] = bernoulli_rng(pr[n]);
}
}
``````
``````set.seed(4507907)

N <- 800
x1 <- rpois(N,10)
x2 <- rbinom(N,1,0.5)
x1x2 <- x1*x2

b <- c(-3,2,2,1)

odds <- exp(b)*(1+ b*x1 + b*x2 + b*x1x2)
p <- odds/(1+odds)

y <- rbinom(N,1,p)
table(y)

X = cbind(x1,x2,x1x2)

dat = list(N=N, X=X, y=y, ncol=ncol(X))

fit <- stan("linOdds.stan",

``````

Tested it with a simulated data. However, many divergent transitions.
When I truncated the priors to be non-negative (see below).

``````parameters{
real a;
vector<lower=0>[ncol] b;
}
``````

No divergent transitions, posterior looks okay.

``````Inference for Stan model: linOdds.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
a    -3.05    0.03  0.85 -4.67 -3.61 -3.07 -2.48 -1.42   770 1.00
b  3.38    0.26  5.56  0.34  1.19  2.20  3.94 11.95   474 1.01
b  5.05    1.02 19.95  0.11  0.96  2.23  4.54 21.46   382 1.01
b  1.57    0.06  1.85  0.06  0.49  1.05  2.00  6.02  1115 1.00
``````

However, the pairs plots show some interesting shapes Also there could be negative effects in real data examples.

The issue is the coefficients can be negative in real data.
So I guess my question is could any one suggest a reparameterisation?
To avoid divergent transitions?

Do I need to do something with the target distribution?

Sorry for the slow response. Have you tried modeling on the logit scale? That is, using the log odds instead of odds. Stan has a special bernoulli_logit distribution that might be helpful.

Hi Jonah,

May be I am not understanding something, but if I model it on logit scale it would be simillar to logistic regression.

I found the soultion as suggested by: Haitao Chu et.al here https://insights.ovid.com/pubmed?pmid=21228700

There needs to be two more constrains

``````  real<lower=-1> b1;
real<lower=-1> b2;
real<lower=-(1+b1+b2)> b3;
``````

Thank you again.

Sreenath

Glad you found the solution. I suggested the log-odds because fitting a Bernoulli model is often more stable numerically using bernoulli_logit than bernoulli, but that may not make sense here (only had a chance to take a brief look). I haven’t had a chance to look into the paper you linked to but if it provides a solution that results in good sampling then that’s great.