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.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
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
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