I coded the sigmoid model in Stan as follows (good fit, no divergent transitions):
data {
int<lower=1> N; // number of rows
vector<lower=0>[N] y; // outcome
vector<lower=0>[N] t; // time in days
}
parameters {
real<lower=0> alpha;
real<lower=0> beta;
real<lower=0> M;
real<lower=0> sigma;
}
model {
vector[N] mu;
mu = M ./ (1 + exp(-beta * (t - alpha)));
alpha ~ normal(100, 100); // number of days at which expected # of counts is halfway the maximum
beta ~ normal(0.5, 5); // growth parameter
M ~ normal(1000, 10000); // maximal number of cases
y ~ normal(mu, sigma);
}
However, a multilevel model seems appropriate because so many countries have so few cases at this point (thankfully). To that end, I tried rewriting this as a multilevel model as follows:
data {
int<lower=1> N; // number of rows
int<lower=1> J; // number of locations
int<lower=1,upper=J> location[N]; // the locations themselves
vector<lower=0>[N] y; // outcome
vector<lower=0>[N] t; // time in days
}
parameters {
vector<lower=0>[J] alpha;
vector<lower=0>[J] beta;
vector<lower=0>[J] M;
real<lower=0> M_mu; // average maximal number of cases
real<lower=0> M_sigma;
real<lower=0> alpha_mu; // average number of days at which expected # of counts is halfway the maximum
real<lower=0> alpha_sigma;
real<lower=0> beta_mu; // average growth parameter value
real<lower=0> beta_sigma;
real<lower=0> sigma;
}
model {
vector[N] mu;
for ( i in 1:N ) {
mu = M[location[i]] ./ (1 + exp(-beta[location[i]] * (t - alpha[location[i]])));
}
M ~ normal(M_mu, M_sigma);
alpha ~ normal(alpha_mu, alpha_sigma); // number of days at which expected # of counts is halfway the maximum
beta ~ normal(beta_mu, beta_sigma); // growth parameter
// hyperparameters
M_mu ~ normal(1000, 10000);
alpha_mu ~ normal(100, 100);
beta_mu ~ normal(0.5, 5);
y ~ normal(mu, sigma);
}
However, this model crashes during sampling as it uses up >16GB of ram. That surprises me since there are only 17k rows in this dataset. Any ideas for troubleshooting this?
If the size of the data isn’t a problem, it could be an issue with the indexing. What is the value of J? Since you need to be able to allocate 3 vectors of length J in memory. You could also try re-running the model with a subset of the data (i.e. 1700 rows instead) to see if the model works fine. If the model still fails, that would indicate an error with the syntax (or how the data is being passed).
Also, this may be unrelated to your issue, but I don’t think this loop is doing what you want:
vector[N] mu;
for ( i in 1:N ) {
mu = M[location[i]] ./ (1 + exp(-beta[location[i]] * (t - alpha[location[i]])));
}
Because you don’t have an index on mu (or t), you’re effectively just overwriting muN times, and then the observed data is fitted using the last iteration. You could also use the inv_logit function for better resistance to over-/under-flow:
vector[N] mu;
for ( i in 1:N ) {
mu[i] = M[location[i]] * inv_logit(beta[location[i]] * (t[i] - alpha[location[i]]));
}
N = ~17k and J = 254
Fixing the mu rewrite issue also fixes the memory issue, although I’m not sure why!
It’s still taking forever to run, but I’m trying additional speedups like vectorizing: M[location] .* inv_logit(beta[location] .* (t - alpha[location]))
I think I can take it from here :) Thank you so much for your help!
To speed things up you should have a look into using non-centered parameterisations. Hierarchical parameters can have some tricky posterior geometry to explore which results in divergences. These can instead be parameterised using manipulations of the standard normal distribution. This will often also have the benefit of being faster, since Stan can explore a standard normal distribution very efficiently. More information (and examples) are here in the Stan manual. However these can also perform worse with large sample sizes, so your mileage may vary.
You could also look into the map_rect functionality for parallelising the model. This would likely provide the biggest speed-ups, but requires re-writing the model a little. More info in the Map-Reduce chapter in the manual
Thank you, this is great. I’ve tried de-centering everything as follows but it’s still running very slowly. I’ll keep tinkering including using map_rect
data {
int<lower=1> N; // number of rows
int<lower=1> J; // number of locations
int<lower=1,upper=J> location[N]; // the locations themselves
vector<lower=0>[N] y; // outcome
vector<lower=0>[N] t; // time in days
}
parameters {
vector<lower=0>[J] alpha_raw;
vector<lower=0>[J] beta_raw;
vector<lower=0>[J] M_raw;
real<lower=0> M_mu_raw;
real<lower=0> M_sigma_raw;
real<lower=0> alpha_mu_raw;
real<lower=0> alpha_sigma_raw;
real<lower=0> beta_mu_raw;
real<lower=0> beta_sigma_raw;
real<lower=0> sigma;
}
transformed parameters{
real<lower=0> M_mu; // average maximal number of cases
real<lower=0> alpha_mu; // average number of days at which expected # of counts is halfway the maximum
real<lower=0> beta_mu; // average growth parameter value
real<lower=0> M_sigma;
real<lower=0> alpha_sigma;
real<lower=0> beta_sigma;
vector<lower=0>[J] alpha;
vector<lower=0>[J] beta;
vector<lower=0>[J] M;
M_mu = 1000 + M_mu_raw * 25000;
alpha_mu = 100 + alpha_mu_raw * 50;
beta_mu = 0.3 + beta_mu_raw * 2.5;
M_sigma = M_sigma_raw * 25000;
alpha_sigma = alpha_sigma_raw * 50;
beta_sigma = beta_sigma_raw * 2.5;
alpha = alpha_raw * alpha_sigma + alpha_mu;
beta = beta_raw * beta_sigma + beta_mu;
M = M_raw * M_sigma + M_mu;
}
model {
vector[N] mu;
// use inv_logit for better computational stability, an vectorize for speed
mu = M[location] .* inv_logit(beta[location] .* (t - alpha[location]));
M_raw ~ std_normal();
alpha_raw ~ std_normal();
beta_raw ~ std_normal();
// hyperparameters
M_mu_raw ~ std_normal();
alpha_mu_raw ~ std_normal();
beta_mu_raw ~ std_normal();
M_sigma_raw ~ std_normal();
alpha_sigma_raw ~ std_normal();
beta_sigma_raw ~ std_normal();
y ~ normal(mu, sigma);
}
I was trying to do something similar a few days ago with a hierarchical Gompertz model implemented through brms - see thread here: Hierarchical Gompertz model
I gave up in the end as no matter what I tried I could not quite make it work. I think there might be an identifiability problem with the equation - the equation has 3 fixed effect parameters and one variable. But, I don’t know quite so much about identifiability so its just a guess.
Hey jroon. I don’t think it’s an identifiability issue as I can recover all the parameters (I’m sure constraining everything to positive helps). When I fit the model to a subset of a thousand rows of the dataset as suggested by @andrjohns I get very sensible results that plot over the data just fine. The decentered model I posted does fit (slowly) without divergent transitions – once the kaggle markdown is finished I’ll edit this post to include a link
It runs and fits the data well when I subset to just a few countries (no divergences) but now it’s giving weird results for the full dataset! I’m going to try different prior values – I think the ones I chose were unrealistically large and wide–and will take your advice on selecting ~10 countries
The blog notes that certain values of the parameter bring out the pathology - perhaps in a multi-level model the chance to hit a pathological value is higher because you now have as many potentially pathological variables as you do countries? Perhaps @martinmodrak who wrote the blog has some thoughts here?
Sigmoids are challenging to fit. In this case I would guess the biggest problem would be that the upper plateau parameter is almost not constrained by data at all… (this is also problem for any predictions from such a model). Finally the data is messy and has all sorts of reporting biases and non-stationarity due to evolving reactions of the governments and public, so a sigmoid might be a problematic model.
With epidemic modelling in this time I would also ask you to be very careful about communicating the results. There is a lot of noise in the public communication channels that can easily drown out the voices of the best experts so please make sure you are not adding to this noise.