[Gmm-commits] r168 - in pkg: . SEIR SEIR/R SEIR/man SEIR/src SEIR/vignettes

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Apr 14 05:41:35 CEST 2020


Author: chaussep
Date: 2020-04-14 05:41:34 +0200 (Tue, 14 Apr 2020)
New Revision: 168

Added:
   pkg/SEIR/
   pkg/SEIR/DESCRIPTION
   pkg/SEIR/NAMESPACE
   pkg/SEIR/R/
   pkg/SEIR/R/seir.R
   pkg/SEIR/man/
   pkg/SEIR/man/plot.Rd
   pkg/SEIR/man/print.Rd
   pkg/SEIR/man/seir.Rd
   pkg/SEIR/src/
   pkg/SEIR/src/Makevars
   pkg/SEIR/src/SEIR.h
   pkg/SEIR/src/solveseif.f
   pkg/SEIR/src/src.c
   pkg/SEIR/vignettes/
   pkg/SEIR/vignettes/seir.Rnw
   pkg/SEIR/vignettes/seir.bib
   pkg/SEIR/vignettes/seir.pdf
Log:
nothing to do with GMM, just a package to simulate the spread of a virus

Added: pkg/SEIR/DESCRIPTION
===================================================================
--- pkg/SEIR/DESCRIPTION	                        (rev 0)
+++ pkg/SEIR/DESCRIPTION	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,13 @@
+Package: SEIR
+Version: 0.1-0
+Date: 2020-04-13
+Title:  Susceptible, Exposed, Infectious, Removed Model
+Author: Pierre Chausse <pchausse at uwaterloo.ca>
+Maintainer: Pierre Chausse <pchausse at uwaterloo.ca>
+Description: This is a simple model for the spread of a virus. It is the solution to a system of differential equations with plot methods to visualize the solution. 
+Depends: R (>= 3.0.0)
+Imports: stats
+Suggests: knitr
+License: GPL (>= 2)
+NeedsCompilation: yes
+VignetteBuilder: knitr

Added: pkg/SEIR/NAMESPACE
===================================================================
--- pkg/SEIR/NAMESPACE	                        (rev 0)
+++ pkg/SEIR/NAMESPACE	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,15 @@
+useDynLib(SEIR, .registration = TRUE, .fixes="F_")
+
+importFrom("stats", "spline")
+importFrom("utils", "tail")
+importFrom("graphics", "plot", "lines", "legend")
+
+export(SEIR, solveSEIR, plot.seir, print.seir)
+
+S3method(plot, seir)
+S3method(print, seir)
+S3method("[", seir)
+
+
+
+

