Longitudinal sampling with SEEPS

library(SEEPS)
set.seed(1945)

In this vingette, we address how to perform longitudinal sampling. In most applications, data is not taken from a single time point, but rather taken over investigations and multiple years.

We’ll look at how to implement each of these components in SEEPS.

Simulation

The basic simulation function gen_transmission_history_exponential_constant generates a transmission history for a single population. However, the act of sampling a population removes a set of individuals from a population. Usually, we assume that sampling an individual would effectively remove them from the population (from the large-population view). To properly account for this, we need to use the more detailed API for designing and implementing the simulation.

The simulation API

The lower level simulation API is concentrated on the step function. Users seeking to develop new simulation methods should concentrate on using this function.

The step function takes in a state at integer time t with a list of parameters, and returns the state at time t+1. This setup is designed to be comfortable to users familiar with dynamical systems which often express \(x(t_{n+1}) = f(x(t_n); \theta)\).

The wrap_parameters function is provided to encapsulate simulation parameters into a named list. Given the list of parameters, the initialize method generates an initial state x0 with introduced one active infection.

SEEPS provides two additional functions to help with simulations that can perform multiple steps. * gen_const_phase performs a fixed number of steps and returns the final state. * gen_exp_phase performs steps until the the number of active cases is at least the minimum_population (a parameter) and is at least 90% of the maximum_population_target (another parameter).

Removing individuals from the active population

Observing the simulation to draw a sample generally removes the individual from the active population. SEEPS provides the remove_samples function to perform this operation. remove_samples accepts a state and a vector of sample indices to remove from the simulation. A new state is returned with the samples removed.

state <- SEEPS::remove_samples(state, c(1, 2, 3))

A simple example

We’ll start with a simple example. We’ll simulate a small population with a a target maximum population size of 20 individuals and sample 5 individuals every 12 months, starting when the population is at least 18 individuals (90% of 20), with \(R_0 = 3\).

states_sampled <- list()
samples_all <- list()

samples_per_year <- 5
# Discover each contact with 80% probability
contact_tracing_parameter <- 0.8

sim_parameters <- SEEPS::wrap_parameters(
    minimum_population = 6,
    offspring_rate_fn = SEEPS::get_biphasic_HIV_rate(params = list("R0" = 3)),
    maximum_population_target = 20,
    total_steps = 0, # We are going to manually control this so the value will not be used
    spike_root = FALSE
)

state <- SEEPS::initialize(sim_parameters) # Setup
state <- SEEPS::gen_exp_phase(state, sim_parameters) # Run until we have at least 90% of the maximum population

# Take a first sample at the end of the exponential phase
sample_taken <- SEEPS::contact_traced_uniform_restarts_ids(
    state$active, state$parents, samples_per_year, contact_tracing_parameter
)
# Record the sampling
samples_all[[1]] <- list("curr_step" = state$curr_step, "samples" = sample_taken$samples)
# Remove from the population
state <- SEEPS::remove_samples(state, sample_taken$samples)

# Now simulate forward
for (time_pointer in 1:3) {
    state <- SEEPS::gen_const_phase(state, sim_parameters, num_steps = 12)

    sample_taken <- SEEPS::contact_traced_uniform_restarts_ids(
        state$active, state$parents, samples_per_year, contact_tracing_parameter
    )
    # Record the sampling
    samples_all[[time_pointer + 1]] <- list("curr_step" = state$curr_step, "samples" = sample_taken$samples)
    # Remove from the population
    state <- SEEPS::remove_samples(state, sample_taken$samples)
}

# Now that we have collected all of the sampling data, we can flatten it into a single list

sample_times_ <- list()
sample_ids <- list()
ptr <- 1
for (i in seq_along(samples_all)) {
    sample_times_[[ptr]] <- samples_all[[i]]$curr_step
    sample_ids[[ptr]] <- samples_all[[i]]$samples
    ptr <- ptr + 1
}

# Convert to

genealogy <- SEEPS::reduce_transmission_history_bpb2(
    samples = sample_ids,
    parents = state$parents,
    current_step = sample_times_
)

# Convert to a phylogeny

phylogeny <- SEEPS::geneology_to_phylogeny_bpb(
    transmission_history = genealogy$parents,
    infection_times = genealogy$transmission_times,
    sample_times = genealogy$sample_times,
    a = 5, b = 5,
    leaf_sample_ids = genealogy$transformed_sample_indices
)

# Convert to a newick string for plotting
newick_tree_string <- SEEPS::phylogeny_to_newick(phylogeny$phylogeny, mode = "mean", label_mode = "abs")

# Plot the tree
library(ape)
tree <- drop.tip(read.tree(text = newick_tree_string), "")
plot(tree, cex = 0.5)
axis(1, las = 1, col.axis = "black", col.ticks = "black")

