From b1944996b45a1622d6a04a160e6601ad18be4093 Mon Sep 17 00:00:00 2001 From: Zachary Susswein Date: Fri, 29 Nov 2024 11:49:41 -0500 Subject: [PATCH 1/6] Add a day-of-week effect to simulated dataset The current timeseries doesn't have a day of week effect, which makes it hard to test how the day of week effect implementation is performing. This modifies the simulated timeseries to add a large-ish day of week reporting effect. In implementing this change, I converted the timeseries over to using a `data-raw` script for reproducibility. I moved the original Gostic "Practical considerations" dataset into an RDS object in `data-raw`. This change makes the dataset tracked and available for re-generation of the package dataset, but doesn't expose it to or ship it with the package. It's not 100% reproducible because I don't implement the git clone of the "Practical considerations" repo, but I think it gets us a fair but further than before. I also clean up the package-provided `stochastic_sir_rt` datset. I rename the `obs_incidence` col to `obs_cases`, add a `obs_cases_no_dow` for testing, and rename `date` to `reference_date` in anticipation of a future `report_date`. Closes #90 --- .pre-commit-config.yaml | 2 +- data-raw/gostic2020.rds | Bin 0 -> 9024 bytes data-raw/stochastic_sir_rt.R | 43 +++++++++++++++++++++++++++++++++++ data/stochastic_sir_rt.rda | Bin 10678 -> 5596 bytes 4 files changed, 44 insertions(+), 1 deletion(-) create mode 100644 data-raw/gostic2020.rds create mode 100644 data-raw/stochastic_sir_rt.R diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index e3b32e5..c7a0fcc 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -45,7 +45,7 @@ repos: entry: Cannot commit .Rhistory, .RData, .Rds or .rds. language: fail files: '\.(Rhistory|RData|Rds|rds)$' - exclude: ^tests/testthat/data/ + exclude: '(^tests/testthat/data/)|(^data-raw/)' ##### # Secrets - repo: https://github.com/Yelp/detect-secrets diff --git a/data-raw/gostic2020.rds b/data-raw/gostic2020.rds new file mode 100644 index 0000000000000000000000000000000000000000..afc144b78d94a6eaabdb2020d64d3f3d45dcb796 GIT binary patch literal 9024 zcmai%XEYlO*#4E4N3E8|s?nNltVdPNBoC@KEkbRz_oz)s)zjJ(ReQI!S5?gjwW>xa z6??=iAwozbB=2+nAOGjP=e<7s?sKj$_qo5^=emTG|DyYUVmQc;=r-D-5F1-(SE3F= zwl@T#lqrRA3P0aJ;0Tet%bb@##;4Es;1g?}{^Nx58?Jxp>FW)d>OKbctcq!49B~S~ zB1xA~a#80-a0p5mb%Oa_8O7nfy5HGI=>(&9I?)lz#*bM8J08j}b>9lBOdGF7@3cHz z5qtIcW<9Wik+9E`+>T|5WyT06@7*`N#?{0QyFMF->Sb)Y=|f-8$K3>gvG>p|_U5cfH3?$r)cVkl zJS0XeV^nXOqud-*4+Eu-xW;rOs!YfJ8XgOeE9t#^q@c_ib491OaZM@UmH-1XZexuj zfWGjm_S!9Fwux7NC%kfcO0TcMr_3@Dd$LAP;aYRl zKX1&S{V$mB)}VTCB0tdofTv!d|1Z>w+@UERJBa@OXa%xc4k}Ag-T?3X6+BYR%LsWj z?s$cl{*b{c6SRJd#LId$TRiFld?f9d7TBdB+$-Q*l7Zbpy|%~%x#hd0gYV5(H0YoG zFvz?V2oKG2N#~WDM+Kl>`0Jg$6U_uE+kattphFWIIaGpvp~z-n?LlH)6jbJN$MY}3 z3*z&FwwcT1e3!QNd-Fza`l0>g%;=~>7t8H1z(r=!5N_;KX2_>Tmp5H{v+Fiu`vs3O zTMzX{tVxM3Mm4flmyqnbvr4hrK8MWgEKw``sQ8I#*obD>j**p zepxL@EO?;yNTFCj2ttYe$LZXAm<+W{lBWAnNMryIXj~&|K(qjD9f+r1S*S_)lG5g!yjt~ zuIlRAi1EP@_xU~riEV6Jk0?JdeGsEA%h#2oJM^tme4~Xw6W?b%@_dML?XAn0^&P&f zFoU5v>rD9zt&v?^M#?Z*e@L{`<*T4H-#GX~$`c;3jkD{S&pN7IT;O7STMA;LF|{rp zmEt6{D5vzf)oaNbgQEa<<`|^G>n>Hjfc+>=+nubn8Y(Ajr|Rk|uVf!zH76{v;Hr+? zyR}JvPS~*p^O>?pU;H=yrt!zjk8EW6hAhRO?-yRJg^6l>SK@8I)AwG0IiEvsZQ%v8 zZJkUkhA0-`)O|JWXudB;tf{TyYUoTxFSmGX$d*_W*^^oNwp45zHKVr>%$XbWm|hMf z8XFj6M0tFq=c}#(ux# ztQzpwy2hzXv!mBb7GD#Txe}d9qRro<8rb7YWNyI2J4q%IiL_l^ys$AhtSr>1_!=JmOpH2i+XlA0@QzA?p9BZAs)M^pi)&0of*!RZj4tYV zai@tip@xqABoeN-;cdu(15cM10!Yml9B@EJ&0yQGOppIwGH(ASW-;(NY|%RjH(I%c4E|6b=i}?V&uvpO*V#YI_sza( zA7R&ckgv8_oqoG;A?WA&rc1oSlng&|x@w(2#G=A)_ozAzJXJlY8uI;>4x3kp2HVRr zi0zD)lJ$_$_(>_i5Pb4|-`*=-dpW71SiBdVcoe}IBw;jmR|fflZsYQZ zjF{=(_Q?A@GbXWHtycXkwYs0p$}&3G*lxRFFxqcePO-YWt^Bg4rN_&o3OEO6`Nmyt zbbZn=0#=ed^h1NP`4kR(hLc)v)L&}jJe7}+%+#u%_(a`&Y7$h1?3np%8*C@L0OUEy zN@oA@>1=z#1&kfJ%`~YQ;ZdzlFaw{;bHKNrEuKs=MC@j$dqse)<%E6n(y605TH9ug z$JO5FSUa$)gfJ3N3Sp_Lzef7D_MU-Ye5JGWUHH9CisttHfzO(=mao9~UEtTN(AUPy z{#llDk$E$-PJh<+?(|lZLSL(I!0C-V47*u)cPf zlL{Vr(&kzb*fMFHM^abVOTl_d@XKUH3W+pm?hf7Hm-*Ko)jP~L_byhJqrCQYHwUKT zPCoe6TjsS5I@P)K34F^Vi0XViEqWl9V|>qf;&`{`c)&n)LAofa5<++6&QCPfqI|yD z=((Hs604;$YnSFddNXHh^4``_Adsr7y-|rtzqPIV;^00>>KY-g&|`iDD+uGh+kk zz4b~1#EuHWG53p42f5S&J=&FnUkKzsN^4KiY8!5r3mniBkoPE4lv%MMFmQ+3)S*XL zqB>TRKNC5T_M`z8Vt7Y94+SxY`Q(56dGL5etijM}%mrom6=Zu~qS-ytc?15^!YJ|7 zC6pty0{4XigoPkz%@m=ooQQ>W*)GebkV{GSo|g{~s&lE`4ICo?HowWa+o`0fxiCTB zYK2N;|1#&CG(W`zt+?=zGC@1SEkE)DCH0NFkJL%O26F6QE{mXK&Qyhr_3or`TA4in znFHmYrgiJOuU>`?IRUwDOZD=o5CpCfc%bW3Hi>vukLkCYnXe@aZt}>i&&-l8uPB#+ zE}unRETmCuWI%S_STAjd-bHk5w96u+SGU=q0!97!9*6MO&76*J!U;qHa*e?D#EZtDcG=+u2EtVKE=7upIpT@3z7kXO zPow3Pp6XtEsAI8DxI_VHbZFbhNs$HVinrsA$j}4*{74%SgJ{;>eQf_GQd#{i?m@zz zJCxgEC9emd?`+u|y-V30mt(xjOAthbMh-RD+)n@S*80$KfUPRXLA5#07-29x8yxxN zA?Y5Qy{@;0bbtI-dxHSNLIw7$!3aC|NuX0CJLKcB9kB3S(yHkk@s^i6u*c|yoN2Dy zHA+_EzRQ_-UZK%v3(2hM5z<4?wtjb+jXE0YUI8CGnTRw-kt5BcC53hY!-hP5J5NN)|R5{QblN=#)Hu8G4u9|MuX0t6v zF*4CiuID(@!Tt=DY>#(Bwx#vEUrP(U>tG0T`vPuRgl|YwU_ZeQUr_I`L3b&Im2dhj z&#Y3SW1LFhb33=r&c9Qi3r8c|xA@r#9%7#L#}rZaP(#Wiy>Q0XfjL&;n<&}f{DLW? zhF@BVBRHe8QG;UZIU${tfh_mLuz@}8>nHC;k<;^%!X3hQEHvjbKt*#GX>$@Wetolh zgyXi2HZtkR`0#fMq6e<{AS0$*D)2-(;=#Y&5kkdzPl{3Tk}Jb$GT4Kdml$PA3Db=pWw|gl&jfszf(mZgiT)w{0sh&E?1!6?GpsTc}=mI3L+{XDa z@$?o*$%J5|Rj#xo8`+{*4mo0TRl*s)wo#tmeveF&C`XR9wwOIGf!Gl~wIx%j znp5cU0bjPqw1r!jUE}4nTQ@H6ZTt}`FR2wJXA=ki9!k!5YLY^*b8@NYsc=l^ZS%WP z5qo!NSwr4hk*8vVe`&=krAKR!cR&$Ey+1@6b?SyaNkRV_e162_ijkF{t>fNkIdOu` zYo>r~29-S$@%t|FNoPrcv!X-ULYN2CTb3Bjz;jA+Ek-sFdj8b5RVsNutpejUWoye} zdy6OX_wTNZcE+~@icab9hP#RiA)lJ*ZYk%siWB=3!zfxbAI2E-tStP~`K zkZw>I9|$S#XxL6%Y%yFK=^7MF;R-X34z!GZw$yr#O@m_mR;-qMZ4xN!o|wb)_VyC& z-NU;)10b)OcF6`F;#c>9%gFBb)t@(r55P8;E$!$7F5;{{W(6-7h>NybVk$_aB z3&A`l!~K{)OXlhz%j6w=RX8Jeg8SSZd`e3mL(dArE3TX5*D) zB_rZC-t08nAmN4xaeCw>v0a2qhzNU~JN_Yg{(;dL#U^?F%uIff=dSP1sl$8nT`*}k zqM{S{EHyXUBH?^s#|hkOI2L>177!*MKbUyVN{l-6mDgzloSZK^r3uS!vk((h{xHZ_ ziY2LpG~`YP4HnL8k6Gc8&aIgVW_n{v1xYIKSGoR1H-e^oa-`G!h62y&qm&Vq+qQO^#7j49Ft$sJo1_(V(6b61p&9leWx7ztvttM1oOEv)E| zRm~qSEs>v8deYnW)Hk)N~iJRiGetR_G@{CLj z5?n34Dg7|>^xq&vG-_B-27<1R!>N>1T#@#w(jITdMr|;$thBkQphE+}Ex3{Vy5AG( zGomad&k2UJ-V-ov{F^7ujDzH>$)K-{PQnYr+=4sQe`FL;O{FDhB0?3Pg4 zU)rBHved7CS9-$i7%dBPzG22lm7mOF>`|DMVZ2P&fZZl2GE$t^;0%kb$G>@*j|w#W zcn8vBeHqCEN5d=o1eod>A0bYaI5ovd4r6~_kk&An2oh8~ooAmcRVAXQ0IK`5QyZ#L zhIa-|iGd3>KnGm_MujrKs?KbQ6Lf@cq7?%V_+0RZo2Q*V<07a%K0p^Vos&GMLbou2{La z!_qnxI-m@i46KK+!4b)W}5rw`;t5%HUP>OsTJ6UsIWJZcfl z5bLo6SGE<2W-yxaIAwYa$Jac~f{Xf-gCIEBEB{jKQ((+ zlf`DvaXrWzS#i*)j1;%1Z=KlZ5CpeG9C*@6+7KjWEgcj36K8I*9Qb9Fa0|FhaalzR z+Dq*mc|vytvR0P5Q0w*m&A7&(Su}gPb54%^PYTQjiI{=N1$tn3rnkxOdB;kez>lg8 zBot>p7UvG6M(CvU8@B%#P;v{yFj3S8KDsMAjte=BXh}!8sjN0bgtiGcf{OD8mVfic zM;x`?E@}S*St9*5f~*ZaQo7iQSziBpucM$_n_cZno*?IPD?&JJUcnsEhx(vwwGw^?*_^0#su)4?QaPlenE{^bMV9J0yz%)X}@nT zT)5Ee`csjvN2*P{nnQ~zk-;c!yr}Q^Kq*_aY&VoqE{0(Jl!5Do>i(n_+63*Id=t~nYFM3YF>8J0I zQ@@rR_O6^^49=B?U=m~-Qsw%wKHyQ;XmVG`bc%B;A^qWOo^j7-XLr?88NmkquF-Y5 z!z_K{N>dZF-E*_al>4N6%*OgdD`J#$y3TtyQ!|cqu!ai7?MV)a``_0b>xPjl=*7-~ z*%$Rr)+1d*BVC^+BE&!qVXo2+2Ql-DgXeypK7uPZClXNI17&EioNkyhwkKDp#b-R8 z){RqBe{63ZfY`#cw8(?DmvOt|vwqDc+!3y+0oX6TAJ3j+Yv0fs_B*rLi3`Z)Q!!B2 zXcItpnOT<_d-%}c531uVmS*Niix%aL;>DS^1-aEqG1Eq-|iTqiyL!;yI$0qr6 zPT#|}b>pY!+_UBs=b1PbzV{bLhAQ)q#ely~O}%*ZM6Hwcknog6{Y~koivF1YjtbJr zT31u*3D{8h4_qaFXM=Ah!qyFP#?@fBFHy{dw=6w@Wvh0inZe@ry;wc=undT-Pw374 z1_363px{tU2T&R2hH@?$ZE|zg#`+BX!^2lADC6e+_t~#fdw$Z2<2a+O(+wp?@{c=~ zte787D$F*l?WQ!D&m)8^c6c+54Ad#PnXcMgm;Msm(#qE}H6mG*|Jf`j7(1Q#?9yH* z`9+t#_>hU}w>>L)!Owuf_PtKomJ(d_&aOfM+IVz)<;cxHyhW*d{`Lx%&HwDJQdEBd zV*NjxX5sHalktZBbxB4RY|u(7Iu|zovB+Y`+xs~HU^_RlTaYK0304SO!1Fm+m4w5e z6nW<ztUdItql;nT~Ay7y5s85_p|QJySXJ6;{B< z%Y5~`W+MAG!q6#i_yu^)z9?|R3gIU@RJ}-!_DmP_4(1O^kq?T=V$3op?8QoisDk@R z$Ee#>*!$*%s{4s4w2R519Ob;gi?Kf(?(%PA{CJis%i-|ghUxRO6ynw+l1O>EseVBG8))!P0{c2+vU#wfb=NVAb`)j{5^y}Kra=~JV>C|SsAKHy+rWLXnHgEq_$Bq)QZrw7je`>P* zLUE~8$L2W9euC~{HIsaA_Fj{4${NfY%6;AdL^sQFCuN($RJ zbV_c#EtoV{k<`PH1A@Fr29n8hHF%X;^QqYwyGL(n7TP6Pn|Z8KBhXqg4DG-082UBb z^?U4+R*i@K8=SE1j41%o2_M)cgP>FI6LA7X_#9yVSKbtD4REl{}q|?V4ZC!$fzKdc)5T9S_2mdJMB=$ZH+liuf8I zdsKy*#d}NFM#Cp}B(BWIHU$5j{T*>SHhYW165HX_%2HUuSARtPa)54rSlme zboRx@3|=b6(zRZ$hwOULXg-fHp&)ffEM3EN7Q@bD{+nZ1regb{}D7)$;0Y~dDo}b;_CZ-n8)CT1N4LV+;`o=V& z+z!6;Kiny`4~V3dSPtvjUmm~J;>}$zcXloXXiGXZ1oH_jzy@jg37O&*8ai_q+DWUL zZx8i{c5Y-BTCM6p|2)~yPfw(n|6w8;I?QJ7bDKz4cO6#&i&~uA(cEu=(GfbQ*Z7>x z|IEFmKn^GgK*9gM&9@R(=qD3G-HOWv5W3)+NlC;aY6!93)s_1b=AXvx?dRkVBA^Rv zZLM6~$I)#5TaAp;mY-BAb}arJ)c{W^I}V6o{p=x*>l1!K(Co)>+$3z(3i`4vBWh?y zs%CSB42kh2d;zNb!s%T8398v-+DC}j#?`p&^lYH>eayuRdU_rnmz#lTWbD%2%#7lgBSA+B0juX0zjdhWzF2+Dwi6TXMb$X5t7b z6Byi53VN)&eKsM<3jMYH#_r&-$y{&xnl&YhqdJqEQ)_GAY?3{*`_1{~ffF<{r_03y zVFT4W%Gyeen>N=w<}wh5#9~+t#>~WcH++{5 zda+;5Q=<9|3J?x<=s77~HOChN`Q-dMbO8R8lcz;_$SMpDuR*Qs(?QUO9hqD;eE_5|LgHaKe+%_BLbI%h3}1^%r`YLkg%??>4sV<6Im>%FR|S51DO)VCc9mD1>hSTwk(xVm`HwXL57%HoOyc=|X@|(XxgyA^tl=Md& z5OIN%5bOd}{cs*$AQ;zw6tN?MLtiX50*{aOBqh?0YIWnaqW6`2>qLCnxr?g^rAc2e z+(OXj1<=jnX%-4z!2#eN(DXu94ut&^$@V3x_HQ42V#65#I+qk@Z)J`V2K)g-QFf^k8l z5(OY95gRle_0j=lcbMI1ZgvwfPOc*Z%t_X>b5w^$^K~gAoa$Tu0SIeGBHWfaWAr=y}*Y-#n^;g6iB54Qt>-F}`;zBzut?oh*Y`A-@E z`xC2S4F~yEN+k~tyJNae1fr=gCOHB2;2xtY|wbb^Hc>*uUyUvjJU!Sedzucz+I)LgYm-4vB*0o0jiV{ad z^s_gWF_+s0=NBXfq$wh6 zhdl{JpE5xYV?j&=K47yJ@#>Bix+#W%tlZeju4MnP3tzsK%5zY}6N$Uc!-ptny+Gy2 z_{029Z_m{_;J^$~%egcBgt^H>L?rvxL^H?O3$3@rKj4qsG?*F|&gad>_pp43E$%?_ w8d;gjlWMO)pEJ8N=+gCvGr*%_de<634TWNH`--GWQY0}aC52J&@~ literal 0 HcmV?d00001 diff --git a/data-raw/stochastic_sir_rt.R b/data-raw/stochastic_sir_rt.R new file mode 100644 index 0000000..1d67935 --- /dev/null +++ b/data-raw/stochastic_sir_rt.R @@ -0,0 +1,43 @@ +# Start with the original Gostic, 2020 dataset +gostic2020 <- readRDS(file.path("data-raw", "gostic2020.rds")) + +# Add a date col for package use, abitrarily starting from Jan 1, 2023 +gostic2020["reference_date"] <- gostic2020[["time"]] + as.Date("2023-01-01") - 1 + +# Add an outcome col with a day of week effect + +# Use a sine curve to make effect large and obvious +dow_effect <- 0.25 * cos(seq(from = 0, to = pi, length.out = 7)) +# Get dims to match +dow_col <- rep(dow_effect, ceiling(nrow(gostic2020) / 7))[seq_along(gostic2020)] +gostic2020["dow"] <- dow_col + +# Apply dow to log expected cases +gostic2020["incidence_dow"] <- exp(log(gostic2020[["incidence"]]) + dow_col) + +# Apply outcome noise +withr::with_seed(12345, { + gostic2020["obs_cases"] <- as.integer(rnbinom( + nrow(gostic2020), + mu = gostic2020[["incidence_dow"]], + size = 10 + )) + gostic2020["obs_cases_no_dow"] <- as.integer(rnbinom( + nrow(gostic2020), + mu = gostic2020[["incidence"]], + size = 10 + )) +}) + +# Subset to cols of interest +cols <- c( + "reference_date", + "true_rt", + "dow", + "incidence", + "obs_cases", + "obs_cases_no_dow" +) + +stochastic_sir_rt <- gostic2020[, cols] +usethis::use_data(stochastic_sir_rt, overwrite = TRUE) diff --git a/data/stochastic_sir_rt.rda b/data/stochastic_sir_rt.rda index 526660facb829b7b8335b251ef5d07b1517de289..443792422ba05e1477ba5d18ea047922c075ebd3 100644 GIT binary patch literal 5596 zcmbW5RaX=MqeUq}N`=m0zVc2TP;!EXn%eAzp_IR^5*YFeN)Nev>yr``o2+l zHZA_xe4w?&;<}Ih%FQ$#24^@_S*`q=6i4bw04o9xS38QSY!!4IWZ^_DxdcE2Z5nHyt`Jb}Bp?@bg1IO_tBB{QqZ$AzD{dv=$c=`^8jxI-A&@dk z(nXGWy4a5(+8oOAa7XxJA?!v+*Nzsfm3}RsrwIi*=`6Ym0;Ce-j#JWfw4|Ugm+}m{ z3>oFz@)K8%eIYw(wMP;-ocSt8MjahVfxWAVfMY+VQRo2S%$(}U3!$vSId8w0S>{^7 z2CXWfiP>7x`dek?m5-?pYR7REb3nvG#abT6aRnUVuZ&n#cJkMRtF^>b+(aOrm4OH( z>WSk*8Dio@D&@rhD}5ptEDb;>9A7EM(TtGeM_5u~mtfpz2|Gbqvvci)m?-n~uRtQv zX)197RBEZmWhINW%EjenWsZxq<*r8=+B9VZC@wB^f&vdrA$*mmMF}QX0HCoV z6%nlJA?2{b61Bt%0ckC1Y4(7y-8^;?-L3C|6hYl4-GyIGkYpbJx@nKus}XT{UB4T)~q))(Sbt9 z@rQS_?;9m(Awr_4*yk$l^zc%u?PY@y@|DM@OO4(wAj}C>h*)0hy`3xsRImcUR9dxE z#l2j*3^vxb(=Ezipb&5;U;xyh#|GhNA!Ncbd*xv2ScNJ_Y9tgS)aiSk1dm}m{~0vU z>m=wt&P;o8j}y%VtF5bO*;}5l+PHZRqUOxQv#D7>G44Ek;&A3?@bT}bdk?oMu!~En zs|Z1%7v<}wqYTj^U&22qkxd(=m&HIGh@y250_|gugN}5{*UOKx1+C(5`ib7)Hso3y z0?rGPn67SWZ^MdBtf(C(XCLb?1&+yYXTIGMDSxvuEMRW$_(C<4t3vjuen(L~_-R%a z7>n`kJ30S<+#N~hs9U$bi{}y){QIg(8x_g_T>s_0oE#ANTCATLimzFd4a^~;E5%j~ z7VyPMkPQ$Oui`_CjT0)`16)ynxyA9Obd2sRetC7uo~Z0HEC;&e#P?T^Tdgpj;mShgB@#tU`XJ(nPcsFiP-TflQDfNj_67l)c2uY#g!$P zYzffo&h*7eVqG3TLxT&X~{cl^0rQy4hJ3*;@e7DX0#7?ed2=6KprphDPw8B%jfG_%_ z4f0x)=PrpPAcqSgd)9R%mfPSQ55&)y_070lX>;k-V9%u@|Ii0pF97u~*qOAO zP4yezhv8})oBnPh+rl0cAxzYU-fQP_3TUp*x_WY>#l(@tz3I;gOtG`2n~2cWs7C!A z(CX6LBW#2ciN@7szPxE`fJ66X!RJA!BRiTmfa#^gOQQO8F#%sShmm4KY= zc>jEv2O8j>U+J=P)x@)O8>xYa;R2Tj%TMba;QeKP&_>*Z36|PNIBPw)-kqMxOX@#OIx}!>D&|d0WcQ~4ajn0T3bL{nmxB*CE?`QQTk~v ztzYD)!JY%hL(%Xm!_MQ8Ahm>~xW~JgI%qxD1jlzwkA(|n1J3lHwo-v6aV+uS&tw`r zMFK9Lz4>?0;#-rHcr4Q(VWID30r9g_v!x{Yfw9WnX9(jkBJB0lULJ7)b*_d)enBP} z6EFAh4;9RwGJL0ZS12X(Q=0k=C}$8^JS0WuX!-1lA%7P{0kp z`hJ-+ToVEuiM7+xw%s?Gmju-V%k<=3A$pigQQ{Lj&TPS62OqH;w6{zDi2159-r_;N zwy;rotMv&|GdZwHrxkWYc_Fc1x^fv}iia9-PmR^p7GN2vx{$OIt1>lv_le_2C+r&A zo}o6aw2^k3A0mcpUzWZw?|*jYVNQpHAYdc5xjwGMhqTYKQ)hg}dn@XyZYkIX%pWt_ zXFlOC>U2aR_N+|M!pIcFB4JU@O|Y1?mHS{`{9XBjl>HUzG+butjekGEar@q4dfIn% zO|!0iw)_j`%98WN_OzMpaqVV7Sem~sRzz4&pm)E_LFvmJa$N5u)1V7Gej7BqChBTO zf!!h9uriImVEGWq&VC}v^^v(>>Q>Vd-TLz&0pTniaF`yqu(wK14W33b$?qj00tA2y zc{qpC$Ii)0td!usY~9F!m)~C7(Mp|8%Hw41*;Ou{a0UcWi7o9Utk(A}#49Ra zq19+((YWK73@su=&ZZd;hTXVHU4DEMMJ*5h*LYblZMb99g>A*DxGVu*+1xZJ4%GXi z)7DuH%XVlTCOY;H@DqGWO`N3yd;Mcd3m3&&KQ0Pav=TZMf{NRPNcF6cG@kNt-+(W? z=Jz;FT@`GZxdDg5f0ZJ;A=$Bs>F@i8lqwk}sKSZSgc3~0lN5NB`#is?GV_lN-OMd8 zh(&s-bn3-GIKRAY%^A8|z;J~J?zY6~NI15Z{TvpS@<0dZ`B%*vzm9cxO^!z+6xEJJ zMpL3^RhvY`MO9Zl6vbXs)@R#Q(x`sg>gCN~%VE`2jrKH3i>c|I77a4%ZcBYd&IY3a za60UqKw^`nZrz+n)C_;I3bWWx2Koi;3@r>R0Tz1_Fe{<10fv7y4pRx&0=WoK)-w${Yqlv|o3vi;&02b*6uEym(Xu^ZK9g$cbA|HdyJ7faQj|F(fa zB9rdsQ@<@OnxQ~8A{xw&nQ82+?}`eDk>3_1xdy}Eor((6{P%W5-St{BCi}8BLkeys zYtR_6Fi>Psu)h=1V~V)wt-8#LT%N1HFI}r$_R-8wHon_?99I}LxQZDrb4j>S!G00=nPR7bYhFV> zh^~g+{^eKXs@)k)b8xo}n!?o)WK5y;^j5s*^4UbZ{ zmTfuRi{J+tKJ#T?;D3wAJ_yiiz41D$Tj99ZQxongf{KjMWGP1oL_XU`Cw!YG=Waio z`Bc2A`_{aC?y6MWD%YFVwv!U{agJ1CN~oHL5bx5Wbw=`FKNH&5@7t$ zwP`p<@Gkg3Oe^FiVegTkd@PWbFMVGN2YsC3q{kI72evWD-~iK+txWXxw8l4!((20|c`aed4s_xA?tka8= zvjYG$C_|{-FE~G6=nwgss9j5}A9r0+$HHNTFX^Q6N55XOawmQri@~`p^6GtSEHb7b)l&eJvEA| zmcgXHf8?E0-=71QYnszVliG7gW}qswUljP8vY}Qw6dtiAt2JGo1MtCRLoMz7Dmb_D z)ZGn$%(L1eeBLeI#4Bk}b4H`h$-r#$h=`0B%c345jDm~4-r{O1oh||;mXYQ5R9I_$ zX_CwNkK=~VPqKznD zl>?OjNT{iYHQ!qbQZn{s@n3w#Q2}qsy(_kPl2QE7IsYf`wM3$fY^JmM(<@L3$TE=K zri~a@_$zjQxdh+-sFP^SU@Y6^^*cM(w%5(9cwe${6<+J{jDlXs8#oI-MLKn{j|wP!^f{Oh!~Wk`KzWbbg@IRKdnJ196hQxSqr7- zxf}UK#J0$eSA~1uomP5sEYTjR!oKP+TE$6wh;cnMHWGlp)}Bt3Yy0f;!a9ecBznxh zi#B(ps!)gs7uP(FE*bUaHu^-|x&L*d5%y{C zJ@Ozei-$L>5EQEJmys`DGE{_g6Q$s>FVimg`s>Sw1m4<}LI$c`B_?0XZnv#p;@A7=_m^V?jyZn{sLrw47W6+kfG&5|p>Hb*i)CiWe1Fa#x~rYy@4OUxw_E>ao!gVmqj4$U-+ViN z+$bW#punZpVHxfsBEsj6pSOt{VpUyIPEW*86z)~v^9`1whRAB%>-*g`FcZj-RZ0~6 z!}Ui_H&-J;)DhGU!R9OPH`*P=$x>ak$IelaD3lnv~?^LN6Vw5Pec_=hF_+~`Zm z#%)Y%VhJV4!>ih%I!6u2JSBworgk7jaueRx9WArf;-)f#pBkVOLB2*x7|2%h{!TibP z`MW2BmDRFY=AAdmbY}!~OwE99dMD;yj>&-ITeigfwieHMk0F-~lXtu>eRp(rRb}9V zzK=+HJaYQis6+~O{oh-er=gGc>MAL=6M$R6DV(LLwAFr}LMP3EnIz1&bcW$f zUUP3IZ1_dxaRTu}73H#HBRDz!3Oyzs7`f{-L-h&=-1_*!R-C-d+JyT{k=M)49-X{& zWq*71!|G)r*4B2$Y4h}`!dXr6mgB!4rjJ#AT}7DefaybEtDUk4tVdx$gf z3pl$si8ZH>aAs!u(-K=Ze$GHm0d5iy5A4p@-2Dw{o(lHx3%!}?T)$}~R9vEAv#-tL zA?r?wo`CRqLE1GQsmjjh27Y)MoBa;eM(=K)zb5_KTt{Fv|LjgT z;P62AGwRQ$Z6#A^QcIc#(v0ux=vjl$j&2S+m4-h4Hb#5%mD(CT7yfYIS@!AF})5(0_J?ZBHhf3a8= zNjd9#=I>reCjYBA$CwX4nS04cwfE4va<6=9)^FS2EO%GAf#N+xDp?!p%^F@uhZns) S{m%CC&O=E|H|4)??*9M`G|gH7 literal 10678 zcmb_?MNB0O%q=#!%i!)h$i>~A3tSuq8{A=VcXxMpcXxM(!C{b#yTbtgw|JXBdCA+n zrpd`^+D&^lZ7QZ?!OJG9O#{+S5i~@CAiDMY`hVFT=@#*ISDp~`gatuyGr(rP70nw0 zf|aBn!kNrpO#-dQM?J8ccxv1dElWRvHW1Lchv#P6^q={&@s}AQ};=+)_m#kSrZGpyIGBp%fuW3O0NqX}l~+ zK`~zuh*40gNOw=QQJSl$c$S3)R1^m++bR;O{VOZZwaPf7$t{ZjJXM_8t8Ex5EW(!1 zKp0yV089n|3*i=nmZcf&1OOV@Hw7%U57rC?PLhaH12F$;RMb=WYDY#`_bYP0Cy0(V)7$^1Ee7O zN@sCrR_IF;pfDgxQ=Uv)k+cA$%C(>tOcBgXFH2e;EvqQUZkdo?y?hX(KS}1=M_^;SAefd$qDwqCL zu;N7JRklhr=S z0#%~y|7X8@P1nEf`Es~r^Jw3?@)Ivw{Mz{%SUPj@-1b;k(B`mMHEmw<=omx&sWb=d z$}ez#ynmy^cCaK@9g0x=&M7FV;91#$lhU$w(!NU6FUz?HllGgxBr1X_06G@YafPv< zy_zb84JKS?k@|a$94E7GgaWk;0;VBmO|BjnqU8KBZ)#^n_jwj1$}A>!JL-zKvaFA z_BkdvS_Bv|3i1TBLJY{0P>S%(1ylBNzslwO*fpw-WHoHvirZSGiJ#CIQFT}d2p;?^ zX(O5t5ah^8B?W~dbf;Je8Uz%MO!SduEPp$kw$(VdS0j545>_msHG^Z}H%aPoY6%zg zQV`$n(~O3-2!0bDtkQ8+GF5NGelD)*qp7$SIA^w`q*u4%R|S-(1Z16$l?D3ax58Hr|+UmqUB7G;+itSvS#oqG7`e77_d5Q z2$fQLE2@@t7UVPu@QMI_hfE17a_r`%`QvRY&5K6mB_`}lsa&yVJN5)dA#C^XqdT=s0L6vpo9 z^n^2P=>2gVn`*%m%4bI5ze%!xrEQoacbH#fX$^=TS+kYWrH&r{j2ypJi+?6D{1xit zh-N#Q{%jIOzq7lG#FDaNn~T6U+YK88eP}m|3l3J)-ao69X9AInIo#9 zZ48x9o>j}1zH-3msabapS$Q359|Ib%WK~xn#`?x^;$<5V`E4OTmxo*elsS8=EW7dj zR&YZG@0aAjd=yr7`lt%r5Y^pVAWnsif zFKR4teUX6@?0%Wc7gR`}5ly;Hb|wg2f7CFD1V2&u-U+-wybi>1F7WmK)lUwuDWoEZ zSSZVHk6kyNJ!8h#!T1QqGm(o1PknGvLnTz~g;Aa&%J8V_jptRY@`puE5y;5LKDcM0 z^yW{!fwlT3@qa7pQMK?%3?Qb*P(=0gw<}cK5A*Eg6o4c+Up-wGi7>ox^_(1?4lEff zZ}1Z#x2giVe6X;J9J-VYPHXaoxy5M4)KeQG7rF~`q{`C;J}Yz~bek1X^({RxMT>ue zV*A|LNK?ck?_q6zhcHQ!?dEH3WMbVztXhcmtj)q5b?(O-#VWr=%{m$g$K!#Tz*irN zr*SO96%&ois_U(6vE8~=ki4q&-1_w&Hf&{QWkDQEly-QrFu7}$2eFsve8P&MNns|S zI2%LvB>gr*;4nA5K|@`UR`P9iBCqk&#<(O$%Sb@g(7l|1(NK-!18)|7p?9$H2gRpa z^qkv9%jkQ$uX@AKt)aIx2aL>EcUqhVs@c_Vf8`7w1X`!mE*&q}on%Mi3)Rs?Z{PUL z$h+7ihzHBBfPv=bTk)8K{6ux@WDpKXKr|!nu-0duBq+HwD!iB!;R`+s&sHs{?5wQ*(Bx*Sq(07FS2ugPjx3t^%O^K6==CViwxehfPjs`WE00fFVRr(< zL3Yx3J}L>H^pHHC`BxNY5*r?FQNBr%1rj5J+pYej*hRae1dArjxGX&~Q_3#&`_>OH19TyDR(4 z`WF2S>?Uv6C6Jm-Cp5G}H*(Ye5r}$|SvQRzh9m9@Z-Z0h3NHXDGGQW5s8DkY`^N4^ z*w;{|cBo>7#xwVxx;@LAm~Z$al%GClt~!gx8s~ntqlA|ql+D`P!)e++B0I`SJ+Uqn zQZzqn0i^kY%iC2)VS$A=BRJ18p_68>|4k6$OA0^FCl&b>)nuvW5Yce6hBRW^Jc?{^ z8oN8lU$|&+7XO_Bf)1If2g84gymu>Rhm|y$P@ymhbYpbW%>MDS}q+P(#5uq z6c6uNv%2zGYBC`e2h%P`NpQQ;*0Mp!RtNqg&RELB-Yw~6v&b4MB&MNaZr%>XBoFUG zhL~z#9>g~k&n?LdsKY{Eab4KcT=fA2!Se)bLYSX>3#x~eq{IffJLopC2`ZN8=AJ@B zgzqC;Mc3=AtWudHPfSmUUW+vcQ-qZCq@20i<$H-Iw>sZh3(r(k>osTl&S-azi;lUhat}SGD`x38B29;(xv_ z&>CNG)(2{U)2)A3ZVJ#PvH%O&D{bxG7iiyOHutrIG{c{CnMho?HNjPZD6EEwL<7QgK ziLf8gz7rK90UN9jyl=3kxs<}Z)Wv*MD1mUnNy8=r#r(;<&ZvZfFj>w#%MTd>$sf!7 zdUj0}8_MmT|?W!X_O<5E8HRwZ>FWL^Mj zN65i${%@OgG9Oqsn+q(OA&HDO$lw)r4}=!;E_bHai+1rYM`hMEpTi6(v-wB#T_?3L zo1{NJTzcdkUWZTf*GoP@vf~%7F3<}``y7Yl&GK!*LyhL$vZE*EcJU|L;HGXRB>pEQ zqcDdL?46z`@Lg#`o+vc)i@T6%rW!Z$BgpPT`9yU!G0`q*Z9{qehW_nz@X!2m6XDU4 za+c2}T5sqhm;I+*r@~57j~C1p`aLMcSJeT_JL1LHVlw-r)a=K~#&uJAKZ;sjhg`PuFh3{1PnsUROV6MT98WiOH1V$95 zj4ZP7!k;_czN8tOzNYN1F1>JYW9H(a_~A|+^OX3&b}3hV<54LD|7WaA*wI~_gtEI` zD+~DugHU@;b5!PkV`G+=Ty594$G~u8CFAHU-Kl2F#MsZ_zE;d>j^rc0U|HK``?1^< zYr{0XK$%&o?ahvZRGjl!uB=o^hx=;i7(}0u`^dG0Fb;|oR)`WC>2u~dQ}K>=IGvY^ z?QQ+FtQR-h5N`{!M(sP#;t*S%{>9d2&(+})K~>J%GxnasnN9rD%vmMl zS)b$67KtNYM-W$r>CIrspmP)82IH9);l(gi$01P{&>oEb$s&XXu%O3BqG!aB_yZ^V znPeXU=8HN>aVCQ2Wa{wp$pOe!C`AaGBO>5Ei8o7VYe!i=D(2{TA$F0*HIRzW@|3Wj z^)?}YbePKO^9qDOFxBcEHnix|r!a=G7|hcRgUD?>HJ&vlwwU~rY|ps`=BNzbMPeLj zV1|~}jri!;XSNbrNp5#W@oNN>A3yZ+<-g87)e_aUD^haC5iT88(Nv%oRW~+TtdAEgX;)T%QW|rvlmj1r3I)@VB~Z;*h-e z(9KyTF0u*PFdbgf2uly+Nil>9usYl??IKk9)RZ_~Pcg|R>=d)mlfn`2F|!RYpomW~ zb$Zz6z?qicuA@G)v4t7yDh><{p{ZdQT*Nz3r|wR}Z@v;vJnS>75occVNJ|AmW=fr= zavTd}@wMg22Lr16Otk}Y=ZH+)%r(sUNqU!&BZ;!WWqtBQuJ-el(5)qpPmX;6xXmy9 z=cbm>8C~Q=QS)zr3Wk{>9-Atk+~Li2=!T{FNfY`GFBpcopnuCLjyL57zd#LQM{3&Q z6*-13@s0`iCsoj5n@G?UrgbEd!y<$533O6+;DP-V-?wn>w}*$j6M}@qU1j7ve4@V{$k0|d%lj^FGzjiIEizB+nX82vdg3LB4ht#> zQW2&qFwKhRF*sLZw%4=McvfVRK9=(UPy{P1mbzj)T@ue-vzUMmh`~|F_=uWKChst?8E+r*2Pab&iqPUGKp*) zD9f91`b`4To}NJZK95jA7#K*b{Z{gDtK8n6f*V&9ox#i2X+1WcRv*9t^L+j-k;}lo z(%GAak00VHDoB3p1urs?<RRAyI1$-4Ggb z0YkSR2NpLQ-~y#KEE!v#qD`$vbiDqhw-j2cjHv#ft5nK^?qcI8gTFN(HSm*hti`sH9~AXcIE z?M-m3!VeujG{^2sg7Pwi++p22c} zW#ZATCl;%Bd4JcxcCVk}jfjQtWZotxTn&z0e?;pGVjD~84zIj^L=BpbT?;;cGOw|y zDEjDv%O2XWZS=TwEq3%!srvpZt;V^P1rc=Eg6b{!Dwm3gX{ z6`Pt~JZX7^%9>)LUwAsx_OGk*A!1RXV)qXm{(LE5~N9R2*fio-MPYSTiFp0V%~b z*Kb~1?UyZ0ZVxLwil*BYcDKbE&y3R{KRKn0Or5$M~(B5zM0+(3Ti z_h26|O>NJ;BL@c~r*vO&P*zw|<>;1*0FUnt7p)u}wkjogSg6t8Ogba;-yJ~V!Oj(Z zuT5t*yawZusP9Z#{|kCd|F7Zyo%Da0cH!LHV>o&sk@T{fSg1!UDyhs>neuoYxcs-) zOZ&3gX`~*mWEO`0GOJ{eJnF`sY%o5er;wJppGVCIZ&x%6{J zLG^ysbFS3rWP@abrWlyZ4FHc7{4ptVK5%SPyh_*kjPPr`gMSDW{iv;9UFbQ9?v+CE zr;w2mWRn;zm%)%?!$#+y)wPTn)+r5E-lwPWpp} z2$+ht6KOZ*ym!k^kUsNm{55V+BK>yQkMoUdfduwdfJSX>v>P*5dHs`BLjSq_+;w$T7GF z;_@zdT=+i3(=seTE%!@tpccueix4%b8Kza?Tpz?4hnqL1&DqS^Cn`o_CqJ#7`}WUTPJpMVR$aBnorF zgmT;>hu3M^_1rZIS`?ZPp}Li$_H!G4^+q!xuVsmiQW(XE6-?Hh^NqUpU_yhxvGbPyIG0M%m->q}WfH9Kd!H`k#X~b9g zL2&sFx7`}NYr;pl2X zh~XKEetoR)BIG9O9G6o(@FEnvuVeub(VG`Ik!1}IUt#m*X>dJZS4{NRQ4=une zevrsb(~?Xj;K%Q{zjY7TAxzVE#Ky--M-9egFpXRtkFGfn#@u_ApQ4VtLcsnCLv>Y(05Js~7xqi{?_N(I7oXBUsoH##;`B%LzV$Qf zR__&5jAlJq+r@qRyke#>G&xIj`q$o6US93`WZweu7dDi;fa`WO#j4*q{es~lZ~5p# z`Ni@UeBroHb?rqrP~IrfU8HwEgG1m5NUm3rsc7*a>tZ!1V&|3;1KM=VSztzI4{1OhA7>t6uRLY|ars=XAOf+t)n%bX z8YrC1zQ+;Im2w1KHrXXVfD(C&O;d`>=2s5dw$GYy%crpShhm!?)X%j~zjy{GVOS~M zVtSVpRJ89Z^j0qZQ=t)_m$;>39{AaOYsGL+mmb6;882RctnW85b@DQUG5h1{+vCHQ zco*UdbI|Z36NOh;Dc1PEl6LveJ~{M;t>_yNDuDS{Cyobe6hkDdiY2RmKef!a^3ajv z->c0L1Ix`XZ*{TSDnWhGJu^Q8WX&C(+P9cT`6EPIPmKS0!`uDT+L+V)1LSv$Pk4KX z&mHrt50bxN`%f#HQW7jbFILqh>y55jjCyaiEM}+oURs6d;BK1ZKat?1B_cSmJ(8Ew zJhYZy0}e)085+2>1sg5(H_q3T6D058pyc2pHE8J#wd9oF{7?k)QrC)cCj^f>3XG<3 zp>F|1%C4isCOo(w(KPXnSY7RN)466xaH@PoDp6dH-GUi+`toO~V#*Sb9!Gy+wk2kP z5m&pTW%QkUozs0yG4-Eaez|c4R#V=ZeUml6YP6&e{V==rbeCRz4BpRYdgRoosW=o} zy?!;%v3cmy1^FfvRyN2Pw67Y;r;zG@-NOZr^WI*cTIO- zDs)H}P(WBt_Eu5ZcoXJMZAgh|v3L1i{2uM8R9KJ7!O~;-%e@Oejw}#bHH@lsr!3(4 z_>q+M7`qgPi3^#9r$(3xGlf%UQchT71FSc{En-LubEuow)0=kP^M}WFt~{iXM~Hlk z5~{3t_GeR*C8c|F@8X{|xyK69#jsaeTC!(q$B8Robzgd<2;>bCbTQPBJe>E|?8zMf zS7A}3bqWvHzkB(wzV((xtfInEQzf4%a~G~T9iqf0M{qM*+&PfFVC_Z#DHkqakz$bA zmXAcrtN=`WvQ(g-OaXJdygR{REl#sEpLa2FPRS`_B^3h?olrde1%6A;M6{(M9qWA3 zv~k`Ob^8cwJV5K;t+^ak&Ue-r15viy#$<_*i`NJfv-(T2!YBRFYwO-$LCvp-J&ibK*pO1dbrGn&^BCp_Q~MKZ!{g6m?pb zbepEVd7!rujz6_Gep$v~hf??l@x`fAjR40Fiim5~=8I~P!@J!YU$BOq@0})|*uc~T zQj>)-GGZ=b5)*T;YOuM1U$M)e$Z2l*oq54!#?R8ZG58uAtmw@SdgBq zQ;-G+eT0#nN%f7wY$0gw*lHgAo4;=6bM3+B?}dI~9Nu36 zBne7oqunE~YF`^fl0R0D9{ozYlhOnIDM1OwF`kc@2YsMB_097Snh1yOQP=5}n9;sZ z@7C^!9JlKYE2_^>b)QF*-GUlxzodOrKLI%gIZ>sqEozwA~G-FZQ*B(yHq&+;5frIp9zU^x)tKmjqy->de zRsTU6G&vECr;_Ob0z-leMib%e*zels*QptJSrROG%H|i)R z9519vutkKLK;3;v1><+~awThBUa-E5xtANb9_MPKqhmdmRCyp1uCivx}4G=;BQf`y@j^JG2_=2LYbl z{lWq>SPQAqy$4S|X3zgR!hQ=;_XzRNR*m>tN_FFT7@~{6gPYhvaG-iGiwzE|EHxZi0WVfN%(34e7VS1O86C2`f?*+YQ%A68+oi1H3 zNzrMKg#^iwqW#KXE_J- zz7ka7`&GLWiVB0YA)!tOO|b>c~!Nh1et$a)Llwz!0vB}S`+sGEhWSFIY~Z!Pm@ z9D`S0U`y~VY~w!L*TXVCI^;axq#SFq_c z{f!m+g%L6ovNPqA%k8YE*&n}k7z$Ggbbo9*q7f?eYWo;JzmqaD5~Kb99vEQ1l?;Ub zKxDU~_)sP#EXMDMF@L)^$h6R>B!l50LLxfWb3W3>QlsmPG4l8O=Xl7;>MVv|~nG^BC z*6m}jD(PM~r>)cH)@f#jORdS!%?C|c`+w3bEy?|&2)AYR@x>KT{*pySX(^1wbm^g% zNS72tk?0&^TM2MvcDs0=X54ErCNSe0R@ii z1W%yQIW*%4mC6E{**(7<9x>J%BGrljIKW6@A@|VYSy8ZxmYcn;IE7pa06vVTqrTV zD09r#><-oil0~Ah^SW)tjOM_${_J+TKY&lf4o|P%wXOcuF!#{wrhVR#YUXeK_Xs1> zKCL+f4BL*??hfJpcY~~`0A;DI_S4H#o0Q(8|s;>zJWMB zi&hF(Y%%M5_+LMAy2JK>wpml#F=_%Zs>hQdjo%8}!{MsiZ;Jxd;ICj>jBX9-Op<=W z^wahZGHlaiK+|D|v#OM-pD%%m3OG=Ps8{H{6w+`*f@uatMk1pa`ZrBcI81jng7TeaAckAO1v`bs!g?Ow zwZ5sI9dXnQYR^3$Nrz-)Fd~~!%J@+Rr7Pr{_xO2RW;O~UuX+V`*h|elKAqjXZNeWZ45Uobtb7Feicg{ zMu^efaAxFrl1CLj>CA4wGDeu_?TP3D25UVAU`d2sL~D0jlyhU6tS5@K8mvvJQSSpj zdZMgyYe;?%FG=Pdj2V<;c_vc9Go%4nAi?gf+9yA4@P1yJYLy!caw5aR7+%6HS(Q$t z+jqN@u+KTD5QJEAD{&F?o?yJ{(iBs2Ae5)J9z@4`aiV)JKdrr@-xFy$b{Wr|(r$rk zR`{jceuK4-4Hx=9eKr-QOs+YoJvc~(EK+9xa8kU^P?2*%w3II!VV604v)e7Tk7(I) gtRct@s(!@D&_zbn;Y7f%s1Ay3IvKg#`fvXK17f_o$^ZZW From 094835a5ccc35cd5dacee4cecfbb2372cbc0ebba Mon Sep 17 00:00:00 2001 From: Zachary Susswein Date: Fri, 29 Nov 2024 12:12:25 -0500 Subject: [PATCH 2/6] Update vignette to use new col names --- vignettes/RtGam.Rmd | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/vignettes/RtGam.Rmd b/vignettes/RtGam.Rmd index b40198e..dc77be1 100644 --- a/vignettes/RtGam.Rmd +++ b/vignettes/RtGam.Rmd @@ -38,14 +38,14 @@ library(ggplot2) The `RtGam` package requires a time series of observed incident case counts. For demonstration purposes, the package includes a simulated epidemic dataset, `stochastic_sir_rt`, adapted from @gostic2020practical. -In this example dataset, `obs_incidence` is the vector of observed case counts, and `date` is the corresponding vector of dates. +In this example dataset, `obs_cases` is the vector of observed case counts, and `reference_date` is the corresponding vector of dates. The dataset has a timeseries of 299 days that we'll subset to 60 days to make plotting easier. ```{r plot data} data <- stochastic_sir_rt[41:100, ] ggplot(data) + - geom_point(aes(date, obs_incidence)) + + geom_point(aes(reference_date, obs_cases)) + theme_bw() + labs(x = "Reference date", y = "Observed cases") ``` @@ -53,7 +53,7 @@ ggplot(data) + Because this dataset was simulated, we know the true $R_t$ that generated the epidemic timeseries. ```{r plot true rt} ggplot(data) + - geom_point(aes(date, true_rt)) + + geom_point(aes(reference_date, true_rt)) + theme_bw() + labs(x = "Reference date", y = "True Rt") + scale_y_continuous( @@ -71,8 +71,8 @@ The model can automatically select hyperparameters like the smoothing basis dime ```{r fit model} fit <- RtGam( - cases = data[["obs_incidence"]], - reference_date = data[["date"]] + cases = data[["obs_cases"]], + reference_date = data[["reference_date"]] ) print(fit) ``` @@ -84,8 +84,8 @@ See `smooth_dim_heuristic()` and `penalty_dim_heuristic()` for more on how these ```{r fit model again} fit_non_adaptive <- RtGam( - cases = data[["obs_incidence"]], - reference_date = data[["date"]], + cases = data[["obs_cases"]], + reference_date = data[["reference_date"]], m = 1 ) print(fit_non_adaptive) @@ -147,7 +147,7 @@ plot( mean_delay = 0, gi_pmf = sir_gt_pmf ) + - geom_point(aes(date, true_rt), data = data) + geom_point(aes(reference_date, true_rt), data = data) ``` Future development effort will aim to improve model performance on challenging datasets like this example. From 7b00d8084c4d7d9cd62a959669a1fa1cc40992e4 Mon Sep 17 00:00:00 2001 From: Zachary Susswein Date: Fri, 29 Nov 2024 12:13:42 -0500 Subject: [PATCH 3/6] Bump NEWS --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index c15aa43..235fc84 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# RtGam v0.3.0 + +## Features +* Add day of week effect to simulated timeseries (#91) + # RtGam v0.2.0 This release introduces new functionality for working with RtGam models, including S3 methods for prediction and plotting. It also includes bug fixes, improved handling of input data, and enhanced developer tooling to streamline contributions and maintenance. From 21375d78f0bd9feb169bfa515cc9b19bfc5afa55 Mon Sep 17 00:00:00 2001 From: Zachary Susswein Date: Fri, 29 Nov 2024 12:35:05 -0500 Subject: [PATCH 4/6] Clean up examples and docs --- R/data.R | 54 +++++++++++++++++++++------------------- R/plot.R | 4 +-- README.md | 2 +- man/plot.RtGam.Rd | 4 +-- man/sir_gt_pmf.Rd | 5 ++++ man/stochastic_sir_rt.Rd | 50 ++++++++++++++++++------------------- 6 files changed, 63 insertions(+), 56 deletions(-) diff --git a/R/data.R b/R/data.R index 24e70ad..c7150cb 100644 --- a/R/data.R +++ b/R/data.R @@ -1,9 +1,12 @@ +# nolint start: line_length_linter #' Synthetic dataset of stochastic SIR system with known Rt #' -#' A dataset from Gostic, Katelyn M., et al. "Practical considerations for -#' measuring the effective reproductive number, Rt." PLoS Computational Biology -#' 16.12 (2020): e1008409. The data are simulated from a stochastic SEIR -#' compartmental model. +#' A simulated dataset derived from Gostic, Katelyn M., et al. "Practical +#' considerations for measuring the effective reproductive number, Rt." +#' PLoS Computational Biology 16.12 (2020): e1008409. The data are simulated +#' from a stochastic SEIR compartmental model. The original timeseries and Rt +#' are available in the `incidence` and `true_rt` columns and with other +#' added or modified to increase noise and add a day-of-week effect. #' #' This synthetic dataset has a number of desirable properties: #' @@ -25,33 +28,28 @@ #' incidence and/or observed cases is often undesirably low for testing. Many #' empirical datasets are much noisier. As a result, models built with these #' settings in mind can perform poorly on this dataset or fail to converge. We -#' manually add observation noise with `rnbinom(299, mu = -#' stochastic_sir_rt[["obs_cases"]], size = 10)` and the random seed 123456 and -#' store it in the `obs_incidence` column. +#' manually add observation noise and a day-of-week reporting effect. #' #' @name stochastic_sir_rt -#' @format `stochastic_sir_rt` A data frame with 301 rows and 12 columns: +#' @format `stochastic_sir_rt` A data frame with 299 rows and 6 columns: #' \describe{ -#' \item{time}{Timestep of the discrete-time stochastic SEIR simulation} -#' \item{date}{Added from the original Gostic, 2020 dataset. A date -#' corresponding to the assigned `time`. Arbitrarily starts on January 1st, -#' 2023.} -#' \item{S, E, I, R}{The realized state of the stochastic SEIR system} -#' \item{dS, dEI, DIR}{The stochastic transition between compartments} -#' \item{incidence}{The true incidence in the `I` compartment at time t} -#' \item{obs_cases}{The observed number of cases at time t from -#' forward-convolved incidence.} -#' \item{obs_incidence}{Added from the original Gostic, 2020 dataset. The -#' `incidence` column with added negative-binomial observation noise. -#' Created with `set.seed(123456)` and the call -#' `rnbinom(299, mu = [["incidence"]], size = 10)` Useful for -#' testing.} -#' \item{true_r0}{The initial R0 of the system (i.e., 2)} -#' \item{true_rt}{The known, true Rt of the epidemic system} +#' \item{reference_date}{The date the cases were observed.} +#' \item{true_rt}{The known, true Rt of the epidemic system.} +#' \item{dow}{The magnitude of the day-of-week effect added to the log of `incidence`.} +#' \item{incidence}{The true number of cases occurring on the `date` in the +#' simulated system.} +#' \item{obs_cases}{The observed number of cases on `date` after the day-of-week +#' reporting effect and observation noise have been applied to `incidence`.} +#' \item{obs_cases_no_dow}{The observed number of cases on `date` after observation +#' noise has been applied to `incidence`. It does not include a day of week effect.} #' } #' @source -#' # nolint -#' # nolint +#' +#' +#' +#' @references Gostic, Katelyn M., et al. "Practical considerations for measuring the +#' effective reproductive number, R t." PLoS computational biology 16.12 +#' (2020): e1008409. "stochastic_sir_rt" #' Generation interval corresponding to the sample `stochastic_sir_rt` dataset @@ -73,4 +71,8 @@ #' @name sir_gt_pmf #' @format `sir_gt_pmf` A numeric vector of length 26 that sums to one within #' numerical tolerance +#' @references Gostic, Katelyn M., et al. "Practical considerations for measuring the +#' effective reproductive number, R t." PLoS computational biology 16.12 +#' (2020): e1008409. +# nolint end: line_length_linter "sir_gt_pmf" diff --git a/R/plot.R b/R/plot.R index e8b6fea..b807585 100644 --- a/R/plot.R +++ b/R/plot.R @@ -13,8 +13,8 @@ #' @returns A ggplot2 object #' @examples #' fit <- RtGam( -#' stochastic_sir_rt[["obs_incidence"]], -#' stochastic_sir_rt[["date"]] +#' stochastic_sir_rt[["obs_cases"]], +#' stochastic_sir_rt[["reference_date"]] #' ) #' # Plot draws from the fitted model against the data #' plot(fit) diff --git a/README.md b/README.md index a4c6b90..ca8ab8f 100644 --- a/README.md +++ b/README.md @@ -39,7 +39,7 @@ This snippet uses the simulated dataset to demonstrate how to use the main model library(RtGam) fit <- RtGam( - cases = stochastic_sir_rt[["obs_incidence"]], + cases = stochastic_sir_rt[["obs_cases"]], # Randomly chosen date reference_date = stochastic_sir_rt[["date"]] ) diff --git a/man/plot.RtGam.Rd b/man/plot.RtGam.Rd index f6972ec..aa1140e 100644 --- a/man/plot.RtGam.Rd +++ b/man/plot.RtGam.Rd @@ -29,8 +29,8 @@ draws. } \examples{ fit <- RtGam( - stochastic_sir_rt[["obs_incidence"]], - stochastic_sir_rt[["date"]] + stochastic_sir_rt[["obs_cases"]], + stochastic_sir_rt[["reference_date"]] ) # Plot draws from the fitted model against the data plot(fit) diff --git a/man/sir_gt_pmf.Rd b/man/sir_gt_pmf.Rd index a58324d..67b9146 100644 --- a/man/sir_gt_pmf.Rd +++ b/man/sir_gt_pmf.Rd @@ -27,4 +27,9 @@ distribution. We produce the PMF using https://doi.org/10.1101/2024.01.12.24301247 for more information on double-censoring biases and corrections. } +\references{ +Gostic, Katelyn M., et al. "Practical considerations for measuring the +effective reproductive number, R t." PLoS computational biology 16.12 +(2020): e1008409. +} \keyword{datasets} diff --git a/man/stochastic_sir_rt.Rd b/man/stochastic_sir_rt.Rd index 0f49cfb..08724d8 100644 --- a/man/stochastic_sir_rt.Rd +++ b/man/stochastic_sir_rt.Rd @@ -5,38 +5,34 @@ \alias{stochastic_sir_rt} \title{Synthetic dataset of stochastic SIR system with known Rt} \format{ -\code{stochastic_sir_rt} A data frame with 301 rows and 12 columns: +\code{stochastic_sir_rt} A data frame with 299 rows and 6 columns: \describe{ -\item{time}{Timestep of the discrete-time stochastic SEIR simulation} -\item{date}{Added from the original Gostic, 2020 dataset. A date -corresponding to the assigned \code{time}. Arbitrarily starts on January 1st, -2023.} -\item{S, E, I, R}{The realized state of the stochastic SEIR system} -\item{dS, dEI, DIR}{The stochastic transition between compartments} -\item{incidence}{The true incidence in the \code{I} compartment at time t} -\item{obs_cases}{The observed number of cases at time t from -forward-convolved incidence.} -\item{obs_incidence}{Added from the original Gostic, 2020 dataset. The -\code{incidence} column with added negative-binomial observation noise. -Created with \code{set.seed(123456)} and the call -\verb{rnbinom(299, mu = [["incidence"]], size = 10)} Useful for -testing.} -\item{true_r0}{The initial R0 of the system (i.e., 2)} -\item{true_rt}{The known, true Rt of the epidemic system} +\item{reference_date}{The date the cases were observed.} +\item{true_rt}{The known, true Rt of the epidemic system.} +\item{dow}{The magnitude of the day-of-week effect added to the log of \code{incidence}.} +\item{incidence}{The true number of cases occurring on the \code{date} in the +simulated system.} +\item{obs_cases}{The observed number of cases on \code{date} after the day-of-week +reporting effect and observation noise have been applied to \code{incidence}.} +\item{obs_cases_no_dow}{The observed number of cases on \code{date} after observation +noise has been applied to \code{incidence}. It does not include a day of week effect.} } } \source{ -\url{https://github.com/cobeylab/Rt_estimation/tree/d9d8977ba8492ac1a3b8287d2f470b313bfb9f1d} # nolint -\url{https://github.com/CDCgov/cfa-epinow2-pipeline/pull/17} # nolint +\url{https://github.com/cobeylab/Rt_estimation/tree/d9d8977ba8492ac1a3b8287d2f470b313bfb9f1d} +\url{https://github.com/CDCgov/cfa-epinow2-pipeline/pull/91} +\url{https://github.com/CDCgov/cfa-epinow2-pipeline/pull/17} } \usage{ stochastic_sir_rt } \description{ -A dataset from Gostic, Katelyn M., et al. "Practical considerations for -measuring the effective reproductive number, Rt." PLoS Computational Biology -16.12 (2020): e1008409. The data are simulated from a stochastic SEIR -compartmental model. +A simulated dataset derived from Gostic, Katelyn M., et al. "Practical +considerations for measuring the effective reproductive number, Rt." +PLoS Computational Biology 16.12 (2020): e1008409. The data are simulated +from a stochastic SEIR compartmental model. The original timeseries and Rt +are available in the \code{incidence} and \code{true_rt} columns and with other +added or modified to increase noise and add a day-of-week effect. } \details{ This synthetic dataset has a number of desirable properties: @@ -58,7 +54,11 @@ In practice, we've found that the amount of observation noise in the incidence and/or observed cases is often undesirably low for testing. Many empirical datasets are much noisier. As a result, models built with these settings in mind can perform poorly on this dataset or fail to converge. We -manually add observation noise with \code{rnbinom(299, mu = stochastic_sir_rt[["obs_cases"]], size = 10)} and the random seed 123456 and -store it in the \code{obs_incidence} column. +manually add observation noise and a day-of-week reporting effect. +} +\references{ +Gostic, Katelyn M., et al. "Practical considerations for measuring the +effective reproductive number, R t." PLoS computational biology 16.12 +(2020): e1008409. } \keyword{datasets} From 064b4e5b1f1b3ce9bd88c17a1a9d8f1b94f52466 Mon Sep 17 00:00:00 2001 From: Zachary Susswein Date: Fri, 29 Nov 2024 13:14:45 -0500 Subject: [PATCH 5/6] Fix silent bug in day of week Was silently producing and broadcasting day of week vector of length 13 out to 299. This bug was leading to a phasing of the day-of-week effect and making it impossible to fit a meaningful effect. --- data-raw/stochastic_sir_rt.R | 5 +++-- data/stochastic_sir_rt.rda | Bin 5596 -> 5589 bytes 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/data-raw/stochastic_sir_rt.R b/data-raw/stochastic_sir_rt.R index 1d67935..8ef6ea4 100644 --- a/data-raw/stochastic_sir_rt.R +++ b/data-raw/stochastic_sir_rt.R @@ -7,9 +7,10 @@ gostic2020["reference_date"] <- gostic2020[["time"]] + as.Date("2023-01-01") - 1 # Add an outcome col with a day of week effect # Use a sine curve to make effect large and obvious -dow_effect <- 0.25 * cos(seq(from = 0, to = pi, length.out = 7)) +dow_effect <- 0.5 * cos(seq(from = 0, to = pi, length.out = 7)) +nrep <- ceiling(nrow(gostic2020) / 7) # Get dims to match -dow_col <- rep(dow_effect, ceiling(nrow(gostic2020) / 7))[seq_along(gostic2020)] +dow_col <- rep(dow_effect, nrep)[seq_len(nrow(gostic2020))] gostic2020["dow"] <- dow_col # Apply dow to log expected cases diff --git a/data/stochastic_sir_rt.rda b/data/stochastic_sir_rt.rda index 443792422ba05e1477ba5d18ea047922c075ebd3..8b19ab1aede4d897e6de3e77bd508838adf53dfb 100644 GIT binary patch literal 5589 zcmbVQS5OlGlLib`N(7NA66pv?Aasz9R6|EVM0yELia?MOB{ZqhJCYE3FH!{*2n0io zbVC>Ey^4DNFE?{@4>xn)%g*lX!+ta0?9R7pR<3ejRZD&oOT}FadlGv4@W1~TmKeh5 zug>3-n{kpbkSOfix(Hs#z$@b5{pvMa^ESl9MKxeKV9Xv;S`XEV{pkh<7Z;YN5^zBV zJa{Ia+eb;fLGc(o85E+WajM84FSO;1A^bMM1xaX?VzQw^D+Po}eF;XmTN6bb94aLc zA+!}{LRK)A%_f1=p$StV4Cu$^(D9!dC>g{Fl%iDA?|pze37~}POa|eA8AP$$ra`?4 z+`Ryjm4H5RV}tZ7CaGjiYC;av#?!b`p+cow1t@~Ic2sovx^xcIXq`g{HEEzGh%%8* zoyZX=ZKi%}V(%XI)_@c}NnIXq0A(z!fWzH$LU_5w-GixqYYZgMvC}~&wFnd@1g;MV zB@l*5KVi#2Iwe8f;HntRp#3$6N1ebrsjZQ&tT8!WbmP9g9jmK5yOfdvdjV7xYJffg z43uxd`3XG*yhCfOsXP-(0O@S4ailbcYW!!;KUxsC)CspO1{0uLLo_B*3DCoDz>1?P z?J}Zqq@x2AxJ%^n=Rn=~bKq>CtsoP1f;5H#LpVVcu`4k_b&5_Q>~w69F?P5LiINIC z#L@NWt;*mj1af4+VT6LJLEw(L8e`N_f+-bns4zU&ksq<81<4&t={6>ra4&2f)d#Y%0#4c4UjXRXU%0bM2y9WOU<}yh=dB`0{$C*bj=j~_ zF~29H{Bghik8&>l`|jEC6>)L@<6k1nTZS)EGP2Qk@BD1Psw7nw3wRea|Ss zg-$AXlsQ54=(eeUB4>L^VEZ0%dOQ-8@?k83(L<(~5Gh$j#nzM@C7Jh3NDqaO`?+*W zw?kUqn^}^lffku4o1zZ~;cmB*6%SQkt+1JZtZ%7^KC^U*%0GU4qb_oJ=pP+2>J}ur zbLp=MPjcjyv)i0kf|~Mh$!LYoZnfLc<=YC)=@;2~L(1}bcEa7c-uK=WeRdIC*|^>o zq8|g2khDA&{_KqkrQrr~=`pTsy# zWqwws6ByzhXbqzAJr{xvKUUL*T{$%uX62a>TMT1qJXRX`JK3me@eW5X>mZSa&(Vzy2qwPSD|2gNcE^sCrCVO z_GpDxbx)eMrlYPdHwgRVQKy58{Nc%Cx(=p5Z?_s}1A%FL^bN@cFDdc<<;E7v#m8dW-X=1ZXjO)s^_moZx&6-fgO*A>HXMFp?i@% zp_teEKDN`VjDPH2@YPUUJ)m)B#l73A4@&7m2RDQrdhpGU#v`ul#5UITAL{JAw5^{z z=2fhlt$^y*$u#JDaylo){(N?I`_ww}yBq2+XqeR(nfCIo z6sx3jpLMvVV1KSHvD?qkXIZABkoJQ_sWK8i#r9sSDs@xhfoyU{PG2w;l<&P{hG?P!o<<=P~WpEk$lDWZDS1I z2(-K_>~_ad+)4!9&DNELsK$D_c>G8e(K>UjdQi6seV@|3X#h~7Tv~P>9K&&4t2Wl1 zkO{9Z4eY7A9V!CnibHapYp6i?^2Kzhzy{Ac)|l&Fx`VxO)ioJ!8dy3f0u-5qr5^1q zMMyA5Olq_~U|&-RI`8Fhc{`*rPesEE83*k!S!pxs4D&{&J}p0}vy%J6UmLr|OA@V2 zo)V!VQAH=Avh*l$r7jvc1QNg=RGy$aGlAqc%8x4ow5^tdJ%IwUJyNk$&&mMGxYudO z*dBCdF57I8{#7zEi!@X3Md&M?wg$DkkJm@Pw@@JY1eSHR+0JecDZXN^EGRh65PS9Y zMwG#0%egyaD9`GAYoZ!&no2L>YRddXyu@l#L@93vMA$8irS1sX|oL8a8MyXf5;+0dI z*=)nG+HfWJ;0f%Z>NbbTg9+0N-jzy)tWcYwlQXNFDW!Sp+^~DEQLFF z3MF5AvA|q_Xh?xgM);Z#Er3oco3+2)R1`XG3a$Gof)d(3!zWJSJ{1|xxaSSuHa4c8 zc#zwyXL0?*bWbl@!)!{hMhExAY(M`O_W+o&E9~r9V^A+2nAYL_(S^CInJDoVb3iMI z?YZRv$7htLnT_reo35z-v42bL#orq59@u6xL$l%2YX{bx4bI0*PX_Lzh2+Z951PCm zwW~CWWf`7{6nh*qSe@7lE6|xYCXXpe*Und&KL-OXz1m_1iz*$`?_;N23twGrp2;Zm zEt)jX)@Ajm!Qgs-gqp%8zv`~M>U6JenqB$+9&&)$apl6F2(ZkuR>)Cb57!^DJo-|D z!&198WgpA72@L1=f_V^q-e<}@djWotep!x)UXR9F=k(LQa>Zo66|EshsWU_W2S}TC zESC3`pOf0D#Af}xGcT)eB3$-5cifyp|6>XF5cD$buXM@DXHOTq#g2kRLB3t_^qTL& zeu|fz$$vinW4W}S|8gB3ds@8OlB8)Pyqp1`VSB{ke8Pe&8CZ(%`#M~y6yS5^_Kn9l z(H=kzooJUlY{g4*wq;xz(?zYVOuCEBY@Kd^p_crMnBBNZHP|x- z9MQO?lV7qmN^2w0@pZNGE01V&k^n|LjXO80#v}0PC?_G-$j} zTQQnVw_#h+aXo4343iw4IXzWfi4kRJ506PNtg8K!Q<9lUx=>To?=>9cuS%jay&Nd_ zGlx8Gw#?1op!U{7hCxdWYf>?uC{|SxZm*1gqj5jrcj(8Z2yDkYhtzUM_@@|m(V*ae z_ueCP69fgvjjH-Nh}Ekn6z?`>*~+?&jJ0H_WrVkUA%w70D?S%;g*DCVu8SG#xvc%$ zH=h?2bhmTe?Lb{xxYEVLGm<)OD0EJ7cdD1do1F7$E>7E+FPZMjTU}A@mHCMgJ#?oM ze*MgN2D#&lLsiipw&B+Z_I5-&1NYyS{%mJBkKsGBTIv@`4=Q}V{pLUEr2X{cZMcj! zc7HLtb6@)x{!vURXM1yQl=172uDZ8-<4582V6%Xu(Yt4^ zvS7jAMnn4=6`onI*<&X8dWM?Yd>p^bXlO4RvMV%wUq~9ccQetur2ovC-`f$-?^}Of#nJym+4A3+WzJwz*LLSUs3&ey2{T^|A;UIP_!Pg6*Z^8uc{404vkAi* z6^`GTke7T8a#q4t-7)L`ECS`~&!TT?cRafhadII433D2HQg8+Vm7h+{&CqIXqxPk= zRmG2yBM;-sexmtoaN$l)J;UqY=IV0+4=?%iFeAj@hBRgpTr>hC@2FI`few#G9Os4X z3m;~J+zG+du2QcOhI|)E_sCBX0?OH&vj*p6&QDyFg2cXel?11&OB*5HT#+-sqTNK$ zz=AMD-omRlZ+L8tMnuO{amD+;IzM^37cZ=Ha`1pLdn`p8M}eo>$Ki3@uctP#`5M8z z^(2eoHSNlumM45?S^-1yGGnu?<^dP82R?;NKHOAfluf3$2Zy^m)#wWay3;(};e_al z>)*Y@Tf5IH)OD`}0Ng-ppB$K^j?kFy`;5=SnC;Ygt8W1B$6t=rFGalO4i@e`M}Yk; z>}aeo$^|dQ7skz-u#{GR#NHqY?6bT~N_b}seux}F60Uv`j@(ThjeiO7@Dt+liS z`Ao5ICdPD}UC&LR;os%Eqjz4t_&8FB*Yoe<0eemoXDq1eq;QeJ_(ry&gTUL4v$2U~ zhwogpet*EpU|lS_C_A1dChTYzeNWbKEtH69utsauZu-XG?5v9U(CXe?B>YY`%_%t^ zOYjv@3F=YVnBo*HeO_rk<>(V%;e^-^_EU@lt~#E3Ep15 ze3D+e@jSo0s2ivPlvEuIQDO9WDiXY0dY|zo;@ZPCdalJKR&&M1jLC$^WsXvJY)B{6 zCIvLFuJQ2`-@_+YJI;2|o`0Y4zDH~p`YE(94r)uV)SqSt#mtSmFH-hkC$;Koo^`XM z%V`gc=`w2e<(pF6HLS$wP{pqH6t=i1upRZH1jYwncH%9uTy1p=^o%R9WlSC6295=TzGg+m=wQtdrpiOeUb@ve4HU*2) zi9#d0gbR;(vHK>3`}2R=)fuA{vL7;rGYe}{?s%H}qf2pp^B%C-*{wAKKZP5rahOpy z^1XV+Ku>=SIBa}?nABCZ$<9NR>J z!`^W3FIdyv@q4%YM;>LN?o~|nKzfORed7b9cKr_jY}bvkN& zAoD{x`#O$;oHjF6A7YJu&Fqg{=b9-QdtlGZz$jwXgH{rTeM6k17ffg@^7W|>*_J-;Ul_Z6e7pAWoD=oICP zlNNzVWI4?BFI<~(xH?7lk$33jT#{FwnfG_11_!n4Y!6Pb&d!q6=O0GK3SmN7ZO6WP zo-m<&gP+%WIi|#nRr6TRUykrOi%{p2?M?yl&WmETd@+_E=#eou`ER^L%ChTMKFS`I zi?0+kyAnhz-Y1foHygH=^S&UcN49(fvij=E>(7A3`4&Dnrvv`rDf;Q^1xg?btq|vt zv&$=tOdwK;5Z_}MG%l5axTjXdjW&esf)mfRD5q5n@s;vf>l$i*3H9MgM+XgpueqD2lRGetCKIn3F>6xQ(e-#g-50;5)aA2~c=4is=|t1C9|TnbNpJH56zr^Q_BL_L7pc zHI^@s|IEx*w;*0(7c;oP#@zP$Zi46Us4ald9d+>`DQP4Ak^)I*`rRZqfp^P`mNIM6 z>P{+G)cTTx<`fiN9uiEWF5vL?S^8TUc7(JM1bqgcq6>)(xw9jO9sO1wI(>ku42*0| z(Z?}w8GV-7ygu}<0t0Iveb7(j!k4zFhqRe2VekdTj`Q)kq3*jdd8LQ}Ldd6yJ<5rN z>RcB!`zdAS?dIkXDv$IbY9^ujdT+l9<++(D?W!&@%~-2WLc9a`6bAp=MVetl)>o)4 zH8Yo4ESgqs_z6v4YCaDGt?w_&5(?G>5)4qEb+!Xr**R!R*F&F*G--}k?XY_&T+>;f zHU`iJ8`IyNa>=i^(-uD5%!=Gm3KzMaXs@i0g=BeCugS+xD<<;940uR;4Dg%HI&TRV zyyC#b34Xm4Y53;dHN5+@^`j%Gyu$|RFl&2*&nbuhlamQdyD_C9Vqr$deNZ)(5e=$H Rj{eKyt0`-k{GT`Oe*x>*zli_< literal 5596 zcmbW5RaX=MqeUq}N`=m0zVc2TP;!EXn%eAzp_IR^5*YFeN)Nev>yr``o2+l zHZA_xe4w?&;<}Ih%FQ$#24^@_S*`q=6i4bw04o9xS38QSY!!4IWZ^_DxdcE2Z5nHyt`Jb}Bp?@bg1IO_tBB{QqZ$AzD{dv=$c=`^8jxI-A&@dk z(nXGWy4a5(+8oOAa7XxJA?!v+*Nzsfm3}RsrwIi*=`6Ym0;Ce-j#JWfw4|Ugm+}m{ z3>oFz@)K8%eIYw(wMP;-ocSt8MjahVfxWAVfMY+VQRo2S%$(}U3!$vSId8w0S>{^7 z2CXWfiP>7x`dek?m5-?pYR7REb3nvG#abT6aRnUVuZ&n#cJkMRtF^>b+(aOrm4OH( z>WSk*8Dio@D&@rhD}5ptEDb;>9A7EM(TtGeM_5u~mtfpz2|Gbqvvci)m?-n~uRtQv zX)197RBEZmWhINW%EjenWsZxq<*r8=+B9VZC@wB^f&vdrA$*mmMF}QX0HCoV z6%nlJA?2{b61Bt%0ckC1Y4(7y-8^;?-L3C|6hYl4-GyIGkYpbJx@nKus}XT{UB4T)~q))(Sbt9 z@rQS_?;9m(Awr_4*yk$l^zc%u?PY@y@|DM@OO4(wAj}C>h*)0hy`3xsRImcUR9dxE z#l2j*3^vxb(=Ezipb&5;U;xyh#|GhNA!Ncbd*xv2ScNJ_Y9tgS)aiSk1dm}m{~0vU z>m=wt&P;o8j}y%VtF5bO*;}5l+PHZRqUOxQv#D7>G44Ek;&A3?@bT}bdk?oMu!~En zs|Z1%7v<}wqYTj^U&22qkxd(=m&HIGh@y250_|gugN}5{*UOKx1+C(5`ib7)Hso3y z0?rGPn67SWZ^MdBtf(C(XCLb?1&+yYXTIGMDSxvuEMRW$_(C<4t3vjuen(L~_-R%a z7>n`kJ30S<+#N~hs9U$bi{}y){QIg(8x_g_T>s_0oE#ANTCATLimzFd4a^~;E5%j~ z7VyPMkPQ$Oui`_CjT0)`16)ynxyA9Obd2sRetC7uo~Z0HEC;&e#P?T^Tdgpj;mShgB@#tU`XJ(nPcsFiP-TflQDfNj_67l)c2uY#g!$P zYzffo&h*7eVqG3TLxT&X~{cl^0rQy4hJ3*;@e7DX0#7?ed2=6KprphDPw8B%jfG_%_ z4f0x)=PrpPAcqSgd)9R%mfPSQ55&)y_070lX>;k-V9%u@|Ii0pF97u~*qOAO zP4yezhv8})oBnPh+rl0cAxzYU-fQP_3TUp*x_WY>#l(@tz3I;gOtG`2n~2cWs7C!A z(CX6LBW#2ciN@7szPxE`fJ66X!RJA!BRiTmfa#^gOQQO8F#%sShmm4KY= zc>jEv2O8j>U+J=P)x@)O8>xYa;R2Tj%TMba;QeKP&_>*Z36|PNIBPw)-kqMxOX@#OIx}!>D&|d0WcQ~4ajn0T3bL{nmxB*CE?`QQTk~v ztzYD)!JY%hL(%Xm!_MQ8Ahm>~xW~JgI%qxD1jlzwkA(|n1J3lHwo-v6aV+uS&tw`r zMFK9Lz4>?0;#-rHcr4Q(VWID30r9g_v!x{Yfw9WnX9(jkBJB0lULJ7)b*_d)enBP} z6EFAh4;9RwGJL0ZS12X(Q=0k=C}$8^JS0WuX!-1lA%7P{0kp z`hJ-+ToVEuiM7+xw%s?Gmju-V%k<=3A$pigQQ{Lj&TPS62OqH;w6{zDi2159-r_;N zwy;rotMv&|GdZwHrxkWYc_Fc1x^fv}iia9-PmR^p7GN2vx{$OIt1>lv_le_2C+r&A zo}o6aw2^k3A0mcpUzWZw?|*jYVNQpHAYdc5xjwGMhqTYKQ)hg}dn@XyZYkIX%pWt_ zXFlOC>U2aR_N+|M!pIcFB4JU@O|Y1?mHS{`{9XBjl>HUzG+butjekGEar@q4dfIn% zO|!0iw)_j`%98WN_OzMpaqVV7Sem~sRzz4&pm)E_LFvmJa$N5u)1V7Gej7BqChBTO zf!!h9uriImVEGWq&VC}v^^v(>>Q>Vd-TLz&0pTniaF`yqu(wK14W33b$?qj00tA2y zc{qpC$Ii)0td!usY~9F!m)~C7(Mp|8%Hw41*;Ou{a0UcWi7o9Utk(A}#49Ra zq19+((YWK73@su=&ZZd;hTXVHU4DEMMJ*5h*LYblZMb99g>A*DxGVu*+1xZJ4%GXi z)7DuH%XVlTCOY;H@DqGWO`N3yd;Mcd3m3&&KQ0Pav=TZMf{NRPNcF6cG@kNt-+(W? z=Jz;FT@`GZxdDg5f0ZJ;A=$Bs>F@i8lqwk}sKSZSgc3~0lN5NB`#is?GV_lN-OMd8 zh(&s-bn3-GIKRAY%^A8|z;J~J?zY6~NI15Z{TvpS@<0dZ`B%*vzm9cxO^!z+6xEJJ zMpL3^RhvY`MO9Zl6vbXs)@R#Q(x`sg>gCN~%VE`2jrKH3i>c|I77a4%ZcBYd&IY3a za60UqKw^`nZrz+n)C_;I3bWWx2Koi;3@r>R0Tz1_Fe{<10fv7y4pRx&0=WoK)-w${Yqlv|o3vi;&02b*6uEym(Xu^ZK9g$cbA|HdyJ7faQj|F(fa zB9rdsQ@<@OnxQ~8A{xw&nQ82+?}`eDk>3_1xdy}Eor((6{P%W5-St{BCi}8BLkeys zYtR_6Fi>Psu)h=1V~V)wt-8#LT%N1HFI}r$_R-8wHon_?99I}LxQZDrb4j>S!G00=nPR7bYhFV> zh^~g+{^eKXs@)k)b8xo}n!?o)WK5y;^j5s*^4UbZ{ zmTfuRi{J+tKJ#T?;D3wAJ_yiiz41D$Tj99ZQxongf{KjMWGP1oL_XU`Cw!YG=Waio z`Bc2A`_{aC?y6MWD%YFVwv!U{agJ1CN~oHL5bx5Wbw=`FKNH&5@7t$ zwP`p<@Gkg3Oe^FiVegTkd@PWbFMVGN2YsC3q{kI72evWD-~iK+txWXxw8l4!((20|c`aed4s_xA?tka8= zvjYG$C_|{-FE~G6=nwgss9j5}A9r0+$HHNTFX^Q6N55XOawmQri@~`p^6GtSEHb7b)l&eJvEA| zmcgXHf8?E0-=71QYnszVliG7gW}qswUljP8vY}Qw6dtiAt2JGo1MtCRLoMz7Dmb_D z)ZGn$%(L1eeBLeI#4Bk}b4H`h$-r#$h=`0B%c345jDm~4-r{O1oh||;mXYQ5R9I_$ zX_CwNkK=~VPqKznD zl>?OjNT{iYHQ!qbQZn{s@n3w#Q2}qsy(_kPl2QE7IsYf`wM3$fY^JmM(<@L3$TE=K zri~a@_$zjQxdh+-sFP^SU@Y6^^*cM(w%5(9cwe${6<+J{jDlXs8#oI-MLKn{j|wP!^f{Oh!~Wk`KzWbbg@IRKdnJ196hQxSqr7- zxf}UK#J0$eSA~1uomP5sEYTjR!oKP+TE$6wh;cnMHWGlp)}Bt3Yy0f;!a9ecBznxh zi#B(ps!)gs7uP(FE*bUaHu^-|x&L*d5%y{C zJ@Ozei-$L>5EQEJmys`DGE{_g6Q$s>FVimg`s>Sw1m4<}LI$c`B_?0XZnv#p;@A7=_m^V?jyZn{sLrw47W6+kfG&5|p>Hb*i)CiWe1Fa#x~rYy@4OUxw_E>ao!gVmqj4$U-+ViN z+$bW#punZpVHxfsBEsj6pSOt{VpUyIPEW*86z)~v^9`1whRAB%>-*g`FcZj-RZ0~6 z!}Ui_H&-J;)DhGU!R9OPH`*P=$x>ak$IelaD3lnv~?^LN6Vw5Pec_=hF_+~`Zm z#%)Y%VhJV4!>ih%I!6u2JSBworgk7jaueRx9WArf;-)f#pBkVOLB2*x7|2%h{!TibP z`MW2BmDRFY=AAdmbY}!~OwE99dMD;yj>&-ITeigfwieHMk0F-~lXtu>eRp(rRb}9V zzK=+HJaYQis6+~O{oh-er=gGc>MAL=6M$R6DV(LLwAFr}LMP3EnIz1&bcW$f zUUP3IZ1_dxaRTu}73H#HBRDz!3Oyzs7`f{-L-h&=-1_*!R-C-d+JyT{k=M)49-X{& zWq*71!|G)r*4B2$Y4h}`!dXr6mgB!4rjJ#AT}7DefaybEtDUk4tVdx$gf z3pl$si8ZH>aAs!u(-K=Ze$GHm0d5iy5A4p@-2Dw{o(lHx3%!}?T)$}~R9vEAv#-tL zA?r?wo`CRqLE1GQsmjjh27Y)MoBa;eM(=K)zb5_KTt{Fv|LjgT z;P62AGwRQ$Z6#A^QcIc#(v0ux=vjl$j&2S+m4-h4Hb#5%mD(CT7yfYIS@!AF})5(0_J?ZBHhf3a8= zNjd9#=I>reCjYBA$CwX4nS04cwfE4va<6=9)^FS2EO%GAf#N+xDp?!p%^F@uhZns) S{m%CC&O=E|H|4)??*9M`G|gH7 From 90f5f700e686fcd473deca09ab3db43fa7c0b7fe Mon Sep 17 00:00:00 2001 From: Zachary Susswein Date: Mon, 2 Dec 2024 14:17:07 -0500 Subject: [PATCH 6/6] One more pass for language --- R/data.R | 27 +++++++++++++-------------- man/stochastic_sir_rt.Rd | 27 +++++++++++++-------------- 2 files changed, 26 insertions(+), 28 deletions(-) diff --git a/R/data.R b/R/data.R index c7150cb..ad73592 100644 --- a/R/data.R +++ b/R/data.R @@ -5,15 +5,14 @@ #' considerations for measuring the effective reproductive number, Rt." #' PLoS Computational Biology 16.12 (2020): e1008409. The data are simulated #' from a stochastic SEIR compartmental model. The original timeseries and Rt -#' are available in the `incidence` and `true_rt` columns and with other -#' added or modified to increase noise and add a day-of-week effect. +#' are available in the `incidence` and `true_rt` columns and with additional +#' columns added or modified to increase noise and add a day-of-week effect. #' #' This synthetic dataset has a number of desirable properties: #' #' 1. The force of infection changes depending on the Rt, allowing for sudden #' changes in the Rt. This allows for modeling of sudden changes in infection -#' dynamics, which might otherwise be difficult to capture. Rt estimation -#' framework +#' dynamics, which might otherwise be difficult to capture. #' #' 2. The realized Rt is known at each timepoint #' @@ -24,24 +23,24 @@ #' frameworks, providing practical guidance on how to use this dataset to #' evaluate Rt estimates. #' -#' In practice, we've found that the amount of observation noise in the -#' incidence and/or observed cases is often undesirably low for testing. Many -#' empirical datasets are much noisier. As a result, models built with these -#' settings in mind can perform poorly on this dataset or fail to converge. We -#' manually add observation noise and a day-of-week reporting effect. +#' In practice, we've found that the amount of observation noise is often +#' undesirably low for testing. Many empirical datasets are much noisier. As a +#' result, models built with these realistically noisy settings in mind can +#' perform poorly on this dataset or fail to converge. To better reflect realistic +#' settings, we manually add observation noise and a day-of-week reporting effect. #' #' @name stochastic_sir_rt #' @format `stochastic_sir_rt` A data frame with 299 rows and 6 columns: #' \describe{ -#' \item{reference_date}{The date the cases were observed.} +#' \item{reference_date}{The date cases were observed.} #' \item{true_rt}{The known, true Rt of the epidemic system.} #' \item{dow}{The magnitude of the day-of-week effect added to the log of `incidence`.} -#' \item{incidence}{The true number of cases occurring on the `date` in the -#' simulated system.} +#' \item{true_cases}{The true number of cases occurring on the `date` in the +#' simulated system before observation noise or day-of-week effects.} #' \item{obs_cases}{The observed number of cases on `date` after the day-of-week -#' reporting effect and observation noise have been applied to `incidence`.} +#' reporting effect and observation noise have been applied to `true_cases`.} #' \item{obs_cases_no_dow}{The observed number of cases on `date` after observation -#' noise has been applied to `incidence`. It does not include a day of week effect.} +#' noise has been applied to `true_cases`. It does not include a day of week effect.} #' } #' @source #' diff --git a/man/stochastic_sir_rt.Rd b/man/stochastic_sir_rt.Rd index 08724d8..595f459 100644 --- a/man/stochastic_sir_rt.Rd +++ b/man/stochastic_sir_rt.Rd @@ -7,15 +7,15 @@ \format{ \code{stochastic_sir_rt} A data frame with 299 rows and 6 columns: \describe{ -\item{reference_date}{The date the cases were observed.} +\item{reference_date}{The date cases were observed.} \item{true_rt}{The known, true Rt of the epidemic system.} \item{dow}{The magnitude of the day-of-week effect added to the log of \code{incidence}.} -\item{incidence}{The true number of cases occurring on the \code{date} in the -simulated system.} +\item{true_cases}{The true number of cases occurring on the \code{date} in the +simulated system before observation noise or day-of-week effects.} \item{obs_cases}{The observed number of cases on \code{date} after the day-of-week -reporting effect and observation noise have been applied to \code{incidence}.} +reporting effect and observation noise have been applied to \code{true_cases}.} \item{obs_cases_no_dow}{The observed number of cases on \code{date} after observation -noise has been applied to \code{incidence}. It does not include a day of week effect.} +noise has been applied to \code{true_cases}. It does not include a day of week effect.} } } \source{ @@ -31,16 +31,15 @@ A simulated dataset derived from Gostic, Katelyn M., et al. "Practical considerations for measuring the effective reproductive number, Rt." PLoS Computational Biology 16.12 (2020): e1008409. The data are simulated from a stochastic SEIR compartmental model. The original timeseries and Rt -are available in the \code{incidence} and \code{true_rt} columns and with other -added or modified to increase noise and add a day-of-week effect. +are available in the \code{incidence} and \code{true_rt} columns and with additional +columns added or modified to increase noise and add a day-of-week effect. } \details{ This synthetic dataset has a number of desirable properties: \enumerate{ \item The force of infection changes depending on the Rt, allowing for sudden changes in the Rt. This allows for modeling of sudden changes in infection -dynamics, which might otherwise be difficult to capture. Rt estimation -framework +dynamics, which might otherwise be difficult to capture. \item The realized Rt is known at each timepoint \item The dataset incorporates a simple generation interval and a reporting delay. @@ -50,11 +49,11 @@ Gostic et al. benchmark the performance of a number of Rt estimation frameworks, providing practical guidance on how to use this dataset to evaluate Rt estimates. -In practice, we've found that the amount of observation noise in the -incidence and/or observed cases is often undesirably low for testing. Many -empirical datasets are much noisier. As a result, models built with these -settings in mind can perform poorly on this dataset or fail to converge. We -manually add observation noise and a day-of-week reporting effect. +In practice, we've found that the amount of observation noise is often +undesirably low for testing. Many empirical datasets are much noisier. As a +result, models built with these realistically noisy settings in mind can +perform poorly on this dataset or fail to converge. To better reflect realistic +settings, we manually add observation noise and a day-of-week reporting effect. } \references{ Gostic, Katelyn M., et al. "Practical considerations for measuring the