Added: pkg/SEIR/R/seir.R
===================================================================
--- pkg/SEIR/R/seir.R	                        (rev 0)
+++ pkg/SEIR/R/seir.R	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,124 @@
+SEIR <- function(t, y, y0, c=0, sigma=1/5.2, gamma=1/18,
+                 r0=matrix(c(0,20,70,84,90,3,2.6,1.9,1,.5),ncol=2),
+                 type=c("Lin", "Const"))    
+{
+    type <- match.arg(type)
+    type <- ifelse(type=="Lin", 0, 1)
+    par <- c(c, sum(y0), sigma, gamma)
+    res <- .Fortran(F_sysseir, as.double(t),
+                    as.double(y), as.double(par), as.integer(nrow(r0)),
+                    as.double(r0[,2]), as.double(r0[,1]), as.integer(type),
+                    RZero=double(1), dy=double(4))$dy
+}
+
+solveSEIR <- function(h=1e-2, T=180,
+                      c=0, y0=c(11e6, 40, 800, 0),
+                      sigma=1/5.2, gamma=1/18,
+                      r0=matrix(c(0,20,70,84,90,3,2.6,1.9,1,.5),ncol=2),
+                      type=c("Lin", "Const"))
+{
+    type <- match.arg(type)
+    typeF <- ifelse(type=="Lin", 0, 1)
+    n <- floor(T/h)
+    res <- .Fortran(F_solveseir, as.integer(n), as.double(0),
+                    as.double(T), as.double(c), as.double(sigma),
+                    as.double(gamma), as.double(y0),
+                    as.integer(nrow(r0)),
+                    as.double(r0[,2]), as.double(r0[,1]), as.integer(typeF),
+                    RZero=double(n+1), t=double(n+1), y=double((n+1)*4),
+                    h=double(1))
+    ans <- list(y=matrix(res$y,ncol=4), t=res$t, h=res$h,
+                RZero=res$RZero)
+    colnames(ans$y) = c("S","E","I","R")
+    class(ans) <- "seir"
+    ans$Pop <- sum(y0)
+    ans$h <= res$h
+    ans
+}
+
+plot.seir <- function(x, which=c("acase","tcase", "S", "E", "I", "R", "both","rzero"),
+                      add=FALSE, bothAdd=list(), ...)
+{
+    which <- match.arg(which)
+    x$y <- x$y/1e3
+    acase <- rowSums(x$y[,2:3])
+    tcase <- rowSums(x$y[,2:4])
+    if (which == "acase")
+    {
+        if (!add)
+            {
+                plot(x=x$t, y=acase, type="l", xlab="Days",
+                     ylab = "Thousands",
+                     lwd=2,
+                     main="SEIR Model: Number of Active Cases", ...)
+            } else {
+                lines(x=x$t, y=acase, ...)
+            }
+    } else if (which == "tcase") {
+        if (!add)
+        {
+            plot(x=x$t, y=tcase, type="l", xlab="Days",
+                 ylab = "Thousands",
+                 lwd=2, main="SEIR Model: Number of Cases", ...)
+        } else {
+            lines(x=x$t, y=tcase, ...)
+        }
+    } else if (which %in% c("S","E","I","R")) {
+        w <- which(c("S","E","I","R") == which)
+        wn <- c("Susceptibles", "Exposed", "Infectious", "Removed")[w]
+        main <- paste("SEIR Model: Number of ", wn, sep="")
+        y <- x$y[,w]
+        if (!add)
+        {
+            plot(x=x$t, y=y, type="l", xlab="Days",
+                 ylab =  "Thousands", lwd=2, main=main, ...)
+        } else {
+            lines(x=x$t, y=y, ...)
+        }        
+    } else if (which == "both") {
+        ylim <- range(c(tcase,acase))
+        if (!add)
+        {
+            plot(x=x$t, y=acase, type="l", xlab="Days",
+                 ylab =  "Thousands", ylim=ylim, lwd=2,
+                 main = "SEIR Model: Total and Active Cases", ...)
+            lines(x=x$t, y=tcase, col=2, lwd=2) 
+            legend("topleft", c("Active Cases", "Cumulative"),
+                   col=1:2, lty=1, lwd=2, bty='n')
+        } else {
+            lines(x=x$t, y=acase, ...)
+            bothAdd$x <- x$t
+            bothAdd$y <- acase
+            do.call(lines, bothAdd)
+        }
+    } else {
+        plot(x=x$t, y=x$RZero,
+             main="R-Zero function\n(number infected by one infectious)",
+             xlab="Days", ylab="infected/infectious", type='l', lwd=2, ...)
+    }
+    invisible()
+}
+
+print.seir <- function(x, digits=2, ...)
+{
+    cat("SEIR model for the spread of a virus\n")
+    cat("************************************\n")
+    cat("Initial condition\n")
+    cat("\tTotal Cases (thousands): ",
+        formatC(round(sum(x$y[1,2:4])/1000,digits), format='f'),"\n",sep="")
+    cat("\tActive Cases (thousands): ",
+        formatC(round(sum(x$y[1,2:3])/1000,digits), format='f'),"\n",sep="")
+    cat("Condition after ",
+        formatC(round(tail(x$t,1),0), format='d')," Days\n",sep="")
+    cat("\tTotal Cases (thousands): ",
+        formatC(round(sum(tail(x$y[,2:4],1))/1000,digits), format='f'),"\n",sep="")
+    cat("\tActive Cases (thousands): ",
+        formatC(round(sum(tail(x$y[,2:3],1))/1000,digits), format='f'),"\n",sep="")
+}
+
+"[.seir" <- function(x, i)
+{
+    ans <- sapply(1:4, function(j) spline(x$t,x$y[,j],xout=i)$y)
+    dimnames(ans) <- list(i, colnames(x$y))
+    ans
+}

