neděle 29. března 2020

Co bylo dřív a jak je to s exponencielou

Je jasné, že růst počtu potvrzených případů může mít dvě příčiny: že je nakaženo víc lidí a/nebo že se víc testuje. Kdyby se zítra provedlo deset milionů testů, počet nakažených by vyskočil astronomicky, nicméně každý, kdo je nakažený, by se izoloval a epidemie by skončila dříve. Skutečnou hrozbou nejsou pozitivně testovaní, ale neodhalení přenašeči. V tomto článku se budu zabývat otázkou, do jaké míry počty potvrzených případů reflektují velikost této hrozby a zda a jak ji lze predikovat.
Budu uvažovat zjednodušeně takto: v populaci I nakažených, o kterých nevíme a v některých případech to o sobě něvědí ani oni. Počet nakažených v čase t má střední hodnotu m(t). Každý den pocítí p procent nakažených symptomy a jde se nechat otestovat. Kroně nich k testování přijde A lidí, kteří nakaženi nejsou. Předpokládejme, že je test spolehlivý a že pokud z něho někdo vyjde pozitivní, bude se testovat dalších k lidí z jeho okolí. Za těchto předpokladů je počet testů provedených za den roven T=pI+A+kpI, zatímco denní přírůstek pozitivních je P=pI+qpI, kde q je poměr nakažených mezi lidmi z okolí. V čase t tedy platí

T(t)=pm(t)+a+kpm(t)+e(t)=a+bn(t)+e(t),     P(t)=pm(t)+qkpm(t)+f(t)=n(t)+f(t), 

kde je střední hodnota A, n(t) je násobek trendu (samotný ho z rovnic nedostaneme, ale to nevadí,  zajímá nás hlavně jeho tvar), b je konstanta a e,f jsou náhodné chyby. V souladu s obecným přesvědčením, že u nás trend není exponenciální, budeme ho hledat mezi polynomy třetího řádu (více už by bylo moc, ne vzhledem k možnému růstu, ale čím větší řád bychom použili, tím více by se "trend" přizpůsoboval náhodným fluktuacím). Výsledky lze vidět na obrázcích
Denní přírůstek pozitivních (10-28.3.)

Denní počet testů  (10-28.3.)

Co se týče čísel, trend vychází n(t)=-4+17.6t-1.4t2+0.08t3, střední hodnota počtu zdravých, co se přišli sami testovat, je a=132 +/- 25 a intervalová předpověď na následuící dny je

        95% dolní   střední  95% horní
29.3.         316       406        612
30.3.         372       464        688
31.3.         434       530        772
1.4.          505       603        865
2.4.          584       685        966

Nyní se vraťme k otázce, jakým způsobem počet testů ovlivňuje či neovlivňuje počty pozitivních. To, že tato čísla spolu souvisí, je jasné i z grafů i z korelace odchylek od trendu, která je rovna 0.81. Spíš je tady otázkou, co bylo dříve, zda.. znáte to. Více testů může být proto, že přišlo více nakažených, nebo proto, že se jich zkrátka udělalo víc, následkem čehož přibylo pozitivních (například poslední den našich dat - v sobotu - vidíme na obou grafech ptopad), Protože však víme, že testů cca. desetkrát více než nakažených - viz obrázek
Procento pozitivních výsledků

