Introduction

This a quick start for R FlywayNet package that allows to infer the structure of migratory bird flyway networks given observed counts.

The objective is just to illustrate how it is possible to use the package.

For a description of the package, see the package documentation.

Load the package

library(FlywayNet)

Define a toy migration structure with 5 sites

The function generate_toy_migration allows to define a toy example of a migratory bird flyway network: 100 birds, 5 sites with some transitions forbidden. A time interval, called horizon, of 20 weeks (time steps) is considered. The function creates a migration structure (see package documentation) with observed counts every week for each site.

migr <- generate_toy_migration()
print(migr)
#> $site_name
#> [1] "departure" "site1"     "site2"     "site3"     "arrival"  
#> 
#> $link_knowledge
#>       [,1]  [,2]  [,3]  [,4]  [,5]
#> [1,] FALSE  TRUE FALSE FALSE FALSE
#> [2,] FALSE FALSE  TRUE  TRUE FALSE
#> [3,] FALSE FALSE FALSE FALSE  TRUE
#> [4,] FALSE FALSE FALSE FALSE  TRUE
#> [5,] FALSE FALSE FALSE FALSE FALSE
#> 
#> $flight_duration
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    1    0    0    0
#> [2,]    0    0    1    2    0
#> [3,]    0    0    0    0    2
#> [4,]    0    0    0    0    1
#> [5,]    0    0    0    0    0
#> 
#> $initial_state
#> [1] 100   0   0   0   0
#> 
#> $horizon
#> [1] 20
#> 
#> $death_probability
#> [1] 0.1 0.2 0.2 0.1 1.0
#> 
#> $transition_law_type
#> [1] "multinomial"
#> 
#> $transition_law_param
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0  0.9 0.00 0.00  0.0
#> [2,]    0  0.0 0.55 0.25  0.0
#> [3,]    0  0.0 0.00 0.00  0.8
#> [4,]    0  0.0 0.00 0.00  0.9
#> [5,]    0  0.0 0.00 0.00  0.0
#> 
#> $sojourn_law_type
#> [1] "Poisson"
#> 
#> $sojourn_law_param
#> [1] 3 4 2 4 5
#> 
#> $observation_law_type
#> [1] "Poisson"
#> 
#> $observation_law_param
#> NULL
#> 
#> $observation
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> [1,]   NA   91   58   30   11    6    2    2    0     0    NA    NA    NA    NA
#> [2,]   NA    3   18   49   51   45   46   33   20    12    11     5     7     6
#> [3,]   NA    0    0    0    3   12   20    5    8    22     7     4     2     5
#> [4,]   NA    0    0    0    0    1    1    0   11     9     4     5     7     7
#> [5,]   NA    0    0    0    0    0    0    0    0    16    25    15    37    39
#> [6,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA    NA
#> [7,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA    NA
#>      [,15] [,16] [,17] [,18] [,19] [,20] [,21]
#> [1,]    NA    NA    NA    NA    NA    NA    NA
#> [2,]     3     0     0    NA    NA    NA    NA
#> [3,]     2     0     0     1     2     1     0
#> [4,]    10     8     4     1     2     0     0
#> [5,]    44    49    41    49    47    58    52
#> [6,]    NA    NA    NA    NA    NA    NA    NA
#> [7,]    NA    NA    NA    NA    NA    NA    NA
#> 
#> attr(,"class")
#> [1] "migration"

It is possible to check the validity of the structure. This is useful when the structure is created using the new_migration function that allows to set all attributes. The function validate_migration reports identified problems (nothing if OK).

validate_migration( migr )
#> list()

As it is possible to set sites link knowledge ( see migr$link_knowledge ), it is possible to visualize the possible network.

plot( migr )

#> IGRAPH 8a22443 DNW- 5 5 -- 
#> + attr: name (v/c), TRUE (v/c), weight (e/n)
#> + edges from 8a22443 (vertex names):
#> [1] departure->site1   site1    ->site2   site1    ->site3   site2    ->arrival
#> [5] site3    ->arrival

It is also possible to visualize the observed count (migr$observation attribute).

plot_observedcounts( migr, migr$observation )

Estimate parameters with ABC

Consider that we want to estimate both transition and sojourn duration parameters
considering their respective laws set in migr \(transition_law_type and migr\)sojourn_law_type attributes. We first try the ABC method.

set.seed(123)  # just to be able to reproduce results
# We set verbose to FALSE because in this vignette the display
# is much more verbose than normal execution.
migr_ABC <- estimate_migration_ABC( migr, verbose=FALSE )  # takes several seconds, 32 steps required
print(migr_ABC$estimation_method$output$transition_law_param)
#>      [,1] [,2]      [,3]      [,4] [,5]
#> [1,]    0  0.9 0.0000000 0.0000000  0.0
#> [2,]    0  0.0 0.4819864 0.3180136  0.0
#> [3,]    0  0.0 0.0000000 0.0000000  0.8
#> [4,]    0  0.0 0.0000000 0.0000000  0.9
#> [5,]    0  0.0 0.0000000 0.0000000  0.0
print(migr_ABC$estimation_method$output$sojourn_law_param)
#> [1] 1.3959284 1.0175416 0.6463718 0.2047319 5.0000000

Note that in real application, it would be necessary to refine arguments of the function.

Estimate parameters with MCEM

Lets now try the other available MCEM method. For this method, it is necessary to estimate parameters with different initial values. Here we just consider two initial values for illustrative purpose but note that real application more initial values have to be considered.

set.seed(123)  # just to be able to reproduce results
migr_MCEM1 <- estimate_migration_MCEM( migr,
    log_transitions = TRUE, log_sojourns = TRUE, log_loglikelihood = TRUE) 
