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 a 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 P (10-28.3.)

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

Co se týče čísel, trend vychází n(t)=-4+17.6t-1.4t²+0.08t³, 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 R 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