- pravda bude asi spíš to druhé. Zkusme tedy (už bez teorie, to bych s tím strávil den) "zkorigovat" počty pozitivních poměrem počtu testů a trendu: Q=P/(T/(a-bn)). Následující obrázky obsahují porovnání P a Q a jejich trendů.
Srovnání P a Q
Srovnání trendů P a Q
Jakkoli korigovaný trend vypadá o něco méně rozháraně, a jakkoli kvadratická chyba odhadu korigovaného trendu vyšla dvakrát menší než u nekorigovaného, ze statistického hlediska jsou trendy skoro stejné (máme jen 2x19 pozorování!). 
Dále zkoumejme hypotetickou situaci, kdy se v našem zjednodušeném světě dodatečně provede S testů, ať už proto, že se zvětšílo k (testoval by se širší okruh), nebo se třeba testovalo namátkou mezi pravděpodobnými přenašeči. Je jasné, že v obou případech se zvedne i počet pozitivních, a to o rS, kde r je poměr nakažených v nově testované skupině. Protože se dá předpokládat, že toto r bude menší než q (vzdálené okolí bude mít pravděpodobnost nákazy menší než to nejbližší), a protože mezi pozitivními jsou i ti, co přišli s příznaky (s pravděpodobností pozitivity I/(I+A), tedy blízkou jedné), bude mít zmiňované navýšení za následek snižování R - procenta pozitivních výsledků testů. Když se však podíváme na graf výše, nic takového u nás nenastává. Asi by bylo lepší, kdyby se dodatečně testovalo víc, na druhé straně tato stabilita spolu s faktem, že počty testů a počty pozitivních vykazují velmi podobný trend, nasvědčuje tomu, že přírůstky pozitivně testovaných poskytují dobrý odhad vývoje trendu epidemie.
Než skončíme, vraťme se ještě k otázce trendu. Na sítích se sice šíří optimistické přesvědčení, že "jsme pod exponenciálou", obávám se však, že je předčansé. Toto přesvědčení ravděpodobně vzniklo na základě chybného odhadu z kumulovaných počtů, tj. z regrese log(P(1)+..+P(t))=a+bt, který vyjde takto:
 
Problém ale je, že tento odhad může být značně vychýlený, protože zde neplatí předpoklad nezávislosti reziduí (není žádný důvod, aby log(P(1)+P(2))-(a+2b) bylo nezávislé na log(P(1))-(a+b)). 
Pokud budeme trend odhadovat z logaritmů přírůstků (které v exponencíále také rostou exponenciálně), dostaneme poněkud jiný obrázek.

Zelená křivka znázorňuje exponencíální trend s počátkem odhadu 10. března (kdy se zavřely školy), modrý trend je odhadován z dat od 5. března (odkdy už jsou přírůstky nenulové). Jasně zde vidíme, jak je zde odhad ošemetný: zelená exponenciela "vyhladila" data skoro stejně jako polynom. Navíc i kdybychm nade vší pochybnost dokázali, že byla křivka až doteď polynomiální, kde můžeme vzít jistotu, že se epidemie nezdačne chovat odteď jinak? Na definitivní optimismus je evidentně ještě dost brzo.

Poznámka: Odhad dvourovnicového lze proést "standardně" pomocí vážených nejmenších čtverců: vážených proto, že je variance náhodných chyb úměrná trendu (to lze okamžitě vidět z grafu reziduí) - kdybychom toto zanedbali a použili "obyčejné" nejmenší čtverce, ztratili bychom cenné informace ze začátku řady, kdy byly počty ještě malé.

Gretl skript:

# nahraďte svým. Musí mít sloupce t - čas (16.březen=43900), P, T - kumulativní počty pozitivních / testovaných
open /home/martin/Documents/s/covid/data/cr.ods



setobs 1 1 --special-time-series




diff P

diff T





smpl 39 62




logs d_P

series tau = t-43900



ols l_d_P const tau

exptrend = exp($coeff[1]+tau*$coeff[2])



logs P

ols l_P const tau
cexptrend = exp($coeff[1]+tau*$coeff[2])
gcexp <- gnuplot P cexptrend --time-series --with-lines





smpl 44 62




ols l_d_P const tau

exp2trend = exp($coeff[1]+tau*$coeff[2])
gexp <- gnuplot d_P exptrend exp2trend --time-series --with-lines



series pomer = d_P / d_T

gr <- gnuplot pomer --time-series --with-lines
gexp <- gnuplot d_P exptrend --time-series --with-lines



scalar cnst=0

scalar lin=0
scalar quad=0
scalar cub=0
scalar a = 0
scalar b = 0



mle loglik = - (d_P - tr )^2 / ( 2 * (d_P(-1)+1)) \

             - (d_T - b * tr - a)^2 / ( 2 * (d_T(-1)+10))   
    series tr = cnst + lin*tau + quad*tau^2 + cub* tau^3 
    params cnst lin quad cub a b
end mle  



trendp =  cnst + lin*tau + quad*tau^2 + cub* tau^3 

