diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2e6adc6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.Rhistory +.Rproj.user diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..e69de29 diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..e69de29 diff --git a/R/getFinalSize.R b/R/getFinalSize.R index 3d8e57d..b27884f 100644 --- a/R/getFinalSize.R +++ b/R/getFinalSize.R @@ -1,5 +1,6 @@ library(deSolve) +#' @export getFinalSizeAnalytic <- function(Rinit, Iinit, Vinit, N, R0, a, eps, q) { if (sum(Iinit) == 0) Iinit <- N / sum(N) @@ -24,6 +25,7 @@ getFinalSizeAnalytic <- function(Rinit, Iinit, Vinit, N, R0, a, eps, q) { opt$par + Iinit + Rinit } +#' @export exposure.SIR <- function(Time, state, Pars) { with(as.list(c(Time, state, Pars)), { @@ -48,6 +50,7 @@ exposure.SIR <- function(Time, state, Pars) { }) } +#' @export rescale.R0 <- function(beta, gam, pop.p, N, R0.value) { NGM.unscaled <- N * pop.p * beta * 1 / gam dom.eigen <- as.numeric(eigen(NGM.unscaled)$values[1]) @@ -57,6 +60,7 @@ rescale.R0 <- function(beta, gam, pop.p, N, R0.value) { return(scaling.factor) } +#' @export sim.exposure.SIR <- function(Rinit, Iinit, Vinit, tm, N, R0, gam, a, eps, q) { Ntot <- sum(N) @@ -96,7 +100,7 @@ sim.exposure.SIR <- function(Rinit, Iinit, Vinit, tm, N, R0, gam, a, eps, q) { simulation } - +#' @export getFinalSize <- function(vacTime, vacPortion, popSize, R0, recoveryRate, contactRatio, contactWithinGroup, suscRatio) { # vacTime: time after first case at which all vaccinations are delivered @@ -125,6 +129,7 @@ getFinalSize <- function(vacTime, vacPortion, popSize, R0, recoveryRate, } +#' @export finalSizeExample <- getFinalSize(vacTime = 100, vacPortion = c(0.1, 0.1), popSize = c(800000, 200000), R0 = 1.5, recoveryRate = 1 / 7, contactRatio = 1.7, contactWithinGroup = c(0.4, 0.4), suscRatio = 1) diff --git a/man/exposure.SIR.Rd b/man/exposure.SIR.Rd new file mode 100644 index 0000000..209b93f --- /dev/null +++ b/man/exposure.SIR.Rd @@ -0,0 +1,82 @@ +\name{exposure.SIR} +\alias{exposure.SIR} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ +A Capitalized Title (ideally limited to 65 characters) +} +\description{ +%% ~~ A concise (1-5 lines) description of what the function does. ~~ +} +\usage{ +exposure.SIR(Time, state, Pars) +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{Time}{ +%% ~~Describe \code{Time} here~~ +} + \item{state}{ +%% ~~Describe \code{state} here~~ +} + \item{Pars}{ +%% ~~Describe \code{Pars} here~~ +} +} +\details{ +%% ~~ If necessary, more details than the description above ~~ +} +\value{ +%% ~Describe the value returned +%% If it is a LIST, use +%% \item{comp1 }{Description of 'comp1'} +%% \item{comp2 }{Description of 'comp2'} +%% ... +} +\references{ +%% ~put references to the literature/web site here ~ +} +\author{ +%% ~~who you are~~ +} +\note{ +%% ~~further notes~~ +} + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +\seealso{ +%% ~~objects to See Also as \code{\link{help}}, ~~~ +} +\examples{ +##---- Should be DIRECTLY executable !! ---- +##-- ==> Define data, use random, +##-- or standard data sets, see data(). + +## The function is currently defined as +function (Time, state, Pars) +{ + with(as.list(c(Time, state, Pars)), { + S <- c(S1, S2) + I <- c(I1, I2) + R <- c(R1, R2) + N <- c(N1, N2) + beta <- (1 - epsilon) * outer(activities, activities)/sum(c(N1, + N2) * activities) + epsilon * activities/c(N1, N2) * + diag(2) + beta <- beta/scaling.factor + dS <- -(beta \%*\% I) * S + dI <- (beta \%*\% I) * S - gam * I + dR <- gam * I + return(list(c(dS, dI, dR))) + }) + } +} +% Add one or more standard keywords, see file 'KEYWORDS' in the +% R documentation directory (show via RShowDoc("KEYWORDS")): +% \keyword{ ~kwd1 } +% \keyword{ ~kwd2 } +% Use only one keyword per line. +% For non-standard keywords, use \concept instead of \keyword: +% \concept{ ~cpt1 } +% \concept{ ~cpt2 } +% Use only one concept per line. diff --git a/man/finalSizeExample.Rd b/man/finalSizeExample.Rd new file mode 100644 index 0000000..8b9b832 --- /dev/null +++ b/man/finalSizeExample.Rd @@ -0,0 +1,28 @@ +\name{finalSizeExample} +\alias{finalSizeExample} +\docType{data} +\title{ +A Capitalized Title for the Data Set +} +\description{ +%% ~~ A concise (1-5 lines) description of the dataset. ~~ +} +\usage{data("finalSizeExample")} +\format{ + The format is: + num [1:2] 233224 102557 +} +\details{ +%% ~~ If necessary, more details than the __description__ above ~~ +} +\source{ +%% ~~ reference to a publication or URL from which the data were obtained ~~ +} +\references{ +%% ~~ possibly secondary sources and usages ~~ +} +\examples{ +data(finalSizeExample) +## maybe str(finalSizeExample) ; plot(finalSizeExample) ... +} +\keyword{datasets} diff --git a/man/getFinalSize.Rd b/man/getFinalSize.Rd new file mode 100644 index 0000000..3c34554 --- /dev/null +++ b/man/getFinalSize.Rd @@ -0,0 +1,97 @@ +\name{getFinalSize} +\alias{getFinalSize} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ +A Capitalized Title (ideally limited to 65 characters) +} +\description{ +%% ~~ A concise (1-5 lines) description of what the function does. ~~ +} +\usage{ +getFinalSize(vacTime, vacPortion, popSize, R0, recoveryRate, contactRatio, contactWithinGroup, suscRatio) +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{vacTime}{ +%% ~~Describe \code{vacTime} here~~ +} + \item{vacPortion}{ +%% ~~Describe \code{vacPortion} here~~ +} + \item{popSize}{ +%% ~~Describe \code{popSize} here~~ +} + \item{R0}{ +%% ~~Describe \code{R0} here~~ +} + \item{recoveryRate}{ +%% ~~Describe \code{recoveryRate} here~~ +} + \item{contactRatio}{ +%% ~~Describe \code{contactRatio} here~~ +} + \item{contactWithinGroup}{ +%% ~~Describe \code{contactWithinGroup} here~~ +} + \item{suscRatio}{ +%% ~~Describe \code{suscRatio} here~~ +} +} +\details{ +%% ~~ If necessary, more details than the description above ~~ +} +\value{ +%% ~Describe the value returned +%% If it is a LIST, use +%% \item{comp1 }{Description of 'comp1'} +%% \item{comp2 }{Description of 'comp2'} +%% ... +} +\references{ +%% ~put references to the literature/web site here ~ +} +\author{ +%% ~~who you are~~ +} +\note{ +%% ~~further notes~~ +} + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +\seealso{ +%% ~~objects to See Also as \code{\link{help}}, ~~~ +} +\examples{ +##---- Should be DIRECTLY executable !! ---- +##-- ==> Define data, use random, +##-- or standard data sets, see data(). + +## The function is currently defined as +function (vacTime, vacPortion, popSize, R0, recoveryRate, contactRatio, + contactWithinGroup, suscRatio) +{ + Isim1 <- c(0, 0) + Rsim1 <- c(0, 0) + if (vacTime > 0) { + sim1 <- sim.exposure.SIR(Rinit = c(0, 0), Iinit = c(0, + 0), Vinit = c(0, 0), tm = vacTime, N = popSize, R0 = R0, + gam = recoveryRate, a = c(1, contactRatio), eps = contactWithinGroup, + q = c(1, suscRatio)) + Isim1 <- as.numeric(sim1[nrow(sim1), c("I1", "I2")]) + Rsim1 <- as.numeric(sim1[nrow(sim1), c("R1", "R2")]) + } + getFinalSizeAnalytic(Rinit = Rsim1, Iinit = Isim1, Vinit = popSize * + vacPortion, N = popSize, R0 = R0, a = c(1, contactRatio), + eps = contactWithinGroup, q = c(1, suscRatio)) + } +} +% Add one or more standard keywords, see file 'KEYWORDS' in the +% R documentation directory (show via RShowDoc("KEYWORDS")): +% \keyword{ ~kwd1 } +% \keyword{ ~kwd2 } +% Use only one keyword per line. +% For non-standard keywords, use \concept instead of \keyword: +% \concept{ ~cpt1 } +% \concept{ ~cpt2 } +% Use only one concept per line. diff --git a/man/getFinalSizeAnalytic.Rd b/man/getFinalSizeAnalytic.Rd new file mode 100644 index 0000000..1bdd802 --- /dev/null +++ b/man/getFinalSizeAnalytic.Rd @@ -0,0 +1,100 @@ +\name{getFinalSizeAnalytic} +\alias{getFinalSizeAnalytic} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ +A Capitalized Title (ideally limited to 65 characters) +} +\description{ +%% ~~ A concise (1-5 lines) description of what the function does. ~~ +} +\usage{ +getFinalSizeAnalytic(Rinit, Iinit, Vinit, N, R0, a, eps, q) +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{Rinit}{ +%% ~~Describe \code{Rinit} here~~ +} + \item{Iinit}{ +%% ~~Describe \code{Iinit} here~~ +} + \item{Vinit}{ +%% ~~Describe \code{Vinit} here~~ +} + \item{N}{ +%% ~~Describe \code{N} here~~ +} + \item{R0}{ +%% ~~Describe \code{R0} here~~ +} + \item{a}{ +%% ~~Describe \code{a} here~~ +} + \item{eps}{ +%% ~~Describe \code{eps} here~~ +} + \item{q}{ +%% ~~Describe \code{q} here~~ +} +} +\details{ +%% ~~ If necessary, more details than the description above ~~ +} +\value{ +%% ~Describe the value returned +%% If it is a LIST, use +%% \item{comp1 }{Description of 'comp1'} +%% \item{comp2 }{Description of 'comp2'} +%% ... +} +\references{ +%% ~put references to the literature/web site here ~ +} +\author{ +%% ~~who you are~~ +} +\note{ +%% ~~further notes~~ +} + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +\seealso{ +%% ~~objects to See Also as \code{\link{help}}, ~~~ +} +\examples{ +##---- Should be DIRECTLY executable !! ---- +##-- ==> Define data, use random, +##-- or standard data sets, see data(). + +## The function is currently defined as +function (Rinit, Iinit, Vinit, N, R0, a, eps, q) +{ + if (sum(Iinit) == 0) + Iinit <- N/sum(N) + Sinit <- N - Rinit - Iinit - Vinit + fn <- (1 - eps) * a * N + f <- fn/sum(fn) + cij <- diag(eps) + outer((1 - eps), f) + R0i <- R0/eigen(a * q * cij)$values[1] * a * q + Zrhs <- function(Z) c(Sinit * (1 - exp(-R0i * cij \%*\% ((Z + + Iinit)/N)))) + optfn <- function(x) ifelse(all(x > 0), max(abs(x - Zrhs(x))), + Inf) + optVal <- Inf + while (optVal > 1) { + opt <- optim((0.01 + 0.98 * runif(length(N))) * N, optfn) + optVal <- opt$value + } + opt$par + Iinit + Rinit + } +} +% Add one or more standard keywords, see file 'KEYWORDS' in the +% R documentation directory (show via RShowDoc("KEYWORDS")): +% \keyword{ ~kwd1 } +% \keyword{ ~kwd2 } +% Use only one keyword per line. +% For non-standard keywords, use \concept instead of \keyword: +% \concept{ ~cpt1 } +% \concept{ ~cpt2 } +% Use only one concept per line. diff --git a/man/rescale.R0.Rd b/man/rescale.R0.Rd new file mode 100644 index 0000000..096008d --- /dev/null +++ b/man/rescale.R0.Rd @@ -0,0 +1,79 @@ +\name{rescale.R0} +\alias{rescale.R0} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ +A Capitalized Title (ideally limited to 65 characters) +} +\description{ +%% ~~ A concise (1-5 lines) description of what the function does. ~~ +} +\usage{ +rescale.R0(beta, gam, pop.p, N, R0.value) +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{beta}{ +%% ~~Describe \code{beta} here~~ +} + \item{gam}{ +%% ~~Describe \code{gam} here~~ +} + \item{pop.p}{ +%% ~~Describe \code{pop.p} here~~ +} + \item{N}{ +%% ~~Describe \code{N} here~~ +} + \item{R0.value}{ +%% ~~Describe \code{R0.value} here~~ +} +} +\details{ +%% ~~ If necessary, more details than the description above ~~ +} +\value{ +%% ~Describe the value returned +%% If it is a LIST, use +%% \item{comp1 }{Description of 'comp1'} +%% \item{comp2 }{Description of 'comp2'} +%% ... +} +\references{ +%% ~put references to the literature/web site here ~ +} +\author{ +%% ~~who you are~~ +} +\note{ +%% ~~further notes~~ +} + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +\seealso{ +%% ~~objects to See Also as \code{\link{help}}, ~~~ +} +\examples{ +##---- Should be DIRECTLY executable !! ---- +##-- ==> Define data, use random, +##-- or standard data sets, see data(). + +## The function is currently defined as +function (beta, gam, pop.p, N, R0.value) +{ + NGM.unscaled <- N * pop.p * beta * 1/gam + dom.eigen <- as.numeric(eigen(NGM.unscaled)$values[1]) + NGM.scaled <- NGM.unscaled/dom.eigen * R0.value + scaling.factor <- beta/(NGM.scaled/(N * pop.p) * gam) + return(scaling.factor) + } +} +% Add one or more standard keywords, see file 'KEYWORDS' in the +% R documentation directory (show via RShowDoc("KEYWORDS")): +% \keyword{ ~kwd1 } +% \keyword{ ~kwd2 } +% Use only one keyword per line. +% For non-standard keywords, use \concept instead of \keyword: +% \concept{ ~cpt1 } +% \concept{ ~cpt2 } +% Use only one concept per line. diff --git a/man/sim.exposure.SIR.Rd b/man/sim.exposure.SIR.Rd new file mode 100644 index 0000000..b49e8a2 --- /dev/null +++ b/man/sim.exposure.SIR.Rd @@ -0,0 +1,109 @@ +\name{sim.exposure.SIR} +\alias{sim.exposure.SIR} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ +A Capitalized Title (ideally limited to 65 characters) +} +\description{ +%% ~~ A concise (1-5 lines) description of what the function does. ~~ +} +\usage{ +sim.exposure.SIR(Rinit, Iinit, Vinit, tm, N, R0, gam, a, eps, q) +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{Rinit}{ +%% ~~Describe \code{Rinit} here~~ +} + \item{Iinit}{ +%% ~~Describe \code{Iinit} here~~ +} + \item{Vinit}{ +%% ~~Describe \code{Vinit} here~~ +} + \item{tm}{ +%% ~~Describe \code{tm} here~~ +} + \item{N}{ +%% ~~Describe \code{N} here~~ +} + \item{R0}{ +%% ~~Describe \code{R0} here~~ +} + \item{gam}{ +%% ~~Describe \code{gam} here~~ +} + \item{a}{ +%% ~~Describe \code{a} here~~ +} + \item{eps}{ +%% ~~Describe \code{eps} here~~ +} + \item{q}{ +%% ~~Describe \code{q} here~~ +} +} +\details{ +%% ~~ If necessary, more details than the description above ~~ +} +\value{ +%% ~Describe the value returned +%% If it is a LIST, use +%% \item{comp1 }{Description of 'comp1'} +%% \item{comp2 }{Description of 'comp2'} +%% ... +} +\references{ +%% ~put references to the literature/web site here ~ +} +\author{ +%% ~~who you are~~ +} +\note{ +%% ~~further notes~~ +} + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +\seealso{ +%% ~~objects to See Also as \code{\link{help}}, ~~~ +} +\examples{ +##---- Should be DIRECTLY executable !! ---- +##-- ==> Define data, use random, +##-- or standard data sets, see data(). + +## The function is currently defined as +function (Rinit, Iinit, Vinit, tm, N, R0, gam, a, eps, q) +{ + Ntot <- sum(N) + beta <- (1 - eps) * outer(a, a)/sum(N * a) + eps * a/(N) * + diag(2) + scaling.factor <- rescale.R0(beta = beta, gam = gam, pop.p = N/Ntot, + N = Ntot, R0.value = R0) + pars <- list(N = Ntot, N1 = N[1], N2 = N[2], activities = a, + epsilon = eps, gam = gam, scaling.factor = scaling.factor) + times <- c(0, tm) + S.init <- c(S1 = N[1] - Rinit[1] - Iinit[1] - Vinit[1], S2 = N[2] - + Rinit[2] - Iinit[2] - Vinit[2]) + I.init <- c(I1 = Iinit[1], I2 = Iinit[2]) + R.init <- c(R1 = Rinit[1], R2 = Rinit[2]) + if (sum(I.init) == 0) { + I.init <- c(I1 = N[1]/Ntot, I2 = N[2]/Ntot) + S.init <- S.init - I.init + } + y0 <- c(S.init, I.init, R.init) + simulation <- as.data.frame(ode(y0, times, exposure.SIR, + pars)) + simulation + } +} +% Add one or more standard keywords, see file 'KEYWORDS' in the +% R documentation directory (show via RShowDoc("KEYWORDS")): +% \keyword{ ~kwd1 } +% \keyword{ ~kwd2 } +% Use only one keyword per line. +% For non-standard keywords, use \concept instead of \keyword: +% \concept{ ~cpt1 } +% \concept{ ~cpt2 } +% Use only one concept per line.