Introduction

In FlywayNet: Estimate Migration Birds Flyways Network with HSMM

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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)

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 )

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

plot( migr )

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)
print(migr_ABC$estimation_method$output$sojourn_law_param)

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) 
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) 

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)
print(migr_MCEM$estimation_method$output$sojourn_law_param)

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)
print(migr_traj$estimation_method$output$sojourn_law_param)

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 )
L_MCEM <- get_paramlikelihood( migr_MCEM, migr$observation, n=10000)
print( L_MCEM )
L_traj <- get_paramlikelihood( migr_traj, migr$observation, n=10000)
print( L_traj )

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 )

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')

FlywayNet documentation built on March 18, 2022, 7:21 p.m.