trendt = b * trendp + a
up = trendp-d_P
ut = trendt-d_T
corr up ut

uhat = up/sqrt(trendp)
stu=sqrt(var(uhat))

print stu


g_P <- gnuplot d_P trendp --time-series --with-lines

g_T <- gnuplot d_T trendt --time-series --with-lines



mle loglik = - (d_P - tr)^2 / ( 2 * (d_P(-1)+1))

    series tr = cnst + lin*tau + quad*tau^2 + cub* tau^3 
    params cnst lin quad cub 
end mle    
trend =  cnst + lin*tau + quad*tau^2 + cub* tau^3 
series uhls = (trend - d_P) / sqrt(d_P(-1)+1)





d_Q = d_P/(d_T/trendt)




mle loglik = - (d_Q - tr)^2 / ( 2 * (d_Q(-1)+1))

    series tr = cnst + lin*tau + quad*tau^2 + cub* tau^3 
    params cnst lin quad cub 
end mle    
trendq =  cnst + lin*tau + quad*tau^2 + cub* tau^3 
series uhls = (trend - d_P) / sqrt(d_P(-1)+1)



g_cor <- gnuplot d_P d_Q --time-series --with-lines

g_trends <- gnuplot trendp trendq --time-series --with-lines


středa 25. března 2020

Jaktože v Itálii počty neklesají?

V hojně (i mnou) sdíleném článku je řečeno, že už pět dní po zavedení "Kladiva" (lockdownu = radikálního omezení sociálních kontaktů) by měly počty nově nakažených začít klesat, což je v článku dokumentováno jednak daty z WuChanu a jednak výpočtem z epidemického simulátoru. V Itálii už je to více než dva týdny a furt nic (tedy alespon nic statisticky významného). Někde musí bejt chyba. První na ráně je samozřejmě simulátor (sám se zabývám matematickými modely, znám jejich omezení...). Pak to může býz zpoždění hlášení. Oficiálně potvrzené údaje v Číně tam začaly klesat taky až po čtrnácti dnech, těch pět dnů se týkalo zpětně dovozeného počtu skutečných případů podle toho, jak pacienti zpětně hlásili příznaky. (Nicméně se vtírá podezření, jestli tento údaj nebyl tak trochu "postrčen", aby argument pro karanténu byl pádnější - klesat to začíná PŘESNĚ pět dnů po lockdownu, přičemž pět dnů je nejčastěji udávaná inkubační doba). Další možnost je, že karanténa v itálii není tak efektivní jako karanténa v Číně (což je asi taky pravda). Osobně jsem přemýšlel o roli neodhalených přenašečů: pokud je takový v (doma zavřené) rodině, může členy nakazit po pěti dnech, u nich pak bude trvat dalších pět dnů, než se nemoc projeví.
Všechny tyto argumenty jsem už na síti slyšel. Otázka je, jak to kvantifikovat, abychom věděli, kdy to teda začne klesat? Souhlasím s René Levinským, že je potřeba mít model. Otázka, kde ho vzít. Podle mé (velice omezené, uznávám) znalosti existuje spoustu simulátorů, tak či onak vycházejících ze SIR modelu - co je ale u těchhle modelů obecně problém, je kalibrace. Pokud ale chci model kalibrovat (statisticky odhadovat jeho parametry) musí být relativně jednoduchý (ve složitých vznikají nelinearity, které zabraňují korektně určovat míru nejistoty, která je kruciálně důležitá - v těch simulátorech nic takového není).
Prozatím se myslím politické rozhodování děje clkem ad hoc: dokud počty nakažených nezačnou nade vší pochybnost klesat, bude se držet lockdown. Co bude dál, to je myslím dost těžko předjímat. Co je ale myslím jasné teď je, že aby byl "lockdown" co nejkratší, měl se podpořit vyhledáváním asymptomatických přenašečů a jejich přisnější izolací (což je myslím to, co v Asii dělají).
Já - ačkoli je to možná blbost - zatím budu zoušet "statisticky korektně" chytit vývoje akumulovaných hodnot, třeba to k něčemu bude, až budou nějaké první výsledky, dám vědět... Vlastně nic nového: "vědci" pořád něco zkouší...

