\documentclass[12pt,a4paper]{article}
\usepackage[latin1]{inputenc}
\usepackage[none]{hyphenat}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{graphicx}
\usepackage{float}
\usepackage{amsthm}
\usepackage{framed}
\usepackage{cite}
\usepackage{url}
\usepackage{color}
\usepackage{multirow}
\usepackage{lscape}

%\definecolor{gray}{RGB}{160,160,160}
\definecolor{gray}{RGB}{0,0,0}
\renewcommand{\baselinestretch}{1.5}


\newcommand{\Bin}{\text{Binom}}
\newcommand{\Mun}{\text{Multinom}}



\newcommand{\rM}{\alpha_M}
\newcommand{\rH}{\alpha_H}
\newcommand{\rI}{\alpha_I}






\author{Shubin M., Virtanen M., Toikkanen S., Lyytik\"ainen O., Auranen K.}
\title{Estimating the burden of  A(H1N1)pdm09 influenza in Finland during two seasons}
\begin{document}

Epidemiology and Infection

\begin{large}
Estimating the burden of  A(H1N1)pdm09 influenza in Finland during two seasons
\end{large}

M. Shubin, M. Virtanen, S. Toikkanen, O. Lyytik\"ainen, K. Auranen.

\begin{large}
Supplementary Material
\end{large}


\paragraph*{Model representation.}
 
In the main article, the complete-data likelihood for a given stratum (age group$\times$region) in one season is presented as a product of binomials:

\begin{footnotesize}
\begin{align*}
&P(M, m, H, h, I, i \vert  S; p, s, g, \rM, \rI, \rH)   
= \\
&\text{Binom}(M+H+I \vert S;, p) \times
\text{Binom}(H+I \vert M+H+I; s) \times
\text{Binom}(I \vert H+I; g)  \\
&\text{Binom}(m \vert M; \rM) \times
\text{Binom}(h \vert H; \rH) \times
\text{Binom}(i \vert I; \rI). 
\end{align*}
\end{footnotesize}
This likelihood is equivalent to a multinomial model with 7 entries (6 free entries) (Figure S1):

\newcommand{\sm}{\!\! - \!\!} % minus sign with no spacing

\begin{footnotesize}
\begin{align*}
&P(M, m, H, h, I, i \vert  S; p, s, g, \rM, \rI, \rH)   
=  \Mun( U,\      m,\     M\sm m,\     h,\     H \sm h,\     i,\     I \sm i  
   \vert \\
   &S; \ 1 \sm p  ,\  p(1 \sm s) \rM,\    p(1 \sm s)(1 \sm \rM),\     ps(1 \sm g) \rH,\     ps(1\sm g)(1 \sm \rH),\   psg\rI,\     psg(1\sm\rI ) ).
\end{align*}
\end{footnotesize}
Here $U=S-M-H-I$ stands for uninfected individuals. The equivalence  of	the two model formulations can be proved using the rule for combining the marginals of the multinomial:

\begin{align*}&\Mun \left(a_1, a_2, a_3, \dots ,a_n  \Bigg\vert \sum_{i=1}^n a_i; p_1, p_2, p_3, \dots ,p_n\right)  = \\
%&\dfrac{ (a_1+a_2+a_3 \dots + a_n)! }{a_1!a_2!a_3! \dots  a_n!}p_1^{a_1}p_2^{a_2}p_3^{a_3}\dots
%p_n^{a_n} = \\
%&\dfrac{ (a_1+a_2+a_3 \dots + a_n)! }{(a_1+a_2)!a_3! \dots  a_n!}(p_1+p_2)^{a_1+a_2}p_3^{a_3}\dots
%p_n^{a_n} 
%\dfrac{ (a_1+a_2)! }{a_1!a_2!}
%\left(\frac{p_1}{p_1+p_2}\right)^{a_1}
%\left(\frac{p_2}{p_1+p_2}\right)^{a_2}
%= \\
&\Mun \left(a_1+ a_2, a_3 \dots ,a_n  \Bigg\vert \sum_{i=1}^n a_i; p_1+ p_2, p_3, \dots, p_n\right) 
\Bin \left(a_1 \Big\vert a_1+a_2 ;\dfrac{p_1}{p_1+p_2} \right)
\end{align*}


It can be clearly seen from  the multinomial formulation that the model has 3 observations and 6 parameters for a single stratum. The parameters define the outcome probabilities in a product form, which explains the hyperbolic posterior dependence between pairs of parameters (see Figure S2). However, the use of the Bayesian  model  with the assumption about the equality  of the parameters between different strata  helps ``borrowing strength" acrossstrata sharing the same age in the parameter estimation, in addition to the impact of the prior distributions.  


[figure S1 about here]

\paragraph*{The numbers of prevented cases. }

The expected number of mild infections for a given `age group$\times$region' stratum in the second season is

$$M^{(2)}=S^{(2)}p^{(2)}(1-s)=N(1-p^{(1)})(1-v)p^{(2)}(1-s).$$
The hypothetical number of infections in the absence of vaccination (assuming the vaccination had no effect of attack rate per susceptible in the second season $p^{(2)}$)  can be estimated as

$$M^{*(2)}=              N(1-p^{(1)})     p^{(2)}(1-s).$$
Then the number of prevented cases is

$$M^{*(2)}-M^{(2)}=
M^{(2)} \left( \dfrac{M^{*(2)}}{M^{(2)}} -1 \right) =
M^{(2)} \left( \dfrac{1}{1-v} -1 \right) =
M^{(2)} \dfrac{v}{1-v}
$$
The number of prevented hospitalized and IC infections was estimated in the similar way.


\paragraph*{Computational methods.}

The  posterior distribution of the model parameters was explored using a Gibbs sampler [1]. Due to the fact that the complete-data likelihood is a product of binomials the full conditional distribution of each of the parameters is a Beta distribution and the full conditional distribution of the other model unknowns (numbers of infections) is a Binomial distribution. 
For example, the full conditional for probability $s$ (hospitalization/infection ratio) is

$$
p(s_{a} \vert *) = \text{Beta} \left(s_{a}  \Bigg\vert\sum_{r}H_{a,r}+ \sum_{r}I_{a,r}  + \alpha , \sum_{r}M_{a,r}  + \beta \right)
$$
here $\alpha$ and $\beta$ are the parameters of the prior Beta distribution for $s$.
The number of unobserved hospitalized cases $H -h$ has the following distribution:

$$
p(H_{a,r}-h_{a,r} \vert *) =  \text{Binom} \left( H_{a,r}-h_{a,r} \Big\vert N_{a,r}-M_{a,r}-h_{a,r}-I_{a,r} ; \dfrac{ ps(1-\alpha_{Ha})}{ps(1-\alpha_{Ha}) + (1-p)} \right)  
$$

We run the algorithm for $10^7$ iterations  recording values at every $10^\text{th}$ iteration. The first 1\% samples were discarded as  ``burn-out" iterations. The Gibbs sampler was implemented in Python 2.7.1 with the NumPy external library for sampling from random distributions and the Python Image Library (PIL) for visualization. 



\paragraph*{Posterior uncertainty.}

Figure S2 presents an example of the joint posterior distribution of  parameters $p$, $s$, $\rM$, $g$ in one stratum (Helsinki region, age 15-20 years).  There was a clear posterior dependence among the parameters $p$, $s$ and $\alpha_M$.
Nevertheless, the peak area of the posterior was still identifiable.  Similar patterns were evident in all strata.

[figure S2 about here]

\paragraph*{Sensitivity to the prior distributions.}

Table S1 presents the impact of different prior choices on the estimated number of A(H1N1)pdm09 infections in the two seasons. 

[Table S1 about here]

\paragraph*{References}
$ $

 (1)	\textbf{ Geman S, Geman D}. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. \textit{IEEE transactions on pattern analysis and machine intelligence} 1984; \textbf{6}: 721-741.

 
\newpage 
 
Tables and figures
 
 
 

%\bi
%\includegraphics[width=400pt]{schemaS1.PNG}
%\ei

\textbf{Fig S1.} The model as a directed acyclic graph. The dotted lines present deterministic relations and the solid lines present stochastic relations. The rectangles and circles represent known and unknown parameters, respectively.

%\bi
%\includegraphics[width=400pt]{schemaS2.PNG}
%\ei

\textbf{Fig S2.}  Joint posterior distribution of parameters $p$, $s$, $\alpha_M$, $g$ in one stratum (Helsinki region, age 15-20 years)


\textbf{Table S1.}
The impact of prior distributions on the estimated number of all influenza A(H1N1)pdm09 infection in the two seasons. 
The results are obtained be changing one prior at a time. SD = standard deviation.  

\begin{scriptsize}

\begin{tabular}{llrrrrlrl}
Parameters& Prior & Mode & Mean & SD &  Total & & Severe&  \\
\hline
Base-case analysis & & & & &  258000 & & 3000&\\
\\
Incidence $p$ & Beta(2,4)  & 0.25 & 0.33 & 0.17 & base case & \\
 & Uniform(0,1)            & -    & 0.50 & 0.28 &  -8300 & (-3\%) & 0 \\
 & Beta(3,6)              & 0.28 & 0.33& 0.14 &  +12000 & (+5\%) & 0\\
\\
Severity $s$  & Beta(1.33,34)  & 0.01 & 0.04 & 0.032 &  base case &  \\
 & Uniform(0,1)            & -    & 0.5 & 0.28 &  0 &  & +50 &(+2\%)\\
 & Beta(1.13,15)               & 0.01 & 0.07& 0.06 &  0 & & +30 &(+1\%)\\
 &  Beta(1,50)              & 0 & 0.02& 0.02 &  +3700 & (+1\%)& -60 &(-2\%) \\
 & Beta(1,99)             & 0 & 0.01& 0.01 &  +8300 & (+3\%)& -140 &(-5\%)  \\
\\
IC/hospitalization ratio $q$ & Beta(5.3,40)  & 0.1 & 0.11 & 0.047 &  base case & \\
 & Uniform(0,1)            & -    & 0.5 & 0.28 &   0 & & +150 & (+5\%)\\
  & Beta(1,9)            &  0   & 0.1 & 0.09 &   0 & & +210 & (+7\%)\\
  & Beta(10,90)            & 0.09  & 0.1 & 0.029 &   0 & & -40 & (-4\%)\\
\\
Ascertainment  & Beta(1.33,34)  & 0.01 & 0.04 & 0.032 &  base case & \\
probability $\rM$  & Uniform(0,1)            & -    & 0.5 & 0.28 &  -234000 & (-90\%)& -290 & (-10\%) \\
of the mild cases & Beta(1.13,15)               & 0.01 & 0.07 & 0.06 &  -120000 & (-45\%) & -50 & (-2\%)  \\
 & Beta(1,50)               & 0 & 0.02& 0.02 &  +240000 & (+90\%) & +30 & (+1\%)  \\
 & Beta(1,99)               & 0 & 0.01& 0.01 &  +700000 & (+270\%) & +50 &(+3\%)\\
\\
Ascertainment  & Beta(9.6,3.9)  & 0.75 & 0.71 & 0.119 &  base case & \\
probability $\rH$  & Uniform(0,1)            & -    & 0.5 & 0.28 &  0 & & +180 & (+6\%) \\
of the hospitalized cases & Beta(50,1)      & 1 & 0.98&0.02 &  0 &  & -700 & (-24\%)  \\
& Beta(99,1)      & 1 & 0.99& 0.01  &  0 &   & -700 & (-24\%)  \\
\end{tabular}

\end{scriptsize}



\end{document}