library(dplyr) # Data manipulation
library(ggplot2) # Plotting
library(tidyr) # Tidy data
library(sf) # Simple features for R (spatial objects)
library(sit) # Analyse MRR data from SIT experiments
Printing the sit
object gives a very short description
of its contents:
sit_prototype
#> ── sit ─────────────────────────────────────────────────────────────────────────
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Release events: 1 (areal) + 3 (point)
#> Release sites : 2
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Adult traps : 21 (sit)
#> Egg traps : 7 (control) + 14 (sit)
#> Trap types : 21 (BGS) + 21 (OVT)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Surveys : 3486 (adult) + 462 (egg)
Use summary()
to get a more detailed summary of each of
the sit
components:
summary(sit_prototype)
#> ── Releases ────────────────────────────────────────────────────────────────────
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Number of point releases: 1
#> Number of areal releases: 1
#>
#> Release events:
#> id type site_id date colour n geometry
#> 1 1 point 1 2019-11-24 23:00:00 yellow 10000 POINT (0 0)
#> 2 2 point 1 2019-11-30 23:00:00 red 10000 POINT (0 0)
#> 3 3 point 1 2019-12-12 23:00:00 blue 10000 POINT (0 0)
#> 4 4 areal 2 2019-12-20 23:00:00 pink 10000 POINT EMPTY
#>
#> Release sites:
#> Simple feature collection with 2 features and 1 field (with 1 geometry empty)
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: 0 xmax: 0 ymax: 0
#> Projected CRS: +proj=laea +x_0=0 +y_0=0 +lon_0=16.41472 +lat_0=48.23528
#> id geometry
#> 1 1 POINT (0 0)
#> 2 2 POINT EMPTY
#>
#> ── Traps ───────────────────────────────────────────────────────────────────────
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Number of traps in each group:
#> Area Type Coords n
#> 1 control OVT FALSE 7
#> 2 sit OVT FALSE 14
#> 3 sit BGS TRUE 21
#>
#> ── Adult surveys (total captures) ──────────────────────────────────────────────
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> , , sit
#>
#> female male
#> yellow NA 58
#> red NA 79
#> blue NA 86
#> pink NA 90
#> wild 226 299
#>
#>
#> ── Egg surveys (total captures) ────────────────────────────────────────────────
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Sterile Fertile
#> control 3407 18432
#> sit 4488 6808
Use plot()
to get a very basic image of the experimental
spatial set up.
plot(sit_prototype)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
You can make more professional and even interactive maps by
leveraging other R packages of your choice. This is beyond the scope of
sit
, but I’ll provide a few examples here using the package tmap
. To
learn more on mapping, see the chapter Making maps with
R in Robin Lovelace’s excellent book Geocomputation with
R.
library(tmap)
tmap_cmode <- tmap_mode("view")
#> tmap mode set to interactive viewing
tm_basemap("CartoDB") +
tm_shape(
sit_traps(sit_prototype, area = "sit", stage = "adult")
) +
tm_dots() +
tm_shape(
sit_revents(sit_prototype, type = "point") %>%
st_as_sf()
) +
tm_bubbles(
col = "colour",
jitter = 1 # Jitter a bit to see the overlapping release points
)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
tmap_mode(tmap_cmode)
#> tmap mode set to plotting
Note how we retrieved only adult
traps from the
sit
area, and only the point
releases, out of
all available traps and releases. We will see how that worked in the
following section.
Extraction functions allow retrieving the data contained in a
sit
object, while filtering specific groups as needed.
Here we describe the list of arguments of each extraction function, with a few examples. This is very much borrowed from the help pages, but focused on the extraction method only.
sit_trap_types()
x
|
A sit object.
|
This one is simple: it only takes a sit
object and
returns the list of trap types. E.g.
sit_trap_types(sit_prototype)
#> id name label stage description
#> 1 1 Ovitrap OVT egg NA
#> 2 2 BG-Sentinel BGS adult NA
#> 3 3 Human Landing Catch HLC adult NA
The result is a data.frame
itself. Thus, if you really
need to filter some observations, you can put your R-skills to use:
sit_traps()
x
|
A sit object.
|
type
|
Character vector. Labels of requested trap types. By default,
sit_trap_types()$label .
|
stage
|
Character vector. Either adult , egg or both
(default).
|
area
|
Character vector. Either control , sit or both
(default).
|
Retrieve the table of ovitraps:
sit_traps(sit_prototype, type = "OVT")
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Simple feature collection with 21 features and 8 fields (with 21 geometries empty)
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
#> Projected CRS: +proj=laea +x_0=0 +y_0=0 +lon_0=16.41472 +lat_0=48.23528
#> First 10 features:
#> id code area label type_name type_label type_stage type_description
#> 1 22 CS01 sit NA Ovitrap OVT egg NA
#> 2 23 CS02 sit NA Ovitrap OVT egg NA
#> 3 24 CS03 sit NA Ovitrap OVT egg NA
#> 4 25 CS04 sit NA Ovitrap OVT egg NA
#> 5 26 CS05 sit NA Ovitrap OVT egg NA
#> 6 27 CS06 sit NA Ovitrap OVT egg NA
#> 7 28 CS07 sit NA Ovitrap OVT egg NA
#> 8 29 CS08 sit NA Ovitrap OVT egg NA
#> 9 30 CS09 sit NA Ovitrap OVT egg NA
#> 10 31 CS10 sit NA Ovitrap OVT egg NA
#> geometry
#> 1 POINT EMPTY
#> 2 POINT EMPTY
#> 3 POINT EMPTY
#> 4 POINT EMPTY
#> 5 POINT EMPTY
#> 6 POINT EMPTY
#> 7 POINT EMPTY
#> 8 POINT EMPTY
#> 9 POINT EMPTY
#> 10 POINT EMPTY
Note that this coincides with the table of traps for the
stage = 'egg'
.
However, traps in the control
area are only a
subset:
sit_traps(sit_prototype, area = "control")
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Simple feature collection with 7 features and 8 fields (with 7 geometries empty)
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
#> Projected CRS: +proj=laea +x_0=0 +y_0=0 +lon_0=16.41472 +lat_0=48.23528
#> id code area label type_name type_label type_stage type_description
#> 1 36 BL01 control NA Ovitrap OVT egg NA
#> 2 37 BL02 control NA Ovitrap OVT egg NA
#> 3 38 BL03 control NA Ovitrap OVT egg NA
#> 4 39 BL04 control NA Ovitrap OVT egg NA
#> 5 40 BL05 control NA Ovitrap OVT egg NA
#> 6 41 BL06 control NA Ovitrap OVT egg NA
#> 7 42 BL07 control NA Ovitrap OVT egg NA
#> geometry
#> 1 POINT EMPTY
#> 2 POINT EMPTY
#> 3 POINT EMPTY
#> 4 POINT EMPTY
#> 5 POINT EMPTY
#> 6 POINT EMPTY
#> 7 POINT EMPTY
sit_revents()
x
|
A sit object.
|
type
|
Character vector. Which types of release events to
retrieve. Either point or areal or both
(default).
|
sit_revents(sit_prototype, type = 'point')
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> id type site_id date colour n geometry
#> 1 1 point 1 2019-11-24 23:00:00 yellow 10000 POINT (0 0)
#> 2 2 point 1 2019-11-30 23:00:00 red 10000 POINT (0 0)
#> 3 3 point 1 2019-12-12 23:00:00 blue 10000 POINT (0 0)
sit_revents(sit_prototype, type = 'areal')
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> id type site_id date colour n geometry
#> 1 4 areal 2 2019-12-20 23:00:00 pink 10000 POINT EMPTY
sit_revents(sit_prototype)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> id type site_id date colour n geometry
#> 1 1 point 1 2019-11-24 23:00:00 yellow 10000 POINT (0 0)
#> 2 2 point 1 2019-11-30 23:00:00 red 10000 POINT (0 0)
#> 3 3 point 1 2019-12-12 23:00:00 blue 10000 POINT (0 0)
#> 4 4 areal 2 2019-12-20 23:00:00 pink 10000 POINT EMPTY
sit_adult_surveys()
and
sit_egg_surveys()
Extract survey data, possibly filtering results.
Arguments:
x
|
A sit object
|
area
|
Character. Either control , sit or both
(default).
|
trap_type
|
Character. Any of sit_trap_types(x)$label . Typically,
OVT for sit_egg_surveys() and BGS
or HLC for sit_adult_surveys() .
|
following_releases
|
A sit_release object with a subset of release events or
missing (default) for all release events. Use in combination with
following_days to filter surveys of the target populations
within a given number of days after release. Note that counts of wild
populations will always be included in the results with
a value of NA in pop_col .
|
following_days
|
Integer or missing (default). Number of days after releases to return,
if following_releases is not missing.
|
species
|
Character or NULL (default). Filter surveys for a particular species, or ignore the species if NULL. |
Retrieve all the egg surveys in the control area:
sit_egg_surveys(sit_prototype, area = "control")
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> # A tibble: 154 × 6
#> id trap_code datetime_start datetime_end fertile n
#> <int> <chr> <dttm> <dttm> <lgl> <dbl>
#> 1 309 BL01 2019-11-21 23:00:00 2019-11-28 23:00:00 TRUE 41
#> 2 310 BL01 2019-11-21 23:00:00 2019-11-28 23:00:00 FALSE 37
#> 3 311 BL02 2019-11-21 23:00:00 2019-11-28 23:00:00 TRUE 102
#> 4 312 BL02 2019-11-21 23:00:00 2019-11-28 23:00:00 FALSE 20
#> 5 313 BL03 2019-11-21 23:00:00 2019-11-28 23:00:00 TRUE 204
#> 6 314 BL03 2019-11-21 23:00:00 2019-11-28 23:00:00 FALSE 18
#> 7 315 BL04 2019-11-21 23:00:00 2019-11-28 23:00:00 TRUE 354
#> 8 316 BL04 2019-11-21 23:00:00 2019-11-28 23:00:00 FALSE 54
#> 9 317 BL05 2019-11-21 23:00:00 2019-11-28 23:00:00 TRUE 31
#> 10 318 BL05 2019-11-21 23:00:00 2019-11-28 23:00:00 FALSE 10
#> # ℹ 144 more rows
Retrieve the adult surveys from the day that followed the areal release. Only counts from the target release, or the wild population, please.
sit_adult_surveys(
sit_prototype,
following_releases = sit_revents(sit_prototype, type = 'areal'),
following_days = 1
)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> # A tibble: 63 × 7
#> id trap_code pop_col datetime_start datetime_end sex n
#> <int> <chr> <chr> <dttm> <dttm> <chr> <dbl>
#> 1 2227 1 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 female 0
#> 2 2228 1 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 male 0
#> 3 2232 1 pink 2019-12-20 23:00:00 2019-12-21 23:00:00 male 0
#> 4 2233 2 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 female 1
#> 5 2234 2 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 male 2
#> 6 2238 2 pink 2019-12-20 23:00:00 2019-12-21 23:00:00 male 0
#> 7 2239 3 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 female 0
#> 8 2240 3 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 male 0
#> 9 2244 3 pink 2019-12-20 23:00:00 2019-12-21 23:00:00 male 0
#> 10 2245 4 wild 2019-12-20 23:00:00 2019-12-21 23:00:00 female 1
#> # ℹ 53 more rows
Dispersal of sterile males in space and time.
## Retrieve age since release date of each survey
sm_ages <- survey_ages(
surveys = sit_adult_surveys(sit_prototype),
releases = sit_revents(sit_prototype, type = 'point')
)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
## Summarise counts by trap and age-group (by 5 days)
sm_summary <-
sm_ages %>%
mutate(
age_group = cut(age, 5*(0:6))
) %>%
filter(!is.na(age_group)) %>% # remove ages 30+ (max is 36)
group_by(trap_code, age_group) %>%
summarise(n = sum(n), .groups = 'drop')
## Join this data to the geospatial information of the traps
sm_sf <- sit_traps(sit_prototype, area = "sit", stage = "adult") %>%
left_join(sm_summary, by = c(code = 'trap_code'))
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
tm_shape(sm_sf) +
tm_bubbles(size = 'n') +
tm_facets(by = 'age_group')
The competitiveness γ of the sterile male individuals is defined as their relative capacity to mate with a wild female, compared to a wild male.
Thus, in a homogeneously mixed population with Ms sterile males and Mw wild males, the probability that a mating occurs with a sterile individual is $$ p_s = \frac{\gamma M_s}{M_w + \gamma M_s} = \frac{\gamma R_{sw}}{1 + \gamma R_{sw}}, $$ where Rsw = Ms/Mw is the sterile-wild male ratio. At a given sterile-wild male ratio Rsw = Ms/Mw we observe a fertility rate Hs in the field.
Assuming a residual fertility rate Hrs for sterile males and a natural fertility rate Hw for wild males, the observed fertility rate Hs in the field is:
Solving the above equation for γ we obtain the Fried index, which estimates the competitiveness based on observed fertility rates and the sterile-wild male ratio: where Hw is the percentage egg hatch in the control area (i.e. the natural fertility); Hs is the percentage egg hatch in the release area (observed fertility under a sterile-wild male ratio Rsw) and Hrs is the residual fertility of sterile males.
This requires the estimation of:
The sterile-wild males ratio Rsw. Estimated from the observed ratios in adult surveys in the following 7 days after areal releases by default.
The natural fertility Hw. Estimated from egg-hatching rates in the control area.
The observed fertility in the release area Hs. Estimated from egg-hatching rates in the SIT area.
The residual fertility of sterile males Hrs. This depends on the level of radiation used for sterilisation, and the influx of females from beyond the SIT area. This cannot be estimated from SIT data and needs to be fixed by the user. Using Hrs = 0 reduces to the so-called CIS index.
sit_competitiveness(sit_prototype, residual_fertility = 0.05)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Sterile male competitiveness ────────────────────────────────────────────────
#> ℹ Estimated value: 0.6
#>
#>
#> Component Value
#> ------------------------------- ------
#> Sterile-wild male ratio 1.13
#> Natural fertility 0.84
#> Observed fertility in SIT area 0.52
#> Sterile fertility (assumed) 0.05
If there were more than one areal release, we could compute the
estimated competitiveness for each one separately by specifying
following_release
. We can also control the number of days
following the release to be used for the calculation of the sterile-wild
male ratio.
sit_competitiveness(
sit_prototype,
following_releases = sit_revents(sit_prototype, type = "areal")[1,],
following_days = 10,
residual_fertility = 0.05
)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Sterile male competitiveness ────────────────────────────────────────────────
#> ℹ Estimated value: 1.18
#>
#>
#> Component Value
#> ------------------------------- ------
#> Sterile-wild male ratio 0.58
#> Natural fertility 0.84
#> Observed fertility in SIT area 0.52
#> Sterile fertility (assumed) 0.05
We can be shocked by the difference in the estimate after simply widening the observation window to 10 days rather than 7, and legitimately ask which time-span is more appropriate.
Notice that the only component of the calculation that changed sensibly is the sterile-wild male ratio, which essentially halved.
Let’s see how the relationship between captured sterile and wild males evolve in time. We can look at the ratio, but the difference is a more natural scale.
areal_release <- sit_revents(sit_prototype, type = "areal")
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
areal_surveys <- sit_adult_surveys(
sit_prototype,
area = "sit",
following_releases = areal_release,
following_days = 10
)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
sterile_wild_male_ratio(areal_surveys) %>%
mutate(
age = as.numeric(datetime_end - areal_release$date)
) %>%
group_by(age) %>%
summarise(
sterile = sum(sterile),
wild = sum(wild),
swdiff = sterile - wild
) %>%
ggplot(aes(age, swdiff)) +
geom_point() +
geom_smooth(method = 'loess', formula = 'y ~ x') +
labs(x = "Age", y = "Sterile - wild difference")
Starting from day 8, the number of captured sterile males drops consistently with respect to the captured wild males, which explains the difference in the estimate of competitiveness.
Since the calculation pools observations assuming a constant ratio, we can conclude that 7 days was an appropriate time-span for these data.
Since the dispersion evolves in time, typically increasing with time since release, we compute the MDT at each day j and, potentially, for the population corresponding to the release event k.
where nt is the total number of traps, di is the distance from trap i to the release site of event k and nijk is the number of sterile males surveyed at trap i, in day j after release event k.
Furthermore, a set of weights wik for
each trap and release event can be used, in order to account for a
inhomogeneous arrangement of the traps, as explained below. This is
activated by the argument spatial_adjustment
in
sit_mdt()
and sit_flight_range()
. Without
spatial adjustment, the weights are set as wik = 1
and the MDT reduces to the so-called Mean Dispersal Distance (MDD).
Depending on the needs and hypotheses, the calculation could also aggregate observations by days, populations or both.
If most traps were concentrated far away from the release points, then MDT would be artificially inflated as a consequence. A way of correcting for this is to weight the counts by the inverse relative local trap density at each trap. However, a spatial density estimate of traps would be too noisy with only a couple dozen traps and suffer from considerable edge effects. Thus, we assume isotropy (neglect the direction) and consider only the radial density from the release site.
Trap density with respect to the distance from the release point increases linearly with the distance for a uniform distribution of traps.
The relative density is the ratio between the observed density and the expected density under uniformity (straight line). While the ωi is the inverse of the ratio at the distance of the corresponding trap (Figure below).
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
This approach also suffers from edge effect. However, the trick to solve it is much simpler and the results more robust. The figure above shows the uncorrected density in grey and the adjusted density as a black curve with points at the observed distances.
The only issue would be for trap patterns that are strongly concentrated in specific directions at different distances. However, this seems unlikely.
With this spatial adjustment (performed by default), we can assess the MDT by each released population (point releases only), and we can also plot the histogram of travelled distances weighted by the number of observations and the spatial adjustment:
(mdt_prototype <- sit_mdt(sit_prototype, by = "population"))
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Mean Distance Travelled ─────────────────────────────────────────────────────
#>
#>
#> population mdt
#> ----------- -------
#> blue 110.09
#> red 63.48
#> yellow 69.35
plot(mdt_prototype)
The histogram excludes the most extreme 2.5 % of observations on each side, to improve visualisation. The vertical lines display the mean distances travelled by population. They appear shifted towards 0 due to a large number of adjusted distances near 0 that are not shown.
Contrast with the results without the spatial adjustment. The estimated distances are shorter, due to a over-represented number of traps at shorter distances.
sit_mdt(sit_prototype, by = "population", spatial_adjustment = FALSE)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Mean Distance Travelled ─────────────────────────────────────────────────────
#>
#>
#> population mdt
#> ----------- ------
#> blue 66.97
#> red 52.44
#> yellow 58.67
We can also consider the MDT by age (pooling across populations).
sit_mdt(sit_prototype, by = "age")
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Mean Distance Travelled ─────────────────────────────────────────────────────
#>
#>
#> age mdt
#> ---- -------
#> 1 88.40
#> 2 51.75
#> 3 57.35
#> 4 68.11
#> 5 109.28
#> 6 79.75
#> 7 154.34
#> 8 79.37
#> 9 104.39
#> 10 128.51
#> 11 128.51
#> 12 37.38
#> 13 330.73
#> 14 128.51
#> 15 76.27
There is no obvious increase of MDT with time. It could be one particular population which is distorting the results, but we don’t have enough observations to split the MDT by age and population.
On the other hand, estimates get noisier for larger ages since observations are more scarce. If we restricted the analysis to the first 10 days of age, the pattern is more apparent.
sit_mdt(sit_prototype, by = "age") %>%
filter(age <= 10) %>%
ggplot(aes(age, mdt)) +
geom_point() +
geom_smooth(method = "lm", formula = 'y ~ x')
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
This is an alternative measure of the dispersal capacity of sterile males based on quantiles rather than averages.
The method regresses the log-distances (+1 to avoid divergence) from the release
point as a function of the cumulative (trap-density adjusted) catches.
This relationship is modelled linearly, and the expected distances at
the 50 and 90% adjusted catches with respect to the total are defined as
the corresponding flight ranges FR50
and
FR90
.
The function sit_flight_range()
uses the levels 50 and
90 by default and allows making a basic plot of the underlying model for
verification.
(fr_prototype <- sit_flight_range(sit_prototype))
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Flight range ────────────────────────────────────────────────────────────────
#>
#>
#> population level cum_adj_catch FR
#> ----------- ------ -------------- -------
#> blue 50 15.93 73.54
#> blue 90 28.67 240.58
#> red 50 11.63 51.89
#> red 90 20.94 165.90
#> yellow 50 9.14 64.86
#> yellow 90 16.45 172.20
Moreover, the plot
method returns the underlying data
invisibly, which allows making a more customised figure. E.g.
fr_plot %>%
ggplot(aes(cx, logd, colour = population)) +
geom_point(show.legend = FALSE) +
geom_smooth(method = "lm", formula = y ~ x, se = FALSE, show.legend = FALSE) +
labs(x = "Cumulative captures", y = "Log-distance (log-m)") +
scale_colour_manual(values = unique(fr_plot$population))
The function also allows defining any number of quantiles and gives results either disaggregated by population (default) or all populations pooled together.
sit_flight_range(sit_prototype, levels = c(25, 50, 75), pool = TRUE)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Flight range ────────────────────────────────────────────────────────────────
#>
#>
#> population level cum_adj_catch FR
#> ----------- ------ -------------- -------
#> pooled 25 18.35 29.31
#> pooled 50 36.70 61.46
#> pooled 75 55.04 127.70
Assuming that mosquitoes fly randomly (Brownian motion), their dispersion follow Fick’s first law of diffusion, which postulates that flux goes from regions of high concentration to regions of low concentration, with a magnitude that is proportional to the concentration gradient (spatial derivative). The proportionality constant is the diffusivity D and depend of the species’ behaviour for that location, season and weather conditions.
Under this model, the Mean Squared Displacement (MSD) of mosquitoes from their release point at time t is MSD(t) = 4Dt.
From which we can solve for the Diffusion coefficient by using the distances from the release point at which sterile males were captured to estimate the MSD at a given age.
The function sit_diffusion()
computes the MSD of the
individuals (by population, by default, or possibly pooling populations)
at the different observed ages, and derives D as the regression slope. The
resulting object can be plotted to see the underlying regression
model.
(diff_prototype <- sit_diffusion(sit_prototype))
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Diffusion coefficient ───────────────────────────────────────────────────────
#>
#>
#> population Dest
#> ----------- -------
#> blue 861.09
#> red 139.38
#> yellow 173.77
We assume that for a released population of sterile males there is a certain probability of daily survival (πPDS), which represents the probability that an individual (or the proportion of individuals that) makes it to the next day.
The proportion pj of sterile males still alive at day j decreases geometrically, starting from some initial survival p0 at day 0:
pj = p0 πPDSj.
Assuming that each day j we capture a fraction c of the released individuals that remain alive, we have that the expected number of surveyed individuals at day j is
E[yj] = c p0 πPDSj.
In log-scale, log (E[yj]) = α + β ⋅ j, where α = log (c p0) and β = log (πPDS).
This leads to estimating PDS by regressing the log-number of catches (+ 1, to avoid divergences) on the number of days since release (age).
The methodology was taken from Kay and Muir (1998), but originally developed in Gillies (1961). There is no appropriate statistical model other than the observation the the rate of decrease in the caught individuals is approximately constant, which leads to the regression approach.
The recapture θ and survival S rates are computed as
where α and β are the regression coefficients of the log-catches against age, N is the number of individuals released and d is the number of days after release (age in days). We use d = 1.
The function sit_survival()
returns a table with the
Probability of Daily Survival (PDS), the Average Life
Expectancy ALE
, Recapture Rate
RR
, Survival Rate SR
), by population
(by default, unless pool = TRUE
)
(surv_prototype <- sit_survival(sit_prototype))
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Survival ────────────────────────────────────────────────────────────────────
#>
#>
#> population N_released PDS ALE RRx1e3 SR
#> ----------- ----------- ----- ----- ------- -----
#> blue 10000 0.81 4.83 1.96 0.81
#> red 10000 0.81 4.83 1.76 0.82
#> yellow 10000 0.81 4.83 1.38 0.82
An estimation of the density, or equivalently, the size of the wild mosquito population is necessary for the determination of the appropriate quantity of sterile males to be released for controlling the population effectively.
A simple estimate is obtained as follows (Thompson 2012, Ch. 18). Let the total captures at a day t be the sum of mt marked and nt wild mosquitoes. Assuming that the proportion of marked individuals in the sample equals that in the population of size P: $$ \frac{m_t}{m_t + n_t} = \frac{M_t}{M_t + P}, $$ where Mt = R Sat is the number of marked individuals captured at time t, with R the number of released adults, S the daily survival rate and at the number of days since release (age). I.e., the number of marked individuals at time t is the remaining number from those released that survived for at days.
The Lincoln Index (a.k.a. the Petersen estimator) has been used as a simple estimate of the wild male population size, assuming that the survival rate of an individual remains constant. Here we use a modified version that corrects for small samples and compensates for daily survival. Pt = R Sat (nt + 1)/(mt + 1).
The values of R, nt, mt and t can be gathered from the adult surveys data. The calculation requires the estimation of the survival rate S.
This provides multiple estimates of the population size, one for each age t and for each released population, which can be averaged to compute a final estimate based on all the available data.
The calculation needs to be split by released population. First, because survival rates could be different among populations, but also because the number nt of captured wild males at a given date can be compared against the number of captured sterile males mt at different ages. Still, the function allows pooling data across population as if they were the same.
The function returns a point estimate and the standard deviation.
sit_wild_size(sit_prototype)
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#> Upgrading dm object created with dm <= 1.0.3.
#>
#> ── Wild male population size ───────────────────────────────────────────────────
#>
#>
#> population p_est p_sd
#> ----------- -------- --------
#> blue 4981.57 5420.85
#> red 3471.37 1919.16
#> yellow 6527.88 4248.72