L_{m,h}}\varphi_{k,-1}\big{(}\sqrt{v}\lambda,z\big{)}d\mu(z)=\sum_{b=0}^{\lfloor
p/2\rfloor}\frac{\operatorname{Tr}_{m,h}(L_{z}^{2b}f)}{(4\pi
v)^{b}b!}+\operatorname{Tr}^{\mathrm{c}}_{m,h}(f,v)+\\\
&+\operatorname{Tr}^{\mathrm{cc}}_{m,h}(f,v)-\sum_{\ell\in\Gamma\backslash\operatorname{Iso}(V)}\frac{\varepsilon_{\ell}i^{k}}{\sqrt{N}}a(\Theta_{k,\ell},m,h,v)\sum_{l=0}^{p}c_{\ell}(0,l)\phi_{0}(k-l,T)+O\big{(}\tfrac{1}{T}\big{)}.\end{split}$

###### Proof.

For $\lambda\in L_{m,h}$ we set
$\mathcal{H}_{\lambda,T}:=\mathcal{H}\setminus\big{(}B_{\epsilon}(\ell_{\lambda})\cup
B_{\epsilon}(\ell_{-\lambda})\big{)}$ for $\epsilon=e^{-2\pi T}$. The usual
unfolding argument, the fact that the stabilizer of $\lambda\in L_{m,h}$ is
trivial by Lemma 2.5, and the decay from Remark 3.15 allows us to replace, as
in Lemma 5.2 of [BF2], the left hand side by

$v^{\frac{1-k}{2}}\\!\sum_{\lambda\in\Gamma\backslash
L_{m,h}}\\!\int_{H_{T}}f(z)\varphi_{k,-1}\big{(}\sqrt{v}\lambda,z\big{)}d\mu(z)=v^{\frac{1-k}{2}}\\!\\!\sum_{\lambda\in\Gamma\backslash
L_{m,h}}\\!\int_{H_{\lambda,T}}\\!f(z)\varphi_{k,-1}\big{(}\sqrt{v}\lambda,z\big{)}d\mu(z)+O\big{(}\tfrac{1}{T}\big{)}$

(the error term here is, in fact, much better). The argument from the proof of
Proposition 4.9 expresses the term associated with $\lambda\in\Gamma\backslash
L_{m,h}$ as

(4.17) $\sum_{b=0}^{\lfloor p/2\rfloor}\frac{1}{(4\pi
v)^{b}b!}\int_{c(\lambda)\cap
Y_{T}}\\!(L_{z}^{2b}f)(z)(\lambda,Z(z))^{k-1-2b}dz+\sum_{\nu=0}^{p}v^{\frac{1-k}{2}}\int_{\partial\mathcal{H}_{\lambda,T}}\\!(L_{z}^{\nu}f)(z)\varphi_{k-\nu-1,\nu}\big{(}\sqrt{v}\lambda,z\big{)}dz.$

Next, $\partial\mathcal{H}_{\lambda,T}=\partial
B_{\epsilon}(\ell_{\lambda})\cup\partial B_{\epsilon}(\ell_{-\lambda})$ for
this $\epsilon$, both with the opposite orientation. In the integral along the
first part, as $\sigma_{\ell}^{-1}$ takes $\partial
B_{\epsilon}(\ell_{\lambda})$ to $\mathbb{R}+iT$, Equation (2.14) expresses
the corresponding term in Equation (4.17) as

$-\sum_{\nu=0}^{p}\int_{-\infty}^{\infty}v^{\frac{1-k}{2}}(L_{z}^{\nu}f\mid_{2k-2\nu}\sigma_{\ell_{\lambda}})(x_{\ell_{\lambda}}+iT)\varphi_{k-\nu-1,\nu}\bigg{(}\sqrt{\tfrac{mv}{N}}\big{(}\begin{smallmatrix}1&-2r_{\lambda}\\\
0&-1\end{smallmatrix}\big{)},x_{\ell_{\lambda}}+iT\bigg{)}dx_{\ell_{\lambda}}.$

Equation (3.19) thus implies that the integral along $\partial
B_{\epsilon}(\ell_{-\lambda})$ is evaluated in the same manner, with $\lambda$
replaced by $-\lambda$, and multiplied by $(-1)^{k}$.

Applying Equation (3.18), noting that the pairings with
$Z(x_{\ell_{\lambda}}+iT)$ and $Z^{\perp}(x_{\ell_{\lambda}}+iT)$ are
$2\sqrt{mv}(x_{\ell_{\lambda}}+iT-r_{\lambda})$ and
$\frac{2\sqrt{mv}}{T}(x_{\ell_{\lambda}}-r_{\lambda})$ respectively, we
substitute $\eta:=2\sqrt{2\pi mv}$ and
$\xi:=\frac{\eta}{T}(x_{\ell_{\lambda}}-r_{\lambda})$ and use Equations (4.5)
and (3.24) and the definition of the Fourier transform to present the latter
expression as

$\displaystyle-$ $\displaystyle
v^{\frac{1-k}{2}}\sum_{\nu=0}^{p}\int_{-\infty}^{\infty}\sum_{n\in\mathbb{Z}}\sum_{l=\nu}^{p}\frac{(-1)^{\nu}l!c_{\ell_{\lambda}}(n,l)}{(l-\nu)!T^{l-\nu}}e^{-2\pi
nT/\alpha_{\ell_{\lambda}}}\mathbf{e}\Big{(}\tfrac{nx_{\ell_{\lambda}}}{\alpha_{\ell_{\lambda}}}\Big{)}\frac{T^{k-\nu-1}}{(2\pi)^{k/2}}(\xi+i\eta)^{k-\nu-1}h_{\nu}(\xi)dx_{\ell_{\lambda}}$
$\displaystyle=$
$\displaystyle-v^{\frac{1-k}{2}}\sum_{n\in\mathbb{Z}}\sum_{l=0}^{p}\frac{l!c_{\ell_{\lambda}}(n,l)T^{k-l-1}}{(2\pi)^{k/2}}e^{-2\pi
nT/\alpha_{\ell_{\lambda}}}\int_{-\infty}^{\infty}\mathbf{e}\Big{(}\tfrac{nx_{\ell_{\lambda}}}{\alpha_{\ell_{\lambda}}}\Big{)}\mathbf{g}_{k-l,l}(\xi;\eta)dx_{\ell_{\lambda}}$
$\displaystyle=$
$\displaystyle-v^{\frac{1-k}{2}}\sum_{n\in\mathbb{Z}}\sum_{l=0}^{p}\frac{l!c_{\ell_{\lambda}}(n,l)T^{k-l}}{(2\pi)^{k/2}\eta}e^{-2\pi
nT/\alpha_{\ell_{\lambda}}}\mathbf{e}\Big{(}\tfrac{nr_{\lambda}}{\alpha_{\ell_{\lambda}}}\Big{)}\widehat{\mathbf{g}_{k-l,l}}\Big{(}-\tfrac{nT}{\alpha_{\ell_{\lambda}}\eta};\eta\Big{)}.$

Propositions 3.22 and 3.25 show, via Equation (4.4), that the contribution of
fixed $n$ and $l$ to the latter expression is $\frac{i^{k}\eta^{k-1}}{(2\pi
v)^{(k-1)/2}}c_{\ell_{\lambda}}(n,l)\mathbf{e}\big{(}\frac{nr_{\lambda}}{\alpha_{\ell_{\lambda}}}\big{)}=i^{k}\big{(}2\sqrt{m}\big{)}^{k-1}c_{\ell_{\lambda}}(n,l)\mathbf{e}\big{(}\frac{nr_{\lambda}}{\alpha_{\ell_{\lambda}}}\big{)}$
times

