From f7b1bc36224942340984c30f1039c00077131f93 Mon Sep 17 00:00:00 2001 From: ben18785 Date: Thu, 3 Oct 2024 17:54:42 +0100 Subject: [PATCH] added in real covid-19 case data + vignette --- DESCRIPTION | 1 + R/data.R | 14 ++ data/covid_colombia_cases_deaths_mobility.rda | Bin 0 -> 13002 bytes .../articles/fitting_real_covid19_data.Rmd | 136 ++++++++++++++++++ ...ng_synthetic_data_including_covariates.Rmd | 3 +- .../fitting_synthetic_data_using_epidp.Rmd | 7 +- 6 files changed, 154 insertions(+), 7 deletions(-) create mode 100644 R/data.R create mode 100644 data/covid_colombia_cases_deaths_mobility.rda create mode 100644 vignettes/articles/fitting_real_covid19_data.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index 439208c..b7fafc8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -39,3 +39,4 @@ Suggests: VignetteBuilder: knitr URL: https://ben18785.github.io/epidp/ Config/testthat/edition: 3 +LazyData: true diff --git a/R/data.R b/R/data.R new file mode 100644 index 0000000..104b511 --- /dev/null +++ b/R/data.R @@ -0,0 +1,14 @@ +#' Daily COVID-19 case data for Bogota and Medellin and corresponding Google mobility data +#' +#' @format ## `covid_colombia_cases_deaths_mobility` +#' A data frame with 1348 rows and 10 columns: +#' \describe{ +#' \item{date}{Daily} +#' \item{city}{Bogota or Medellin} +#' \item{cases}{Daily counts of COVID-19 cases} +#' \item{deaths}{Daily counts of COVID-19 deaths} +#' ... +#' } +#' @source +#' @source +"covid_colombia_cases_deaths_mobility" diff --git a/data/covid_colombia_cases_deaths_mobility.rda b/data/covid_colombia_cases_deaths_mobility.rda new file mode 100644 index 0000000000000000000000000000000000000000..863b3d2a89e2105aa0bed3513bac74a810efaff2 GIT binary patch literal 13002 zcmd^jV{j!v)a8qviSuF`6Wiv5FSc!KV%s)fY}>Y-%)~RXCYYUXx4uQ~R_))d+WX`5 zt=p%&@2T$UzSUwnR=ljD+SF>=qiopB0Mm@`fB$cG!0rbC@ID!R0YDrG?F1IGfrhChi}Bnv9K)h$fBm0a0f2WFNA9)MF9XHK*9i#e+>YdV?iNs0U@by z$sWC^oTf}^RS;|jS88Q1m4%0IeZjF;2+NlFk$r%7ObrB@v4flq?R#PZ14eA9K_|M03nrxBcX{84D$a! z4gSBGFiGo=SHsI6xOxJ8!W=oLa9JlFtR`1A+eS6!TgGoKB0u&PN+mbTvj=7# zh(v7W2$Yq9k=Phy;2=uu=TDMIdBZxIFM=W^LGKgU0`aEk5Uk8UA|>IW;VbcSFcO2q zV4Wl1$Y79Qgl|!7d`8I~RhfVf0ul`x>HlnE(2YfaFjggEb|F&M@K64N!E7Eh0jwA>U2B4vXyFm0)|BWqbebtIW7q}0+D zkb=p*%$nk7!7}xvDeOyCd8H<=Q1(YEu3!w@o1Wtb?KR zkJxYwLUnoC%hAdW;Qf3F7HM+Wcdj^*gXt>e*I3VG@iHF`#ByTxz}~OmM;vaux=Hwx zei+T}J9oJ@l|T3nk1Rd+eYyExxsmjv=b7u+KSX^P0mxH33@Mru&5BR&i&bLf?kdOA zpV`93LD!8u4M1SICKgFY4DtJX*YXQtRc&3B=qSy*=gIpw8kMfAlf7obdTxy;zx|9i zYqji?WW4f+qBuW7AQ$0-C{l;o<(;{RCV&IbJuDH;d_&X-)a61scOVgEBoH zz8cAhb>WkcUzUITwDC84AA~aG8(&WvKm8I#%=@*M-sC$zLxo9F3QE4u*i>bm9T4nn zeVL$+L4{5ooZRB?{alwI2}x*ODj59Si`(;PQ+TUDVVc7Xz z^w{otV*Goli&EYo9pFmxMpSY2>vimbkN<@#<@|a&=Gd3b0&V~N$c zm^*0n2EBHdIDqNvjf;IqSf9U;;o}2i8{__Oq8#RmW4G_wg8dUtj**fIYCnO=@n8+)3-|3@7_3@u?1sPu5eH9X$G}~y`TO1 zzz*Ic-_?|iqq!H9vfhq2gr9pAx7^)n?pO~Q(p6%0G~f4}J4Tz~hF{y0cJ45Hwy!&v z!>$t*FFbcQcOQA8KEGB}yco88(%{&PW`Jaai1dTdE=qU%p}Oi@|Q*Jj#pDx@M5?$7H+slZ|(;o`1jQ7u7OsAd_>cKt&qX zM}g+?{p?e|rthZnr4>@wT@3g0Bj1|lSoM6wzWw1#iuK5QyZaW04yLxv! zwEd%gO1a@`XSQ9Dz|SiNf|a@QKb)@~RmGKCZ;Zr7sK003ubp3QR(qvOzklZ;Y<(v=8K+bmyNz= zjq^NHx6|DcqunmJbw4)duZG+^9v8LmnvE$vH5M>T1S(Prv#>c3L*E(yF7z z^+Ybk&3N6KrcDHvC_P#_V5S{uozgC|TG0BGt%Vv<`aU~wW*ugSzM~GQ(>ZHO;ZY+b zXzl7xZ-%_}$3d}!ELVsV+>c&g);ywIP&o`NpX935OW zd1pp&w1OsjJM@rLgp$IuWP*(|!Shg`h;f5z??O_T}dKGmtggdd9)-{^>t^MX08do3E^%vJkm$ttUf>VKf(tb7|&&Q<=16A2lO;n2Q&hNH<1A z_i=vJka*LI{YJPGz{?Kc^UyxN&HY2U)s(xf7>i8SiXmLJ;wRfTHs9yxW zNwC(*u=`r3Mt`6jU}&xj;!vU}yj@Pt?Qod5>WIX1RQSeXBD;y2B5~G`wxn(+myzsY zEfNho77mb#{nSHGoK#kFde}^-4iq9%`VvDR-%sHfj03PV-aR)f<^e+LUBGAv$Bgg5R+776=F5J0^!I$s_UHp0_^bv_yf!n4-g zY*=85w_g-!a~WzZR+R{=5=Z(45__8|=f=En;k&VZT}9`=d=4PYcF8G(i9n z4%0H0yR|(eD9qZad2%9(^+bZ4-o>s4q2^>Fx{?!tuE|}L{)I~7J~HLkLj+XQt&$`z zNkYB0YP}?lvRx}CPc$%js9{cxr|Y+FIq7GHDSjmA4`{_->yd>;oRcAst4OV%FTX~o)t4%GOZVR$2O(t3MW4`=)Ljo zVxLY+vp3;NxDSJOtm_}Pdd2NwSord!qCe`aXa3deWlD=piB@F0)QdILgzj; zeuY=@yR~c9ULW#wx4p_te&x0-#H%=}hX`7$s(+gz-n0wtcBd|?)(b$DDI!N3UOb>|FU?!Lj9+otN?5;r$xj@vQ z-fp?Q4H_c5MHA|2bS(NvlYurCVj3$ zPKoAp)#SMr$R_CFWvYc*jL{7RNwN0Ou5&A>p7jN0TuLDJKm`aE1ze8@Z%;TJ#%62L zE8|=@91C?DJgS4hB2%(L{gO_fCZyw{c(nx5b8<&2q?okq7_rc-N$WsHb!ccYju?E= z@fg|P2|2g=py&{_BtRHYFJ-C_0Jmy82vg@)oDQ@E*n3+^+|2T50r3nZG71mMv=fj0 z5P2^H!mCnS$>cTZm~1n3@RiS$Wuh5aAo1(lDQde3POIC$SKRJj@p?cH{7aE|6Z9uM zs9Vb8^-O0*qie(95f;5GHEJ(17!gQ;gF91;h<7((m9#8Q2zh?PXD}5>fhahOJ5P@d zkr))X!VZY+PzbELut8BOGo$~q!bULO`ZnAnewc{6xPj~2%3gEqat8IR)a*! zpa^@$i7&sH4s?jJPg`fH6SgZA@c zAFfs_Q)#Io@mV@$QW{LN>WsO2i}6#(FJc}B)KFTamA3#xXnTN2GI+6r4D7d2egT;+ zo#*@9tS7(6o(jhY!A@rFNKT{>RWR%UMJZuAU!0>0zoBAWFKX646|;CX2!AY09=3o_ z>b?S`4N1QW3i^$y1><7o%%JNP9zX|q3u&Jm0mTSI7I;q@X* zDYgbr)hH5Up00$Y6mMuEY*uSR0myV%k(r{q&SX69PH%bdceB%hv3Lv{x{P1S34=QM z!-41qCYcb(+>Be#ood0VIY6plg5NI)L`VJQT(Z90bjJlr9LEhcgZ&x) zheZ$-JV6buFt07PETY+K!US$Ch&`#P*^hL_veV~Ga@E0sm3)H)HqaI&ERAjT+v^v- z3-ei7l#REfq~7;a8GN=G{2wJvPf1UEDvyfGpugq(d|nPV=Q^ z^qz{_w<VJ7Qs&p$-^t8!S!kx7Tt zx+5jHc%6vmJenyjIW`BzESI_?$CHinNwpJaVFyC1z(8HXMty-ZIB-kr2u97wBQvZs zw$@`R^WZNjeMJ(@0KSYA273dELxZ?KmTgBU_bAX}gg!~DKS8)MTtQ)*4_wCsR%siu zrchbG_dGe%e)lw|M!LcpS;wK0Ku1!cPqf6#Ff0Grj8b#Tw)w(x+_ME2HA)@pG>ruv zK0sM&@N(8h0n<@EY8!@=S&lPx900!^FB`}+@jwrQv9qdXCfIa|P>Y1nn$Nae?*93g zaY}kG&@dqaev}7xSCr(mQB77h!c}>^DDcoNVq_VhMl0ED;3thw0Mvn!Ly?xuA+=+0 z?{W@;lYsF$7(y3I{E>E=WqjIVa;dyn^JsX-Il3pLi)GW$ylCGgmvqiS(Z>Fq zbyPTuYTZvnf#T{RCy6)fB$+6@T zrjQDz6ouoLM-0&uK}x(I#X%y3#_FP7+)mTcakv4I>6df&Ho{iCqfyqEOFib) zSinqLb6-=00cnIQJ&(KN?^b5X+%-)$PR7FlER7GJ!>Q3LsOy+&fiI@GsFuCH`4IQw zwl!h#GJV&`(`VH_QFu{1zNz#253F1QS$DXA5fmst`ty+yhP4|kHH~wUtrjt~>R1mp zC3^2$GKnC?YuMdgP**WIM)lB8zkB$EiCD*)_V;dGXs&~xh2*E$u4}3img)Uor-Z%Y z+39TYF(k9REoNjwEdRQ4Lqw{p%B!t62W+NA@>q)s@+pez?c$Sw4q3dIePf#>+IEBqCcX$Sc<9R1iDz;bUBZ-JHs-27}bQQ7PfqP)kiPl#b?i916g2ltkhU~9> zQJ(CfP#sQ6&6~0Z`A#%r~G?y74#~~Or3PGGMm@ftdZCy7Dh?=s-OI*DUC-? zT)p9MZ*Ip^7x!&BnMs}EkVPVfYn3dr5W!+WGo&go5O3|RYZ#YDb3CYnqnnn=W7?(}LbpK*j`qv%G(N^KO~sm_a6s3f(%ivwifzPu$srAdqkf%mwlUHr5i1 zxbD^yv~@Dp)YRmWNSlklXg}npp)NfJS57%=%t%t1d$1E}E2cJ#YHFCWXxW8}hY#Bm za3^`gv78uEtu$;?6P)#=WE;1|SM|8tu5+ENY9dyTy{dGo=`pWe#H2%h1eCX_8MIfg z9{ZJcs2I@JZLe;VguTPLbz9tT*)c_LSH@)$!?fEBwBr)z%q|`< zuWSHLF(%s2rzE7LMDCcOgb^>K4jsqDh8^J-(m;kEZ<#sWahS!%BaC^ya5gzX>Tn1u$oFs+80py64RGtEFUoe`UT9D!5!; zTFo-YGgM@zVB0q2zR!xi?c|X}dW?vG1RkR{1eJ)Aoi`xYwrc2YX0kU)xl}W%77l>f z+sI(lRf4z?+TrQRiz9Q0GG33KFb+H$=}a5*~GbVD~SFI(`TVwf~y z@vqRI)@(mqNO@ec)TL3;MAh|*$*Dx5Dbg+JvxdOzWgt-zHQFDzOj;6pw9NXcl1`lp z-H4D`qT+O65}Rle%KEV#>K@anQ;5~6(&p&sFnVUfkh0Fi5TnL|;ZNF_DEyZl(~G#P+|sIf>c- z$N2x5{r^hzKP3I9I#;VOM(j&L4?GmiFABlrO|k|sCbGGMp=nS2{hn_+dqmUns+`od zWBy3du|^{d>ucnxY1>X8mgnwx#~ZiozdAnW*S)l_Yy#@0ZRAws za%zBS15bwaSNB%<3P8ENHAL0ki29R@^IOi=_b*T}E5F-xRUMp6(iv_RBs5^?%|kF! z%dqD>41A}(ufsRPGj#P7^!O`tXfQq6>@C|e+}u}z-*!vS!j;Kr(L&*sNinAwDVq56 zV#tXt$j~HqiYJOWJ)UltJ9Tu-eqQ{Tq6d!?B&1NOrmhX6t5;*IH`UZ*F5G$f2|kBR zpLhN0v9){A_34>$(5-Fnb8WVyPcsLB)WLT zL0|d9#XYnkGMfL2939p7#;4!hJA#!k04}o6@%xgBI6U!w(Is^G3-N*A;5}u)sr(_* z{$qT@{ls%;_3GtI?#~k!jf;&hUn<;gHU1^1`@6}8G;*fH#L2pSp zUr_SAcc~Wf^Zhhu(=JUdX~1l%7Nq^d<%;`V?@uW@KHhQJKj9U^$k-rrkTz(@IJm}p zo^|2Qe=j|KZH=e1G0ox0BT>Tom^AWst5&YLYZyj<%@;GFs7SY-`;@GX1D}HNOB0_{_Zr#ZdIZC%R_3EXv>RT&;4&PPyTty`I$XU@Vag5RP&Z`*uQ=DVPpa2!~{bI z!OayG;}n*_h>DPk=VfFyG|UXAlvvDm*^bq6-z*=zRBO9n-f>pVCbrIsys){iSe7pUb4^e=z@JxnNjAMig> zL6w-eib|RYCYu{gf_91qrJ2;wCukpfb#38e)PO+rL;vpWGHrVBzhi@K0j+0J2CmYOE7Jml?VdmTQRbF{2ApL0 zB}nZSAivTzVWU^B4Q=*RSn`w8Z)_+JWOaAD3v>aF2p6F-QGH>FClh6dw<`fwV08e* zJ}d{cSq7AJvMJwO@CvUQa*t z?oj5`bgx_`vb5t@tr?%4HaUumw@Ysv4rqWMYfAIPT02>jxH$xX3?30caY%W|u5$dt zX1lG=@5CsC1rZ}O`=e`;AS@n?1@34A`ocF%7D;jlFIir{ z6xaCkA~(pkjRqb)#H*&^@plZi9X;JQUznuRs>=!s7nZGQl?!jyBw@bHg#P<~LHkb- zA1Q(t&#+n8;QXJ0|D*eV^Ul;au0Yg?C1mm7ai>j%F*U06Y0Ktxw0U!O{b|npqmG#q zpN~NBhzT8w79%G}lFL_6TidNeuUS(*ou5d*uz1Xrbz#teHKT^w<<(EX`NtOy+p!*X zM6XBpWE{GZ2^kk6mLK$Baa5fH`YknCZPnNIZ-i#3an--LF(*Of-i<73A-)!O-2|p; zQ8oeuU1jgW8F$lPpu1oD9GF(08_s)16WZIcXq9geTRQ2nAm0+Wcv*Po{>EfYv-;M3 z?7jcxox0y6$^S{52%uQV-lz$bg2WKn5xq}GpfL$erlnlil1szU;L`Y8`{UxuqpOd^ zzIWs2OKyO0{Pp$&O`#~2ME-L1ByS~Zia5|G+8EM9cRn+&A1sffJ6BHw-%k6e>v1c#MEN7sCm6_e zo3j7({P9@s68!hx?R=O#{w1#+fW)j=4#Io|q4R9_X{O%s?K)JeR~ci}%9B&a2PlBm zDXUdO9pCHz68~-fdH3_)JAc}^QME8lE?tc8id8k4C?2R@nlyc9An)j1R$66-2u?+8 zR1JiLhDSt`J#cBP!#1^*SXRx~8Y&G-^|X&IA7_s&&Bhmz!8Efpnq0D)fX$0^%3^_} zmK8G@ZQ3}#aBBPOnKZY!>7g%nx$wq(~7e*yd$trKM$*06W z?9SiwN~DRz9AQ0s#)EtOG+5h=*gIsxF_uhLBL+L{x{uS;UtjY%QbM5 z&It{ZR=|~-%-&i8jQR@l$Nagh-i%FhqKe~rs9ljGpr67 zHY8g^B6p`Z!C$4K14DVBpX}%i!R7s)mt!{#AHu$}gm3sf!!sWFX%+#BP>pUk{<^nWv-m;UOy2Iq3vA1tn( zxT>t|W3&gz!5Ja?ZABM+tsD8x%YlJrzWvj)Hj45*lGp1U+1ZHFgQxXTo*dL(QV``k zueU4D2EDX*6ofqYutQ(_4b_dtIpNLZw@|^Xx?jwfgl@)3^X!OquaK(-R>pQ_+`nQ4 zQtt+!Po2%)lii@7W`(3?q&$-cV|5~QxV<5*8NIE2CvW0tp+^H<>*3Vk%Br+8-=db? zymcuIs&uDyUbi_$Wkznql3giumz9hzjxWI$4{Szk#J}UWSXOsDCK$Z^CtexEKONJv z1q94p1X(l{jlTD6Z)fL02W2}f!cIS{Q(Ics)WOF_wZ&3zg2d>v4v+I?H-M-w2#$)3 zB9|8eCylJ>w0)~*CQs@$vW;; zNT4yoktO~W>%lVXrTWnpwo}Lp+S&T(NR}x!T;m~bFM7z*$Y!xe<5g-ioduJ$u#^g9 zlC$Y~vlU#7@@Rnu`Nov&Atk~vyTtqa6F-w>%j$e;xuzD)UA0^!*==QDTgwA&_bDD- zzY3C96?|tmMp{(uc(0>?g&&xZ_&}FJP9Zj}4?*nk-w1nMTnyB^{<8ZVL?~pZJihAN zpF=2TQbE0i*hG&fNIdY4R(E(Km86yC%(?>4zNy-mZmo zQt^e!Zr*a(#5Us$P)u;Nchg&_Ug&AnLoLfm9GN&|AeBD6aJEtnA(LB;N`*aEXP3N% zbjt*ecr=~I6R@08HOTkDiQ~c*hhqt`k6brj+RR}id(`$K?RQPw*;SS(u;M0WeK=Yc z+;;QF{deHNSS^=#&Q2eRzk1H9Xj#$V{^F%Egp(cgkh0};QZ78-M^r_rJ5Kd2j!BAE~1=G1}D{*NAP$={#2VkT~_JBG>Y6-yd=PGCJmUJ?; z@6W7DJ4?sIQHd+6LN6UvP6XL$cn-L`9qWp8hSHeh^HuRjl|@jYk905w>(0E@RjUmp zs(b`XJ+teZDVMBP~q6@M5Cazbv#xv*ySW9leLif3>w zklXPNHYb<>hVE;Dpk%ITQ|@y-IWhT$Tbha^0%u86Bcii|tYt$ZJf;@mZUp<~eswdf zcWp*`57xO`F=JZnH>0}b?8CacLS_v3I;K@Ti>d%hVRsOVE>>d$g$@64N9^gCN`V)p z)6vN^N;NwO2NB6QKe~4Ri7-%r0B7+afkGpjYPRm=5(HK>YTq}&twgkkSk>c!(^+B6 z9Chg75U=foVx=b?f`}TYd4#lH(l(n zV~-}b-?Xk~9byShO~9ieWyKN&lJ;~&yUlPrM|Vh+e&~>8;fO|>^gs;*c?o8hR!%^B zKG@UBtKqRC2RBOX0u!;ts45Z|3pXq)EHexd+>Tkq6i}T=y}xz6Y_}{RY{juZF0N_qDl}|QdD1$9dZXTKS1UVU z^Uqe4i>^Fikz;E3mS`{-{C-AW#8pL!=6I7>Eq*xhheaD)i_VTtKg9RD=J}V| zq~2VjZw-Q7p~y4)ngML;!{fU@CiMr$N}n#N4Tht^rot!cuzKFz7Vr%TiKNYRQU+`2 zgTw=T8ZpAYOvc8%I&z6%KZ~g(W<6WfO9=$XMo?Q1FfmRS5Y(f}u*8yvI`t0U&MZ37 zP8+Z4?87TWoKgMjur&FcnmPQFK|IN~GqR6pF;6}Iw|Zix&RDN&t6Y7KO6Lj@4`d;F zhI@krEYteiZ_KiMNcI@Hnms5&;`y$MhS%g8yBF!+-WmcyV9**=kL|aWm7C)PV7ss1I)c5ij$wEGbBAO}dH@2mS<;NYek(kezU%s|oht4a; z6yyJ{;?(re{JJ($)Nx7cD}e(G8l$Il!P&<0OEn#=QH_^9-Rh7WVBN17&O6Kj>N2u~ zp=tJ&9n(X04d9Gv$w8!mRS30w3t;Q^_2UaFh0^XL(t~K*kdT}Pc18M8^G6Pg0f~Z& z^gQPEZZY8myP9bAHeF_axLkkGpK2cTG3|V&N?Op(`vScTE*3_e#tT<#?I{O~8`20jGx*;q=E-8}qtQxip@gmb$ z=<9fG{W7mM=Rxdwwc>sW%m{BlX+ExHkF+$$R1uvMt4Y#i6~`!{#dddN?3u~LOq8+@ z7x=vLSn9r~?wBlnU2h`adS-Z~elS}kFJdC=FsgF}$%;^$05NuunL**1CqV4j@><%7 z-kUjosbMke)iSKVK&3&GK%|CsxjW1vTwPhDNFnS~MiF(rTc{|pukz8f>Wa$qXj4%k z`2L(}k4vl_T-X&PtXu42gV|FL-pDl3J*R3b5eQApZjKz^R!&e`yMPz%{L?yDPiwFU zelF2&w$OcQ|>$QsD&T4m?gpxhl5HuU(-@!G6eAx$GbgiMFUx_H; zUZSBCa8RqUfzQ-fEh|yAS{BT~K)h(w!E*Va?W>oV18&iXxAtL+(t@>Ys|#I!b0vLq ze&Xinc{x?_Kn^o_n0Wh>)6@rKh0(m6Xp{wJcUW~Y2}4&)C#=z}Uq5mOwQq(ZMa(paXY4R=eG$_z@R+;sY*EsroF;;3fs1;zS zI$WUY$$`Il$6gpH2IoW%H52x(AaTu@mZZpa1Mup7loz`3(LfZ-z1a%vt7RaK$MLJ4 z!(*NEPHeIn?+ZxtnE5o3PESrjft7sW$kO{Wvl}B4kjj4XY-k-b%zTejUPUxou+H+X zS3$6`TKnKeX*lCGFZ+CTI#Acljx2$0%U);RLonuo8`)*(C9IN^mPWE7h2NneUnH?5=xG4JM7b}@seXGc{`m;6 z?Jw@m_9cA9XPE>$LZj59$%}_5)p-U7iQ?Ods2PBB&JGqcf0bSHlHpLO#VBZmjbHHL zb_XK`6c+fKIL2-FsKkCpk($xkl|<=i#fa(BM%11-^4(QZo>M!^eAvXL7B*i5Ph;il zcw`lAWteli*I}!IudI&;;e9qrgMqQFNPFb+vSL}F*fWY*JT@Zs5e}A#VFU~WXT?0P zim#3rS1+cFlq;TefwxRrCkld| zM6xPG|0J#t74grpXy31`vEV~x^_rB?L}HV70dS%V>u0ieAZDJ#% + pivot_longer(c(cases, deaths)) %>% + ggplot(aes(x=date, y=value)) + + geom_line() + + facet_grid(vars(name), vars(city), + scales="free") + + xlab("Date") + + ylab("Count") +``` +We first estimate $R_t$ for Bogota using only the case data using optimisation to give us a quick set of estimates. +```{r} +df_bogota <- covid_colombia_cases_deaths_mobility %>% + filter(city=="Bogota") + +# generate serial interval for COVID-19 based on reasonable mean, sd +mean_si <- 6.5 +sd_si <- 4.03 +w <- generate_vector_serial(nrow(df_bogota), mean_si, sd_si) + +# fit using optimisation +fit <- fit_epifilter( + N=nrow(df_bogota), + C=df_bogota$cases, + w=w, + is_sampling=FALSE, + as_vector=FALSE +) + +# plot +R <- fit$par$R +df_bogota %>% + mutate(Rt=R) %>% + select(date, Rt, cases) %>% + filter(date >= as.Date("2020-04-01")) %>% + pivot_longer(-date) %>% + ggplot(aes(x=date, y=value)) + + geom_line() + + scale_color_brewer("R_t", palette = "Dark2") + + ylab("R_t") + + facet_grid(vars(name), scales = "free") +``` +We now fit using a fully Bayesian framework which outputs uncertainty. +```{r} +fit <- fit_epifilter( + N=nrow(df_bogota), + C=df_bogota$cases, + w=w, + is_sampling=TRUE, + iter=50, + chains=1 # as CRAN does not allow multiple cores +) + +# extract posterior quantiles +R_draws <- rstan::extract(fit, "R")[[1]] +lower <- apply(R_draws, 2, function(x) quantile(x, 0.025)) +middle <- apply(R_draws, 2, function(x) quantile(x, 0.5)) +upper <- apply(R_draws, 2, function(x) quantile(x, 0.975)) + +# plot +df_bogota %>% + mutate( + lower=lower, + middle=middle, + upper=upper + ) %>% + select(date, lower, middle, upper) %>% + filter(date >= as.Date("2020-04-01")) %>% + ggplot(aes(x=date)) + + geom_ribbon(aes(ymin=lower, ymax=upper), fill="blue", alpha=0.6) + + geom_line(aes(y=middle), colour="blue") + + geom_hline(yintercept = 1, linetype=2, colour="orange") + + xlab("Date") + + ylab("R_t") +``` + +# Probing the drivers of $R_t$ using mobility data +We assume a relationship between workplace mobility $m_t$ and $R_t$ of the form: + +$$ +\log(R_t) = \alpha + \beta m_t + \epsilon_t, +$$ + +where $\epsilon_t$ represents the components of $R_t$ unrelated to workplace mobility. + +We now fit this model using `epidp`. + +```{r} +X <- tibble( + cons=rep(1, nrow(df_bogota)), + m=df_bogota$workplaces + ) %>% + mutate( + m=scale(m)[, 1] + ) %>% + as.matrix() + +fit <- fit_epifilter_covariates( + N=nrow(df_bogota), + C=df_bogota$cases, + w=w, + X=X, + is_sampling=TRUE, + iter=50, # probably too few iterations + chains=1 # as CRAN does not allow multiple cores; should run with more cores +) + +print(fit, "beta[2]") +``` +This negative association probably is a result of individuals responding to the COVID-19 pandemic conditions or governmental policy. diff --git a/vignettes/articles/fitting_synthetic_data_including_covariates.Rmd b/vignettes/articles/fitting_synthetic_data_including_covariates.Rmd index 1ca9deb..80276cc 100644 --- a/vignettes/articles/fitting_synthetic_data_including_covariates.Rmd +++ b/vignettes/articles/fitting_synthetic_data_including_covariates.Rmd @@ -19,8 +19,7 @@ library(tidyr) options(mc.cores=4) ``` -In this document, we explore how the incorporation of covariate information affects -estimation of $R_t$. +In this document, we explore how the incorporation of covariate information affects estimation of $R_t$. ## Added Gaussian noise around a sinusoidal covariate-driven mean $R_t$ ```{r} diff --git a/vignettes/fitting_synthetic_data_using_epidp.Rmd b/vignettes/fitting_synthetic_data_using_epidp.Rmd index 24c3d64..4cf24bd 100644 --- a/vignettes/fitting_synthetic_data_using_epidp.Rmd +++ b/vignettes/fitting_synthetic_data_using_epidp.Rmd @@ -21,12 +21,10 @@ library(dplyr) library(magrittr) library(purrr) library(tidyr) -options(mc.cores=4) ``` -In this vignette, we show how `epidp` can be used to generate then fit to synthetically -generated infection data. Throughout, we assume a generating process of the form: +In this vignette, we show how `epidp` can be used to generate then fit to synthetically generated infection data. Throughout, we assume a generating process of the form: \begin{equation} I_t \sim \text{Poisson}(R_t \Lambda_t), @@ -149,8 +147,7 @@ fit <- fit_epifilter( print(fit, c("sigma", "R")) ``` -We now overlay the estimated $R_t$ versus the actual values. The estimated $R_t$ values -coincide reasonably with the true values. +We now overlay the estimated $R_t$ versus the actual values. The estimated $R_t$ values coincide reasonably with the true values. ```{r} # extract posterior quantiles R_draws <- rstan::extract(fit, "R")[[1]]