[Pomp-commits] r668 - in pkg/pompExamples: . data inst/data-R inst/doc src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 18 13:50:55 CEST 2012
Author: kingaa
Date: 2012-04-18 13:50:55 +0200 (Wed, 18 Apr 2012)
New Revision: 668
Removed:
pkg/pompExamples/inst/data-R/budmoth-params.rda
pkg/pompExamples/inst/data-R/pertussis-params.rda
Modified:
pkg/pompExamples/.Rbuildignore
pkg/pompExamples/.Rinstignore
pkg/pompExamples/data/pertussis.sim.rda
pkg/pompExamples/inst/data-R/budmoth.sim.R
pkg/pompExamples/inst/data-R/pertussis.sim.R
pkg/pompExamples/inst/doc/budmoth-model-slices.rda
pkg/pompExamples/inst/doc/budmoth-model-true-loglik.rda
pkg/pompExamples/inst/doc/budmoth-model.pdf
pkg/pompExamples/inst/doc/pertussis-model-true-loglik.rda
pkg/pompExamples/inst/doc/pertussis-model.Rnw
pkg/pompExamples/inst/doc/pertussis-model.pdf
pkg/pompExamples/src/pertussis.c
pkg/pompExamples/tests/pertussis.R
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- drop V from pertussis model
- drop finiteness checks in pertussis model
Modified: pkg/pompExamples/.Rbuildignore
===================================================================
--- pkg/pompExamples/.Rbuildignore 2012-04-18 10:01:33 UTC (rev 667)
+++ pkg/pompExamples/.Rbuildignore 2012-04-18 11:50:55 UTC (rev 668)
@@ -3,3 +3,4 @@
inst/doc/(.+?)\.bst$
inst/doc/(.+?)\.R$
inst/doc/(.+?)\.png$
+inst/doc/(.+?)\.pbs$
Modified: pkg/pompExamples/.Rinstignore
===================================================================
--- pkg/pompExamples/.Rinstignore 2012-04-18 10:01:33 UTC (rev 667)
+++ pkg/pompExamples/.Rinstignore 2012-04-18 11:50:55 UTC (rev 668)
@@ -2,3 +2,4 @@
inst/data-R/Makefile
inst/doc/(.+?)\.bst$
inst/doc/(.+?)\.rda$
+inst/doc/(.+?)\.pbs$
Modified: pkg/pompExamples/data/pertussis.sim.rda
===================================================================
(Binary files differ)
Deleted: pkg/pompExamples/inst/data-R/budmoth-params.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/data-R/budmoth.sim.R
===================================================================
--- pkg/pompExamples/inst/data-R/budmoth.sim.R 2012-04-18 10:01:33 UTC (rev 667)
+++ pkg/pompExamples/inst/data-R/budmoth.sim.R 2012-04-18 11:50:55 UTC (rev 668)
@@ -1,7 +1,17 @@
require(pompExamples)
-load("budmoth-params.rda")
+paramcsv <- "
+,alpha,sig.alpha,gam,lambda,sig.lambda,g,delta,a,sig.a,w,beta0,beta1,u,sigQobs,sigNobs,sigSobs,Q.0,N.0,S.0,seed
+tri,0.5,0.1,50,22,0.25,0.08,10,1.7,0.1,0.15,0,35,0.9,0.03,0.5,0.1,0.96,0.02,0.22,1691699385
+para1,0.5,0.1,50,22,0.25,0.08,0.5,1.7,0.1,0.15,0,35,0.9,0.03,0.5,0.1,0.96,0.02,0.22,1116757478
+food,0.5,0.1,20,5,0.25,0.02,10,1,0.1,0,0,35,0.9,0.03,0.5,0.1,0.96,0.02,0.22,1054866677
+para2,0.5,0.1,50,10,5,0.08,0.5,1.7,1,0.15,0,35,0.9,0.03,0.5,0.1,0.96,0.02,0.22,1361101458
+"
+tc <- textConnection(paramcsv)
+params <- read.csv(tc,row.names=1)
+close(tc)
+
datasets <- rownames(params)
seeds <- params["seed"]
params <- subset(params,select=-seed)
Deleted: pkg/pompExamples/inst/data-R/pertussis-params.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/data-R/pertussis.sim.R
===================================================================
--- pkg/pompExamples/inst/data-R/pertussis.sim.R 2012-04-18 10:01:33 UTC (rev 667)
+++ pkg/pompExamples/inst/data-R/pertussis.sim.R 2012-04-18 11:50:55 UTC (rev 668)
@@ -5,8 +5,22 @@
require(pompExamples)
-load("pertussis-params.rda")
+paramcsv <- "
+,birthrate,deathrate,mean.beta,ampl.beta,imports,sigma,gamma,alpha,alpha.ratio,report.prob,boost.prob,polar.prob,foi.mod,popsize,noise.sigma,tau,S.0,E.0,I.0,R1.0,R2.0,seed
+SEIR.small,0.02,0.02,450,0.15,10,46,26,0,1,0.3,0,0,0,500000,0,0.01,0.0574148031949802,0.0004081763321755,0.00067028956509212,0.941506730907752,0,1831650124
+SEIR.big,0.02,0.02,450,0.15,10,46,26,0,1,0.3,0,0,0,5000000,0,0.01,0.0515635231482973,0.000437143470487014,0.000734641109212043,0.947264692272004,0,908022490
+SEIRS.small,0.02,0.02,150,0.15,10,46,26,0.1,1,0.1,0,0,0,500000,0,0.01,0.157360395940609,0.000837874318852172,0.00124181372794081,0.45913512973054,0.381424786282058,1111340406
+SEIRS.big,0.02,0.02,150,0.15,10,46,26,0.1,1,0.1,0,0,0,5000000,0,0.01,0.157398354546347,0.00132093662562661,0.0022558671035406,0.457185201591761,0.381839640132725,1751228386
+SEIRR.small,0.02,0.02,150,0.15,10,46,26,0.1,1,0.11,0.75,0,0.5,500000,0,0.01,0.128943112158304,0.00068688724266688,0.00114414648269803,0.638074319602244,0.231151534514087,350421545
+SEIRR.big,0.02,0.02,150,0.15,10,46,26,0.1,1,0.11,0.75,0,0.5,5000000,0,0.01,0.127128689912424,0.00126497004491763,0.00216092385991776,0.639879739889535,0.229565676293206,748454784
+full.small,0.02,0.02,150,0.15,10,46,26,0.1,1,0.1,0.75,0.1,0.5,500000,0.01,0.01,0.132553922599906,0.0010539075727066,0.00166100642162314,0.641737544956371,0.222993618449393,581894515
+full.big,0.02,0.02,150,0.15,10,46,26,0.1,1,0.1,0.75,0.1,0.5,5000000,0.01,0.01,0.130980596244438,0.00115096693013597,0.0018994251960431,0.643957103848235,0.222011907781148,301057392
+"
+tc <- textConnection(paramcsv)
+params <- read.csv(tc,row.names=1)
+close(tc)
+
datasets <- rownames(params)
seeds <- params["seed"]
params <- subset(params,select=-seed)
@@ -27,21 +41,21 @@
PACKAGE="pompExamples",
paramnames=c(
"birthrate","deathrate","mean.beta","ampl.beta",
- "imports","sigma","gamma","alpha.1","alpha.2","alpha.ratio",
- "report.prob","boost.prob","polar.prob","vacc.prob",
+ "imports","sigma","gamma","alpha","alpha.ratio",
+ "report.prob","boost.prob","polar.prob",
"foi.mod","noise.sigma","popsize","tau"
),
- statenames=c("S","E","I","R1","R2","V","cases","W","err","simpop"),
+ statenames=c("S","E","I","R1","R2","cases","W","err","simpop"),
zeronames=c("cases","err"),
- ivps=c("S.0","E.0","I.0","R1.0","R2.0","V.0"),
- comp.names=c("S","E","I","R1","R2","V"),
+ ivps=c("S.0","E.0","I.0","R1.0","R2.0"),
+ comp.names=c("S","E","I","R1","R2"),
rmeasure = "negbin_rmeasure",
dmeasure = "negbin_dmeasure",
- logitvar=c("report.prob","boost.prob","polar.prob","vacc.prob"),
+ logitvar=c("report.prob","boost.prob","polar.prob"),
logvar=c(
"birthrate","deathrate","mean.beta","ampl.beta","imports","sigma","gamma",
- "alpha.1","alpha.2","alpha.ratio","foi.mod","noise.sigma","tau",
- "S.0","E.0","I.0","R1.0","R2.0","V.0"
+ "alpha","alpha.ratio","foi.mod","noise.sigma","tau",
+ "S.0","E.0","I.0","R1.0","R2.0"
),
initializer = function (params, t0, statenames, comp.names, ivps, ...) {
states <- numeric(length(statenames))
Modified: pkg/pompExamples/inst/doc/budmoth-model-slices.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/doc/budmoth-model-true-loglik.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/doc/budmoth-model.pdf
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/doc/pertussis-model-true-loglik.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/doc/pertussis-model.Rnw
===================================================================
--- pkg/pompExamples/inst/doc/pertussis-model.Rnw 2012-04-18 10:01:33 UTC (rev 667)
+++ pkg/pompExamples/inst/doc/pertussis-model.Rnw 2012-04-18 11:50:55 UTC (rev 668)
@@ -19,13 +19,6 @@
\setcounter{secnumdepth}{1}
\setcounter{tocdepth}{1}
-%% \begingroup\makeatletter\ifx\SetFigFont\undefined
-%% \gdef\SetFigFont#1#2#3#4#5{
-%% \reset at font\fontsize{#1}{#2pt}
-%% \fontfamily{#3}\fontseries{#4}\fontshape{#5}
-%% \selectfont}
-%% \fi\endgroup
-
\begingroup\makeatletter\ifx\SetFigFont\undefined
\gdef\SetFigFont#1#2{
\reset at font\fontsize{#1}{#2pt}
@@ -79,13 +72,13 @@
We will describe the models, specify how we have used them to generate simulated epidemiological data, lay out the ``bounds of our simulated ignorance'' for the fitting exercises, and do a bit of preliminary investigation of the likelihood surface.
In the simulation studies, we will focus on model selection and hypothesis testing.
This is simplified somewhat by the fact that the candidate models are nested;
-the supermodel in this case is the SVEIRR model.
+the supermodel in this case is the SEIRR model.
\tableofcontents
-\section{The SVEIRR equations}
+\section{The SEIRR equations}
-The SVEIRR system is a simpler alternative to the multiple-compartment model of \citet{Wearing2009}.
+The SEIRR system is a simpler alternative to the multiple-compartment model of \citet{Wearing2009}.
The following key assumptions are made:
\begin{enumerate}
@@ -93,27 +86,24 @@
With probability $\phi$, recovering individuals reactivate full susceptibility and enter the $S$ class.
This captures so-called ``polarising'' immunity.
-\item Vaccinated individuals (in $V$ class) are assumed to be fully protected.
+\item Individuals in the $R_1$ class may lose immunity and become susceptible to immune boosting upon entry to the $R_2$ class, with rate $\alpha_1\,R_1$.
-\item Individuals in the $R_1$ and $V$ classes may lose immunity and become susceptible to immune boosting upon entry to the $R_2$ class, with rates $\alpha_1\,R_1$ and $\alpha_2V$, respectively.
-
\item Within the $R_2$ class, exposure to infection is described by $\epsilon\,\lambda$, where $\epsilon$ modulates the force of infection in this class.
Crucially, those exposed --- with probability $\xi$ --- experience sub-clinical infections that boost immunity, moving affected individuals to the immune $R_1$ class ($R_2\,\rightarrow\,R_1$).
\item A fraction ($1-\xi$) of exposed individuals in the $R_2$ class are assumed to develop full infection, i.e., become fully infectious ($R_2\,\rightarrow\,E$).
-\item In the absence of immune boosting, individuals reactivate full susceptibilibty ($R_2\,\rightarrow\,S$) after a mean period of $1/\alpha_3$ in the $R_2$ class.
+\item In the absence of immune boosting, individuals reactivate full susceptibilibty ($R_2\,\rightarrow\,S$) after a mean period of $1/\alpha_2$ in the $R_2$ class.
\end{enumerate}
-The mean field equations describing this the SVEIRR system are given by:
+The mean field equations describing this the SEIRR system are given by:
\begin{align}
-\dot S&=\nu\,(1-p)\,N-\lambda(t)\,S-\mu\,S+\phi\,\gamma\,I+\alpha_3\,R_2\\
-\dot V&=\nu\,p\,N-(\mu+\alpha_2)\,V\\
+\dot S&=\nu\,\,N-\lambda(t)\,S-\mu\,S+\phi\,\gamma\,I+\alpha_2\,R_2\\
\dot E&=\lambda(t)\,S+(1-\xi)\,\epsilon\,\lambda(t)\,R_2-(\mu+\sigma)\,E\\
\dot I&=\sigma\,E-(\mu+\gamma)\,I\\
\dot R_1&=(1-\phi)\,\gamma\,I-(\mu+\alpha_1)\,R_1+\xi\,\epsilon\,\lambda(t)\,R_2\\
-\dot R_2&=\alpha_1\,R_1-\left(\mu+\alpha_3+\epsilon\,\lambda(t)\right)\,R_2+\alpha_2\,V
+\dot R_2&=\alpha_1\,R_1-\left(\mu+\alpha_2+\epsilon\,\lambda(t)\right)\,R_2
\end{align}
In these equations, the force of infection, $\lambda$, is time dependent and given by
\begin{equation}
@@ -152,18 +142,18 @@
{\color[rgb]{0,0,0}
\put(4861,-3121){\framebox(540,540){R$_2$}}
}%
+ %% {\color[rgb]{0,0,0}
+ %% \put(1621,-3121){\framebox(540,540){V}}
+ %% }%
{\color[rgb]{0,0,0}
- \put(1621,-3121){\framebox(540,540){V}}
- }%
- {\color[rgb]{0,0,0}
\put(1081,-1771){\vector( 1, 0){540}}
- \put(1256,-1681){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$(1-p)\,N\,\nu$}}}}
+ \put(1256,-1681){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$N\,\nu$}}}}
}%
+ %% {\color[rgb]{0,0,0}
+ %% \put(1081,-2851){\vector( 1, 0){540}}
+ %% \put(1306,-2761){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$p\,N\,\nu$}}}}
+ %% }%
{\color[rgb]{0,0,0}
- \put(1081,-2851){\vector( 1, 0){540}}
- \put(1306,-2761){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$p\,N\,\nu$}}}}
- }%
- {\color[rgb]{0,0,0}
\put(2161,-1771){\vector( 1, 0){540}}
\put(2381,-1906){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\lambda$}}}}
}%
@@ -180,10 +170,10 @@
\put(2971,-2761){\vector( 0, 1){720}}
\put(3916,-2716){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$(1-\xi)\,\lambda\,\epsilon$}}}}
}%
-%% {\color[rgb]{0,0,0}
-%% \put(4321,-1771){\vector( 1, 0){540}}
-%% \put(4541,-1906){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\gamma$}}}}
-%% }%
+ %% {\color[rgb]{0,0,0}
+ %% \put(4321,-1771){\vector( 1, 0){540}}
+ %% \put(4541,-1906){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\gamma$}}}}
+ %% }%
{\color[rgb]{0,0,0}
\put(4321,-1771){\vector( 1, 0){540}}
\put(4541,-1906){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$(1-\phi)\,\gamma$}}}}
@@ -198,25 +188,25 @@
\put(4996,-2041){\vector( 0,-1){540}}
\put(4821,-2311){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\alpha_1$}}}}
}%
+ %% {\color[rgb]{0,0,0}
+ %% \put(2161,-2941){\vector( 1, 0){2700}}
+ %% \put(3421,-3076){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\alpha_2$}}}}
+ %% }%
{\color[rgb]{0,0,0}
- \put(2161,-2941){\vector( 1, 0){2700}}
- \put(3421,-3076){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\alpha_2$}}}}
- }%
- {\color[rgb]{0,0,0}
\put(5401,-2851){\line( 1, 0){270}}
\put(5671,-2851){\line( 0, 1){1620}}
\put(5671,-1231){\line(-1, 0){3915}}
\put(1756,-1231){\vector( 0,-1){270}}
- \put(2971,-1186){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\alpha_3$}}}}
+ \put(2971,-1186){\makebox(0,0)[b]{\smash{\SetFigFont{6}{7.2}{$\alpha_2$}}}}
}%
\end{picture}%
}
\end{center}
\caption{
- Diagram of the full SVEIRR model, a simplified version of that studied by \citet{Wearing2009}.
- When $\alpha_{1}=\alpha_{2}=\alpha_{3}=\epsilon=\phi=0$, this reduces to the SEIR model.
- When $\phi=\epsilon=0$ and $\alpha_{1}{\gg}\alpha_{3}$, this reduces to the SEIRS model.
- When $\phi=\epsilon=0$ and $\alpha_{1}=\alpha_{3}$, we obtain an SEIRS system where the waning of immunity is Gamma-distributed with a shape parameter of 2.
+ Diagram of the full SEIRR model, a simplified version of that studied by \citet{Wearing2009}.
+ When $\alpha_1=\epsilon=\phi=0$, this reduces to the SEIR model.
+ When $\phi=\epsilon=0$ and $\alpha_{1}{\gg}\alpha_{2}$, this reduces to the SEIRS model.
+ When $\phi=\epsilon=0$ and $\alpha_{1}=\alpha_{2}$, we obtain an SEIRS system where the waning of immunity is Gamma-distributed with a shape parameter of 2.
}
\label{fig:sveirr}
\end{figure}
@@ -240,7 +230,7 @@
how can we assess the evidence for or against an important role for a particular mechanism?
In the present instance, we are concerned with two mechanisms that have been hypothesized to play an important role in pertussis transmission: waning of immunity and boosting of immunity by natural exposure.
We will assess the strength of evidence for each of these mechanisms using a model selection approach.
-Nested within the full SVEIRR model are several submodels:
+Nested within the full SEIRR model are several submodels:
\begin{enumerate}
\item When the intensity of the white noise $\zeta$ is zero, we lose the environmental stochasticity,
@@ -248,9 +238,9 @@
\item When $\alpha_{1}=\phi=0$, there is no loss of immunity and we have the SEIR system.
-\item When $\eps=\phi=0$, and both $\alpha_{1},\alpha_{3}>0$, we have loss of immunity but no boosting:
+\item When $\eps=\phi=0$, and both $\alpha_{1},\alpha_{2}>0$, we have loss of immunity but no boosting:
the system can be said to be of SEIRS form.
- If in addition $\alpha_{1}=\alpha_{3}$, the waiting time to loss of immunity is Gamma distributed with shape parameter 1.
+ If in addition $\alpha_{1}=\alpha_{2}$, the waiting time to loss of immunity is Gamma distributed with shape parameter 2.
\end{enumerate}
@@ -273,8 +263,7 @@
$\bar\beta$ & mean contact rate & $\bar{\beta}>0$ \\
$b_1$ & amplitude of seasonality & $0\;\leq\;b_1\;\leq\;1$\\
$1/ \alpha_1$ & mean period in R$_1$ & 1~da $\leq 1/\alpha_1\;\leq\;1/\mu$~yr\\
- $1/ \alpha_2$ & mean period in V & 1~da $\leq 1/\alpha_2\;\leq\;1/\mu$~yr\\
- $1/ \alpha_3$ & mean period in R$_2$ & 1~da $\leq 1/\alpha_3\;\leq\;1/\mu$~yr\\
+ $1/ \alpha_2$ & mean period in R$_2$ & 1~da $\leq 1/\alpha_2\;\leq\;1/\mu$~yr\\
$\xi$ & probability of immune boosting & $0\;\leq\;\xi\;\leq\;1$\\
$\phi$ & probability of not developing immunity & $0\;\leq\;\phi\;\leq\;1$\\
$p$ & vaccination probability & $0\;\leq\;p\;\leq\;1$\\
@@ -298,7 +287,7 @@
names(pertussis.sim)
@
There are three submodels (SEIR, SEIRS, SEIRR), representing a model with permanent immunity (SEIR), temporary immunity (SEIRS), and temporary immunity with boosting (SEIRR).
-The model labeled \code{full} is SVEIRR but with polarizing immunity ($\phi>0$) and environmental stochasticity $\sigma_{\zeta}>0$.
+The model labeled \code{full} is SEIRR but with polarizing immunity ($\phi>0$) and environmental stochasticity $\sigma_{\zeta}>0$.
For each of these, there are 2 simulations, one at each of two population sizes (\scinot{5}{5} and \scinot{5}{6}).
The parameter values used in the simulations are given in Table~\ref{tab:sim-params}.
The simulations are performed using an Euler-multinomial approximation to the continuous-time Markov chain implicitly defined by Fig.~\ref{fig:sveirr}.
@@ -316,24 +305,24 @@
"birthrate","deathrate",
"mean.beta","ampl.beta","imports",
"sigma","gamma",
- "alpha.1","alpha.2","alpha.3",
+ "alpha","alpha.ratio",
"report.prob","boost.prob","polar.prob",
- "vacc.prob","foi.mod",
+ "foi.mod",
"popsize","noise.sigma","tau",
- "S.0","E.0","I.0","R1.0","R2.0","V.0"
+ "S.0","E.0","I.0","R1.0","R2.0"
),
"math"] <- c(
"$\\nu$","$\\mu$","$\\hat\\beta$","$b_{1}$","$\\iota$",
"$\\sigma$","$\\gamma$",
- "$\\alpha_{1}$","$\\alpha_{2}$","$\\alpha_{3}$",
+ "$\\alpha_{1}$","$\\alpha_{2}/\\alpha_{1}$",
"$\\rho$","$\\xi$","$\\phi$",
- "$p$","$\\epsilon$",
+ "$\\epsilon$",
"$N$","$\\sigma_{\\zeta}$","$\\tau$",
- "$S(0)/N$","$E(0)/N$","$I(0)/N$","$R_1(0)/N$","$R_2(0)/N$","$V(0)/N$"
+ "$S(0)/N$","$E(0)/N$","$I(0)/N$","$R_1(0)/N$","$R_2(0)/N$"
)
params <- params[c("math","name","SEIR.small","SEIRS.big","SEIRR.small","full.big")]
-params <- subset(params,!name%in%c("popsize","vacc.prob","S.0","E.0","I.0","R1.0","R2.0","V.0"))
+params <- subset(params,!name%in%c("popsize","S.0","E.0","I.0","R1.0","R2.0"))
names(params) <- c("parameter","R name","SEIR","SEIRS","SEIRR","full")
@
@@ -547,8 +536,8 @@
We will suppose that no vaccination occurs ($p=0$) and that all infections result in immunity ($\phi=0$).
In addition, we will initially assume that the process noise is purely demographic: $\sigma_{\zeta}=0$.
Thus, for every model, we must fit the three parameters $\bar{\beta}$, $b_1$, and $\rho$.
-For the SEIR model, these are the only parameters that must be fit, since we assume $\alpha_{1}=\alpha_{3}=\eps=0$.
-For the SEIRS model, we must additionally fit $\alpha_{1}$ and $\alpha_{3}$, though we assume $\eps=0$.
+For the SEIR model, these are the only parameters that must be fit, since we assume $\alpha_{1}=\alpha_{2}=\eps=0$.
+For the SEIRS model, we must additionally fit $\alpha_{1}$ and $\alpha_{2}$, though we assume $\eps=0$.
Finally, in the SEIRR model, we must fit $\eps$ and $\xi$ in addition to all the others, for a total of seven parameters.
This is summarized in Table~\ref{tab:design}.
@@ -563,21 +552,21 @@
\hline\hline
Truth & Hypothesis & Parameters to estimate & Parameters fixed\\
\hline
- SEIRR.small & SEIR & $\bar{\beta}$, $b_1$, $\rho$ & $\alpha_{1}=\alpha_{3}=\eps=\phi=0$ \\
- SEIRR.small & SEIRS & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_3$ & $\eps=\phi=0$\\
- SEIRR.small & SEIRS (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$ & $\alpha_3=\alpha_1$, $\eps=\phi=0$\\
- SEIRR.small & SEIRR & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_3$, $\eps$, $\xi$ & $\phi=0$ \\
- SEIRR.small & SEIRR (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\eps$ & $\alpha_3=\alpha_1$, $\xi=\phi=0$ \\
- full.small & SEIR & $\bar{\beta}$, $b_1$, $\rho$ & $\alpha_{1}=\alpha_{3}=\eps=\phi=0$ \\
- full.small & SEIRS & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_3$ & $\eps=\phi=0$ \\
- full.small & SEIRS (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$ & $\alpha_3=\alpha_1$, $\eps=\phi=0$ \\
- full.small & SEIRR & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_3$, $\eps$, $\xi$ & $\phi=0$ \\
- full.small & SEIRR (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\eps$ & $\alpha_3=\alpha_1$, $\xi=\phi=0$ \\
+ SEIRR.small & SEIR & $\bar{\beta}$, $b_1$, $\rho$ & $\alpha_{1}=\alpha_{2}=\eps=\phi=0$ \\
+ SEIRR.small & SEIRS & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_2$ & $\eps=\phi=0$\\
+ SEIRR.small & SEIRS (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$ & $\alpha_2=\alpha_1$, $\eps=\phi=0$\\
+ SEIRR.small & SEIRR & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_2$, $\eps$, $\xi$ & $\phi=0$ \\
+ SEIRR.small & SEIRR (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\eps$ & $\alpha_2=\alpha_1$, $\xi=\phi=0$ \\
+ full.small & SEIR & $\bar{\beta}$, $b_1$, $\rho$ & $\alpha_{1}=\alpha_{2}=\eps=\phi=0$ \\
+ full.small & SEIRS & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_2$ & $\eps=\phi=0$ \\
+ full.small & SEIRS (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$ & $\alpha_2=\alpha_1$, $\eps=\phi=0$ \\
+ full.small & SEIRR & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\alpha_2$, $\eps$, $\xi$ & $\phi=0$ \\
+ full.small & SEIRR (alt) & $\bar{\beta}$, $b_1$, $\rho$, $\alpha_1$, $\eps$ & $\alpha_2=\alpha_1$, $\xi=\phi=0$ \\
\hline\hline
\end{tabular}
\end{table}
-In the event that it proves overly time-consuming to effectively fit the five parameters of the SEIRS model, we might allow ourselves to assume that $\alpha_{1}=\alpha_{3}$.
+In the event that it proves overly time-consuming to effectively fit the five parameters of the SEIRS model, we might allow ourselves to assume that $\alpha_{1}=\alpha_{2}$.
In any case, it is likely that the relative values of these two parameters will be poorly identified.
Likewise, if it proves really troublesome to fit six or seven parameters in the SEIRR model, we might allow ourselves to assume that $\xi=0$.
This would have the interesting effect of introducing model misspecification into the problem and give us the opportunity to explore and comment on that issue as well.
@@ -586,7 +575,7 @@
At the outset, we do not know where the MLEs for the SEIR and SEIRS models will be beyond the obvious facts that all the parameters are positive and $\rho$ and $\xi$ are probabilities.
Based on other research, we have some reason to expect to uncover identifiability issues between the reporting rate $\rho$ and the other parameters, especially $\bar{\beta}$.
In addition, waning of immunity acts to increase the overall mean prevalance and should thus be reflected in high $\bar{\beta}$ estimates for the SEIR model.
-As mentioned above, we expect identifiability problems to crop up between $\alpha_{1}$ and $\alpha_{3}$ as well.
+As mentioned above, we expect identifiability problems to crop up between $\alpha_{1}$ and $\alpha_{2}$ as well.
For these reasons, it may initially be best for optimization approaches to initialize optimizers over a wide range and Bayesian approaches to use very flat priors.
Once we can form a clearer notion of where the MLEs exist, it may be worthwhile to revisit this issue.
Modified: pkg/pompExamples/inst/doc/pertussis-model.pdf
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/src/pertussis.c
===================================================================
--- pkg/pompExamples/src/pertussis.c 2012-04-18 10:01:33 UTC (rev 667)
+++ pkg/pompExamples/src/pertussis.c 2012-04-18 11:50:55 UTC (rev 668)
@@ -5,14 +5,6 @@
#include <Rdefines.h>
#include <R_ext/Rdynload.h>
-static double expit (double x) {
- return 1.0/(1.0 + exp(-x));
-}
-
-static double logit (double x) {
- return log(x/(1-x));
-}
-
static inline double rgammawn (double sigma, double dt) {
double sigmasq;
sigmasq = sigma*sigma;
@@ -35,39 +27,35 @@
#define IMPORTS (p[parindex[4]]) // annual influx of infectives
#define SIGMA (p[parindex[5]]) // 1/LP
#define GAMMA (p[parindex[6]]) // 1/IP
-#define ALPHA1 (p[parindex[7]]) // 1/mean residence time in R1
-#define ALPHA2 (p[parindex[8]]) // 1/mean residence time in V
-#define ALPHA_RATIO (p[parindex[9]]) // 1/mean residence time in R2
-#define REPORT_PROB (p[parindex[10]]) // reporting probability
-#define BOOST_PROB (p[parindex[11]]) // probability of immune boosting
-#define POLAR_PHI (p[parindex[12]]) // probability of failing to develop immunity
-#define PVAC (p[parindex[13]]) // vaccine coverage (probability)
-#define FOI_MOD (p[parindex[14]]) // force of infection modulator
-#define NOISE_SIGMA (p[parindex[15]]) // white noise intensity
-#define POP (p[parindex[16]]) // population size
-#define TAU (p[parindex[17]]) // reporting SD
+#define ALPHA (p[parindex[7]]) // 1/mean residence time in R1
+#define ALPHA_RATIO (p[parindex[8]]) // 1/mean residence time in R2
+#define REPORT_PROB (p[parindex[9]]) // reporting probability
+#define BOOST_PROB (p[parindex[10]]) // probability of immune boosting
+#define POLAR_PHI (p[parindex[11]]) // probability of failing to develop immunity
+#define FOI_MOD (p[parindex[12]]) // force of infection modulator
+#define NOISE_SIGMA (p[parindex[13]]) // white noise intensity
+#define POP (p[parindex[14]]) // population size
+#define TAU (p[parindex[15]]) // reporting SD
#define SUS (x[stateindex[0]]) // susceptibles
#define EXP (x[stateindex[1]]) // exposed
#define INF (x[stateindex[2]]) // infectious
#define R1 (x[stateindex[3]]) // 1st recovered class
#define R2 (x[stateindex[4]]) // 1st recovered class
-#define VAC (x[stateindex[5]]) // vaccinated
-#define CASE (x[stateindex[6]]) // cumulative cases
-#define W (x[stateindex[7]]) // cumulative noise
-#define ERR (x[stateindex[8]]) // errors
-#define SIMPOP (x[stateindex[9]]) // simulated population size
+#define CASE (x[stateindex[5]]) // cumulative cases
+#define W (x[stateindex[6]]) // cumulative noise
+#define ERR (x[stateindex[7]]) // errors
+#define SIMPOP (x[stateindex[8]]) // simulated population size
#define DSDT (f[stateindex[0]]) // susceptibles
#define DEDT (f[stateindex[1]]) // exposed
#define DIDT (f[stateindex[2]]) // infectious
#define DR1DT (f[stateindex[3]]) // 1st recovered class
#define DR2DT (f[stateindex[4]]) // 1st recovered class
-#define DVDT (f[stateindex[5]]) // vaccinated
-#define DCASEDT (f[stateindex[6]]) // cumulative cases
-#define DWDT (f[stateindex[7]]) // cumulative noise
-#define DERRDT (f[stateindex[8]]) // errors
-#define DSIMPOPDT (f[stateindex[9]]) // simulated population size
+#define DCASEDT (f[stateindex[5]]) // cumulative cases
+#define DWDT (f[stateindex[6]]) // cumulative noise
+#define DERRDT (f[stateindex[7]]) // errors
+#define DSIMPOPDT (f[stateindex[8]]) // simulated population size
#define REPORTS (y[obsindex[0]])
@@ -75,7 +63,7 @@
int *stateindex, int *parindex, int *covindex,
int covdim, double *covar, double t, double dt)
{
- int ntrans = 17;
+ int ntrans = 15;
double rate[ntrans]; // transition rates
// int ntrans = sizeof(rate)/sizeof(rate[0]); // number of transitions
double trans[ntrans]; // transition numbers
@@ -86,50 +74,17 @@
reulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable("pomp","reulermultinom");
- // untransform the parameters
- alpha3 = ALPHA1*ALPHA_RATIO;
-
- if (!(R_FINITE(BIRTHRATE))) return;
- if (!(R_FINITE(DEATHRATE))) return;
- if (!(R_FINITE(MEANBETA))) return;
- if (!(R_FINITE(AMPLBETA))) return;
- if (!(R_FINITE(IMPORTS))) return;
- if (!(R_FINITE(SIGMA))) return;
- if (!(R_FINITE(GAMMA))) return;
- if (!(R_FINITE(ALPHA1))) return;
- if (!(R_FINITE(ALPHA2))) return;
- if (!(R_FINITE(alpha3))) return;
- if (!(R_FINITE(BOOST_PROB))) return;
- if (!(R_FINITE(POLAR_PHI))) return;
- if (!(R_FINITE(PVAC))) return;
- if (!(R_FINITE(FOI_MOD))) return;
- if (!(R_FINITE(NOISE_SIGMA))) return;
-
- if (!(R_FINITE(SUS))) return;
- if (!(R_FINITE(EXP))) return;
- if (!(R_FINITE(INF))) return;
- if (!(R_FINITE(R1))) return;
- if (!(R_FINITE(R2))) return;
- if (!(R_FINITE(VAC))) return;
- if (!(R_FINITE(CASE))) return;
- if (!(R_FINITE(W))) return;
-
+ alpha3 = ALPHA*ALPHA_RATIO;
beta = MEANBETA*(1+AMPLBETA*term_time(t));
dW = rgammawn(NOISE_SIGMA,dt);
- if (!(R_FINITE(dW))) {
- ERR += 1.0e18;
- return;
- }
-
foi = (IMPORTS+INF*beta*dW/dt)/POP; // force of infection
// compute the transition rates
- // rates into S & V
- rate[0] = (1-PVAC)*BIRTHRATE*POP; // birth into susceptible class
- rate[1] = PVAC*BIRTHRATE*POP; // birth into vaccinated class
+ // rates into S
+ rate[0] = BIRTHRATE*POP; // birth into susceptible class
// rates out of S
- rate[2] = foi; // force of infection
- rate[3] = DEATHRATE; // death from susceptible class
+ rate[2] = foi; // force of infection
+ rate[3] = DEATHRATE; // death from susceptible class
// rates out of E
rate[4] = SIGMA; // termination of latent period
rate[5] = DEATHRATE; // death from exposed class
@@ -138,45 +93,32 @@
rate[7] = (1-POLAR_PHI)*GAMMA; // recovery to immune class
rate[8] = DEATHRATE; // death from infectious class
// rates out of R1
- rate[9] = ALPHA1; // waning of immunity
+ rate[9] = ALPHA; // waning of immunity
rate[10] = DEATHRATE; // death from recovered class
// rates out of R2
rate[11] = BOOST_PROB*FOI_MOD*foi; // boosting of immunity
rate[12] = (1-BOOST_PROB)*FOI_MOD*foi; // infection
rate[13] = alpha3; // waning of R2 immunity
rate[14] = DEATHRATE; // death from R2
- // rates out of V
- rate[15] = ALPHA2; // waning of vaccine immunity
- rate[16] = DEATHRATE; // death from vaccinated class
// compute the transition numbers
trans[0] = rpois(rate[0]*dt); // births are Poisson
- trans[1] = rpois(rate[1]*dt);
reulermultinom(2,SUS, &rate[2],dt,&trans[2]);
reulermultinom(2,EXP, &rate[4],dt,&trans[4]);
reulermultinom(3,INF, &rate[6],dt,&trans[6]);
reulermultinom(2,R1, &rate[9],dt,&trans[9]);
reulermultinom(4,R2, &rate[11],dt,&trans[11]);
- reulermultinom(2,VAC,&rate[15],dt,&trans[15]);
// balance the equations
SUS += trans[0] + trans[6] + trans[13] - trans[2] - trans[3];
EXP += trans[2] + trans[12] - trans[4] - trans[5];
INF += trans[4] - trans[6] - trans[7] - trans[8];
R1 += trans[7] + trans[11] - trans[9] - trans[10];
- R2 += trans[9] + trans[15] - trans[11] - trans[12] - trans[13] - trans[14];
- VAC += trans[1] - trans[15] - trans[16];
+ R2 += trans[9] - trans[11] - trans[12] - trans[13] - trans[14];
CASE += trans[4];
W += (NOISE_SIGMA!=0) ? (dW-dt)/NOISE_SIGMA : 0.0; // mean zero, variance = dt
- SIMPOP = SUS + EXP + INF + R1 + R2 + VAC;
- // keep track of errors
- // if (SUS < 0.0) ERR += 1.0e0;
- // if (EXP < 0.0) ERR += 1.0e3;
- // if (INF < 0.0) ERR += 1.0e6;
- // if (R1 < 0.0) ERR += 1.0e9;
- // if (R2 < 0.0) ERR += 1.0e12;
- // if (VAC < 0.0) ERR += 1.0e15;
+ SIMPOP = SUS + EXP + INF + R1 + R2;
ERR += 1.0;
}
@@ -184,20 +126,19 @@
int *stateindex, int *parindex, int *covindex,
int covdim, double *covar, double t)
{
- int ntrans = 17;
+ int ntrans = 15;
double rate[ntrans]; // transition rates
double alpha3;
double beta;
double foi;
- alpha3 = ALPHA1*ALPHA_RATIO;
+ alpha3 = ALPHA*ALPHA_RATIO;
beta = MEANBETA*(1+AMPLBETA*term_time(t));
foi = (IMPORTS+INF*beta)/POP; // force of infection
// compute the transition rates
- // rates into S & V
- rate[0] = (1-PVAC)*BIRTHRATE*POP; // birth into susceptible class
- rate[1] = PVAC*BIRTHRATE*POP; // birth into vaccinated class
+ // rates into S
+ rate[0] = BIRTHRATE*POP; // birth into susceptible class
// rates out of S
rate[2] = foi*SUS; // force of infection
rate[3] = DEATHRATE*SUS; // death from susceptible class
@@ -205,33 +146,29 @@
rate[4] = SIGMA*EXP; // termination of latent period
rate[5] = DEATHRATE*EXP; // death from exposed class
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 668
More information about the pomp-commits
mailing list