$\begin{cases}\frac{\operatorname{He}_{l}(\eta)}{\eta^{l}}\phi_{n}\Big{(}k-l,T;\frac{2\pi}{\alpha_{\ell_{\lambda}}}\Big{)},&n>0\text{
or }n<0\text{ and }l>k,\\\
\frac{\operatorname{He}_{l}(\eta)}{\eta^{l}}\phi_{n}\Big{(}k-l,T;\frac{2\pi}{\alpha_{\ell_{\lambda}}}\Big{)}-(-1)^{k}\sqrt{2\pi}\big{(}\frac{2\pi
n}{\alpha_{\ell_{\lambda}}}\big{)}^{l-k}\frac{l!}{(l-k)!}\frac{J_{l}(\eta)}{\eta^{l}},&n<0\text{
and }l\geq k,\\\
\big{(}\frac{\operatorname{He}_{l}(\eta)}{\eta^{l}}-\frac{\operatorname{He}_{k}(\eta)}{\eta^{k}}\big{)}\phi_{0}(k-l,T),&n=0\text{
and }l\neq k,\\\
\big{(}\frac{\operatorname{He}_{l}(\eta)}{\eta^{l}}-\frac{\operatorname{He}_{k}(\eta)}{\eta^{k}}\big{)}\phi_{0}(k-l,T)-\frac{(-1)^{k}k!}{\eta^{k}}\big{(}I_{k}(\eta)-\tilde{\Omega}_{k}(\eta)\big{)},&n=0\text{
and }l=k,\\\ \end{cases}$

where the error terms from the former proposition go into
$O\big{(}\frac{1}{T}\big{)}$ and the first term with $l=k$ trivially vanishes.

Now, since $Q(\lambda)=m$ for our $\lambda$, Lemma 4.7 shows that the sum over
$n$ and $l$ of the first terms gives $\sum_{b=0}^{\lfloor
p/2\rfloor}\frac{\operatorname{Sing}_{\lambda}(L_{z}^{2b}f)}{(4\pi v)^{b}b!}$.
Hence these terms, the corresponding ones from the integral along $\partial
B_{\epsilon}(\ell_{-\lambda})$, and the first term from Equation (4.17)
combine, via Equation (4.6), to $\sum_{b=0}^{\lfloor
p/2\rfloor}\frac{\operatorname{Tr}_{\lambda}(L_{z}^{2b}f)}{(4\pi v)^{b}b!}$,
and after summing over all $\lambda\in\Gamma\backslash L_{m,h}$ we get the
first asserted term by Equation (4.7). Now fix
$\ell\in\Gamma\backslash\operatorname{Iso}(V)$, and Equation (2.15) implies
that in the sum of the remaining terms over $\lambda\in\Gamma\backslash
L_{m,h}$, we only get contributions to the integral along $\partial
B_{\epsilon}(\ell_{\lambda})$ when $\iota_{\ell}(m,h)=1$, and to the one over
$\partial B_{\epsilon}(\ell_{-\lambda})$ if $\iota_{\ell}(m,-h)=1$, with the
sign $(-1)^{k}$. Lemma 2.6 and Remark 2.7 now imply that for such $\ell$ the
sum of $\mathbf{e}\big{(}\frac{nr_{\lambda}}{\alpha_{\ell_{\lambda}}}\big{)}$
over $\lambda$ with $\ell_{\lambda}=\ell$ gives
$2\varepsilon_{\ell}\sqrt{\frac{m}{N}}\mathbf{e}\big{(}\frac{nr_{\lambda}}{\alpha_{\ell_{\lambda}}}\big{)}$
in case $2\varepsilon_{\ell}\sqrt{\frac{m}{N}}$ divides $n$ and 0 otherwise.
The resulting sums thus produce the remaining required terms by Equations
(4.9), (4.10), and (2.19) and the value of $\eta$. This completes the proof of
the proposition. ∎

We can now state and prove the analog of Propositions 4.8, 4.9, and 4.10 for
$m=0$, again using the expressions from Equations (2.19) and (4.4), and with
the trace from Equation (4.11).

######  Proposition 4.11.

For an isotropic element $h\in D_{L}$ and large $T>1$ we have

$\begin{split}v^{\frac{1-k}{2}}\int_{Y_{T}}&f(z)\sum_{\lambda\in
L_{0,h}}\varphi_{k,-1}\big{(}\sqrt{v}\lambda,z\big{)}d\mu(z)=\sum_{b=0}^{\lfloor
p/2\rfloor}\frac{\operatorname{Tr}_{0,h}(L_{z}^{2b}f)}{(4\pi
v)^{b}b!}+\operatorname{Tr}^{\mathrm{cc}}_{0,h}(f,v)\\\
&-\sum_{\ell\in\Gamma\backslash\operatorname{Iso}(V)}\frac{i^{k}\varepsilon_{\ell}}{\sqrt{N}}a(\Theta_{k,\ell},0,h,v)\sum_{l=0}^{p}c_{\ell}(0,l)\phi_{0}(k-l,T)+O\big{(}\tfrac{1}{T}\big{)}.\end{split}$

###### Proof.

It is easy to see that
$L_{0,h}\setminus\\{0\\}=\bigcup_{\\{\ell\in\operatorname{Iso}(V)\,|\,\iota_{\ell}(0,h)=1\\}}[(L+h)\cap\ell]$,
a union that respects the $\Gamma$-action. The same unfolding argument from
the proofs of Propositions 4.9 and 4.10 allows us to write this part of the
left hand side as

$\sum_{\\{\ell\in\Gamma\backslash\operatorname{Iso}(V)\,|\,\iota_{\ell}(0,h)=1\\}}v^{\frac{1-k}{2}}\int_{\Gamma_{\ell}\backslash\mathcal{H}_{T}}f(z)\sum_{\lambda\in(L+h)\cap\ell}\varphi_{k,-1}\big{(}\sqrt{v}\lambda,z\big{)}d\mu(z),$

and Remark 3.15 implies that replacing
$\Gamma_{\ell}\backslash\mathcal{H}_{T}$ by
$\Gamma_{\ell}\backslash\mathcal{H}_{\ell,T}$, for
$\mathcal{H}_{\ell,T}:=\mathcal{H}_{T}\setminus B_{\epsilon}(\ell)$, produces
an error term that is much smaller that $O\big{(}\frac{1}{T}\big{)}$666Also
here we actually work with representatives for
$\Gamma\backslash\operatorname{Iso}(V)$, but we again allow this abuse of
notation.. As $\ell$ is perpendicular to neither $Z(z)$ nor $Z^{\perp}(z)$ for
any $z\in\mathcal{H}$, when we invoke Lemma 4.6 as the proof of Proposition
4.10, only (negatively oriented) integral along the boundary
$\Gamma_{\ell}\backslash\partial B_{\epsilon}(\ell)$ remains. Now, since
Equation (2.6) shows that the latter maps under $\sigma_{\ell}^{-1}$ onto
$(\mathbb{R}/\alpha_{\ell}\mathbb{Z})+iT$, Equation (2.8) expresses the
summand associated with $\ell$ as