Reporting the entire history of the simulation

One can also use this approach to record when every individual is removed from the population. Rather than recording a set of samples, we record whenever an individual is removed from the active population.

We define an auxiliary function for simulations.

gen_transmission_history_exponential_constant_local <- function(minimum_population, # nolint: object_length_linter
                                                                offspring_rate_fn, maximum_population_target, total_steps,
                                                                spike_root = FALSE) {
    # Top level/user function. Users should start with this function.

    parameters <- SEEPS::wrap_parameters(
        minimum_population, offspring_rate_fn,
        maximum_population_target, total_steps,
        spike_root
    )

    state <- SEEPS::initialize(parameters)
    SEEPS::validate_state(state) # Check that the state is valid
    all_states <- list()
    all_states[[1]] <- state
    state <- gen_exp_phase(state, parameters, all_states = TRUE)
    all_states <- c(all_states, state)

    state <- gen_const_phase(all_states[[length(all_states)]], parameters, num_steps = total_steps, all_states = TRUE)
    all_states <- c(all_states, state)
    return(lapply(all_states, SEEPS::clean_up))
}

We’ll use a different rate function, where the rate is 50 times higher in the first 3 months than in the remainder of the infection. We set the expected length of the each infection to 30 months (2.5 years).

offspring_rate_fn <- get_biphasic_HIV_rate_function(50, 3, 30, params = list("R0" = 1.3))

# maximum_population_target <- 30
maximum_population_target <- 45
simulator_result <- gen_transmission_history_exponential_constant_local(
    minimum_population = 6,
    offspring_rate_fn = offspring_rate_fn,
    # offspring_rate_fn = SEEPS::get_biphasic_HIV_rate(params = list("R0" = 1.3)),
    # offspring_rate_fn = SEEPS::get_biphasic_HIV_rate(params = list("R0" = 1.2)),
    maximum_population_target = maximum_population_target,
    # total_steps = 48
    total_steps = 60
)

removed_list <- list()
sample_times <- list()
birth_list <- list()

ptr <- 1
# Don't include the 0 individual as they are not part of the population
# 0 is the root of the first infection.
removed_list[[ptr]] <- c()
sample_times[[ptr]] <- c()
birth_list[[ptr]] <- c(1)
ptr <- 2

for (i in 2:length(simulator_result)) {
    # Get the active nodes at i-1 and i
    active_i <- simulator_result[[i]]$active
    active_i_1 <- simulator_result[[i - 1]]$active

    # Compare the active nodes at i-1 and i
    # If the active node is in i-1 but not i, add it to a "removed" list at position i
    removed_active <- active_i_1[!(active_i_1 %in% active_i)]
    new_active <- active_i[!(active_i %in% active_i_1)]
    if (length(removed_active) == 0 && length(new_active) == 0) {
        next # Nothing to add to the samples list
    }

    # Add to removed list
    removed_list[[ptr]] <- removed_active
    birth_list[[ptr]] <- new_active
    # Record the sample time
    sample_times[[ptr]] <- simulator_result[[i]]$t_end
    ptr <- ptr + 1
}

# Final step - get the active nodes at the last step
active_final <- simulator_result[[length(simulator_result)]]$active
# Store it
removed_list[[ptr]] <- active_final
# Record the sample time
sample_times[[ptr]] <- simulator_result[[length(simulator_result)]]$t_end