#> -------------- ITER EM MH 1
#> -------------- ITER EM MH 2
#> -------------- ITER EM MH 3
#> -------------- ITER EM MH 4
#> -------------- ITER EM MH 5
#> -------------- ITER EM MH 6
#> -------------- ITER EM MH 7
#> -------------- ITER EM MH 8
#> -------------- ITER EM MH 9
#> -------------- ITER EM MH 10
migr2 <- generate_modified_migration(migr, 
   sojourn_mode=c(rep("unif", 4), "no"), sojourn_domain=c(1,5), transition_mode="unif")
migr_MCEM2 <- estimate_migration_MCEM(migr2, MC_algo="MH",
    log_transitions = TRUE, log_sojourns = TRUE, log_loglikelihood = TRUE) 
#> -------------- ITER EM MH 1
#> -------------- ITER EM MH 2
#> -------------- ITER EM MH 3
#> -------------- ITER EM MH 4
#> -------------- ITER EM MH 5
#> -------------- ITER EM MH 6
#> -------------- ITER EM MH 7
#> -------------- ITER EM MH 8
#> -------------- ITER EM MH 9
#> -------------- ITER EM MH 10

As it is difficult to parametrize the MCEM method, we suggest to let iterate the method then run a post process with the reestimate_migration_from_MCEMruns function that will be able to detect better stopping criteria for each MCEM run.

migr_MCEM <- reestimate_migration_from_MCEMruns(list(migr_MCEM1, migr_MCEM2))  
print(migr_MCEM$estimation_method$output$transition_law_param)
#>      [,1] [,2]     [,3]     [,4] [,5]
#> [1,]    0  0.9 0.000000 0.000000  0.0
#> [2,]    0  0.0 0.515415 0.284585  0.0
#> [3,]    0  0.0 0.000000 0.000000  0.8
#> [4,]    0  0.0 0.000000 0.000000  0.9
#> [5,]    0  0.0 0.000000 0.000000  0.0
print(migr_MCEM$estimation_method$output$sojourn_law_param)
#> [1] 1.7609000 2.5716326 0.9350539 3.0218391 5.0000000

Estimate parameters from trajectories

The package also provide a naive method that estimates parameters from given individual trajectories.

set.seed(123)  # just to be able to reproduce results
traj <- generate_trajectories( migr )
migr_traj <- estimate_migration_from_trajectories(migr, traj)
print(migr_traj$estimation_method$output$transition_law_param)
#>      [,1] [,2]      [,3]      [,4] [,5]
#> [1,]    0  0.9 0.0000000 0.0000000  0.0
#> [2,]    0  0.0 0.4820513 0.3179487  0.0
#> [3,]    0  0.0 0.0000000 0.0000000  0.8
#> [4,]    0  0.0 0.0000000 0.0000000  0.9
#> [5,]    0  0.0 0.0000000 0.0000000  0.0
print(migr_traj$estimation_method$output$sojourn_law_param)
#> [1] 2.910000 4.336842 1.936170 3.709677 5.000000

Compare estimated parameters

To evaluate the results, the log-likelihood of parameters considering observations is computed. The number of simulations required (n argument), depends of the problem. Here to avoid long execution duration, we only required 10.000 simulations.

set.seed(123)
L_ABC <- get_paramlikelihood( migr_ABC, migr$observation, n=10000)
print( L_ABC )
#> [1] -418.1835
L_MCEM <- get_paramlikelihood( migr_MCEM, migr$observation, n=10000)
print( L_MCEM )
#> [1] -402.9847
L_traj <- get_paramlikelihood( migr_traj, migr$observation, n=10000)
print( L_traj )
#> [1] -428.7315

Here the best (maximum) log-likelihood is for MCEM parameters.

Generate trajectories when parameters are set

Considering choosing MCEM parameters found, it is possible to simulate trajectories.

migr_MCEM$transition_law_param <- migr_MCEM$estimation_method$output$transition_law_param
migr_MCEM$sojourn_law_param<- migr_MCEM$estimation_method$output$sojourn_law_param
set.seed(123)  # just to be able to reproduce results
traj <- generate_trajectories( migr_MCEM )
plot_trajectories ( migr_MCEM, traj )

#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> [1,]  100   85   52   25    8    4    0    0    0     0     0     0     0     0
#> [2,]    0    0   14   44   66   69   61   47   30    19     7     3     2     0
#> [3,]    0    0    0    0    0    3   11   13   16    18    11    11     6     4
#> [4,]    0    0    0    0    0    0    1    2    5    11    10    13    15    14
#> [5,]    0    0    0    0    0    0    0    0    0     2     8    13    27    32
#> [6,]    0    1    4    5    5   10   14   17   25    29    34    34    34    36
#> [7,]    0   14   30   26   21   14   13   21   24    21    30    26    16    14
#>      [,15] [,16] [,17] [,18] [,19] [,20] [,21]
#> [1,]     0     0     0     0     0     0     0
#> [2,]     0     0     0     0     0     0     0
#> [3,]     1     1     0     0     0     0     0
#> [4,]     9     8     7     5     4     2     0
#> [5,]    38    39    40    37    32    24    18
#> [6,]    41    46    51    55    63    72    80
#> [7,]    11     6     2     3     1     2     2

We can also get the counts associated to the generated trajectory with get_counts function.

counts <- get_counts( migr_ABC, traj )

Write and read migration structure

Finally, it is possible write and read migrations structures.

# write_migration (migr_ABC, 'migr_MCEM.txt')
# migr_ABC <- read_migration ('migr_MCEM.txt')