$-\sum_{\nu=0}^{p}v^{\frac{1-k}{2}}\int_{\mathbb{R}/\alpha_{\ell}\mathbb{Z}}(L_{z}^{\nu}f\mid_{2k-2\nu}\sigma_{\ell})(x_{\ell}+iT)\sum_{0\neq\xi\in\mathbb{Z}+\omega_{\ell,h}}\varphi_{k-1-\nu,\nu}\Big{(}\sqrt{v}\big{(}\begin{smallmatrix}0&\beta_{\ell}\xi\\\
0&0\end{smallmatrix}\big{)},x_{\ell}+iT\Big{)}dx_{\ell}+O\big{(}\tfrac{1}{T}\big{)}.$

Here after applying Equation (3.18), the pairings with
$Z(x_{\ell_{\lambda}}+iT)$ and $Z^{\perp}(x_{\ell_{\lambda}}+iT)$ are just
$\sqrt{Nv}\beta_{\ell}\xi$ and $\frac{\sqrt{Nv}\beta_{\ell}\xi}{T}$
respectively. Setting $\upsilon:=\frac{\sqrt{2\pi Nv}\beta_{\ell}}{T}$,
Equations (4.5), (3.24), (3.34), and (3.35) and simple Fourier integration
evaluate the main term as

$\displaystyle-\sum_{\nu=0}^{p}v^{\frac{1-k}{2}}\int_{\mathbb{R}/\alpha_{\ell}\mathbb{Z}}\sum_{n\in\mathbb{Z}}\sum_{l=\nu}^{p}\frac{(-1)^{\nu}l!c_{\ell}(n,l)}{(l-\nu)!T^{l-\nu}}e^{-2\pi
nT/\alpha_{\ell}}\mathbf{e}\bigg{(}\frac{nx_{\ell}}{\alpha_{\ell}}\bigg{)}\sum_{0\neq\xi\in\mathbb{Z}+\omega_{\ell,h}}\frac{(T\upsilon\xi)^{k-1-\nu}}{(2\pi)^{k/2}}h_{\nu}(\upsilon\xi)dx_{\ell}$
$\displaystyle=\\!-\sqrt{v}\sum_{l=0}^{p}\frac{l!c_{\ell}(0,l)}{(2\pi
v)^{k/2}}\alpha_{\ell}T^{k-1-l}\\!\\!\\!\sum_{0\neq\xi\in\mathbb{Z}+\omega_{\ell,h}}\\!\\!\\!\mathbf{g}_{k-l,l}(\upsilon\xi;0)\\!=\\!-\sqrt{v}\sum_{l=0}^{p}\frac{l!c_{\ell}(0,l)}{(2\pi
v)^{k/2}}\alpha_{\ell}T^{k-1-l}\mathbf{G}_{k-l,l}(\omega_{\ell,h};\upsilon).$

We apply Proposition 3.26 and substitute the value of $\upsilon$, Equation
(2.7), and Corollary 3.7 for $P_{k}(0)$, and after summing over $\ell$ our
part of the left hand side takes the form

$\displaystyle\sum_{\begin{subarray}{c}\ell\in\Gamma\backslash\operatorname{Iso}(V)\\\
\iota_{\ell}(0,h)=1\end{subarray}}\sum_{l=0}^{p}\frac{l!c_{\ell}(0,l)}{(2\pi
v)^{(k-1)/2}}\beta_{\ell}\varepsilon_{\ell}\big{(}\sqrt{2\pi
Nv}\beta_{\ell}\big{)}^{k-1-l}\bigg{[}P_{l}(0)\Phi_{k-l}(\omega_{\ell,h})+Q_{l}(0)\Xi_{k-l}(\omega_{\ell,h})\bigg{]}$
$\displaystyle+\frac{i^{k}\operatorname{He}_{k}(0)}{(2\pi
v)^{(k-1)/2}}\sum_{\begin{subarray}{c}\ell\in\Gamma\backslash\operatorname{Iso}(V)\\\
\iota_{\ell}(0,h)=1\end{subarray}}\frac{\beta_{\ell}\varepsilon_{\ell}}{\sqrt{2\pi
Nv}\beta_{\ell}}\left(\sum_{\begin{subarray}{c}0\leq l\leq p\\\ l\neq
k\end{subarray}}{c_{\ell}(0,l)}\tfrac{T^{k-l}}{k-l}-c_{\ell}(0,k)\big{(}\log\tfrac{\sqrt{2\pi
Nv}\beta_{\ell}}{T}+C_{k}\big{)}\right)\\!\\!+\\!O\big{(}\tfrac{1}{T}\big{)},$

where the latter terms with $l>k$ can be artificially added from the error
term.

After canceling, writing $\iota_{\ell}(0,h)$ as a multiplier, and evaluating
$P_{l}(0)$ using Corollary 3.7, Equations (4.5) shows that the first terms
combine to the desired combination of the traces from Equation (4.8). Equation
(4.4) the expresses the remaining expression as the main term from Equation
(4.11) plus

$-\frac{i^{k}\operatorname{He}_{k}(0)}{(2\pi
v)^{k/2}}\sum_{\ell\in\operatorname{Iso}(V)}\iota_{\ell}(0,h)\frac{\varepsilon_{\ell}}{\sqrt{N}}\sum_{l=0}^{p}c_{\ell}(0,l)\phi_{n}(k-l,T),$

which is the last asserted term by Equation (2.19). This gives the desired
right hand side when $h\neq 0$, and for $h=0$ there is also the integral
involving $\varphi_{k,-1}\big{(}0,z\big{)}$ from $0\in L_{0,0}$. But Equation
(3.18) evaluates it as $0^{k}h_{-1}(0)=\delta_{k,0}$, and since in weight 0 we
have, by Equation (4.3), the bound

$\int_{Y}^{{\mathrm{reg}}}f(z)d\mu(z)-\int_{Y_{T}}f(z)d\mu(z)=\sum_{\ell\in\Gamma\backslash\operatorname{Iso}(V)}\int_{T}^{\infty}\alpha_{\ell}\sum_{l=0}^{p}\frac{c_{\ell}(0,l)}{y_{\ell}^{l+2}}=O\big{(}\tfrac{1}{T}\big{)},$

we indeed get the remaining term from Equation (4.11). This proves the
proposition. ∎

