library(ramp.xds)Trace Functions
Basic Forcing by Trivial Modules
Overview
Malaria (and other mosquito-transmitted pathogens) can be understood as complex adaptive systems —non-linear interactions among humans, parasites, mosquitoes, and managers that are forced by weather, hydrology, and other factors. Dynamical systems have thus been developed to understand malaria epidemiology, malaria transmission dynamics, mosquito ecology, and malaria control. These models are naturally modular [1]. Given the complexity, it is sometimes useful to isolate a part of the system and reduce inputs from the rest of the system to a trace function.
We have been using the term trace function to describe any function that sets the value of the a dynamical term that has been passed from one of the dynamical components. This makes it possible to study sub-models, to study effects of exposure, forcing by malaria control, weather, or biotic factors.
Trivial Modules
A trace function sets the value of a dynamical term. The trace functions live in the trivial modules in ramp.xds. Five terms can be set up as trace functions:
\(\Lambda(t)\) or
Lambda— the emergence rate of adult, female mosquitoes, per patch is passed from the trivial L module.\(\eta(t)\) or
eta— eggs deposition rate, the number of eggs laid in aquatic habitats passed from the trivial MY module.\(fqZ(t)\) or
fqZ— the net rate of infectious biting by the adult mosquito population, per patch is passed from the trivial MY module.\(\kappa(t)\) or
kappa— the net infectiousness (NI), the probability that a mosquito will become infected after blood feeding on a human in that patch is passed from the trivial X module.\(E(t)\) or
eir— the daily entomological inoculation rate (EIR) for a cohort of humans (or the force of infection in a population or a cohort, FoI) is set up byxds_setup_eir. The EIR is passed to the XH module by a trace function from the EIR object, which is functionally equivalent to an MY module.
The first four terms are passed from trivial modules of the dynamical components. Examples are are documented in their vignettes. Since forcing by the EIR is handled slightly differently, it is the focus of Section 2: EIR forcing (below), which presents some code to illustrate how setup for trace functions works.
Composed Time Series
To build in some flexibility in a way that is easy for users to work with, all trace functions are called as composed time series functions \(x(t)\) that is the product of four elements:
\[ x(t) = \bar x \times F_S(t) \times F_T(t) \times F_K(t) \]
where
\(\bar x\) is a scaling parameter. It is close to the mean value of an unperturbed system. It is the mean over an interval \((t_0, t_1)\) only if \(F_K(t)=1\) and \[\int_{t_0}^{t_1} F_S(t)\; F_T(t)\; dt = t_1-t_0\]
\(F_S(t)>0\) is a seasonal pattern function, normalized over the interval such that \[\int_{t_0}^{t_1} F_S(t) \;dt = t_1 - t_0\]
\(F_T(t)>0\) is a trend pattern function, normalized over the interval such that \[\int_{t_0}^{t_1} F_T(t) \; dt = t_1 - t_0\]
\(F_K(t)>0\) is a shock function
Seasonal Patterns
The seasonal pattern function \(F_S(t,V)\) is always called
F_season.A library of seasonality pattern functions is maintained as part of
ramp.func
Trends
The trend pattern function \(F_T(t,V)\) is always called
F_trend.A library of trend pattern functions is maintained as part of
ramp.func
Shocks
The shock function \(F_K(t,V)\) is always called
F_shock.A library of shock functions is maintained as part of
ramp.func
EIR Forcing
To understand malaria epidemiology in relation to exposure, we can use ramp.xds::xds_setup_eir() to construct a trace function.
WRITE ME.
Two ways of building trace functions are:
User Defined Functions: Build a trace function and set it up:
Use
xds_setup_eir()Use
change_F_season()
Use
ramp.func::make_functionUse
ramp.library::change_F_s
User Defined Functions
Any function can be configured by a user used for a seasonal pattern, a trend, or a shock. Function has have the form \(F(t, V);\) the second variable, \(V,\) is part of the flexible and extensible design.
In some cases, users might want to set up an exogenous forcing variable and then call a trace functions that uses it, so all the components of the composed time series function have the form \(F(t,V).\) In the function calls that evaluate the trace function, the second variable, \(V,\) is retrieved by the function ramp.xds::get_variables so that trace functions can use any variable stored on the xds model object. By default \(V\) is an empty list, and get_variables returns an empty list
If a user defined function does not need \(V,\) the function must define the variable anyway and set its default value to an empty list, like this:
F_eg = function(t, V=list()){1 + sin(2*pi*t/365)}The function uses only \(t,\) ignoring the second argument (even if it’s Santa Claus).
t <- c(0:730)
plot(F_eg(t, "Santa Claus"), type = "l",
xlab = "Time (Days)",
ylab = expression(F[eg]))
Having defined F_eg, we can use it. In this case, we pass it as the F_season argument of xds_setup_eir. The function show_season plots it (by default, for two full cycles).
mod_eg <- xds_setup_eir(F_season = F_eg)
show_season(mod_eg)
change_F_*
change_F_season | change_F_trend| change_F_shock
If we didn’t pass F_season to xds_setup_eir, the default function is F_one, a function that returns the value \(1\).
mod_eg <- xds_setup_eir()
show_season(mod_eg)
We can modify the F_season argument using change_F_season
mod_eg <- change_F_season(F_eg, mod_eg)
show_season(mod_eg)
change_F_season
The function change_F_season dispatches on xds_obj$forced_by, which is set up by the trivial module.
In this case, the trivial L module is the default for xds_setup, and since no value for Lname was passed as an argument, mod_eg set up a trivial L module. No Loptions were passed at setup, so the function F_season was originally set to F_one.
Similarly, a user-defined trend or shock could be added using change_F_trend or change_F_shock and inspected using show_trend or show_shock
make_function
In the function library ramp.func, make_function returns a properly formatted trace function. To use make_function, the user passes a function object F_obj that was constructed by makepar_F_methodname where methodname is one of the built-in function families.
One advantage of these function families is that the functions include shape parameters, and the returned function is normalized (by default) so that it has an average value of \(1,\)
library(ramp.func)
par_spline <- makepar_F_spline(tt = c(0:5)*365, yy = c(.1, 1.4, 1, 1.2, 0.5, 1), X=2)
par_sin <- makepar_F_sin(pw=2)Now we make the functions and examine them and their product:
t5y <- seq(0, 365*5, by = 5)
F_spline <- make_function(par_spline)
F_sin <- make_function(par_sin)
plot(t5y, F_spline(t5y)*F_sin(t5y), type = "l", xlab = "Time (in Days)", ylab =expression(list(F[sin], F[spline])), lwd=2)
lines(t5y, F_sin(t5y), col = grey(0.5))
lines(t5y, F_spline(t5y), col = grey(0.5))
and we can set it up as before. This time, we set the scaling parameter eir to have an annual value of 1 infectious bite, per person, per year.
mod_both <- xds_setup_eir(eir=1/365, F_trend = F_spline, F_season=F_sin)
mod_both <- xds_solve(mod_both, times = t5y)
xds_plot_PR(mod_both)
change_*
change_season | change_trend| change_shock
The change_* functions in ramp.forcing work with F_obj, rather than with functions.
library(ramp.forcing)
par_shock <- makepar_F_sharkbite(D=365, L=365)
mod_new <- change_shock(par_shock, mod_both)
show_shock(mod_new)
mod_new <- xds_solve(mod_new, times = t5y)
xds_plot_PR(mod_both, clrs=grey(0.5))
xds_plot_PR(mod_new, add=T)
I/O
The xds object can be saved to a file using saveRDS and read from a file using readRDS. In some cases, the forcing functions F_season and F_trend and F_shock cause the size of the saved file to be very large.
To save space, saveXDS strips the functions before calling saveRDS.
Alternatively, the functions can be set up as function objects ramp.func. In this case, the function ramp.forcing::readXDS will read the object using readRDS and then rebuild the functions using ramp.forcing::rebuild_forcing_functions().