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.
library(FlywayNet)
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.
<- generate_toy_migration()
migr 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 )
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.
<- estimate_migration_ABC( migr, verbose=FALSE ) # takes several seconds, 32 steps required
migr_ABC 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.
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
<- estimate_migration_MCEM( migr,
migr_MCEM1 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
<- generate_modified_migration(migr,
migr2 sojourn_mode=c(rep("unif", 4), "no"), sojourn_domain=c(1,5), transition_mode="unif")
<- estimate_migration_MCEM(migr2, MC_algo="MH",
migr_MCEM2 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.
<- reestimate_migration_from_MCEMruns(list(migr_MCEM1, migr_MCEM2))
migr_MCEM 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
The package also provide a naive method that estimates parameters from given individual trajectories.
set.seed(123) # just to be able to reproduce results
<- generate_trajectories( migr )
traj <- estimate_migration_from_trajectories(migr, traj)
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
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)
<- get_paramlikelihood( migr_ABC, migr$observation, n=10000)
L_ABC print( L_ABC )
#> [1] -418.1835
<- get_paramlikelihood( migr_MCEM, migr$observation, n=10000)
L_MCEM print( L_MCEM )
#> [1] -402.9847
<- get_paramlikelihood( migr_traj, migr$observation, n=10000)
L_traj print( L_traj )
#> [1] -428.7315
Here the best (maximum) log-likelihood is for MCEM parameters.
Considering choosing MCEM parameters found, it is possible to simulate trajectories.
$transition_law_param <- migr_MCEM$estimation_method$output$transition_law_param
migr_MCEM$sojourn_law_param<- migr_MCEM$estimation_method$output$sojourn_law_param
migr_MCEMset.seed(123) # just to be able to reproduce results
<- generate_trajectories( migr_MCEM )
traj 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.
<- get_counts( migr_ABC, traj ) counts
Finally, it is possible write and read migrations structures.
# write_migration (migr_ABC, 'migr_MCEM.txt')
# migr_ABC <- read_migration ('migr_MCEM.txt')