## References

  * [ANS] Alfes-Neumann, C., Schwagenscheidt, M., Shintani Theta Lifts of Harmonic Maass Forms, to appear in Trans. Amer. Math. Soc., https://doi.org/10.1090/tran/8265.
  * [AE] Alfes, C., Ehlen, S., Twisted Traces of CM values of Harmonic Weak Maass Forms, J. Number Theory, vol 133 issue 6, 1827–1845 (2013).
  * [Bo] Borcherds, R. E., Automorphic Forms with Singularities on Grassmannians, Invent. Math., vol. 132, 491–562 (1998).
  * [BGK] Bringmann, K., Guerzhoy, P., Kane, B., Shintani Lifts and Fractional Derivatives for Harmonic Weak Maass Forms, Adv. Math., vol 255, 641–-671, (2014).
  * [BF1] Bruinier, J. H., Funke, J., On Two Geometric Theta Lifts, Duke Math J., vol 125 no. 1, 45–90 (2004).
  * [BF2] Bruinier, J. H., Funke, J., Traces of CM Values of Modular Functions, J. Reine Angew. Math., vol 594, 1–33 (2006).
  * [BFI] Bruinier, J. H., Funke, J., Imamoğlu, Ö, Regularized Theta Liftings and Periods of Modular Functions, J. reine angew. Math., vol 703, 43–93 (2015).
  * [BFIL] Bruinier, J. H., Funke, J., Imamoğlu, Ö, Li, Y., Modularity of Generating Series of Winding Numbers, Res. Math. Sci., vol 5 paper 8 (2018).
  * [BGHZ] Bruinier, J., van der Geer, G., Harder, G., Zagier, D., The 1-2-3 of Modular Forms, Universitext (2008).
  * [BO] Bruinier, J. H., Ono, K., Heegner Divisors, $L$-Functions, and Harmonic Weak Maass Forms, Ann. of Math., vol 172, 2135–2181 (2010).
  * [DI] Duke, W., Imamoğlu, Ö, A Converse Theorem and the Saito–Kurokawa Lift, Int. Math. Res. Not., vol 1996 issue 7, 357–355 (1996).
  * [DIT] Duke, W., Imamoğlu, Ö, Tóth, Á. Cycle Integrals of the $j$-Function and Mock Modular Forms, Ann. of Math., vol 173 issue 2, 947–981 (2011).
  * [EMOT] Erdélyi, A., Magnus, W., Oberhettinger, F., Tricomi, F., Higher Transcendetal Functions, McGraw–Hill (1953).
  * [FH] Funke, J., Hofmann, E.,The Construction of Green Currents and Singular Theta Lifts for Unitary Groups, to appear in Trans. Amer. Math. Soc., https://doi.org/10.1090/tran/8289 (2020).
  * [GKZ] Gross, B., Kohnen, W., Zagier, D., Heegner Points and Derivatives of $L$-Series, II, Math. Ann., vol 278 no. 1–4, 497–562 (1987).
  * [HZ] Hirzebruch, F., Zagier, D., Intersection Numbers of Curves on Hilbert Modular Surfaces and Modular Forms of Nebentypus, Invent. Math., vol. 36, 57–114 (1976).
  * [JKK1] Jeon, D., Kang, S.-Y., Kim, C. H., Weak Maass–Poincaré Series and Weight 3/2 Mock Modular Forms, J. Number Theory, vol 133 issue 8, 2567–2587 (2013).
  * [JKK2] Jeon, D., Kang, S.-Y., Kim, C. H., Cycle Integrals of a Sesqui-Harmonic Maass Form of Weight Zero, J. Number Theory, vol 141, 92–108 (2013).
  * [K] Kohnen, W., Modular Forms of Half-Integral Weight on $\Gamma_{0}(4)$, Math. Ann., vol 248, 249–-266 (1980).
  * [L] Li, Y., Average CM-values of Higher Green’s Function and Factorization, pre-print, https://arxiv.org/abs/1812.08523 (2018).
  * [LZ] Li, Y., Zemel, S., Shimura Lift of Weakly Holomorphic Modular Forms, Math. Z., vol 290, 37–61 (2018).
  * [MR] Martin, F., Royer, E., Formes Modulaires et Périodes, Formes Modulaires et Transcendance, Sémin. Congr., vol. 12, Soc. Math. France, Paris, 1–-117 (2005).
  * [Sch] Scheithauer, N. R., The Weil Representation of $SL_{2}(\mathbb{Z})$ and Some Applications, Int. Math. Res. Not., no. 8, 1488–1545 (2009).
  * [Sh] Shimura, G., On Modular Forms of Half Integral Weight, Ann. of Math., vol. 97, 440–481 (1973).
  * [Sn] Shintani, T., On the Construction of Holomorphic Cusp Forms of Half-Integral Weight, Nagoya Math. J., vol 58, 83–126 (1975).
  * [Str] Strömberg, F., Weil Representations Associated to Finite Quad-ratic Modules, Math. Z., vol 275 issue 1, 509–527 (2013).
  * [Za] Zagier, D., Traces of Singular Moduli, In _Motives, Polylogarithms and Hodge theory, Part I_ (Irvine, CA, 1998), volume 3 of Int. Press Lect. Ser., pages 211-244. Int. Press, Somervillle, MA (2002).
  * [Ze1] Zemel, S., A $p$-adic Approach to the Weil Representation of Discriminant Forms Arising from Even Lattices, Math. Ann. Québec, vol 39 issue 1, 61–89 (2015).
  * [Ze2] Zemel, S., A Gross–Kohnen–Zagier Type Theorem for Higher-Codimensional Heegner Cycles, Res. Number Theory, vol 1 paper 25, 1–44 (2015).
  * [Ze3] Zemel, S. On Quasi-Modular Forms, Almost Holomorphic Modular Forms, and the Vector-Valued Modular Forms of Shimura, Ramanujan J., vol 37 issue 1, 165–180 (2015).
  * [Ze4] Zemel, S., Regularized Pairings of Meromorphic Modular Forms and Theta Lifts, J. Number Theory, vol 162, 275–311 (2016).
  * [Ze5] Zemel, S., Generalized Riordan Groups and Operators on Polynomials, Linear Algebra Appl., vol 494, 286–308 (2016).
  * [Ze6] Zemel, S., Normalizers of Congruence Groups in $SL_{2}(\mathbb{R})$ and Automorphisms of Lattices, Int. J. Number Theory, vol 13 issue 5, 1275–1300 (2017).
  * [Ze7] Zemel, S., Shimura’s Vector-Valued Modular Forms, Weight Changing Operators, and Laplacians, Res. Number Theory, vol 4 issue 1, paper 8 (2018).

Fachbereich Mathematik, AG 5, Technische Universität Darmstadt,
Schlossgartenstrasse 7, D-64289, Darmstadt, Germany

E-mail address<EMAIL_ADDRESS>
Einstein Institute of Mathematics, the Hebrew University of Jerusalem, Edmund
Safra Campus, Jerusalem 91904, Israel

E-mail address<EMAIL_ADDRESS>


Assessing dengue fever risk in Costa Rica by using climate variables and
machine learning techniques

Luis A. Barboza 1*, Shu-Wei Chou 2, Paola Vásquez3¤, Yury E. García3,4, Juan
G. Calvo 1‡, Hugo G. Hidalgo 5‡, Fabio Sanchez 1

1 Centro de Investigación en Matemática Pura y Aplicada - Escuela de
Matemática, Universidad de Costa Rica, San José, Costa Rica

2 Centro de Investigación en Matemática Pura y Aplicada - Escuela de
Estadística, Universidad de Costa Rica, San José, Costa Rica.

3 Centro de Investigación en Matemática Pura y Aplicada, Universidad de Costa
Rica, San José, Costa Rica.

4 Department of Public Health Sciences, University of California Davis, CA,
USA.

5 Centro de Investigaciones Geofísicas and Escuela de Física, Universidad de
Costa Rica, San José, Costa Rica.

These authors contributed equally to this work.

‡These authors also contributed equally to this work.

¤Current Address: Dept/Program/Center, Institution Name, City, State, Country

†Deceased

¶Membership list can be found in the Acknowledgments section.

*<EMAIL_ADDRESS>
## Abstract

Dengue fever is a vector-borne disease mostly endemic to tropical and
subtropical countries that affect millions every year and is considered a
significant burden for public health. Its geographic distribution makes it
highly sensitive to climate conditions. Here, we explore the effect of climate
variables using the Generalized Additive Model for location, scale, and shape
(GAMLSS) and Random Forest (RF) machine learning algorithms. Using the
reported number of dengue cases, we obtained reliable predictions. The
uncertainty of the predictions was also measured. These predictions will serve
as input to health officials to further improve and optimize the allocation of
resources prior to dengue outbreaks.

