Dev

NAMESPACE

@@ -23,6 +23,7 @@ export(mean_phase_peak)
 export(plot_eirVpr)
 export(plot_eirpr)
 export(pr2Lambda)
+export(pr_ts2eir_history)
 export(ssMYZ)
 export(sse_season)
 export(update_by_ar)

R/index_of_dispersal.R

@@ -1,43 +1,40 @@
 
-#' Compute the Index of Dispersal for a PR orbit
-#' for the \eqn{i^{th}} element of
-#' eirpr$scaling.
+#' Compute the Index of Dispersal for a PR seasonal orbit
+#' for the \eqn{i^{th}} element of `eirpr$scaling`
 #'
+#' @param i the index of the orbit to be evaluated
 #' @param pars an **`xds`** object
 #'
 #' @export
 compute_IoD_pr = function(i, pars){
   pr <- pars$outputs$eirpr$scaling[[i]]$pr
-  var(pr)/mean(pr)
+  stats::var(pr)/mean(pr)
 }
 
 #' Compute the Index of Dispersal for a seasonal pattern
-#' S(t) over one year. This assumes that
+#' \eqn{S(t)} over one year. This assumes that
 #' \deqn{\int_0^{365} S(t)dt = 365}
+#'
 #' If in doubt, use [compute_IoD_F]
 #'
-#' @param pars an **`xds`** object
-#' @param clrs a [character] vector of colors
-#' @param llty a [list]
+#' @param S a function describing a seasonal pattern
 #'
 #' @export
 compute_IoD_S = function(S){
   S2 = function(t, mean){(1-S(t))^2}
-  return(integrate(S2, 0, 365, mean=mean)$val/365)
+  return(stats::integrate(S2, 0, 365, mean=mean)$val/365)
 }
 
 #' Compute the Index of Dispersal for a function
 #' F(t) over one year.
 #'
-#' @param pars an **`xds`** object
-#' @param clrs a [character] vector of colors
-#' @param llty a [list]
+#' @param F a function describing a seasonality function (possibly not normalized)
 #'
 #' @export
 compute_IoD_F = function(F){
-  mean = integrate(F, 0, 365)$val/365
+  mean = stats::integrate(F, 0, 365)$val/365
 
   F2 = function(t, mean){(mean-F(t))^2}
-  var = integrate(F2, 0, 365, mean=mean)$val/365
+  var = stats::integrate(F2, 0, 365, mean=mean)$val/365
   var/mean^2
 }

R/plot-terms.R