Added: pkg/SEIR/man/plot.Rd
===================================================================
--- pkg/SEIR/man/plot.Rd	                        (rev 0)
+++ pkg/SEIR/man/plot.Rd	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,19 @@
+\name{plot}
+\alias{plot.seir}
+\title{The SEIR model plot method}
+\description{
+ It is a plot method for \code{seir} objects.
+}
+\usage{
+\method{plot}{seir}(x, which=c("acase","tcase", "S", "E", "I", "R",
+                    "both","rzero"), add=FALSE, bothAdd=list(), ...)
+}
+
+\arguments{
+  \item{x}{object of class \code{seir}.}
+  \item{which}{A choice of lines to plot.}
+  \item{add}{Should we add the lines to an existing plot?}
+  \item{bothAdd}{A list if argument to pass to the active case plot when
+  add is set to \code{TRUE}}
+  \item{...}{Other arguments to pass to the plot method.}
+}

Added: pkg/SEIR/man/print.Rd
===================================================================
--- pkg/SEIR/man/print.Rd	                        (rev 0)
+++ pkg/SEIR/man/print.Rd	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,15 @@
+\name{print}
+\alias{print.seir}
+\title{The SEIR model print method}
+\description{
+ It is the print method for \code{seir} objects.
+}
+\usage{
+\method{print}{seir}(x, digits=2, ...)
+}
+
+\arguments{
+  \item{x}{Object of class \code{seir}}
+  \item{digits}{Number of digits to print}
+  \item{...}{Arguments to pass to other methods}
+}

Added: pkg/SEIR/man/seir.Rd
===================================================================
--- pkg/SEIR/man/seir.Rd	                        (rev 0)
+++ pkg/SEIR/man/seir.Rd	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,55 @@
+\name{SEIR}
+\alias{SEIR}
+\alias{solveSEIR}
+
+\title{The SEIR model for virus spread}
+\description{
+The function that solves the SEIR system of differential equations.
+}
+\usage{
+
+SEIR(t, y, y0, c=0, sigma=1/5.2, gamma=1/18,
+     r0=matrix(c(0,20,70,84,90,3,2.6,1.9,1,.5),ncol=2),
+     type=c("Lin", "Const"))    
+
+solveSEIR(h=1e-2, T=180, c=0, y0=c(11e6, 40, 800, 0),
+          sigma=1/5.2, gamma=1/18,
+          r0=matrix(c(0,20,70,84,90,3,2.6,1.9,1,.5),ncol=2),
+          type=c("Lin", "Const"))
+}
+\arguments{
+
+  \item{t}{Time value to evaluate the derivatives.}
+  
+  \item{y}{The vector is variable of interest: S, E, I and R.}
+
+  \item{y0}{Initial Values}
+
+  \item{sigma}{Parameter.}
+
+  \item{gamma}{Parameter.}
+
+  \item{c}{Mutation Parameter}
+  
+  \item{h}{Step size for the RK4 method}
+  
+  \item{T}{Final t}
+
+  \item{r0}{Matrix of breaking points for R-Zero. The first column is
+  the time of the breaking points, and the second is the values of the
+  R-zero.}
+
+  \item{type}{Types of R-Zero function. The default is a linear
+    interpolation and the second is a piecewise constant function.}
+}
+
+\value{
+It returns an object of class \code{"seir"}  
+}
+\references{
+Wang, H., Wang, Z., Dong, Y. et al (2020),  Phase-adjusted estimation of
+the number of Coronavirus Disease 2019 cases in Wuhan, China.
+\emph{Cell Discov}, \bold{6}, \bold{10}.
+}
+
+