## Author summary

Dengue fever is a vector-borne disease that is endemic to tropical and
subtropical countries. It is transmitted by mosquitoes and affects
approximately 100 million people every year. Despite the fact that most
infections are mild or asymptomatic, some may cause severe disease known as
dengue shock syndrome or dengue hemorrhagic fever. Using statistical tools to
provide insight to public health officials regarding the risk of dengue in
different cantons of the country gives them the necessary time to prepare and
better allocate resources.

## Introduction

Dengue virus transmission represents a public health challenge for countries
in tropical and subtropical regions worldwide [1]. For the past decades, the
increasing geographical spread of the pathogen and its two main vectors, Aedes
aegypti and Aedes albopictus [2], has led to the development and
implementation of multiple prevention and control measures [3]. However, the
complex interactions and constant variations in population mobility,
socioeconomic, demographic, environmental, and climate factors that modulate
the spatial and temporal distribution of the disease make it difficult to
achieve timely, effective, and sustainable strategies.

In this context, researchers around the world are increasingly working towards
the development of innovative, tailored, and cost-effective tools that enhance
the design of public health policies for vector-borne diseases founded upon
the rapid systematization and analysis of information, as well as an increase
in interdisciplinary collaboration [4]. In these efforts, statistical and
machine learning techniques are increasingly used for public health
surveillance and epidemiological modeling [5]. These methods efficiently
estimate and predict health outcomes from epidemiological, entomological, and
other climate and social determinants surveillance. Through computational
algorithms, this branch of artificial intelligence facilitates the possibility
of integrating scientific knowledge, processing large databases, learning from
past documented reported cases, and ultimately projecting transmission
tendencies to identify and target the most vulnerable at-risk areas. Dengue is
a climate-sensitive disease, where changes in temperature, humidity, and
precipitation affect the biology, behavior, and availability of the mosquito
to reproduce, develop, propagate the virus, and interact with the human host
[6, 7, 8]. The use of satellite imagery and weather monitoring as input data
in machine learning models and other statistical learning approaches have
shown promising results [9, 10] that could effectively predict the relative
risk of dengue transmission.

In Costa Rica, a country of 5,213,362 inhabitants [11], administratively
divided into seven provinces and 82 cantons, there has been endemic
circulation of three of the four dengue virus serotypes (DENV-1, DENV-2,
DENV-3) and more than 393,000 reported cases by the Ministry of Health since
1993 [12]. The country has a tropical climate, which provides ideal conditions
for the mosquito vector to survive, replicate and transmit disease. The
interaction between geographic, oceanic, and atmospheric factors divides Costa
Rica into seven major climatic regions [13], each further subdivided into sub-
regions. Micro-climates separated by short distances make it crucial to
customize dengue transmission risk analysis in a country where health
authorities do not formally use weather information as input in developing
prevention and control activities.

In this article, we showed the results of using two different modeling
approaches, GAMLSS and RF, to forecast the relative risk of dengue infections
in 32 municipalities of Costa Rica. The analysis is a continuation of previous
work in [14], where an initial approach using Generalized Additive Models and
RF allowed us to retrospectively predict the relative risk of dengue for 2017
in five diverse climate municipalities, using as input the information of five
weather stations provided by the National Meteorological Institute [14]. Based
on those results, an initial collaboration with Costa Rican health authorities
allowed identifying 32 municipalities of interest based on their entomological
and epidemiological characteristics. This work also incorporates satellite
information to produce monthly dengue forecasts. We then evaluated the
forecast with available dengue reports from the Ministry of Health.

## Data description

The predictive models applied in this article used the following sources of
information:

Dengue Data:

Data of clinically suspected and confirmed monthly cases of dengue fever in
Costa Rica. The data is collected from all the local administrative areas
(municipalities or cantons) in the country, covering the years 2000-2021, and
provided by the Ministry of Health. The Ministry of Health has a particular
interest in the monthly forecast of 32 specific cantons due to their
particular epidemiological behavior. In order to quantify the relative
incidence of dengue cases at the $i$-th canton compared with the incidence in
the country at time $t$ (monthly basis), we use the relative risk ($RR$):

$\displaystyle
RR_{i,t}=\frac{\frac{\text{Cases}_{i,t}}{\text{Population}_{i,t}}}{\frac{\text{Cases}_{CR,t}}{\text{Population}_{CR,t}}}$

where $\text{Cases}_{i,t}$ ($\text{Population}_{i,t}$) and
$\text{Cases}_{CR,t}$ ($\text{Population}_{CR,t}$) are the number of observed
dengue cases (population size) at canton $i$ and country-level respectively,
at time $t$. The overall behavior of the relative risks, at three specific
months in 2013 (first row), and three specific years in July (second row) is
shown in Figure 1). Most of the studied cantons are located near the Pacific
and Caribbean Oceans, and there is inter-monthly and interannual variation in
the behavior of the relative risk.

Climate Data:

the covariate information consists of:

  * •

Daily Precipitation estimates ($P_{i,t}$) were used to index land surface
rainfall. Data was obtained from the Climate Hazards Group InfraRed
Precipitation with Station data (CHIRPS); see [15]. Due to the spatial high-
resolution nature of this dataset (5̃km by 5 km), we were able to compute
monthly cumulative rainfall estimates for each canton by adding the same
estimate over smaller administrative areas (distritos).

  * •

Weekly ENSO Sea Surface Temperature Anomaly ($S_{i,t}$), also known as SSTA
index. Data was obtained from the Climate Prediction Center (CPC) of the
United States National Oceanographic and Atmospheric Administration (NOAA).

  * •

Normalized Difference Vegetation Index (NVDI) ($N_{i,t}$), an index of the
greenness of vegetation for a 16-day time resolution and 250m spatial
resolution. It was obtained from the Moderate Resolution Imaging
Spectroradiometer (MODIS) satellite, and available through the `MODISTools` R
package (see [16])

  * •

Daytime Land Surface Temperature ($L_{i,t}$) in Kelvin degrees for a 8-day
time resolution and 1km spatial resolution. Obtained using the same resources
as the NDVI covariate.

  * •

Tropical Northern Atlantic Index ($TN_{i,t}$). Anomaly index of the sea-
surface temperature over the eastern tropical North Atlantic Ocean (see [17]).
This index is used, because previous work in the region, such as [18],
suggested that the inclusion of SST information from the Caribbean/Atlantic
improves performance compared to forecasts produced with only Pacific Ocean
ENSO conditions.

Fig 1: Relative Risk ($RR$) over the 32 cantons under study for different
months and years of available data.

## Methods

### Fitting Stage

As a first step to model the relationship among the climate covariates and the
relative risk of dengue fever at a fixed canton, we incorporate the historical
delayed associations between those variables by applying a Distributed Lag
Non-Linear Model (DLNM) framework [19, 20]. It consists of a bi-dimensional
space of functions that specifies a exposure-lag-response function $f\cdot
w(x,l)$ which depends on the predictor $x$ along the time lags $l$ in a
combined way. This combination allows the possibility of specifying a non-
linear and delayed association between the climate covariate and the dengue
incidence. For each covariate we take into account a maximum exposure of 18
months in its lag representation, based on the cross-correlation and wavelet
behavior among the series (see [21]) and a b-spline or linear basis
representation on the variable space. The procedure was performed with the R
package `dlnm` [22].