@@ -88,9 +88,9 @@ eirpr_seasonal_profile = function(ix, pars, clrs){
 #' Draw the orbit for the \eqn{i^{th}} element of
 #' eirpr$scaling.
 #'
+#' @param i the index of the orbit to plot
 #' @param pars an **`xds`** object
-#' @param clrs a [character] vector of colors
-#' @param llty a [list]
+#' @param clr a [character] vector of colors
 #'
 #' @export
 add_orbits = function(i, pars, clr){
@@ -103,9 +103,9 @@ add_orbits = function(i, pars, clr){
 #' eirpr$scaling, and add points at the
 #' minimum and maximum eir and pr
 #'
+#' @param i the index of the orbit to plot
 #' @param pars an **`xds`** object
-#' @param clrs a [character] vector of colors
-#' @param llty a [list]
+#' @param clr a [character] vector of colors
 #'
 #' @export
 add_orbits_px = function(i, pars, clr){

R/pr2eir_history.R

@@ -0,0 +1,38 @@
+
+
+#' @title Reconstruct a history of exposure from a PR time series
+#' @description Construct a function describing the EIR, including
+#' the mean EIR, the seasonal pattern, and the i
+#'
+#'
+#' For set of paired a time series \eqn{X,} compute the
+#' phase of a seasonal pattern for the EIR
+#' @note This utility relies on `xds_cohort`
+#' @param pr_ts the PR observed
+#' @param times the times of the observations
+#' @param model an `xds` model
+#' @return an `xds` object
+#' @export
+pr_ts2eir_history <- function(pr_ts, times, model){
+
+  # First pass: set the mean eir from the mean pr
+  model <- xde_scaling_eir(model, 25)
+  mean_pr <- mean(pr_ts)
+  model$EIRpar$eir <- xde_pr2eir(mean_pr, model, TRUE)
+
+  # First pass: fit the phase and amplitude
+  fit_phase_sin_season(pr_ts, times, model) -> phase0
+  Fs0 <- make_F_sin(phase=phase0)
+  Fs1 <- fit_amplitude_sin_season(Fs0, pr_ts, times, model)
+  model$EIRpar$F_season <- make_function(Fs1)
+
+  # First pass: fit the phase and amplitude
+  fit_phase_sin_season(pr_ts, times, model) -> phase0
+  Fs0 <- make_F_sin(phase=phase0)
+  Fs1 <- fit_amplitude_sin_season(Fs0, pr_ts, times, model)
+  model$EIRpar$F_season <- make_function(Fs1)
+
+
+
+  return(model)
+}

R/season.R

@@ -1,8 +1,8 @@
 #' @title Compute the phase of the peak
-#' @description For a time series \eqn{X,} compute the
+#' @description For a time series of paired values \eqn{c(t,X)} compute the
 #' phase, the time of the year when there is a peak
 #' @param t the times
-#' @param X the data
+#' @param X the variable
 #' @param window days around t
 #' @return a list with the mean peak and the values
 #' @export
@@ -17,6 +17,48 @@ mean_phase_peak = function(t, X, window=170){
   return(list(mean_peak = mean(peak), full = peak))
 }
 
+
+#' @title Fit the phase for a time series
+#' @description For a time series \eqn{X,} compute the
+#' phase of a seasonal pattern for the EIR
+#' @param data the PR observed
+#' @param times the times of the observations
+#' @param model an `xds` model
+#' @param Fp parameters for
+#' @return a list with the mean peak and the values
+#' @export
+fit_phase_sin_season <- function(data, times, model, Fp=NULL){
+  if(is.null(Fp)) Fp <- makepar_F_sin()
+  ph <- seq(0, 360, by = 45)
+  ss = ph*0
+  for(i in 1:length(ph)){
+    Fp$phase <- ph[i]
+    ss[i] <- sse_season(Fp, data, times, model)
+  }
+  ix = which.min(ss)
+  ff <- ss[ix]
+  best <- ph[ix]
+
+  ph1 <- seq(ph[ix]-40, ph[ix]+40, by=10)
+  ss1 = ph1*0
+  for(i in 1:length(ph1)){
+    Fp$phase <- ph1[i]
+    ss1[i] <- sse_season(Fp, data, times, model)
+  }
+  ix1 <- which.min(ss1)
+  ff1 <- ss1[ix1]
+  best1 <- ph1[ix1]
+
+  if(ff1<ff) best=best1
+  ph2 <- seq(best-4, best+4, by=1)
+  ss2 = ph2*0
+  for(i in 1:length(ph2)){
+    Fp$phase <- ph2[i]
+    ss2[i] <- sse_season(Fp, data, times, model)
+  }
+  return(ph2[which.min(ss2)])
+}
+
 #' @title Given data, compute GoF for a seasonal function
 #' @description For a time series c(`times`,`data`),
 #' compute the sum of squared errors for a seasonal
@@ -55,7 +97,7 @@ fit_amplitude_sin_season <- function(Fpar, data, times, model){
   Fpar$phase = Fpar$phase - shift
 
 
-  pars_seas <- optim(c(sqrt(Fpar$floor), log(Fpar$pw)), F_eval,
+  pars_seas <- stats::optim(c(sqrt(Fpar$floor), log(Fpar$pw)), F_eval,
     data = data, times = times, model = model, Fp=Fpar,
     method = "L-BFGS-B"
   )$par
@@ -90,7 +132,7 @@ fit_sin_season <- function(data, times, model){
   Fpar$phase = Fpar$phase - shift
 
 
-  pars_seas <- optim(c(sqrt(Fpar$floor), log(Fpar$pw), 0), F_eval,
+  pars_seas <- stats::optim(c(sqrt(Fpar$floor), log(Fpar$pw), 0), F_eval,
                      data = data, times = times, model = model, Fp=Fpar,
                      method = "L-BFGS-B"
   )$par
@@ -102,43 +144,3 @@ fit_sin_season <- function(data, times, model){
   return(Fpar)
 }
 
-#' @title Fit the phase for a time series
-#' @description For a time series \eqn{X,} compute the
-#' phase of a seasonal pattern for the EIR
-#' @param data the PR observed
-#' @param times the times of the observations
-#' @param model an `xds` model
-#' @param Fp parameters for
-#' @return a list with the mean peak and the values
-#' @export
-fit_phase_sin_season <- function(data, times, model, Fp=NULL){
-  if(is.null(Fp)) Fp <- makepar_F_sin()
-  ph <- seq(0, 360, by = 45)
-  ss = ph*0
-  for(i in 1:length(ph)){
-    Fp$phase <- ph[i]
-    ss[i] <- sse_season(Fp, data, times, model)
-  }
-  ix = which.min(ss)
-  ff <- ss[ix]
-  best <- ph[ix]
-
-  ph1 <- seq(ph[ix]-40, ph[ix]+40, by=10)
-  ss1 = ph1*0
-  for(i in 1:length(ph1)){
-    Fp$phase <- ph1[i]
-    ss1[i] <- sse_season(Fp, data, times, model)
-  }
-  ix1 <- which.min(ss1)
-  ff1 <- ss1[ix1]
-  best1 <- ph1[ix1]
-
-  if(ff1<ff) best=best1
-  ph2 <- seq(best-4, best+4, by=1)
-  ss2 = ph2*0
-  for(i in 1:length(ph2)){
-    Fp$phase <- ph2[i]
-    ss2[i] <- sse_season(Fp, data, times, model)
-  }
-  return(ph2[which.min(ss2)])
-}

_pkgdown.yml

@@ -63,6 +63,7 @@ reference:
   - xde_scaling
   - xde_scaling_eir
   - xde_scaling_lambda
+  - pr_ts2eir_history
 - title: Seasonality
   desc: |
     No methods to set subclass