Я пытаюсь запустить расширение V3 модели COVID Имперского колледжа Лондона (их репозиторий на github находится здесь: https://github.com/ImperialCollegeLondon/covid19model). Единственное, что я изменил, - это добавление данных для четырех крупных канадских провинций к данным из 14 европейских стран, для которых построена модель, и расширение прогноза до 30 дней в будущее (по умолчанию 7, но код разработан, чтобы позволить такое расширение).
Основная функция модели состоит в оценке ожидаемого ежедневного числа погибших на основе фиксированных эффектов страны и принимаемых мер для замедления распространения вируса. Он также производит ожидаемое число погибших при контрфактическом сценарии, когда никакие меры никогда не используются (например, единственное, что влияет на число погибших, - это прогнозируемое общее количество инфекций и уровень смертности от инфекций) в блоке сгенерированных количеств стандартного кода.
Каждый раз, когда я запускаю его, я получаю такие предупреждения:
1: In validityMethod(object) :
The following variables have undefined values: E_deaths0[128,18]. Many subsequent functions will not work correctly.
E_deaths0 относится к ожидаемым смертельным случаям без мер по смягчению последствий на 128-й день, последний день прогноза, в регионе 18 ( 2 июня в Квебе c, в данном случае).
код для генерации этих оценок:
E_deaths0[1, m]= 1e-15 * prediction0[1,m];
for (i in 2:N2){
E_deaths0[i,m] = ifr_noise[m] * dot_product(sub_col(prediction0, 1, m, i-1), tail(f_rev[m], i-1));
}
Где prediction0 [i, m] - это предполагаемое количество заражений в день i в стране m, N2 - общее количество дней для оценки ( 128). ifr_noise [m] - это модификатор смертности от инфекции в стране m, а f_rev [m] - вектор вероятностей смерти в каждый день после заражения. Нет причин, по которым какой-либо из элементов, которые go при вычислении E_deaths0 когда-либо будет NA, и хотя prediction0 оценивается в 0 во многих выборках для поздних дат, а более поздние значения f_rev [m] округляются до 0, это не Похоже, это должно вызвать проблемы, потому что никакого разделения не происходит.
Это происходит каждый раз для разных географических регионов. Даты также различаются, хотя это всегда последние n дней моделирования в данном регионе, где n находится в диапазоне от 1 до 8 или около того. Я предполагаю, что проблема в sub_col(prediction0, 1, m, i-1), потому что f_rev[m] - это данные, поэтому проблем с выборкой нет, и если бы ifr_noise[m] принимал значения NA, то все даты для geography m завершились бы ошибкой. Однако я не понимаю, как это сделать, потому что prediction0 генерируется непосредственно перед E_deaths0 и никогда не вызывает предупреждений.
Когда я извлекаю предсказания E_deaths0 [128,18] или какие бы то ни было оценки имеют проблемы в этом прогоне из соответствия модели , есть столько образцов, сколько я ожидал, и все они находятся в приемлемом диапазоне. Я не уверен, могу ли я доверять этим результатам, когда процесс их генерации выдает случайные предупреждения.
У кого-нибудь есть опыт работы с подобными проблемами или какие-либо идеи? Спасибо.
EDIT: вот код модели:
data {
int <lower=1> M; // number of countries
int <lower=1> P; // number of covariates
int <lower=1> N0; // number of days for which to impute infections
int<lower=1> N[M]; // days of observed data for country m. each entry must be <= N2
int<lower=1> N2; // days of observed data + # of days to forecast
int cases[N2,M]; // reported cases
int deaths[N2, M]; // reported deaths -- the rows with i > N contain -1 and should be ignored
matrix[N2, M] f; // h * s
matrix[N2, P] X[M]; // features matrix
int EpidemicStart[M];
real pop[M];
real SI[N2]; // fixed pre-calculated SI using emprical data from Neil
}
transformed data {
vector[N2] SI_rev; // SI in reverse order
vector[N2] f_rev[M]; // f in reversed order
for(i in 1:N2)
SI_rev[i] = SI[N2-i+1];
for(m in 1:M){
for(i in 1:N2) {
f_rev[m, i] = f[N2-i+1,m];
}
}
}
parameters {
real<lower=0> mu[M]; // intercept for Rt
real<lower=0> alpha_hier[P]; // sudo parameter for the hier term for alpha
real<lower=0> gamma;
vector[M] lockdown;
real<lower=0> kappa;
real<lower=0> y[M];
real<lower=0> phi;
real<lower=0> tau;
real <lower=0> ifr_noise[M];
}
transformed parameters {
vector[P] alpha;
matrix[N2, M] prediction = rep_matrix(0,N2,M);
matrix[N2, M] E_deaths = rep_matrix(0,N2,M);
matrix[N2, M] Rt = rep_matrix(0,N2,M);
matrix[N2, M] Rt_adj = Rt;
{
matrix[N2,M] cumm_sum = rep_matrix(0,N2,M);
for(i in 1:P){
alpha[i] = alpha_hier[i] - ( log(1.05) / 6.0 );
}
for (m in 1:M){
prediction[1:N0,m] = rep_vector(y[m],N0); // learn the number of cases in the first N0 days
cumm_sum[2:N0,m] = cumulative_sum(prediction[2:N0,m]);
Rt[,m] = mu[m] * exp(-X[m] * alpha - X[m][,5] * lockdown[m]);
Rt_adj[1:N0,m] = Rt[1:N0,m];
for (i in (N0+1):N2) {
real convolution = dot_product(sub_col(prediction, 1, m, i-1), tail(SI_rev, i-1));
cumm_sum[i,m] = cumm_sum[i-1,m] + prediction[i-1,m];
Rt_adj[i,m] = ((pop[m]-cumm_sum[i,m]) / pop[m]) * Rt[i,m];
prediction[i, m] = Rt_adj[i,m] * convolution;
}
E_deaths[1, m]= 1e-15 * prediction[1,m];
for (i in 2:N2){
E_deaths[i,m] = ifr_noise[m] * dot_product(sub_col(prediction, 1, m, i-1), tail(f_rev[m], i-1));
}
}
}
}
model {
tau ~ exponential(0.03);
for (m in 1:M){
y[m] ~ exponential(1/tau);
}
gamma ~ normal(0,.2);
lockdown ~ normal(0,gamma);
phi ~ normal(0,5);
kappa ~ normal(0,0.5);
mu ~ normal(3.28, kappa); // citation: https://academic.oup.com/jtm/article/27/2/taaa021/5735319
alpha_hier ~ gamma(.1667,1);
ifr_noise ~ normal(1,0.1);
for(m in 1:M){
deaths[EpidemicStart[m]:N[m], m] ~ neg_binomial_2(E_deaths[EpidemicStart[m]:N[m], m], phi);
}
}
generated quantities {
matrix[N2, M] prediction0 = rep_matrix(0,N2,M);
matrix[N2, M] E_deaths0 = rep_matrix(0,N2,M);
{
matrix[N2,M] cumm_sum0 = rep_matrix(0,N2,M);
for (m in 1:M){
for (i in 2:N0){
cumm_sum0[i,m] = cumm_sum0[i-1,m] + y[m];
}
prediction0[1:N0,m] = rep_vector(y[m],N0);
for (i in (N0+1):N2) {
real convolution0 = dot_product(sub_col(prediction0, 1, m, i-1), tail(SI_rev, i-1));
cumm_sum0[i,m] = cumm_sum0[i-1,m] + prediction0[i-1,m];
prediction0[i, m] = ((pop[m]-cumm_sum0[i,m]) / pop[m]) * mu[m] * convolution0;
}
E_deaths0[1, m]= 1e-15 * prediction0[1,m];
for (i in 2:N2){
E_deaths0[i,m] = ifr_noise[m] * dot_product(sub_col(prediction0, 1, m, i-1), tail(f_rev[m], i-1));
}
}
}
}
Данные представляют собой список (называемый stan_data) с 17 элементами: * M, целое число, равное 18 в моей версии , 14 в модели ICL
* N, вектор из 18 целочисленных значений, указывающих количество дней доступных данных для каждой страны
* смертей, матрица 135x18 наблюдаемых ежедневных смертей. Будущие даты заполняются матрицей -1
* f, 135x18. Каждая запись f [i, m] указывает вероятность смерти в день i с учетом воздействия в день 0, с учетом предполагаемого распределения времени от контакта до смерти и расчетного уровня смертности от инфекции в стране m
* N0, целое число, равное до 6 в обеих версиях. Модель засевает инфекцию, устанавливая равное количество инфекций в дни с 1 по N0 на основе параметра выборки tau
* случаев, матрицы 135x18 выявленных случаев за день в каждой стране.
* SI, вектор 135 удвоен. числа точности. Они представляют вероятность передачи вируса другому человеку каждый день на основе эмпирических данных о серийном интервале инфекций (не опубликовано)
* EpidemicStart, вектор из 18 целых чисел. Статистическая подгонка модели основана на подборе отрицательной биномиальной модели к наблюдаемым количествам смертей в каждой стране с дня = EpidemicStart [m] до даты N [m] (т. Е. со дня начала эпидемии c до последнего дня данных наблюдений). В настоящее время и в версии ICL, и в моей каждая запись EpidemicStart = 31.
* pop, вектор из 18 целых чисел с населением каждой страны или провинции
* N2, целое число, указывающее количество дней для моделирования. В настоящее время установлено значение 135 (с 27 января по 30 дней с последнего дня данных)
* x, вектор из 135 целых чисел, считая от 1 до N2
* P, целое число, указывающее количество вмешательств, которые необходимо сократить. трансмиссия моделируется. На данный момент установлено 6 в обеих версиях модели.