The structure of the models analyzed in this article in terms of dependent
variable and covariates for a fixed location is as follows:

$\displaystyle RR_{t}\sim
f(RR_{t-1},C_{1}P_{t},C_{2}S_{t},C_{3}N_{t},C_{4}L_{t},C_{5}TN_{t},M_{t})$ (1)

where $f$ is a function depending on the method employed, the matrices $C_{i}$
are defined in terms of the DLNM representation and $M_{t}$ is a factor-type
variable describing the monthly fixed effect. Recall that the time $t$ has
month as unit.

The first method that we use for $f$ is the generalized additive model for
location, scale and shape (GAMLSS). It represents a generalization of the GAM
method used in [14]. It is a flexible class of statistical framework where the
location, scale, skewness and kurtosis parameters from the distribution of the
response variable can be modeled as a additive function of covariates [23]. In
this particular case the model can be written as:

$\displaystyle RR_{t}$
$\displaystyle\overset{ind}{\sim}\mathcal{D}(\mu,\sigma,\nu)$ $\displaystyle
g_{1}(\mu)$
$\displaystyle=\beta_{10}+\beta_{11}RR_{t-1}+\beta_{12}C_{1}P_{t}+\beta_{13}C_{2}S_{t}+\beta_{14}C_{3}N_{t}+\beta_{15}C_{4}L_{t}+\beta_{16}C_{5}TN_{t}+\beta_{17}M_{t}$
$\displaystyle g_{2}(\sigma)$ $\displaystyle=\beta_{20}$ (2) $\displaystyle
g_{3}(\nu)$ $\displaystyle=\beta_{30}$

where the response variable is distributed as a three-parameter distribution
$\mathcal{D}$, $\mu$ is the location, $\sigma$ is the scale, $\nu$ is the
parameter related to the skewness of the distribution and $g_{i}$ for
$i=1,2,3$ are link functions.

Due to the fact that the relative risk of dengue is a non-negative skewed
variable with a large frequency of zeros ($16.2\%$), mixed distributions which
have positive domain with a positive probability at zero are appropriate
choices for modeling purposes. In this case, the zero adjusted gamma
distribution (ZAGA) and the zero adjusted inverse Gaussian (ZAIG) are
considered. Results for both choices are similar, and therefore we show our
results for the ZAGA distribution only.

The ZAGA density function is defined by the mixed continuous-discrete
probability density:

$f_{Y}(y)=\left\\{\begin{array}[]{ll}\nu\quad\text{~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}if
}y=0\\\ (1-\nu)f_{W}(y)\text{~{}~{}if~{}~{}}0<y<\infty\end{array}\right.$

for $0\leq y<\infty$, where $W\sim GA(\mu,\sigma)$ is a gamma distribution
with $0<\mu<\infty$, $0<\sigma<\infty$ and $0<\nu<1$, i.e.

$\displaystyle
f_{W}(y|\mu,\sigma)=\frac{1}{(\sigma^{2}\mu)^{1/\sigma^{2}}}\frac{y^{\frac{1}{\sigma^{2}}-1}e^{-y/(\sigma^{2}\mu)}}{\Gamma(1/\sigma^{2})},$

for $y>0$, $\mu>0$ and $\sigma>0$. The advantage of this parametrization is
that $E(W)=\mu$ and $V(W)=\sigma^{2}\mu^{2}$. For the GAMLSS specification,
$ZAGA(\mu,\sigma,\nu)$ defines the log link functions for $\mu$ and $\sigma$,
i.e., $g_{1}(\mu)=\log(\mu)$ and $g_{2}(\sigma)=\log\sigma$; and the logit
link function for $\nu$, i.e. $g_{3}(\nu)=\log\left[\nu/(1-\nu)\right]$.

The second method uses a Random Forest (RF) approach. This method is based on
the construction of bootstrapped ensemble of regression trees and combined
such that the prediction variance can be reduced (see [24] and [25]). One of
the main advantages of this method is the reduced number of tuning parameters
that eases its computational manipulation and stability [24].

The fitting process of the GAMLSS and RF models was performed with the R
packages `gamlss` [26] and `ranger` [27], respectively.

### Prediction Stage

Equation (1) has been fitted over a certain calibration period using any of
the two methodology alternatives, we forecast the relative risk over a testing
period using the past information of climatic covariates and the relative risk
itself. Provided that the response variable in (1) depends on the current
values of the climatic covariates, it is crucial to select an appropriate
method to obtain the climate predictions in the near future, so that it can
supply accurate inputs to our predictive model under both methodologies. Since
the climate covariates used in this study are highly correlated, a suitable
method to describe and predict their interaction is the vector auto-regressive
(VAR) model (see more details in [28]). For each canton, we include the trend
and seasonal factors to fit a VAR model and select the best lag order based on
BIC criterion, over the training period (which is the same period as the one
used for the fitting (1)). Then, we jointly forecast over the testing period
the covariates. Figure 2 shows an illustration of the observed climate
covariates along with its forecast values at Alajuela and Quepos cantons.
These predictions together with the predicted relative risks are used in order
to provide forecasts of the dependent variable over the testing period.
Finally, in order to assess the prediction uncertainty, we apply a non-
parametric bootstrap [29] and construct prediction intervals using the
corresponding forecasts for each bootstrap step.

### Model Comparison

We use two different metrics to compare the predictive performance of each
methodology for each fixed location. The normalized Mean-Squared Error
($NRMSE$):

$\displaystyle
NRMSE=\sqrt{\frac{1}{m\overline{RR}}\sum_{t=1}^{m}(RR_{t}-\widehat{RR}_{t})^{2}}$

where $m$ is the number of months in the testing period, $\overline{RR}$ is
the mean relative risk over the same period and $\widehat{RR}$ is the
estimated relative risk according to any of the two models. The normalized
Interval Score at $\alpha$ level ($NIS_{\alpha}$) is the normalized version of
the Interval Score (see [30] and [31]) so we can compare different locations
regardless the scale of their corresponding relative risk:

$\displaystyle
NIS_{\alpha}=\frac{1}{m\bar{R}R}\sum_{t=1}^{m}\left[(U_{t}-L_{t})+\frac{2}{1-\alpha}(L_{t}-RR_{t})\cdot
1_{RR_{t}<L_{t}}+\frac{2}{1-\alpha}(RR_{t}-U_{t})\cdot
1_{RR_{t}>U_{t}}\right]$

The latter metric is more complete than the former in the sense of evaluating
the predictive capacity of the models when the uncertainty is summarized
through an predictive interval [31].

## Results

We used the dengue and climate data as described above to fit the model in (1)
using both the GAMLSS and RF methodologies. The training period includes
monthly observations for the 32 cantons in the study, from January 2000 to
December 2020. This period was considered related to available satellite data
and epidemiological information. The DLNM basis was chosen to be linear in the
variable and lag space for the TNA index, Land Surface Temperature, and NDVI,
whereas a B-spline basis is assumed for the variable space, and it is linear
for the lag space for precipitation and ENSO. These choices allow an
acceptable balance among the complexity of the models and predictive precision
over all the locations.