# Now record the transmission history
print(removed_list)
#> [[1]]
#> NULL
#> 
#> [[2]]
#> numeric(0)
#> 
#> [[3]]
#> numeric(0)
#> 
#> [[4]]
#> numeric(0)
#> 
#> [[5]]
#> numeric(0)
#> 
#> [[6]]
#> numeric(0)
#> 
#> [[7]]
#> numeric(0)
#> 
#> [[8]]
#> numeric(0)
#> 
#> [[9]]
#> numeric(0)
#> 
#> [[10]]
#> numeric(0)
#> 
#> [[11]]
#> numeric(0)
#> 
#> [[12]]
#> numeric(0)
#> 
#> [[13]]
#> numeric(0)
#> 
#> [[14]]
#> numeric(0)
#> 
#> [[15]]
#> numeric(0)
#> 
#> [[16]]
#> [1] 8
#> 
#> [[17]]
#> [1] 6
#> 
#> [[18]]
#> numeric(0)
#> 
#> [[19]]
#> [1] 10 19
#> 
#> [[20]]
#> [1]  4 17
#> 
#> [[21]]
#> numeric(0)
#> 
#> [[22]]
#> [1] 14
#> 
#> [[23]]
#> [1] 18
#> 
#> [[24]]
#> [1] 11 43
#> 
#> [[25]]
#> [1] 30
#> 
#> [[26]]
#> [1] 20 44
#> 
#> [[27]]
#> numeric(0)
#> 
#> [[28]]
#> [1]  2 39
#> 
#> [[29]]
#> [1] 36 38
#> 
#> [[30]]
#> [1] 41 46
#> 
#> [[31]]
#> [1] 23
#> 
#> [[32]]
#> [1]  1  5 12
#> 
#> [[33]]
#> [1]  3  9 29
#> 
#> [[34]]
#> [1] 47
#> 
#> [[35]]
#> [1]  7 16 22 45
#> 
#> [[36]]
#> [1] 13 24 51
#> 
#> [[37]]
#> [1] 15 26 28
#> 
#> [[38]]
#> [1] 31 40
#> 
#> [[39]]
#> [1] 34 48
#> 
#> [[40]]
#> [1] 25 32 49
#> 
#> [[41]]
#> [1] 21 54
#> 
#> [[42]]
#> [1] 35
#> 
#> [[43]]
#> [1] 27 52
#> 
#> [[44]]
#> [1] 33
#> 
#> [[45]]
#> [1] 42 60
#> 
#> [[46]]
#> [1] 37
#> 
#> [[47]]
#> [1] 56
#> 
#> [[48]]
#> numeric(0)
#> 
#> [[49]]
#> [1] 57
#> 
#> [[50]]
#> [1] 58 62 73
#> 
#> [[51]]
#> [1] 64 74
#> 
#> [[52]]
#> [1] 50
#> 
#> [[53]]
#> numeric(0)
#> 
#> [[54]]
#> [1] 69 84 85
#> 
#> [[55]]
#> [1] 63 79 82
#> 
#> [[56]]
#> numeric(0)
#> 
#> [[57]]
#> [1] 65 66 67 75
#> 
#> [[58]]
#> [1] 83 90
#> 
#> [[59]]
#> [1] 53 61
#> 
#> [[60]]
#> [1] 55
#> 
#> [[61]]
#> [1] 72 87 89
#> 
#> [[62]]
#> [1] 68 77 80 92
#> 
#> [[63]]
#> [1] 59
#> 
#> [[64]]
#> [1] 76
#> 
#> [[65]]
#> numeric(0)
#> 
#> [[66]]
#> [1] 88
#> 
#> [[67]]
#> [1] 78 91
#> 
#> [[68]]
#> [1] 99
#> 
#> [[69]]
#> [1]  70  71 109
#> 
#> [[70]]
#> numeric(0)
#> 
#> [[71]]
#>  [1]  81  86  93  94  95  96  97  98 100 101 102 103 104 105 106 107 108 110 111
#> [20] 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
#> [39] 131 132 133 134 135 136
print(sample_times)
#> [[1]]
#> NULL
#> 
#> [[2]]
#> [1] 2
#> 
#> [[3]]
#> [1] 3
#> 
#> [[4]]
#> [1] 4
#> 
#> [[5]]
#> [1] 5
#> 
#> [[6]]
#> [1] 6
#> 
#> [[7]]
#> [1] 7
#> 
#> [[8]]
#> [1] 8
#> 
#> [[9]]
#> [1] 9
#> 
#> [[10]]
#> [1] 10
#> 
#> [[11]]
#> [1] 11
#> 
#> [[12]]
#> [1] 12
#> 
#> [[13]]
#> [1] 13
#> 
#> [[14]]
#> [1] 14
#> 
#> [[15]]
#> [1] 15
#> 
#> [[16]]
#> [1] 20
#> 
#> [[17]]
#> [1] 21
#> 
#> [[18]]
#> [1] 22
#> 
#> [[19]]
#> [1] 23
#> 
#> [[20]]
#> [1] 24
#> 
#> [[21]]
#> [1] 25
#> 
#> [[22]]
#> [1] 26
#> 
#> [[23]]
#> [1] 27
#> 
#> [[24]]
#> [1] 28
#> 
#> [[25]]
#> [1] 29
#> 
#> [[26]]
#> [1] 30
#> 
#> [[27]]
#> [1] 31
#> 
#> [[28]]
#> [1] 32
#> 
#> [[29]]
#> [1] 33
#> 
#> [[30]]
#> [1] 34
#> 
#> [[31]]
#> [1] 35
#> 
#> [[32]]
#> [1] 36
#> 
#> [[33]]
#> [1] 37
#> 
#> [[34]]
#> [1] 38
#> 
#> [[35]]
#> [1] 39
#> 
#> [[36]]
#> [1] 40
#> 
#> [[37]]
#> [1] 41
#> 
#> [[38]]
#> [1] 42
#> 
#> [[39]]
#> [1] 43
#> 
#> [[40]]
#> [1] 44
#> 
#> [[41]]
#> [1] 45
#> 
#> [[42]]
#> [1] 46
#> 
#> [[43]]
#> [1] 47
#> 
#> [[44]]
#> [1] 48
#> 
#> [[45]]
#> [1] 49
#> 
#> [[46]]
#> [1] 50
#> 
#> [[47]]
#> [1] 51
#> 
#> [[48]]
#> [1] 52
#> 
#> [[49]]
#> [1] 53
#> 
#> [[50]]
#> [1] 54
#> 
#> [[51]]
#> [1] 55
#> 
#> [[52]]
#> [1] 56
#> 
#> [[53]]
#> [1] 57
#> 
#> [[54]]
#> [1] 58
#> 
#> [[55]]
#> [1] 59
#> 
#> [[56]]
#> [1] 60
#> 
#> [[57]]
#> [1] 61
#> 
#> [[58]]
#> [1] 62
#> 
#> [[59]]
#> [1] 63
#> 
#> [[60]]
#> [1] 64
#> 
#> [[61]]
#> [1] 65
#> 
#> [[62]]
#> [1] 66
#> 
#> [[63]]
#> [1] 67
#> 
#> [[64]]
#> [1] 68
#> 
#> [[65]]
#> [1] 69
#> 
#> [[66]]
#> [1] 71
#> 
#> [[67]]
#> [1] 72
#> 
#> [[68]]
#> [1] 73
#> 
#> [[69]]
#> [1] 74
#> 
#> [[70]]
#> [1] 75
#> 
#> [[71]]
#> [1] 75
full_transmission_history <- reduce_transmission_history_mt(
    parents = simulator_result[[length(simulator_result)]]$parents,
    samples = removed_list,
    current_step = sample_times
)
#> [1] "Samples"
#>   [1]   8   6  10  19   4  17  14  18  11  43  30  20  44   2  39  36  38  41
#>  [19]  46  23   1   5  12   3   9  29  47   7  16  22  45  13  24  51  15  26
#>  [37]  28  31  40  34  48  25  32  49  21  54  35  27  52  33  42  60  37  56
#>  [55]  57  58  62  73  64  74  50  69  84  85  63  79  82  65  66  67  75  83
#>  [73]  90  53  61  55  72  87  89  68  77  80  92  59  76  88  78  91  99  70
#>  [91]  71 109  81  86  93  94  95  96  97  98 100 101 102 103 104 105 106 107
#> [109] 108 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
#> [127] 127 128 129 130 131 132 133 134 135 136
#> [1] "Sample times"
#>   [1] 20 21 23 23 24 24 26 27 28 28 29 30 30 32 32 33 33 34 34 35 36 36 36 37 37
#>  [26] 37 38 39 39 39 39 40 40 40 41 41 41 42 42 43 43 44 44 44 45 45 46 47 47 48
#>  [51] 49 49 50 51 53 54 54 54 55 55 56 58 58 58 59 59 59 61 61 61 61 62 62 63 63
#>  [76] 64 65 65 65 66 66 66 66 67 68 71 72 72 73 74 74 74 75 75 75 75 75 75 75 75
#> [101] 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75 75
#> [126] 75 75 75 75 75 75 75 75 75 75 75

