Merge pull request #21 from dd-harp/dev

updated scaling; demography

.gitignore

@@ -10,3 +10,4 @@ docs/*.html
 *.html
 /docs
 docs/reference/ar_compare.html
+*.html

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(F_sse)
+export(UnevenAgingMatrix)
 export(ar_compare)
 export(dts_compute_gof)
 export(dts_init_by_pr)

R/demography.R

@@ -0,0 +1,15 @@
+#' @title Create a cohort matrix
+#' @description Return a matrix `dH`
+#' @param deathrates the death rates for
+#' @param ageMesh a partition on human age
+#' @return the steady states as a named vector
+#' @export
+UnevenAgingMatrix = function(deathrates, ageMesh){
+  N = length(ageMesh)
+  checkIt(deathrates, N+1)
+  ageUp = c(1/ageMesh,0)
+  Aging = diag(-ageUp-deathrates)
+  Aging[cbind(1+1:N, 1:N)] = 1/ageMesh
+  return(Aging)
+}
+

R/scaling.R

@@ -17,10 +17,11 @@ xde_scaling_eir = function(model, N=25){
   for(i in 1:N){
     model$EIRpar$eir <- aEIR[i]/365
     model <- ramp.xds::xds_solve_cohort(model, A=10, da=1)
-    out <- model$outputs$cohort
-    pr_t = tail(out$true_pr, 365); pr[i] = mean(pr_t)
-    ni_t = tail(out$ni, 365);  ni[i]= mean(ni_t)
-    eir_t = tail(out$eir, 365); eir[i] = mean(eir_t)
+    XH <- model$outputs$orbits$XH[[1]]
+    terms <- model$outputs$orbits$terms[[1]]
+    pr_t = tail(XH$true_pr, 365); pr[i] = mean(pr_t)
+    ni_t = tail(XH$ni, 365);  ni[i]= mean(ni_t)
+    eir_t = tail(terms$EIR, 365); eir[i] = mean(eir_t)
     scaling[[i]] = list(aeir = eir_t*365, eir = eir_t, pr = pr_t, ni = ni_t)
   }
 
@@ -61,7 +62,8 @@ xde_scaling_Z = function(model, N=25){
   for(i in 1:N){
     model$MYZpar$aeir = eir[i]
     xde_tmp <- ramp.xds::xde_stable_orbit(model)
-    tmp <- xde_tmp$outputs$stable_orbit
+    browser()
+    tmp <- xde_tmp$outputs$orbits
     H = tmp$XH[[1]]$H
     ni = tmp$terms$ni
     pr = tmp$terms$pr

vignettes/Scaling.R

@@ -5,14 +5,16 @@ library(rootSolve)
 library(ramp.work)
 
 ## ----echo=F-------------------------------------------------------------------
-#devtools::load_all()
+devtools::load_all()
 
-## -----------------------------------------------------------------------------
-F_sin = function(t){(1.01 + sin(2*pi*t/365))}
-F_1 = function(t){0*t + 1/365}
+## ----Fsin---------------------------------------------------------------------
+tt <- seq(0, 730, by=5) 
+p1 <- makepar_F_sin(floor=0.1)
+Fsin <- make_function(p1)
+plot(tt, Fsin(tt), type="l")
 
 ## -----------------------------------------------------------------------------
-xds_setup_cohort(Xname = "SIS", F_season=F_1) -> sis
+xds_setup_cohort(Xname = "SIS", F_season=Fsin) -> sis
 
 ## -----------------------------------------------------------------------------
 xds_solve_cohort(sis) -> sis
@@ -50,7 +52,7 @@ with(sis$outputs$eirpr, points(aeir, pr, pch = 15))
 with(preir_i, points(365*eir, pr, pch = 19, col = "red"))
 
 ## -----------------------------------------------------------------------------
-sis0 <- xds_setup_cohort(Xname = "SIS", F_season = F_sin)
+sis0 <- xds_setup_cohort(Xname = "SIS", F_season = Fsin)
 xde_scaling_eir(sis0, 25) -> sis0
 
 ## -----------------------------------------------------------------------------
@@ -66,7 +68,7 @@ with(sis0$outputs$eirpr, lines(scaling[[15]]$aeir, scaling[[15]]$pr, col = clrs[
 with(sis0$outputs$eirpr, lines(scaling[[20]]$aeir, scaling[[20]]$pr, col = clrs[20]))
 
 ## ----eval=F-------------------------------------------------------------------
-#  sip = xds_setup_cohort(Xname = "SIP", F_season=F_sin)
+#  sip = xds_setup_cohort(Xname = "SIP", F_season=Fsin)
 #  sip$Xpar[[1]]$eta = 1/40
 #  xde_scaling_eir(sip, 25) -> sip
 

vignettes/Scaling.Rmd

@@ -15,43 +15,48 @@ Load the required packages:
 
 ```{r}
 library(ramp.xds)
-library(deSolve)
-library(rootSolve)
 library(ramp.work)
+library(deSolve)
 ```
 
 ```{r, echo=F}
 #devtools::load_all()
 ```
+An important question in malaria is the relationship between various metrics, especially the relationship between the annual entomological inoculation rate (aEIR) and the parasite rate (PR). To facilitate analysis of malaria data, we have developed some functions to compute scaling relationships.  
 
-# Scaling
 
-This is a dog.
 
 ## xde_scaling 
 
-The function `xde_scaling()` defines the relationship between the EIR and the PR, and it outputs stable orbits for each value of aEIR in a mesh from $10^{-1}$ up to $10^{3}$ The code is in `mob_library/Work`
-Do we need two versions?
+The function `xde_scaling()` defines the relationship between the EIR and the PR. It analyzes stable orbits and outputs the average annual EIR and average annual PR for an even mesh of `log(aEIR)` values running from $10^{-1}$ up to $10^{3}$ The code is in `mob_library/Work`
 
-The cohort trace functions in `ramp.xds` take the form of `F(a, bday=0, scale=1).` 
+To illustrate, we pick a function describing a seasonal pattern using `ramp.xds::make_F_sin`
 
-```{r}
-F_sin = function(t){(1.01 + sin(2*pi*t/365))}
-F_1 = function(t){0*t + 1/365}
+```{r Fsin}
+tt <- seq(0, 730, by=5) 
+p1 <- makepar_F_sin(floor=0.1)
+Fsin <- make_function(p1)
+plot(tt, Fsin(tt), type="l")
 ```
 
+Next, we set up a cohort model: 
+
 ```{r}
-xds_setup_cohort(Xname = "SIS", F_season=F_1) -> sis
+xds_setup_cohort(Xname = "SIS", F_season=Fsin) -> sis
 ```
 
 ```{r}
 xds_solve_cohort(sis) -> sis
 ```
 
+The function `xde_scaling_eir` runs the model over a mesh of `N=25` values: 
+
 ```{r}
 xde_scaling_eir(sis, 25) -> sis
 ```
 
+The results are attached as `sis$outputs$eirpr` 
+
 ```{r}
 plot_eirpr(sis)
 ```
@@ -99,76 +104,3 @@ with(preir_i, points(365*eir, pr, pch = 19, col = "red"))
 ```
 
 
-## Scaling {.tabset}
-
-### Seasonality
-
-```{r}
-sis0 <- xds_setup_cohort(Xname = "SIS", F_season = F_sin)
-xde_scaling_eir(sis0, 25) -> sis0
-```
-
-
-```{r}
-clrs = turbo(25)
-with(sis$outputs$eirpr, plot(aeir, pr, type = "l", log = "x", xaxt= "n", xlab = "aEIR", ylab = "PR"))
-axis(1, 10^(-1:3), c(0.1, 1, 10, 100, 1000))
-lines(sis$outputs$eirpr$aeir, sis0$outputs$eirpr$pr, col = "tomato", lwd=2) 
-
-with(sis0$outputs$eirpr, points(aeir, pr, col = clrs))
-with(sis0$outputs$eirpr, lines(scaling[[5]]$aeir, scaling[[5]]$pr, col = clrs[5]))
-with(sis0$outputs$eirpr, lines(scaling[[10]]$aeir, scaling[[10]]$pr, col = clrs[10]))
-with(sis0$outputs$eirpr, lines(scaling[[15]]$aeir, scaling[[15]]$pr, col = clrs[15]))
-with(sis0$outputs$eirpr, lines(scaling[[20]]$aeir, scaling[[20]]$pr, col = clrs[20]))
-```
-
-### Drug taking 
-
-```{r, eval=F}
-sip = xds_setup_cohort(Xname = "SIP", F_season=F_sin)
-sip$Xpar[[1]]$eta = 1/40
-xde_scaling_eir(sip, 25) -> sip
-```
-
-
-```{r, eval=F}
-sip1 = setup_exposure_nb(sip, 1/50)
-xde_scaling_eir(sip1, 25) -> sip1
-```
-
-```{r, eval=F}
-with(sis$outputs$eirpr, plot(aeir, pr, type = "l", log = "x", xaxt= "n", xlab = "aEIR", ylab = "PR"))
-axis(1, 10^(-1:3), c(0.1, 1, 10, 100, 1000))
-with(sip$outputs$eirpr, lines(aeir, pr, col = "darkorange"))
-with(sip1$outputs$eirpr, lines(aeir, pr, col = "brown"))
-``` 
-
-### Environmental Heterogeneity 
-
-```{r}
-sis4 <- setup_exposure_nb(sis, 1/50)
-xde_scaling_eir(sis4, 25) -> sis4
-```
-
-```{r}
-with(sis$outputs$eirpr, plot(aeir, pr, type = "l", log = "x", xaxt= "n", xlab = "aEIR", ylab = "PR"))
-axis(1, 10^(-1:3), c(0.1, 1, 10, 100, 1000))
-#with(sis2$outputs$eir, lines(aeir, pr, col = "blue"))
-#with(sis3$outputs$eir, lines(aeir, pr, col = "purple"))
-with(sis4$outputs$eir, lines(aeir, pr, col = "darkblue"))
-```
-
-### Travel
-
-```{r}
-sis5 <- setup_travel_static(sis, delta = 1/5/365)
-xde_scaling_eir(sis5, 25) -> sis5
-```
-
-
-```{r}
-with(sis$outputs$eirpr, plot(aeir, pr, type = "l", log = "x", xaxt= "n", xlab = "aEIR", ylab = "PR"))
-axis(1, 10^(-1:3), c(0.1, 1, 10, 100, 1000))
-with(sis5$outputs$eir, lines(aeir, pr, col = "darkgreen")) 
-```
-