Once the transformed covariates in (1) are determined, we fitted both
methodologies over the 32 locations individually. Then adjusted a VAR model
using the climate information for each location over the training period to
obtain predicted values of the covariates. Hence, we can predict the relative
risk of dengue for the first three months of 2021. Figure 2 shows the fitted
values and predictions of the climate covariates for Alajuela and Quepos. We
observed that the general features like trend and seasonality of the multiple
time series are well-captured for the testing period.

Fig 2: Observed climate covariates and forecast values at two specific
cantons: Alajuela (left panels) and Quepos (right panels). Black line:
observed climate covariates, red line: forecast values, and red shaded areas:
$95\%$ confidence regions.

Due to the auto-regressive nature of the model in (1), we used the predicted
value of $RR_{t}$ as a covariate in the prediction of $RR_{t+1}$. Once the
predicted relative risks over the first three months of 2021 are computed, we
compare the observed and the predicted values with the $NRMSE$ and $NIS_{.95}$
metrics and show the behavior of the best six cantons and worst three cantons
according to the latter metric, regardless of the method employed. The
comparison is shown for the training period in Figure 3

Fig 3: Comparison over the fitting period. Upper six panels: best cantons
according to NIS metric. Lower three panels: worst cantons according to NIS
metric. Black line: observed $RR$, red line: estimated $RR$.

and the testing period in Figure 4. Moreover, there is not a significant
difference among the methods in terms of the $NIS_{95}$ metric, as can be seen
in S1 Fig.. S1 Table. shows which model is chosen for each location and their
respective metrics, where in general, there is not a model that is predominant
over all the locations.

Together with the expected values of the relative risks in Figure 4, we
computed predicted uncertainties at $95\%$-level using a blocked non-
parametric bootstrap [29] with 100 replicates and a block size of six months,
over the testing period only.

Fig 4: Forecast comparison over the testing period. Upper six panels: best
cantons according to NIS metric. Lower three panels: worst cantons according
to NIS metric.

In general terms, the fitting performs well in the training period, except for
some extreme observations that the model does not capture closely, for
example, in the cantons of Quepos and Osa. The monthly trend over the testing
period is well captured in cantons like Orotina and Desamparados, and it is
partially captured as Alajuela and San Jose. The uncertainty is wide enough to
conserve the trend information while it covers most of the observed values
throughout the best six cantons. The worst-behaved cantons do not necessarily
have a significant uncertainty compared with the well-behaved ones, while the
bias is present in all the cases, where perhaps under-reporting of observed
cases can be a possible cause of such difference. Forecast comparison over the
testing period. Upper six panels: best cantons according to NIS metric. Lower
three panels: worst cantons according to NIS metric.

## Discussion

Climate conditions in Costa Rica make it ideal for the proliferation of
mosquitoes and, simultaneously, a challenging territory for the
epidemiological surveillance of dengue. The importance of climate information
regarding the incidence of dengue fever has been well established in previous
studies [7, 6].

The understanding of the ecology of the Ae. aegypti mosquito for the
development of climate-based early warning systems for dengue prediction is
crucial. In Costa Rica, dengue dynamics change geographically and temporally
due to microclimates. Traditionally, health officials carry out localized
studies and explore different methodologies to determine what works well
depending on geographic location. Furthermore, this is an opportunity for
health authorities to have a tool that provides situational insight into the
risk of dengue infection in the country.

Here, we implemented a Generalized Additive Model and a Random Forest Model
using mean temperature, precipitation, relative humidity, and weekly ENSO SSTA
to predict the relative risk of dengue in five of the 32 riskiest cantons in
this study. Although models showed a successful performance in the
predictions, their efficacy and capacity were limited by the availability and
quality of meteorological data. In addition, knowing the climatic covariates
was a necessary condition to predict dengue behavior. In this article, these
limitations were overcome using satellite information and proposing a VAR
model to obtain climate predictions, which are incorporated in the models to
forecast the relative risk.

Models were trained using data from 2000 to 2020. Results showed that the VAR
model performs well-predicting climatic variables, letting prospective
analyses be carried out instead of retrospective studies. The GAMLSS and the
Random Forest model successfully predict relative dengue risk in the first
three months of 2021, corresponding to the testing period. The GAMLSS
flexibility allows capturing the dynamic in cantons with a low number of
cases.

Dengue outbreaks affect the health system and the economy of a country. A
person with mild symptoms of dengue can be disabled on average for seven days
[de2013guide], which generates a reduction in the workforce affecting the
productivity of companies. The costs of the Costa Rican Social Security Fund
in care for users with dengue were estimated at $20.3 million for 2013, which
includes care for hospitalized patients, medical consultations, and
disabilities. The Ministry of Health has estimated at $6.5 million the
investment in preventive campaigns and combat actions in this same period
[32].

Dengue disease surveillance and vector surveillance are supposed to detect
outbreaks early, reducing the social and economic impact. The cost of
preventive measures has been reported to be less than the cost of treating an
outbreak [33]. Therefore, early prediction tools are valuable because they
allow taking preventive measures and directing and optimizing resources,
particularly in countries like Costa Rica, where economic resources are
limited.

The predictive capacity of machine learning models and the greater use of non-
traditional information in the development of monitoring efforts allows more
comprehensive surveillance and early identification of possible outbreaks, as
well as opens the possibility to diversify the inclusion of the community in
the identification of priority areas in collaboration with local governments
and authorities.

Public health authorities must be prepared to face the future challenges that
vector-borne diseases pose. External factors such as climate change and the
emergence of new infectious diseases highlight the need for a preventive
approach instead of a reactive one as the basis for the formulation of health
strategies in the country.

## Supporting information

#### S1 Fig.

Comparison of the distribution of NIS metric among methods.

#### S1 Table.

Best model for each canton

Canton | $NRMSE$ | $NIS_{95}$ | Best Model
---|---|---|---
Alajuela | 0.2218 | 2.2789 | RF
Alajuelita | 0.2650 | 13.9639 | GAMLSS
Atenas | 1.7615 | 16.6668 | GAMLSS
Cañas | 0.4572 | 8.8197 | GAMLSS
Carrillo | 0.7468 | 13.2383 | GAMLSS
Corredores | 7.4754 | 32.9063 | RF
Desamparados | 0.0333 | 0.9103 | RF
Esparza | 0.5269 | 13.5993 | RF
Garabito | 32.7929 | 151.7930 | GAMLSS
Golfito | 0.7657 | 13.0531 | RF
Guacimo | 1.9182 | 8.1629 | RF
La Cruz | 9.1306 | 127.2412 | RF
Liberia | 2.2955 | 57.4666 | RF
Limon | 2.1313 | 18.4764 | RF
Matina | 3.9168 | 31.3834 | GAMLSS
Montes de Oro | 18.9354 | 162.9258 | RF
Nicoya | 5.3693 | 47.4291 | GAMLSS
Orotina | 0.0860 | 0.6048 | RF
Osa | 2.4703 | 2.8810 | GAMLSS
Parrita | 8.0959 | 33.1650 | GAMLSS
Perez Zeledón | 1.7181 | 14.7112 | GAMLSS
Pococí | 0.6451 | 10.2410 | RF
Puntarenas | 1.7821 | 22.3798 | GAMLSS
Quepos | 34.3754 | 179.4324 | GAMLSS
San Jose | 0.0195 | 2.5938 | RF
Santa Ana | 0.6331 | 31.8002 | RF
SantaCruz | 0.7159 | 13.4637 | GAMLSS
Sarapiquí | 0.5572 | 1.7992 | RF
Siquirres | 0.8202 | 6.6540 | RF
Talamanca | 10.9225 | 26.2873 | RF
Turrialba | 0.9400 | 7.0685 | GAMLSS