Now we can visualize the entire transmission history.

string <- SEEPS::phylogeny_to_newick(full_transmission_history$geneology, mode = "mean", label_mode = "abs")
tree <- read.tree(text = string)
plot(ladderize(tree))
axis(1, las = 1, col.axis = "black", col.ticks = "black")

We can also plot out the number of new infections per time step.

births_ <- lapply(birth_list, length)
deaths <- lapply(removed_list, length)
# births_ <- births[1:length(births) - 1]  # Don't report the final step
deaths_ <- deaths[1:length(deaths) - 1] # Don't report the final step
sample_times_ <- sample_times[1:length(sample_times) - 1] # Don't report the final step
deaths_ <- unlist(deaths_)
births_ <- unlist(births_)
print(length(births_))
#> [1] 70
print(length(deaths_))
#> [1] 70
plot(seq_along(births_), births_, col = "black", type = "S") # label="Births"
lines(seq_along(deaths_), deaths_, col = "red", type = "S") # label="Deaths"
legend(39, y = NULL, legend = c("Births", "Deaths"), col = c("black", "red"), lty = 1:2, cex = 0.8)


print(length(births_))
#> [1] 70
print(length(deaths_))
#> [1] 70

deltas <- births_ - deaths_
print(length(deltas))
#> [1] 70
plot(1:length(deltas), cumsum(deltas), ylab = "Number of active infections", xlab = "Time index [steps]", ylim = c(-1, maximum_population_target + 2), type = "S")
# Draw a horizontal line at the maximum population target
abline(h = maximum_population_target, col = "red", lty = 2)