Add BDD(r) sampling to LBDP model

Author: Horopter

R/lbdp.R

@@ -10,6 +10,7 @@
 ##' @param lambda per capita birth rate
 ##' @param mu per capita death rate
 ##' @param psi per capita sampling rate
+##' @param r probability that a sampled lineage is removed (must be between 0 and 1; r=0 is non-destructive; r=1 is destructive)
 ##' @param n0 population size at time t0
 ##' @param time final time
 ##' @param t0 initial time
@@ -27,13 +28,15 @@ NULL
 ##' @export
 runLBDP <- function (
   time,  t0 = 0,
-  lambda = 2, mu = 1, psi = 1,
+  lambda = 2, mu = 1, psi = 1, r = 0,
   n0 = 5
 ) {
+  if (!is.numeric(r) || length(r) != 1L || !is.finite(r) || r < 0 || r > 1)
+    pStop(sQuote("r")," must be between 0 and 1.")
   n0 <- round(n0)
   if (n0 < 0)
     pStop(sQuote("n0")," must be a nonnegative integer.")
-  params <- c(lambda=lambda,mu=mu,psi=psi)
+  params <- c(lambda=lambda,mu=mu,psi=psi,r=r)
   ivps <- c(n0=n0)
   x <- .Call(P_makeLBDP,params,ivps,t0)
   .Call(P_runLBDP,x,time) |>
@@ -44,9 +47,13 @@ runLBDP <- function (
 ##' @inheritParams simulate
 ##' @export
 continueLBDP <- function (
-  object, time, lambda = NA, mu = NA, psi = NA
+  object, time, lambda = NA, mu = NA, psi = NA, r = NA
 ) {
-  params <- c(lambda=lambda,mu=mu,psi=psi)
+  if (!isTRUE(is.na(r))) {
+    if (!is.numeric(r) || length(r) != 1L || !is.finite(r) || r < 0 || r > 1)
+      pStop(sQuote("r")," must be between 0 and 1.")
+  }
+  params <- c(lambda=lambda,mu=mu,psi=psi,r=r)
   x <- .Call(P_reviveLBDP,object,params)
   .Call(P_runLBDP,x,time) |>
     structure(model="LBDP",class=c("gpsim","gpgen"))

R/lbdp_exact.R

@@ -5,12 +5,17 @@
 ##' @return \code{lbdp_exact} returns the log likelihood of the genealogy.
 ##' Note that the time since the most recent sample is informative.
 ##' @param x genealogy in \pkg{phylopomp} format (i.e., an object that inherits from \sQuote{gpgen}).
+##' @param r probability that a sampled lineage is removed (must be between 0 and 1; must be 0 for \code{lbdp_exact})
 ##' @references
 ##' \Stadler2010
 ##'
 ##' \King2024
 ##' @export
-lbdp_exact <- function (x, lambda, mu, psi, n0 = 1) {
+lbdp_exact <- function (x, lambda, mu, psi, r = 0, n0 = 1) {
+  if (!is.numeric(r) || length(r) != 1L || !is.finite(r) || r < 0 || r > 1)
+    pStop(sQuote("r")," must be between 0 and 1.")
+  if (!isTRUE(all.equal(r, 0)))
+    pStop(sQuote("r")," must be 0 for ",sQuote("lbdp_exact"),".")
   x |>
     lineages(prune=TRUE,obscure=TRUE) |>
     encode_data() -> data

R/lbdp_pomp.R

@@ -5,8 +5,12 @@
 ##' @importFrom pomp pomp onestep covariate_table
 ##' @inheritParams lbdp_exact
 ##' @export
-lbdp_pomp <- function (x, lambda, mu, psi, n0 = 1, t0 = 0)
+lbdp_pomp <- function (x, lambda, mu, psi, r = 0, n0 = 1, t0 = 0)
 {
+  if (!is.numeric(r) || length(r) != 1L || !is.finite(r) || r < 0 || r > 1)
+    pStop(sQuote("r")," must be between 0 and 1.")
+  if (!isTRUE(all.equal(r, 0)))
+    pStop(sQuote("r")," must be 0 for ",sQuote("lbdp_pomp"),".")
   x |> gendat() -> gi
   n0 <- round(n0)
   if (n0 < 0)

man/lbdp.Rd

@@ -16,10 +16,11 @@ typedef struct {
   double lambda;
   double mu;
   double psi;
+  double r;
   int n0;
 } lbdp_parameters_t;
 
-using lbdp_proc_t = popul_proc_t<lbdp_state_t,lbdp_parameters_t,3>;
+using lbdp_proc_t = popul_proc_t<lbdp_state_t,lbdp_parameters_t,4>;
 using lbdp_genealogy_t = master_t<lbdp_proc_t,1>;
 
 template<>
@@ -29,6 +30,7 @@ std::string lbdp_proc_t::yaml (std::string tab) const {
     + YAML_PARAM(lambda)
     + YAML_PARAM(mu)
     + YAML_PARAM(psi)
+    + YAML_PARAM(r)
     + YAML_PARAM(n0);
   std::string s = tab + "state:\n"
     + YAML_STATE(n);
@@ -41,6 +43,9 @@ void lbdp_proc_t::update_params (double *p, int n) {
   PARAM_SET(lambda);
   PARAM_SET(mu);
   PARAM_SET(psi);
+  PARAM_SET(r);
+  if (!R_FINITE(params.r) || params.r < 0 || params.r > 1)
+    err("r must be between 0 and 1.");
   if (m != n) err("wrong number of parameters!");
 }
 
@@ -57,7 +62,8 @@ double lbdp_proc_t::event_rates (double *rate, int n) const {
   double total = 0;
   RATE_CALC(params.lambda * state.n);
   RATE_CALC(params.mu * state.n);
-  RATE_CALC(params.psi * state.n);
+  RATE_CALC(params.psi * params.r * state.n);
+  RATE_CALC(params.psi * (1 - params.r) * state.n);
   if (m != n) err("wrong number of events!");
   return total;
 }
@@ -78,6 +84,9 @@ void lbdp_genealogy_t::jump (int event) {
     state.n -= 1; death();
     break;
   case 2:
+    state.n -= 1; sample_death();
+    break;
+  case 3:
     sample();
     break;
   default:                      // #nocov

src/lbdp.cc

@@ -0,0 +1,31 @@
+options(digits=3)
+suppressPackageStartupMessages({
+  library(phylopomp)
+})
+
+check_sampling <- function(r_value, expected_sat) {
+  set.seed(10101)
+  runLBDP(time=1,lambda=0,mu=0,psi=10,r=r_value,n0=5) |>
+    gendat() -> gi
+  sample_idx <- gi$nodetype == 1L
+  if (!any(sample_idx))
+    stop("no sample nodes observed for r=", r_value)
+  if (!all(gi$saturation[sample_idx] == expected_sat))
+    stop("unexpected sample saturation for r=", r_value)
+}
+
+check_sampling(0, 1)
+check_sampling(1, 0)
+
+check_bad_r <- function(r_value) {
+  ok <- FALSE
+  tryCatch(
+    runLBDP(time=1,lambda=0,mu=0,psi=1,r=r_value,n0=1),
+    error = function(e) ok <<- TRUE
+  )
+  if (!ok)
+    stop("expected error for r=", r_value)
+}
+
+check_bad_r(-0.1)
+check_bad_r(1.1)

tests/lbdp_bddr.R

@@ -9,7 +9,8 @@
 ##' @param lambda per capita birth rate
 ##' @param mu per capita death rate
 ##' @param psi per capita sampling rate
-##' @param n0 initial population size
+##' @param r probability that a sampled lineage is removed (must be between 0 and 1; r=0 is non-destructive; r=1 is destructive)
+##' @param n0 population size at time t0
 ##' @inheritParams sir
 ##' @return \code{runLBDP} and \code{continueLBDP} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{LBDP}.
 ##'
@@ -19,9 +20,9 @@ NULL
 ##' @export
 runLBDP <- function (
   time, t0 = 0,
-  lambda = 1.3, mu = 1, psi = 1, n0 = 10
+  lambda = 2, mu = 1, psi = 1, r = 0, n0 = 5
 ) {
-  params <- c(lambda=lambda,mu=mu,psi=psi)
+  params <- c(lambda=lambda,mu=mu,psi=psi,r=r)
   ivps <- c(n0=n0)
   x <- .Call(P_makeLBDP,params,ivps,t0)
   .Call(P_runLBDP,x,time) |>
@@ -32,10 +33,10 @@ runLBDP <- function (
 ##' @export
 continueLBDP <- function (
   object, time,
-  lambda = NA, mu = NA, psi = NA
+  lambda = NA, mu = NA, psi = NA, r = NA
 ) {
   params <- c(
-    lambda=lambda,mu=mu,psi=psi
+    lambda=lambda,mu=mu,psi=psi,r=r
   )
   x <- .Call(P_reviveLBDP,object,params)
   .Call(P_runLBDP,x,time) |>

yaml/R/lbdp.R

@@ -8,7 +8,7 @@ parameter:
     name: lambda
     description: per capita birth rate
     type: double
-    default: 1.3
+    default: 2
   mu:
     name: mu
     description: per capita death rate
@@ -19,11 +19,16 @@ parameter:
     description: per capita sampling rate
     default: 1.0
     type: double
+  r:
+    name: r
+    description: probability that a sampled lineage is removed (must be between 0 and 1; r=0 is non-destructive; r=1 is destructive)
+    default: 0.0
+    type: double
 ivp:
   n0:
     name: n0
-    description: initial population size
-    default: 10
+    description: population size at time t0
+    default: 5
     type: int
 state:
   'n':
@@ -39,6 +44,9 @@ events:
   death:
     jump: state.n -= 1; death();
     rate: params.mu * state.n
+  sample_death:
+    rate: params.psi * params.r * state.n
+    jump: state.n -= 1; sample_death();
   sample:
-    rate: params.psi * state.n
+    rate: params.psi * (1 - params.r) * state.n
     jump: sample();

yaml/lbdp.yml

@@ -16,10 +16,11 @@ typedef struct {
   double lambda;
   double mu;
   double psi;
+  double r;
   int n0;
 } lbdp_parameters_t;
 
-using lbdp_proc_t = popul_proc_t<lbdp_state_t,lbdp_parameters_t,3>;
+using lbdp_proc_t = popul_proc_t<lbdp_state_t,lbdp_parameters_t,4>;
 using lbdp_genealogy_t = master_t<lbdp_proc_t,1>;
 
 template<>
@@ -29,6 +30,7 @@ std::string lbdp_proc_t::yaml (std::string tab) const {
     + YAML_PARAM(lambda)
     + YAML_PARAM(mu)
     + YAML_PARAM(psi)
+    + YAML_PARAM(r)
     + YAML_PARAM(n0);
   std::string s = tab + "state:\n"
     + YAML_STATE(n);
@@ -41,6 +43,7 @@ void lbdp_proc_t::update_params (double *p, int n) {
   PARAM_SET(lambda);
   PARAM_SET(mu);
   PARAM_SET(psi);
+  PARAM_SET(r);
   if (m != n) err("wrong number of parameters!");
 }
 
@@ -57,7 +60,8 @@ double lbdp_proc_t::event_rates (double *rate, int n) const {
   double total = 0;
   RATE_CALC(params.lambda * state.n);
   RATE_CALC(params.mu * state.n);
-  RATE_CALC(params.psi * state.n);
+  RATE_CALC(params.psi * params.r * state.n);
+  RATE_CALC(params.psi * (1 - params.r) * state.n);
   if (m != n) err("wrong number of events!");
   return total;
 }
@@ -78,6 +82,9 @@ void lbdp_genealogy_t::jump (int event) {
     state.n -= 1; death();
     break;
   case 2:
+    state.n -= 1; sample_death();
+    break;
+  case 3:
     sample();
     break;
   default:                      // #nocov