## Acknowledgments

## References

  * 1. Brady OJ, Gething PW, Bhatt S, Messina JP, Brownstein JS, Hoen AG, et al.  Refining the global spatial limits of dengue virus transmission by evidence-based consensus.  PloS Neglected Tropical Diseases. 2012;6(8):e1760.  doi:10.1371/journal.pntd.0001760.
  * 2. Messina JP, Brady OJ, Scott TW, Zou C, Pigott DM, Duda KA, et al.  Global spread of dengue virus types: mapping the 70 year history.  Trends in Microbiology. 2014;22(3):138–146.
  * 3. Gubler DJ.  Dengue and dengue hemorrhagic fever.  Clinical microbiology reviews. 1998;11(3):480–496.
  * 4. Organization WH, et al.. Global strategy for dengue prevention and control 2012-2020; 2012.
  * 5. Wiemken TL, Kelley RR.  Machine Learning in Epidemiology and Health Outcomes Research.  Annual review of public health. 2019;41:21–36.
  * 6. Ebi KL, Nealon J.  Dengue in a changing climate.  Environmental research. 2016;151:115–123.
  * 7. Naish S, Dale P, Mackenzie JS, McBride J, Mengersen K, Tong S.  Climate change and dengue: a critical and systematic review of quantitative modelling approaches.  BMC infectious diseases. 2014;14(1):1–14.
  * 8. Tran BL, Tseng WC, Chen CC, Liao SY.  Estimating the threshold effects of climate on dengue: A case study of Taiwan.  International Journal of Environmental Research and Public Health. 2020;17(4):1392.
  * 9. Lowe R, Bailey TC, Stephenson DB, Graham RJ, Coelho CA, Carvalho MS, et al.  Spatio-temporal modelling of climate-sensitive disease risk: Towards an early warning system for dengue in Brazil.  Computers & Geosciences. 2011;37(3):371–381.
  * 10. Lowe R, Bailey TC, Stephenson DB, Jupp TE, Graham RJ, Barcellos C, et al.  The development of an early warning system for climate-sensitive disease risk with a focus on dengue epidemics in Southeast Brazil.  Statistics in medicine. 2013;32(5):864–883.
  * 11. INEC. Bienvenido a INSTITUTO NACIONAL DE ESTADISTICA Y CENSOS — INSTITUTO NACIONAL DE ESTADISTICA Y CENSOS; 2022.  Available from: https://www.inec.cr/.
  * 12. de Salud M. Sitio Web del Ministerio de Salud de Costa Rica. Bienvenido; 2022.  Available from: https://www.ministeriodesalud.go.cr/.
  * 13. IMN. Clima en Costa Rica - IMN; 2022.  Available from: https://www.imn.ac.cr/52.
  * 14. Vásquez P, Loría A, Sanchez F, Barboza LA.  Climate-driven Statistical Models as effective predictors of local Dengue incidence in Costa Rica: a Generalized Additive Model and Random Forest Approach.  Revista de Matematica: Teoria y Aplicaciones. 2020;27(1):1–21.
  * 15. Funk C, Peterson P, Landsfeld M, Pedreros D, Verdin J, Shukla S, et al.  The climate hazards infrared precipitation with stations–a new environmental record for monitoring extremes.  Scientific data. 2015;2(1):150066–150066.
  * 16. Tuck SL, Phillips HRP, Hintzen RE, Scharlemann JPW, Purvis A, Hudson LN.  MODISTools – downloading and processing MODIS remotely sensed data in R.  Ecology and Evolution. 2014;4(24):4658–4668.  doi:10.1002/ece3.1273.
  * 17. Enfield DB, Mestas-Nuñez AM, Mayer DA, Cid-Serrano L.  How ubiquitous is the dipole relationship in tropical Atlantic sea surface temperatures?  Journal of Geophysical Research: Oceans. 1999;104(C4):7841–7848.
  * 18. Hidalgo HG, Alfaro EJ, Quesada-Montano B.  Observed (1970–1999) climate variability in Central America using a high-resolution meteorological dataset with implication to climate change studies.  Clim Change. 2017;141(1):13–28.
  * 19. Gasparrini A, Armstrong B, Kenward MG.  Distributed lag non-linear models.  Statistics in Medicine. 2010;29(21):2224–2234.  doi:https://doi.org/10.1002/sim.3940.
  * 20. Gasparrini A.  Modeling exposure–lag–response associations with distributed lag non-linear models.  Statistics in Medicine. 2014;33(5):881–899.  doi:https://doi.org/10.1002/sim.5963.
  * 21. García YE, Barboza LA, Sanchez F, Vásquez P, Calvo JG. Wavelet Analysis of Dengue Incidence and its Correlation with Weather and Vegetation Variables in Costa Rica; 2021.
  * 22. Gasparrini A.  Distributed lag linear and non-linear models in R: the package dlnm.  Journal of Statistical Software. 2011;43(8):1–20.
  * 23. Stasinopoulos D, Rigby R, Heller G, Voudouris V, De Bastiani F.  Flexible regression and smoothing: using GAMLSS in R.  Chapman and Hall/CRC the R Series. Chapman & Hall/CRC; 2017.
  * 24. Breiman L.  Random Forests.  Machine learning. 2001;45(1):5–32.
  * 25. Hastie T, Tibshirani R, Friedman J.  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  2nd ed. Springer Series in Statistics. Springer; 2009.  Available from: http://www.amazon.com/The-Elements-Statistical-Learning-Prediction/dp/0387848576.
  * 26. Rigby RA, Stasinopoulos DM.  Generalized additive models for location, scale and shape,(with discussion).  Applied Statistics. 2005;54:507–554.
  * 27. Wright MN, Ziegler A.  ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R.  Journal of Statistical Software. 2017;77(1):1–17.  doi:10.18637/jss.v077.i01.
  * 28. Tsay RS.  Multivariate Time Series Analysis with R and Financial Applications; 2015.
  * 29. Efron B, Tibshirani R.  An Introduction to the Bootstrap.  Monographs on Statistics and Applied Probability. Chapman & Hall, CRC; 1993.
  * 30. Winkler RL, Murphy AH.  “Good” Probability Assessors.  Journal of Applied Meteorology and Climatology. 1968;7(5):751–758.  doi:10.1175/1520-0450(1968)007¡0751:PA¿2.0.CO;2.
  * 31. Gneiting T, Raftery AE.  Strictly Proper Scoring Rules, Prediction, and Estimation.  Journal of the American Statistical Association. 2007;102(477):359–378.  doi:10.1198/016214506000001437.
  * 32. de la República PG. Plan General de la Emergencia: Decreto Nº 39526-MP-S “Estado de Emergencia por la Proliferación del Vector del Dengue, Chikungunya y el Zika”; 2016.  Available from: http://www.pgrweb.go.cr/scij/Busqueda/Normativa/Normas/nrm_texto_completo.aspx?param1=NRTC&nValor1=1&nValor2=82834&nValor3=106102&strTipM=TC.
  * 33. Stahl HC, Butenschoen VM, Tran HT, Gozzer E, Skewes R, Mahendradhata Y, et al.  Cost of dengue outbreaks: literature review and country case studies.  BMC public health. 2013;13(1):1–11.