sobota 21. března 2020

Jak to teda bude s tou epidemií?

Vynáří se spousta modelů vývoje epidemie, ale zdá se, že zlatým standardem je SIR model. Detaily zde nebo třeba česky zde.
SIR model je jednoduchý a srozumitelný. Pro ty, kterým se nechce klikat na odkazy, se ho pokusím shrnout.
Každá nemoc má určitou schopnost nakazit ostatní (to je to často zmiňované reprodukční číslo R0). Ze začátku je populace z velké většiny nenakažená. Ti, kteří se nakazí, jsou chvilku infekční, pak nemoc skončí a oni ji už dostat znovu nemonou (jedni protože umřeli, druhé protože získali imunitu). Systém se vyvíjí v čase a výsledkem je vždy to, že nemoc postupně vymizí, otázkou je však jak. Pokud R0<1, bude se počet infikovaných postupně snižovat, pokud však R0>1, pak vypukne epidemie, která se zastaví až v okamžiku, kdy nemocí projde a imunní se stane dostatečně velká část populace na to, aby se v ní už virus dále nešířil. Toto procento závisí na R0 a u koronaviru je to patrně oněch Johnsonových 60-70%. Cílem opatření, která jsou v současnosti aplikována, je R0 co nejvíce uměle snížit - ideální by bylo dosáhnout hodnoty pod jedničkou, každé snížení je ale dobré, protože snižuje maximální počet současně infikovaných, o které se pak zdravotnictví může lépe postarat.
Jakkoli je SIR jen velmi jednoduchý model, je intuitivní a - alespoň pro agregátní data - celkem sedí: stačí se podívat na křivky ve Wikipedii a v Jižní Koreji (která snad data nefalšuje),
V tomto kontextu by mě ale zajímala jedna věc: kam v dané situaci naše vlády míří. Je jasné, že se pokoušejí R0 ze současnbé hodnoty 2 co nejvíc stlačit, když už ne pod jedničku, tak alespoň co nejblíže jedničce. Ale i kdyby se to podařilo pod tu jedničku, tak tu je - oproti modelu - určité zpoždění, takže počet infikovaných bude ještě nějakou dobu růst, a pak začne klesat.
Otázka je, co pak. Nákaza nikdy zcela nevymizí (i v modelu je nulová až v nekonečnu), takže jakmile opatření poleví a R0 se vrátí na dvojku, jsme zas na začátku (března) jen s tím rozdílem, že budeme mít pár desítek tisic imunizovaných k dobru.
Máme tedy před sebou následující možnosti
1. karanténa dokud se nenajde lék/vakcína (odhadováno na rok až dva).
2. přerušovaná karanténa (karanténu přerušíme, až když bude infikovaných dost málo, pak si chvíli vydechneme, a jakmile se to zas rozjede, zase nastolíme karanténu, a tak znovu a znovu, až bude dost imunizovaných
Tohle není příliš optimistické. Ale nepropadejme panice - máme ještě jednu možnost v rukávu
3. Testování - trasování: Až pomine to nejhorší a karanténa skončí, budeme tentokrát situaci skutečně hlídat: jakmile někdo bude mít příznaky, nemilosrdně on, jeho kontakty, a asi i jejich kontakty do karantény, kterou tentokrát všichni budou dodržovat (hygienici se to snaží dělat už teď, ale vzhledem k počtu už to začíná být nemožné). Jestli se to podaří, bude záležet hlavně na nás.
Pak je tu ještě další otázka: až republiku republiku nějakým zázrakem vyčistíme, co s hranicemi, aby "to" sem zas odněkud nepřišlo? Ale to už asi předbíhám.
A taky jsem se od popisu modelu dostal pěkně daleko. Ale tak to u nás "vědců" chodí: zezačátku si uvědumujeme, že naše uvažování je zjednodušené, a pak ho pozvolna začneme považovat za realitu.
Čás ukáže, do jaké míry to bude tak, jak jsem naznačil.
Navíc tu máme ještě jednu možnost
5. Zázrak (pak bychom možná na ten Mariánskej sloup tak nehudrali).