Added: pkg/SEIR/src/Makevars
===================================================================
--- pkg/SEIR/src/Makevars	                        (rev 0)
+++ pkg/SEIR/src/Makevars	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1 @@
+PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lgomp
\ No newline at end of file

Added: pkg/SEIR/src/SEIR.h
===================================================================
--- pkg/SEIR/src/SEIR.h	                        (rev 0)
+++ pkg/SEIR/src/SEIR.h	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,15 @@
+#ifndef R_SEIR_H
+#define R_SEIR_H
+
+#include <R_ext/RS.h>
+
+
+void F77_SUB(solveseir)(int *n, double *c, 
+			double *sigma, double *gamma, 
+			double *rzero, double *y0,
+			double *x, double *y, double *h);
+		 
+void F77_SUB(sysseir)(double *x, double *y, double *par,
+		   double *dy);
+
+#endif		     

Added: pkg/SEIR/src/solveseif.f
===================================================================
--- pkg/SEIR/src/solveseif.f	                        (rev 0)
+++ pkg/SEIR/src/solveseif.f	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,84 @@
+      subroutine rk4(f, parf, nparf, nb, n, k, a, b, y0,
+     *     rzerob, atx, type, rzero, x, y, h)
+      integer n, i, k, nparf, nb, type
+      double precision a, b, y0(k), h, x(n+1), y(n+1,k)
+      double precision f1(k), f2(k), f3(k), f4(k), parf(nparf)
+      double precision rzerob(nb), atx(nb), rzero(n+1), tmp
+      external f
+ 
+      h = (b-a)/n
+      y(1,:) = y0
+      x(1) = a
+      rzero(1) = rzerob(1)
+      do i=2,(n+1)
+         x(i) = a+h*(i-1)
+         call f(x(i-1), y(i-1,:),parf,nb,rzerob, atx, type,
+     *        tmp, f1) 
+         call f(x(i-1)+h/2,y(i-1,:)+h*f1/2, parf, nb, rzerob, atx, type,
+     *        tmp,f2)
+         call f(x(i-1)+h/2,y(i-1,:)+h*f2/2, parf, nb, rzerob, atx, type,
+     *        tmp, f3)
+         call f(x(i), y(i-1,:)+h*f3, parf, nb, rzerob, atx, type,
+     *        rzero(i), f4)
+         y(i,:) = y(i-1,:) + h*(f1+2*f2+2*f3+f4)/6
+      end do
+      end
+
+c     y = {S, E, I, R}
+c     par = {c, pop, sigma, gamma, nb, rzerob(nb), atx(nb)}    
+      subroutine sysseir(x, y, par, nb, rzerob, atx, type, rzero, dy)
+      integer nb, type
+      double precision par(4), x, y(4), dy(4)
+      double precision pop, sigma, gamma, rzero, beta
+      double precision rzerob(nb), atx(nb)
+      external R0
+      c = par(1)
+      pop = par(2)
+      sigma = par(3)
+      gamma = par(4)     
+      call R0(x,rzerob, atx, nb, type, rzero)
+      beta = rzero*gamma     
+      dy(1) = -beta*y(1)*y(3)/pop
+      dy(2) = beta*y(1)*y(3)/pop-sigma*y(2)
+      dy(3) = sigma*y(2)-gamma*y(3)+c*y(4)*y(3)/pop
+      dy(4) = gamma*y(3)-c*y(4)*y(3)/pop
+      end      
+
+      subroutine solveseir(n, a, b, c, sigma, gamma, y0,
+     *     nb, rzerob, atx, type, rzero, x, y, h)
+      integer n, nb, type
+      double precision c, sigma, gamma, y0(4), x(n+1)
+      double precision y(n+1,4), par(4), h, a, b
+      double precision rzerob(nb), atx(nb), rzero(n+1)
+      external sysseir
+      par(1) = c
+      par(2) = sum(y0)
+      par(3) = sigma
+      par(4) = gamma
+      call rk4(sysseir, par, 4, nb, n, 4, a, b, y0,
+     *     rzerob, atx, type, rzero, x, y, h)
+      end      
+      
+      subroutine R0(x, rzerob, atx, nb, type, v)
+      integer nb, i, type
+      double precision x, v, rzerob(nb), atx(nb), m
+
+      if (nb == 1) then
+         v = rzerob(1)
+      else
+         do i=2,nb
+            if (x < atx(i) .and. x >= atx(i-1)) then
+               if (type == 1) then
+                  v=rzerob(i-1)
+               else
+                  m = (rzerob(i)-rzerob(i-1))/(atx(i)-atx(i-1))
+                  v = rzerob(i-1)+m*(x-atx(i-1))
+               end if
+            end if
+         end do
+         if (x >= atx(nb)) then
+            v = rzerob(nb)
+         end if
+      end if
+      end
+    

