SEIR анализ на данните за Covid-19 Индия

В тези уравнения β може да се приеме като процент на контакт между заразените и податливите лица, докато γ е просто процент на възстановяване или средно време, необходимо на заразените лица да се възстановят. ае инкубационен период. Горното уравнение също има условия за смъртност и раждаемост, но в този анализ няма да ги разглеждаме.

Как да използвате тази информация за данни за Covid-19 India?

В Индия с каквито и данни да разполагаме, повечето от заразените хора отнемат около 15–25 дни, за да се възстановят, така че за тази статия γможе да бъде 1/25. Инкубационният период на короната е около 5,2 дни (спорно, но в повечето статии е 5,2), така че за тази статия aе 1/5, но не знаем β.Ще използваме напасване на кривата, за да определим този параметър β. Scipy (пакет в Python) има някакъв инструмент за напасване на крива, а също и за решаване на диференциално уравнение.

Създаване на функция за решаване на диференциално уравнение

# Total population, N.
N = 1210568111
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 1, 0
E0=0
# Everyone else, S0, is susceptible to infection initially.
S0 = N — I0 — R0
# Contact rate, beta, and mean recovery rate(rec),incubation period(gamma),
beta, gamma,rec = 5,5,25 
# A grid of time points (in days)
t = df_train[‘days’]
x_test = df_test[‘days’]
# The SIR model differential equations.
def deriv(t, y, N,gamma,rec,beta):
 S, E, I,R = y
 dSdt = -beta * S * I / N
 dEdt = (beta * S * I / N) — ((1/gamma) * E)
 dIdt = ((1/gamma) * E)-((1/rec)*I)
 dRdt = (1/rec) * I
 return dSdt,dEdt, dIdt, dRdt

В горния кодов фрагмент deriv е името на функцията, която ще генерира необходимите диференциални уравнения. Сега тези диференциални уравнения ще бъдат решени с друга функция.

def fit_odeint(t, beta):
 y0 = S0, E0, I0,R0
 ret = solve_ivp(deriv,[0,t.max()],y0,args=(N,gamma,rec,beta),t_eval=np.arange(t.max()+1),method=’Radau’)
 S, E, I,R = ret.y
 if(I.shape[0]==t.shape[0]):
   if((S.sum()*E.sum()*I.sum()*R.sum())>0):
       res=np.concatenate(((I+R),R),axis=0)
       return(np.array(res))
   else:
       return(np.zeros((2*t.shape[0],)))
 else:
 return(np.zeros((2*t.shape[0],)))

Тази функция използва функцията solve_inp на пакета scipy. Той решава система от диференциални уравнения, която в този случай е дадена от функция deriv. Сега всички вие може би сте забелязали, че тази функция fit_odeint има още един параметър бета (коефициент на контакт), чиято стойност е неизвестна и ние ще намерим тази стойност с помощта на напасване на кривата.

popt, pcov = optimize.curve_fit(fit_odeint,t,y,maxfev=500000000,p0=[beta])

Резултати

Сега е време да видим резултата от цялата тази огромна работа, която свършихме.

Данните от миналата седмица се считат за тестови данни, а всички останали като данни за влака.

Прогнозата за данните от теста не е толкова добра(R2=0,83)но защо е така? Това ли е защото не сме взели под внимание ефекта отблокирането? Нека да проверим това.

Как да вземем предвид ефекта на блокиране в модела SEIR?

Поради блокирането процентът на контактите ще намалее. Има няколко опции за функция за разпадане. Източник - https://github.com/SwissTPH/openmalaria/wiki/ModelDecayFunctions

Използвах функцията на хълм с K=2 и използвах подходяща крива за намиране на L (скорост на разпадане).

Модифицирана функция:

def deriv_1(t, y, N,gamma,rec,L,k):
 S, E, I,R = y
 beta=0.27/((1 + (t/L)**k))
 dSdt = -beta * S * I / N
 dEdt = (beta * S * I / N) — ((1/gamma) * E)
 dIdt = ((1/gamma) * E)-((1/rec)*I)
 dRdt = (1/rec) * I
 #print(dSdt)
 #print(S,R)
 return dSdt,dEdt, dIdt, dRdt

Заключването беше наложено на 25 март 2020 г., така че до този ден се използваше постоянна скорост на контакт и след това беше приложена тази функция на хълм. Тези 0,27 се намират след известна итерация.

Установено е, че L е голям, така че честотата на контакт не намалява много. Възможно е да са необходими повече данни, за да се знае скоростта на намаляване.Отбележете, че всички данни след като блокирането е използвано за намиране на L. За разлика от предишния пример, той не беше разделен на тестване и тренировка, само за да се увеличи наборът от данни за тренировка, така че човек да може да придобие по-добро разбиране на данните.

Нека сега се опитаме да предскажем бъдещето

Ключови изводи:

1› Според този модел ще отнеме около 210 дни повече от 31/01/2020 (първата седмица на август), за да се стабилизира кривата. В този анализ намаляването на процента на контакт е много по-малко, но в действително състояние поради блокиране и други мерки честотата на контакта ще намалее по-бързо и така тази продължителност може да бъде много по-кратка. Можете също така да забележите, че стойността на потвърдения случай е доста висока. Това е така, защото процентът на контакт е амортизацията не е много забележим в модела.

2› Увеличаването на броя на болничните легла и изолирането на всички заразени лица ще намали значително процента на контакт и това ще помогне за стабилизиране на кривата.