Added: pkg/SEIR/src/src.c
===================================================================
--- pkg/SEIR/src/src.c	                        (rev 0)
+++ pkg/SEIR/src/src.c	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,18 @@
+#include <stdlib.h>
+#include <R_ext/Rdynload.h>
+#include "SEIR.h"
+
+static const R_FortranMethodDef fortranMethods[] = {
+  {"sysseir", (DL_FUNC) &F77_SUB(sysseir), 9},
+  {"solveseir", (DL_FUNC) &F77_SUB(solveseir), 15},
+  {NULL, NULL, 0}
+};
+
+void R_init_SEIR(DllInfo *dll)
+   {
+     R_registerRoutines(dll,
+			NULL, NULL, 
+			fortranMethods, NULL);
+     R_useDynamicSymbols(dll, FALSE);
+   }
+

Added: pkg/SEIR/vignettes/seir.Rnw
===================================================================
--- pkg/SEIR/vignettes/seir.Rnw	                        (rev 0)
+++ pkg/SEIR/vignettes/seir.Rnw	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,237 @@
+\documentclass[11pt,letterpaper]{article}
+\usepackage{amsthm}
+
+\usepackage[hmargin=2cm,vmargin=2.5cm]{geometry}
+\newtheorem{theorem}{Theorem}
+\newtheorem{col}{Corollary}
+\newtheorem{lem}{Lemma}
+\usepackage[utf8x]{inputenc}
+\newtheorem{ass}{Assumption}
+\usepackage{amsmath}
+\usepackage{verbatim}
+\usepackage[round]{natbib}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage{hyperref}
+\hypersetup{
+  colorlinks,
+  citecolor=black,
+  filecolor=black,
+  linkcolor=black,
+  urlcolor=black
+}
+
+\bibliographystyle{plainnat}
+
+
+\author{Pierre Chauss\'e}
+\title{The Susceptible, Exposed, Infectious, Removed (SEIR) Model}
+\begin{document}
+
+\maketitle
+
+\newcommand{\E}{\mathrm{E}}
+\newcommand{\diag}{\mathrm{diag}}
+\newcommand{\Prob}{\mathrm{Pr}}
+\newcommand{\Var}{\mathrm{Var}}
+\newcommand{\Vect}{\mathrm{Vec}}
+\newcommand{\Cov}{\mathrm{Cov}}
+\newcommand{\conP}{\overset{p}{\to}}
+\newcommand{\conD}{\overset{d}{\to}}
+\newcommand\R{ \mathbb{R} }
+\newcommand\N{ \mathbb{N} }
+\newcommand\C{ \mathbb{C} }
+\newcommand\rv{{\cal R}}
+\newcommand\Q{\mathbb{Q}}
+\newcommand\PR{{\cal R}}
+\newcommand\T{{\cal T}}
+\newcommand\Hi{{\cal H}}
+\newcommand\La{{\cal L}}
+\newcommand\plim{plim}
+\renewcommand{\epsilon}{\varepsilon}
+
+\abstract{This vignette presents how to use the SEIR package}
+%\VignetteIndexEntry{SEIR Model}
+%\VignetteKeywords{virus, spread}
+%\VignettePackage{SEIR}
+%\VignetteEngine{knitr::knitr}
+<<echo=FALSE>>=
+library(knitr)
+opts_chunk$set(size='footnotesize')
+@ 
+
+\section{The SEIR Model}\label{sec:single}
+
+This package solves the SEIR model for the spread of a virus. In the
+example, the parameters of \cite{wang-all20} are used for the case of
+Wuhan, China. The model is available in Matlab on Peter Forsyth
+webpage
+(\href{https://cs.uwaterloo.ca/~paforsyt/SEIR.html}{https://cs.uwaterloo.ca/~paforsyt/SEIR.html}). In
+fact, the R package is based on the work of Peter. You are invited to
+read his page for more details about the model. The system of
+differencial equations is:
+\begin{eqnarray*}
+\dot{S}(t) &=& -\beta(t) \frac{S(t)I(t)}{N}\\
+\dot{E}(t) &=& \beta(t) \frac{S(t)I(t)}{N} - \sigma E(t)\\
+\dot{I}(t) &=& \sigma E(t) - \gamma I(t) +c\frac{R(t)I(t)}{N}\\
+\dot{R}(t) &=& \gamma I(t)- c \frac{R(t)I(t)}{N}
+\end{eqnarray*}
+where $\beta(t)=\gamma RZero(t)$ and $RZero(t)$ is the number of
+individuals infected by an infectious person. It is a function of time
+because the transmission if function of the social distancing policy
+in place.
+
+\section{The Package}
+
+The main function is solveSEIR. It creates an object of class ``seir''. The main arguments are (See Peter's Matlab code for the source):
+\begin{itemize}
+\item h: The time interval used in the Runge-Kutta method.
+\item T: The number of days to solve.
+\item c: Mutation parameter (0 when no mutation is present)  
+\item sigma: Set to 1/5.2 for the Wuhan case.
+\item gamma: Set to 1/18 for the Wuhan case.
+\item y0: The initial value  $\{S(0), E(0), I(0), R(0)\}$.
+\item r0: The matrix of breaks for RZero(t).  See below.
+\item type: Type of RZero(t). See below.
+\end{itemize}
+
+In this document, we use the default parameter used by Peter
+Forsyth. The important factor is $RZero(t)$, which is determined by
+the arguments tyee an r0. The default type is ``LIN'' and the default r0 is:
+
+<<>>=
+matrix(c(0, 20, 70, 
+    84, 90, 3, 2.6, 1.9, 1, 0.5), ncol = 2)
+@ 
+
+It means that the value of RZero(t) changes at t=20, 70, 84 and 90
+(days). The type argument determines how it changes. By default it
+looks like the following:
+\begin{center}
+  \begin{minipage}{.7\textwidth}
+<<echo=FALSE, fig.height=4.5>>=
+library(SEIR) # we first load the package
+s <- solveSEIR()
+plot(s, "rzero")
+@     
+  \end{minipage}
+\end{center}
+Therefore, $RZero(t)$ is a linear interpolation between the different
+values. The other option is to have a constant $RZero(t)$ between
+periods as in the following graph.
+\begin{center}
+  \begin{minipage}{.7\textwidth}
+<<echo=FALSE, fig.height=4.5>>=
+s <- solveSEIR( type="C")
+plot(s, "rzero")
+@     
+\end{minipage}
+\end{center}
+The first $RZero(t)$ is the default option. If we run the function
+with all default values we obtain (the print method summarizes the result):
+
+<<>>=
+Sol <- solveSEIR()
+Sol
+@ 
+
+The object contains the following elements:
+
+<<>>=
+names(Sol)
+@ 
+
+The element y is an $nT\times4$ matrix, where $nT$ is the number of
+grid points (the closest integer of T/h), with the $i^{th}$ row being
+the solution $\{S(t_i), E(t_i), I(t_i), R(t_i)\}$. The element t is
+the time grid and $RZero$ is an $nT$ vector with the $i^{th}$ element
+being $RZero(t_i)$. The "["  method can be used to get the solution
+  at different points in time:
+<<>>=
+Sol[c(10,50,100,150)]  
+@   
+
+It is possible to use those elements to plot the solution, but the
+easiest ways is to use the plot method. The plot method gives you
+options to plot any of the curve:
+
+\begin{center}
+  \begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "acase")
+@     
+\end{minipage}
+\begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "tcase")
+@     
+\end{minipage}
+\end{center}
+
+\begin{center}
+\begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "both")
+@     
+\end{minipage}
+\begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "rzero")
+@     
+\end{minipage}
+\end{center}
+
+You can also plot them one by one:
+
+\begin{center}
+\begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "E")
+@     
+\end{minipage}
+\begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "I")
+@     
+\end{minipage}
+\end{center}
+
+As a last example, we can reproduce the second wave case produced by
+Peter Forsyth on his website. Suppose RZero stays at 0.5 for 30 days
+and then the government start relaxing the social distancing
+restrictions. The following is what he proposes:
+
+\begin{center}
+\begin{minipage}{.7\textwidth}
+<<fig.height=5>>=
+zbreaks <- matrix(c(0,  20,  70,  84,  90, 120, 240, 360,
+                    3, 2.6,  1.9, 1,  0.5, 0.5, 1.75, 0.5), ncol = 2)
+Sol <- solveSEIR(T=450, r0=zbreaks)
+plot(Sol, "rzero")
+@ 
+\end{minipage}
+\end{center}
+
+The new solution is:
+\begin{center}
+  \begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "acase")
+@     
+\end{minipage}
+\begin{minipage}{.4\textwidth}
+<<fig.height=5>>=
+plot(Sol, "both")
+@     
+\end{minipage}
+\end{center}
+
+
+
+
+
+
+\bibliography{seir}
+\end{document} 
+

Added: pkg/SEIR/vignettes/seir.bib
===================================================================
--- pkg/SEIR/vignettes/seir.bib	                        (rev 0)
+++ pkg/SEIR/vignettes/seir.bib	2020-04-14 03:41:34 UTC (rev 168)
@@ -0,0 +1,9 @@
+ at article{wang-all20,
+  AUTHOR={Wang, H. and Wang, Z. and Dong, Y and Chang, R. and Xu, C. and Yu, X. and Zhang, S. and Tsamlag, L. and Shang, M. and Huang, J. and Wang, Y. and Xu, G. and Shen, T. and Zhang, X.},
+  TITLE={Phase-adjusted estimation of the number of Coronavirus Disease 2019 cases in Wuhan, China},
+  JOURNAL={Cell Discov},
+  VOLUME={6},
+  YEAR={2020},
+  Number={10},
+  url={https://doi.org/10.1038/s41421-020-0148-0}
+ } 

Added: pkg/SEIR/vignettes/seir.pdf
===================================================================
(Binary files differ)

Index: pkg/SEIR/vignettes/seir.pdf
===================================================================
--- pkg/SEIR/vignettes/seir.pdf	2020-02-21 14:08:59 UTC (rev 167)
+++ pkg/SEIR/vignettes/seir.pdf	2020-04-14 03:41:34 UTC (rev 168)

Property changes on: pkg/SEIR/vignettes/seir.pdf
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property


More information about the Gmm-commits mailing list