Chapter 13 Special Analyses II: Quasi-Extinction Analysis
A learning experience is one of those things that says, ‘You know that thing you just did? Don’t do that.’
Conservation managers often need to conduct population viability analyses (PVAs), which are demographic projection studies designed to assess the persistence potential of a target population (Beissinger & Westphal 1998). The populations targeted are typically rare and potentially endangered organisms. A quasi-extinction analysis is a type of PVA in which a population is split into control and treatment groups, and both groups are projected with temporal stochasticity. The persistence of the populations across replicates is then compared. This comparison is meant to overcome the limitations of predicting extinction itself, which is notoriously difficult. Instead, quasi-extinction analyses produce relative extinction probabilities, usually expressed as ratios of the proportion of extinct treatment replicates to the proportion of extinct control replicates.
In the simplest type of quasi-extinction analysis, we might have demographic data for a full population that has not been treated experimentally. Alternatively, we might have a population in which some portion of the population has been treated differently than the rest for some number of years, with demographic data collected from both treatment and control groups. In the former case, we might parameterize control and treatment sets of MPMs that use different underlying assumptions. For example, we might assume in the treatment group that we have boosted the pollination rate, or we have increased the survival of established adults via increased surveillance and the use of fences. In the latter case, we have data that stems from the treatment and so do not need to make any different assumptions, but simply parameterize two sets of models corresponding to the control and treatment groups.
Once the models are parameterized, we may project the matrices forward by some amount of time that we view as reasonable for management purposes. Depending on the agency and the conservation plan, this might be anywhere from 5 to 100 years. If we conduct these as stochastic projections replicated 100 or 1000 times, then we can compare the proportion of these replicates that go extinct in both cases. Since we do not know the true extinction rate, we instead estimate the quasi-extinction probabilities for both groups for the time window involved, and what is most meaningful is their ratio. For example, a treatment strategy that lowers the quasi-extinction probability by 20% might be worth incorporating into a management plan.
13.1 Cypripedium candidum example
Here we will set up a quasi-extinction analysis. We will use this analysis to assess whether hand pollination can decrease the extinction risk of an endangered perennial herb, Cypripedium candidum. Let’s utilize the function-based approach, using a life history model previously seen in earlier chapters (see section 5.2 in chapter 5). We’ll start off loading the data.
Next we’ll build the stageframe. Since we will split the population into two groups later, we will assume the same underlying biology and so the same life history model.
stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "V1", "V2", "V3", "V4",
"V5", "V6", "V7", "V8", "V9", "V10", "V11", "V12", "V13", "V14", "V15", "V16",
"V17", "V18", "V19", "V20", "V21", "V22", "V23", "V24", "F1", "F2", "F3",
"F4", "F5", "F6", "F7", "F8", "F9", "F10", "F11", "F12", "F13", "F14", "F15",
"F16", "F17", "F18", "F19", "F20", "F21", "F22", "F23", "F24")
indataset <- c(0, 0, 0, 0, 0, rep(1, 49))
sizevector <- c(0, 0, 0, 0, 0, seq(from = 0, t = 24), seq(from = 1, to = 24))
repvector <- c(0, 0, 0, 0, 0, rep(0, 25), rep(1, 24))
obsvector <- c(0, 0, 0, 0, 0, 0, rep(1, 48))
matvector <- c(0, 0, 0, 0, 0, rep(1, 49))
immvector <- c(0, 1, 1, 1, 1, rep(0, 49))
propvector <- c(1, rep(0, 53))
comments <- c("Dormant seed", "Yr1 protocorm", "Yr2 protocorm", "Yr3 protocorm",
"Seedling", "Veg dorm", "Veg adult 1 stem", "Veg adult 2 stems",
"Veg adult 3 stems", "Veg adult 4 stems", "Veg adult 5 stems",
"Veg adult 6 stems", "Veg adult 7 stems", "Veg adult 8 stems",
"Veg adult 9 stems", "Veg adult 10 stems", "Veg adult 11 stems",
"Veg adult 12 stems", "Veg adult 13 stems", "Veg adult 14 stems",
"Veg adult 15 stems", "Veg adult 16 stems", "Veg adult 17 stems",
"Veg adult 18 stems", "Veg adult 19 stems", "Veg adult 20 stems",
"Veg adult 21 stems", "Veg adult 22 stems", "Veg adult 23 stems",
"Veg adult 24 stems", "Flo adult 1 stem", "Flo adult 2 stems",
"Flo adult 3 stems", "Flo adult 4 stems", "Flo adult 5 stems",
"Flo adult 6 stems", "Flo adult 7 stems", "Flo adult 8 stems",
"Flo adult 9 stems", "Flo adult 10 stems", "Flo adult 11 stems",
"Flo adult 12 stems", "Flo adult 13 stems", "Flo adult 14 stems",
"Flo adult 15 stems", "Flo adult 16 stems", "Flo adult 17 stems",
"Flo adult 18 stems", "Flo adult 19 stems", "Flo adult 20 stems",
"Flo adult 21 stems", "Flo adult 22 stems", "Flo adult 23 stems",
"Flo adult 24 stems")
cypframe_fb <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
propstatus = propvector, immstatus = immvector, indataset = indataset,
comments = comments)
cypframe_fb
> stage size size_b size_c min_age max_age repstatus obsstatus propstatus
> 1 SD 0 NA NA NA NA 0 0 1
> 2 P1 0 NA NA NA NA 0 0 0
> 3 P2 0 NA NA NA NA 0 0 0
> 4 P3 0 NA NA NA NA 0 0 0
> 5 SL 0 NA NA NA NA 0 0 0
> 6 D 0 NA NA NA NA 0 0 0
> 7 V1 1 NA NA NA NA 0 1 0
> 8 V2 2 NA NA NA NA 0 1 0
> 9 V3 3 NA NA NA NA 0 1 0
> 10 V4 4 NA NA NA NA 0 1 0
> 11 V5 5 NA NA NA NA 0 1 0
> 12 V6 6 NA NA NA NA 0 1 0
> 13 V7 7 NA NA NA NA 0 1 0
> 14 V8 8 NA NA NA NA 0 1 0
> 15 V9 9 NA NA NA NA 0 1 0
> 16 V10 10 NA NA NA NA 0 1 0
> 17 V11 11 NA NA NA NA 0 1 0
> 18 V12 12 NA NA NA NA 0 1 0
> 19 V13 13 NA NA NA NA 0 1 0
> 20 V14 14 NA NA NA NA 0 1 0
> 21 V15 15 NA NA NA NA 0 1 0
> 22 V16 16 NA NA NA NA 0 1 0
> 23 V17 17 NA NA NA NA 0 1 0
> 24 V18 18 NA NA NA NA 0 1 0
> 25 V19 19 NA NA NA NA 0 1 0
> 26 V20 20 NA NA NA NA 0 1 0
> 27 V21 21 NA NA NA NA 0 1 0
> 28 V22 22 NA NA NA NA 0 1 0
> 29 V23 23 NA NA NA NA 0 1 0
> 30 V24 24 NA NA NA NA 0 1 0
> 31 F1 1 NA NA NA NA 1 1 0
> 32 F2 2 NA NA NA NA 1 1 0
> 33 F3 3 NA NA NA NA 1 1 0
> 34 F4 4 NA NA NA NA 1 1 0
> 35 F5 5 NA NA NA NA 1 1 0
> 36 F6 6 NA NA NA NA 1 1 0
> 37 F7 7 NA NA NA NA 1 1 0
> 38 F8 8 NA NA NA NA 1 1 0
> 39 F9 9 NA NA NA NA 1 1 0
> 40 F10 10 NA NA NA NA 1 1 0
> 41 F11 11 NA NA NA NA 1 1 0
> 42 F12 12 NA NA NA NA 1 1 0
> 43 F13 13 NA NA NA NA 1 1 0
> 44 F14 14 NA NA NA NA 1 1 0
> 45 F15 15 NA NA NA NA 1 1 0
> 46 F16 16 NA NA NA NA 1 1 0
> 47 F17 17 NA NA NA NA 1 1 0
> 48 F18 18 NA NA NA NA 1 1 0
> 49 F19 19 NA NA NA NA 1 1 0
> 50 F20 20 NA NA NA NA 1 1 0
> 51 F21 21 NA NA NA NA 1 1 0
> 52 F22 22 NA NA NA NA 1 1 0
> 53 F23 23 NA NA NA NA 1 1 0
> 54 F24 24 NA NA NA NA 1 1 0
> immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max
> 1 0 0 0 0.5 -0.5 0.5
> 2 1 0 0 0.5 -0.5 0.5
> 3 1 0 0 0.5 -0.5 0.5
> 4 1 0 0 0.5 -0.5 0.5
> 5 1 0 0 0.5 -0.5 0.5
> 6 0 1 1 0.5 -0.5 0.5
> 7 0 1 1 0.5 0.5 1.5
> 8 0 1 1 0.5 1.5 2.5
> 9 0 1 1 0.5 2.5 3.5
> 10 0 1 1 0.5 3.5 4.5
> 11 0 1 1 0.5 4.5 5.5
> 12 0 1 1 0.5 5.5 6.5
> 13 0 1 1 0.5 6.5 7.5
> 14 0 1 1 0.5 7.5 8.5
> 15 0 1 1 0.5 8.5 9.5
> 16 0 1 1 0.5 9.5 10.5
> 17 0 1 1 0.5 10.5 11.5
> 18 0 1 1 0.5 11.5 12.5
> 19 0 1 1 0.5 12.5 13.5
> 20 0 1 1 0.5 13.5 14.5
> 21 0 1 1 0.5 14.5 15.5
> 22 0 1 1 0.5 15.5 16.5
> 23 0 1 1 0.5 16.5 17.5
> 24 0 1 1 0.5 17.5 18.5
> 25 0 1 1 0.5 18.5 19.5
> 26 0 1 1 0.5 19.5 20.5
> 27 0 1 1 0.5 20.5 21.5
> 28 0 1 1 0.5 21.5 22.5
> 29 0 1 1 0.5 22.5 23.5
> 30 0 1 1 0.5 23.5 24.5
> 31 0 1 1 0.5 0.5 1.5
> 32 0 1 1 0.5 1.5 2.5
> 33 0 1 1 0.5 2.5 3.5
> 34 0 1 1 0.5 3.5 4.5
> 35 0 1 1 0.5 4.5 5.5
> 36 0 1 1 0.5 5.5 6.5
> 37 0 1 1 0.5 6.5 7.5
> 38 0 1 1 0.5 7.5 8.5
> 39 0 1 1 0.5 8.5 9.5
> 40 0 1 1 0.5 9.5 10.5
> 41 0 1 1 0.5 10.5 11.5
> 42 0 1 1 0.5 11.5 12.5
> 43 0 1 1 0.5 12.5 13.5
> 44 0 1 1 0.5 13.5 14.5
> 45 0 1 1 0.5 14.5 15.5
> 46 0 1 1 0.5 15.5 16.5
> 47 0 1 1 0.5 16.5 17.5
> 48 0 1 1 0.5 17.5 18.5
> 49 0 1 1 0.5 18.5 19.5
> 50 0 1 1 0.5 19.5 20.5
> 51 0 1 1 0.5 20.5 21.5
> 52 0 1 1 0.5 21.5 22.5
> 53 0 1 1 0.5 22.5 23.5
> 54 0 1 1 0.5 23.5 24.5
> sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max
> 1 0 1 NA NA NA
> 2 0 1 NA NA NA
> 3 0 1 NA NA NA
> 4 0 1 NA NA NA
> 5 0 1 NA NA NA
> 6 0 1 NA NA NA
> 7 1 1 NA NA NA
> 8 2 1 NA NA NA
> 9 3 1 NA NA NA
> 10 4 1 NA NA NA
> 11 5 1 NA NA NA
> 12 6 1 NA NA NA
> 13 7 1 NA NA NA
> 14 8 1 NA NA NA
> 15 9 1 NA NA NA
> 16 10 1 NA NA NA
> 17 11 1 NA NA NA
> 18 12 1 NA NA NA
> 19 13 1 NA NA NA
> 20 14 1 NA NA NA
> 21 15 1 NA NA NA
> 22 16 1 NA NA NA
> 23 17 1 NA NA NA
> 24 18 1 NA NA NA
> 25 19 1 NA NA NA
> 26 20 1 NA NA NA
> 27 21 1 NA NA NA
> 28 22 1 NA NA NA
> 29 23 1 NA NA NA
> 30 24 1 NA NA NA
> 31 1 1 NA NA NA
> 32 2 1 NA NA NA
> 33 3 1 NA NA NA
> 34 4 1 NA NA NA
> 35 5 1 NA NA NA
> 36 6 1 NA NA NA
> 37 7 1 NA NA NA
> 38 8 1 NA NA NA
> 39 9 1 NA NA NA
> 40 10 1 NA NA NA
> 41 11 1 NA NA NA
> 42 12 1 NA NA NA
> 43 13 1 NA NA NA
> 44 14 1 NA NA NA
> 45 15 1 NA NA NA
> 46 16 1 NA NA NA
> 47 17 1 NA NA NA
> 48 18 1 NA NA NA
> 49 19 1 NA NA NA
> 50 20 1 NA NA NA
> 51 21 1 NA NA NA
> 52 22 1 NA NA NA
> 53 23 1 NA NA NA
> 54 24 1 NA NA NA
> sizebinb_center sizebinb_width binhalfwidthc_raw sizebinc_min sizebinc_max
> 1 NA NA NA NA NA
> 2 NA NA NA NA NA
> 3 NA NA NA NA NA
> 4 NA NA NA NA NA
> 5 NA NA NA NA NA
> 6 NA NA NA NA NA
> 7 NA NA NA NA NA
> 8 NA NA NA NA NA
> 9 NA NA NA NA NA
> 10 NA NA NA NA NA
> 11 NA NA NA NA NA
> 12 NA NA NA NA NA
> 13 NA NA NA NA NA
> 14 NA NA NA NA NA
> 15 NA NA NA NA NA
> 16 NA NA NA NA NA
> 17 NA NA NA NA NA
> 18 NA NA NA NA NA
> 19 NA NA NA NA NA
> 20 NA NA NA NA NA
> 21 NA NA NA NA NA
> 22 NA NA NA NA NA
> 23 NA NA NA NA NA
> 24 NA NA NA NA NA
> 25 NA NA NA NA NA
> 26 NA NA NA NA NA
> 27 NA NA NA NA NA
> 28 NA NA NA NA NA
> 29 NA NA NA NA NA
> 30 NA NA NA NA NA
> 31 NA NA NA NA NA
> 32 NA NA NA NA NA
> 33 NA NA NA NA NA
> 34 NA NA NA NA NA
> 35 NA NA NA NA NA
> 36 NA NA NA NA NA
> 37 NA NA NA NA NA
> 38 NA NA NA NA NA
> 39 NA NA NA NA NA
> 40 NA NA NA NA NA
> 41 NA NA NA NA NA
> 42 NA NA NA NA NA
> 43 NA NA NA NA NA
> 44 NA NA NA NA NA
> 45 NA NA NA NA NA
> 46 NA NA NA NA NA
> 47 NA NA NA NA NA
> 48 NA NA NA NA NA
> 49 NA NA NA NA NA
> 50 NA NA NA NA NA
> 51 NA NA NA NA NA
> 52 NA NA NA NA NA
> 53 NA NA NA NA NA
> 54 NA NA NA NA NA
> sizebinc_center sizebinc_width group comments
> 1 NA NA 0 Dormant seed
> 2 NA NA 0 Yr1 protocorm
> 3 NA NA 0 Yr2 protocorm
> 4 NA NA 0 Yr3 protocorm
> 5 NA NA 0 Seedling
> 6 NA NA 0 Veg dorm
> 7 NA NA 0 Veg adult 1 stem
> 8 NA NA 0 Veg adult 2 stems
> 9 NA NA 0 Veg adult 3 stems
> 10 NA NA 0 Veg adult 4 stems
> 11 NA NA 0 Veg adult 5 stems
> 12 NA NA 0 Veg adult 6 stems
> 13 NA NA 0 Veg adult 7 stems
> 14 NA NA 0 Veg adult 8 stems
> 15 NA NA 0 Veg adult 9 stems
> 16 NA NA 0 Veg adult 10 stems
> 17 NA NA 0 Veg adult 11 stems
> 18 NA NA 0 Veg adult 12 stems
> 19 NA NA 0 Veg adult 13 stems
> 20 NA NA 0 Veg adult 14 stems
> 21 NA NA 0 Veg adult 15 stems
> 22 NA NA 0 Veg adult 16 stems
> 23 NA NA 0 Veg adult 17 stems
> 24 NA NA 0 Veg adult 18 stems
> 25 NA NA 0 Veg adult 19 stems
> 26 NA NA 0 Veg adult 20 stems
> 27 NA NA 0 Veg adult 21 stems
> 28 NA NA 0 Veg adult 22 stems
> 29 NA NA 0 Veg adult 23 stems
> 30 NA NA 0 Veg adult 24 stems
> 31 NA NA 0 Flo adult 1 stem
> 32 NA NA 0 Flo adult 2 stems
> 33 NA NA 0 Flo adult 3 stems
> 34 NA NA 0 Flo adult 4 stems
> 35 NA NA 0 Flo adult 5 stems
> 36 NA NA 0 Flo adult 6 stems
> 37 NA NA 0 Flo adult 7 stems
> 38 NA NA 0 Flo adult 8 stems
> 39 NA NA 0 Flo adult 9 stems
> 40 NA NA 0 Flo adult 10 stems
> 41 NA NA 0 Flo adult 11 stems
> 42 NA NA 0 Flo adult 12 stems
> 43 NA NA 0 Flo adult 13 stems
> 44 NA NA 0 Flo adult 14 stems
> 45 NA NA 0 Flo adult 15 stems
> 46 NA NA 0 Flo adult 16 stems
> 47 NA NA 0 Flo adult 17 stems
> 48 NA NA 0 Flo adult 18 stems
> 49 NA NA 0 Flo adult 19 stems
> 50 NA NA 0 Flo adult 20 stems
> 51 NA NA 0 Flo adult 21 stems
> 52 NA NA 0 Flo adult 22 stems
> 53 NA NA 0 Flo adult 23 stems
> 54 NA NA 0 Flo adult 24 stemsContinuing on, we’ll now standardize the dataset.
cypfb <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
patchidcol = "patch", individcol = "plantid", blocksize = 4,
sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
stageassign = cypframe_fb, stagesize = "sizeadded", NAas0 = TRUE,
age_offset = 4)
summary_hfv(cypfb)
>
> This hfv dataset contains 320 rows, 57 variables, 1 population,
> 3 patches, 74 individuals, and 5 time steps.
> rowid popid patchid individ year2
> Min. : 1.00 Length:320 A: 93 Min. : 164.0 Min. :2004
> 1st Qu.:21.00 Class :character B:154 1st Qu.: 391.0 1st Qu.:2005
> Median :37.50 Mode :character C: 73 Median : 453.0 Median :2006
> Mean :38.45 Mean : 651.5 Mean :2006
> 3rd Qu.:56.00 3rd Qu.: 476.0 3rd Qu.:2007
> Max. :77.00 Max. :1560.0 Max. :2008
> firstseen lastseen obsage obslifespan
> Min. :2004 Min. :2004 Min. :5.000 Min. :0.000
> 1st Qu.:2004 1st Qu.:2009 1st Qu.:6.000 1st Qu.:5.000
> Median :2004 Median :2009 Median :7.000 Median :5.000
> Mean :2004 Mean :2009 Mean :6.853 Mean :4.556
> 3rd Qu.:2004 3rd Qu.:2009 3rd Qu.:8.000 3rd Qu.:5.000
> Max. :2008 Max. :2009 Max. :9.000 Max. :5.000
> sizea1 sizeb1 sizec1 size1added
> Min. :0.000000 Min. : 0.0000 Min. : 0.0 Min. : 0.000
> 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.: 0.0 1st Qu.: 0.000
> Median :0.000000 Median : 0.0000 Median : 1.0 Median : 2.000
> Mean :0.009375 Mean : 0.7469 Mean : 1.9 Mean : 2.656
> 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.: 3.0 3rd Qu.: 4.000
> Max. :1.000000 Max. :18.0000 Max. :13.0 Max. :21.000
> repstra1 repstrb1 repstr1added feca1
> Min. : 0.0000 Min. :0.000000 Min. : 0.0000 Min. :0.0000
> 1st Qu.: 0.0000 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.:0.0000
> Median : 0.0000 Median :0.000000 Median : 0.0000 Median :0.0000
> Mean : 0.7469 Mean :0.009375 Mean : 0.7562 Mean :0.2656
> 3rd Qu.: 1.0000 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.:0.0000
> Max. :18.0000 Max. :1.000000 Max. :18.0000 Max. :7.0000
> fec1added obsstatus1 repstatus1 fecstatus1
> Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
> Median :0.0000 Median :1.0000 Median :0.0000 Median :0.0000
> Mean :0.2656 Mean :0.7469 Mean :0.2875 Mean :0.1344
> 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000
> Max. :7.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
> matstatus1 alive1 stage1 stage1index
> Min. :0.0000 Min. :0.0000 Length:320 Min. : 0.00
> 1st Qu.:1.0000 1st Qu.:1.0000 Class :character 1st Qu.: 6.00
> Median :1.0000 Median :1.0000 Mode :character Median : 8.00
> Mean :0.7688 Mean :0.7688 Mean :14.17
> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:31.00
> Max. :1.0000 Max. :1.0000 Max. :51.00
> sizea2 sizeb2 sizec2 size2added
> Min. :0.000000 Min. : 0.0000 Min. : 0.000 Min. : 0.000
> 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.: 1.000
> Median :0.000000 Median : 0.0000 Median : 2.000 Median : 2.000
> Mean :0.009375 Mean : 0.8969 Mean : 2.416 Mean : 3.322
> 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.: 3.000 3rd Qu.: 4.000
> Max. :1.000000 Max. :18.0000 Max. :13.000 Max. :24.000
> repstra2 repstrb2 repstr2added feca2
> Min. : 0.0000 Min. :0.000000 Min. : 0.0000 Min. :0.0000
> 1st Qu.: 0.0000 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.:0.0000
> Median : 0.0000 Median :0.000000 Median : 0.0000 Median :0.0000
> Mean : 0.8969 Mean :0.009375 Mean : 0.9062 Mean :0.2906
> 3rd Qu.: 1.0000 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.:0.0000
> Max. :18.0000 Max. :1.000000 Max. :18.0000 Max. :7.0000
> fec2added obsstatus2 repstatus2 fecstatus2
> Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.0000
> Median :0.0000 Median :1.0000 Median :0.0000 Median :0.0000
> Mean :0.2906 Mean :0.9531 Mean :0.3688 Mean :0.1562
> 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000
> Max. :7.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
> matstatus2 alive2 stage2 stage2index sizea3
> Min. :1 Min. :1 Length:320 Min. : 6.00 Min. :0.000000
> 1st Qu.:1 1st Qu.:1 Class :character 1st Qu.: 7.00 1st Qu.:0.000000
> Median :1 Median :1 Mode :character Median :10.00 Median :0.000000
> Mean :1 Mean :1 Mean :18.17 Mean :0.009375
> 3rd Qu.:1 3rd Qu.:1 3rd Qu.:32.00 3rd Qu.:0.000000
> Max. :1 Max. :1 Max. :54.00 Max. :1.000000
> sizeb3 sizec3 size3added repstra3
> Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
> 1st Qu.: 0.000 1st Qu.: 1.000 1st Qu.: 1.000 1st Qu.: 0.000
> Median : 0.000 Median : 1.000 Median : 2.000 Median : 0.000
> Mean : 1.069 Mean : 2.209 Mean : 3.288 Mean : 1.069
> 3rd Qu.: 1.000 3rd Qu.: 3.000 3rd Qu.: 4.000 3rd Qu.: 1.000
> Max. :18.000 Max. :13.000 Max. :24.000 Max. :18.000
> repstrb3 repstr3added feca3 fec3added
> Min. :0.000000 Min. : 0.000 Min. :0.0000 Min. :0.0000
> 1st Qu.:0.000000 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:0.0000
> Median :0.000000 Median : 0.000 Median :0.0000 Median :0.0000
> Mean :0.009375 Mean : 1.078 Mean :0.4562 Mean :0.4562
> 3rd Qu.:0.000000 3rd Qu.: 1.000 3rd Qu.:0.0000 3rd Qu.:0.0000
> Max. :1.000000 Max. :18.000 Max. :8.0000 Max. :8.0000
> obsstatus3 repstatus3 fecstatus3 matstatus3 alive3
> Min. :0.0 Min. :0.0 Min. :0.0000 Min. :1 Min. :0.0000
> 1st Qu.:1.0 1st Qu.:0.0 1st Qu.:0.0000 1st Qu.:1 1st Qu.:1.0000
> Median :1.0 Median :0.0 Median :0.0000 Median :1 Median :1.0000
> Mean :0.9 Mean :0.4 Mean :0.2219 Mean :1 Mean :0.9469
> 3rd Qu.:1.0 3rd Qu.:1.0 3rd Qu.:0.0000 3rd Qu.:1 3rd Qu.:1.0000
> Max. :1.0 Max. :1.0 Max. :1.0000 Max. :1 Max. :1.0000
> stage3 stage3index
> Length:320 Min. : 0.00
> Class :character 1st Qu.: 7.00
> Mode :character Median :10.00
> Mean :18.57
> 3rd Qu.:33.00
> Max. :54.00Next let’s explore the variables that we will build our vital rate models from.
hfv_qc(data = cypfb, vitalrates = c("surv", "obs", "size", "repst", "fec"),
size = c("size3added", "size2added", "size1added"))
> Survival:
>
> Data subset has 58 variables and 320 transitions.
>
> Variable alive3 has 0 missing values.
> Variable alive3 is a binomial variable.
>
> Numbers of categories in data subset in possible random variables:
> indiv id: 74
> year2: 5
>
> A total of 6 individuals have only single transitions.
>
>
> Observation status:
>
> Data subset has 58 variables and 303 transitions.
>
> Variable obsstatus3 has 0 missing values.
> Variable obsstatus3 is a binomial variable.
>
> Numbers of categories in data subset in possible random variables:
> indiv id: 70
> year2: 5
>
> A total of 6 individuals have only single transitions.
>
>
> Primary size:
>
> Data subset has 58 variables and 288 transitions.
>
> Variable size3added has 0 missing values.
> Variable size3added appears to be an integer variable.
>
> Variable size3added is fully positive, lacking even 0s.
>
> Overdispersion test:
> Mean size3added is 3.653
> The variance in size3added is 13.41
> The probability of this dispersion level by chance assuming that
> the true mean size3added = variance in size3added,
> and an alternative hypothesis of overdispersion, is 3.721e-138
> Variable size3added is significantly overdispersed.
>
> Zero-inflation and truncation tests:
> Mean lambda in size3added is 0.02592
> The actual number of 0s in size3added is 0
> The expected number of 0s in size3added under the null hypothesis is 7.465
> The probability of this deviation in 0s from expectation by chance is 0.9964
> Variable size3added is not significantly zero-inflated.
>
> Variable size3added does not include 0s, suggesting that a zero-truncated distribution may be warranted.
>
> Numbers of categories in data subset in possible random variables:
> indiv id: 70
> year2: 5
>
> A total of 6 individuals have only single transitions.
>
>
> Reproductive status:
>
> Data subset has 58 variables and 288 transitions.
>
> Variable repstatus3 has 0 missing values.
> Variable repstatus3 is a binomial variable.
>
> Numbers of categories in data subset in possible random variables:
> indiv id: 70
> year2: 5
>
> A total of 6 individuals have only single transitions.
>
>
> Fecundity:
>
> Data subset has 58 variables and 118 transitions.
>
> Variable feca2 has 0 missing values.
> Variable feca2 appears to be an integer variable.
>
> Variable feca2 is fully non-negative.
>
> Overdispersion test:
> Mean feca2 is 0.7881
> The variance in feca2 is 1.536
> The probability of this dispersion level by chance assuming that
> the true mean feca2 = variance in feca2,
> and an alternative hypothesis of overdispersion, is 0.1193
> Dispersion level in feca2 matches expectation.
>
> Zero-inflation and truncation tests:
> Mean lambda in feca2 is 0.4547
> The actual number of 0s in feca2 is 68
> The expected number of 0s in feca2 under the null hypothesis is 53.65
> The probability of this deviation in 0s from expectation by chance is 5.904e-06
> Variable feca2 is significantly zero-inflated.
>
> Numbers of categories in data subset in possible random variables:
> indiv id: 51
> year2: 5
>
> A total of 21 individuals have only single transitions.We have described these models previously, and so encourage the reader to skip back to chapter 5 to read our interpretations.
Next let’s build the vital rate models. Note that we have not yet introduced any differences in the population. So, these models will be used for both treatment and control groups.
cypmodels2 <- modelsearch(cypfb, historical = FALSE, approach = "mixed",
vitalrates = c("surv", "obs", "size", "repst", "fec"), sizedist = "negbin",
size.trunc = TRUE, fecdist = "poisson", fec.zero = TRUE, suite = "full",
size = c("size3added", "size2added", "size1added"), quiet = "partial")
>
> Developing global model of survival probability...
>
> Global model of survival probability developed. Proceeding with model dredge...
>
> Developing global model of observation probability...
>
> Global model of observation probability developed. Proceeding with model dredge...
>
> Developing global model of primary size...
>
> Global model of primary size developed. Proceeding with model dredge...
>
> Developing global model of reproduction probability...
>
> Global model of reproduction probability developed. Proceeding with model dredge...
>
> Developing global model of fecudity...
>
> Global model of fecundity developed. Proceeding with model dredge...
>
> Finished selecting best-fit models.
summary(cypmodels2)
> This LefkoMod object includes 5 linear models.
> Best-fit model criterion used: aicc&k
>
>
>
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ size2added + (1 | year2) + (1 | individ)
> Data: surv.data
> AIC BIC logLik deviance df.resid
> 128.1324 143.2057 -60.0662 120.1324 316
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 1.198379
> year2 (Intercept) 0.008826
> Number of obs: 320, groups: individ, 74; year2, 5
> Fixed Effects:
> (Intercept) size2added
> 2.0352 0.6344
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: obsstatus3 ~ size2added + (1 | year2) + (1 | individ)
> Data: obs.data
> AIC BIC logLik deviance df.resid
> 118.2567 133.1117 -55.1284 110.2567 299
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 1.078e-05
> year2 (Intercept) 8.776e-01
> Number of obs: 303, groups: individ, 70; year2, 5
> Fixed Effects:
> (Intercept) size2added
> 2.4904 0.3134
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Size model:
> Formula: size3added ~ (1 | year2) + (1 | individ)
> Data: size.data
> AIC BIC logLik df.resid
> 1008.2666 1022.9184 -500.1333 284
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.1109
> individ (Intercept) 1.0561
>
> Number of obs: 288 / Conditional model: year2, 5; individ, 70
>
> Dispersion parameter for truncated_nbinom2 family (): 9.34e+09
>
> Fixed Effects:
>
> Conditional model:
> (Intercept)
> 0.5761
>
>
>
> Secondary size model:
> [1] 1
>
>
>
> Tertiary size model:
> [1] 1
>
>
>
> Reproductive status model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: repstatus3 ~ repstatus2 + size2added + (1 | year2) + (1 | individ)
> Data: repst.data
> AIC BIC logLik deviance df.resid
> 333.6176 351.9324 -161.8088 323.6176 283
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.1829
> year2 (Intercept) 0.6250
> Number of obs: 288, groups: individ, 70; year2, 5
> Fixed Effects:
> (Intercept) repstatus2 size2added
> -1.4631 1.6457 0.1715
>
>
>
> Fecundity model:
> Formula: feca2 ~ size2added + (1 | year2) + (1 | individ)
> Zero inflation: ~size2added + (1 | year2) + (1 | individ)
> Data: fec.data
> AIC BIC logLik df.resid
> 248.8609 271.0264 -116.4305 110
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.5760
> individ (Intercept) 0.1639
>
> Zero-inflation model:
> Groups Name Std.Dev.
> year2 (Intercept) 1.642e-06
> individ (Intercept) 3.089e-04
>
> Number of obs: 118 / Conditional model: year2, 5; individ, 51 / Zero-inflation model: year2, 5; individ, 51
>
> Fixed Effects:
>
> Conditional model:
> (Intercept) size2added
> -0.54014 0.06174
>
> Zero-inflation model:
> (Intercept) size2added
> 3.865 -1.574
>
>
> Juvenile survival model:
> [1] 1
>
>
>
> Juvenile observation model:
> [1] 1
>
>
>
> Juvenile size model:
> [1] 1
>
>
>
> Juvenile secondary size model:
> [1] 1
>
>
>
> Juvenile tertiary size model:
> [1] 1
>
>
>
> Juvenile reproduction model:
> [1] 1
>
>
>
> Juvenile maturity model:
> [1] 1
>
>
>
>
>
> Number of models in survival table: 5
>
> Number of models in observation table: 5
>
> Number of models in size table: 5
>
> Number of models in secondary size table: 1
>
> Number of models in tertiary size table: 1
>
> Number of models in reproduction status table: 5
>
> Number of models in fecundity table: 25
>
> Number of models in juvenile survival table: 1
>
> Number of models in juvenile observation table: 1
>
> Number of models in juvenile size table: 1
>
> Number of models in juvenile secondary size table: 1
>
> Number of models in juvenile tertiary size table: 1
>
> Number of models in juvenile reproduction table: 1
>
> Number of models in juvenile maturity table: 1
>
>
>
>
>
> General model parameter names (column 1), and
> specific names used in these models (column 2):
> parameter_names mainparams
> 1 time t year2
> 2 individual individ
> 3 patch patch
> 4 alive in time t+1 surv3
> 5 observed in time t+1 obs3
> 6 sizea in time t+1 size3
> 7 sizeb in time t+1 sizeb3
> 8 sizec in time t+1 sizec3
> 9 reproductive status in time t+1 repst3
> 10 fecundity in time t+1 fec3
> 11 fecundity in time t fec2
> 12 sizea in time t size2
> 13 sizea in time t-1 size1
> 14 sizeb in time t sizeb2
> 15 sizeb in time t-1 sizeb1
> 16 sizec in time t sizec2
> 17 sizec in time t-1 sizec1
> 18 reproductive status in time t repst2
> 19 reproductive status in time t-1 repst1
> 20 maturity status in time t+1 matst3
> 21 maturity status in time t matst2
> 22 age in time t age
> 23 density in time t density
> 24 individual covariate a in time t indcova2
> 25 individual covariate a in time t-1 indcova1
> 26 individual covariate b in time t indcovb2
> 27 individual covariate b in time t-1 indcovb1
> 28 individual covariate c in time t indcovc2
> 29 individual covariate c in time t-1 indcovc1
> 30 stage group in time t group2
> 31 stage group in time t-1 group1
> 32 annual covariate a in time t annucova2
> 33 annual covariate a in time t-1 annucova1
> 34 annual covariate b in time t annucovb2
> 35 annual covariate b in time t-1 annucovb1
> 36 annual covariate c in time t annucovc2
> 37 annual covariate c in time t-1 annucovc1
>
>
>
>
>
> Quality control:
>
> Survival model estimated with 74 individuals and 320 individual transitions.
> Survival model accuracy is 0.947.
> Observation status model estimated with 70 individuals and 303 individual transitions.
> Observation status model accuracy is 0.95.
> Primary size model estimated with 70 individuals and 288 individual transitions.
> Primary size model R-squared is 0.822.
> Secondary size model not estimated.
> Tertiary size model not estimated.
> Reproductive status model estimated with 70 individuals and 288 individual transitions.
> Reproductive status model accuracy is 0.715.
> Fecundity model estimated with 51 individuals and 118 individual transitions.
> Fecundity model R-squared is 0.562.
> Juvenile survival model not estimated.
> Juvenile observation status model not estimated.
> Juvenile primary size model not estimated.
> Juvenile secondary size model not estimated.
> Juvenile tertiary size model not estimated.
> Juvenile reproductive status model not estimated.
> Juvenile maturity status model not estimated.Now we will finally introduce some differences in the population. Here, we will build two supplement tables, one for the control and one for the treatment group. We will assume that our treatment group is a hand-pollination treatment, in which volunteers will go into the population and pollinate these plants. Our preliminary data suggests that we might be able to triple the fruiting rate from 16% to 48% using this approach.
First the control supplement.
seeds_per_fruit <- 15000 # Mean number of seeds per fruit
germ_c <- 0.15
sl_mult <- 0.7 # Survival to 1st year
cypsupp2_control <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL", "D",
"V1", "V2", "V3", "SD", "P1"),
stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "SL", "rep",
"rep"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA),
givenrate = c(0.08, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, sl_mult, sl_mult, sl_mult, sl_mult,
0.5 * seeds_per_fruit * germ_c, 0.5 * seeds_per_fruit * germ_c),
type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "R", "R"),
stageframe = cypframe_fb, historical = FALSE)
cypsupp2_control
> stage3 stage2 stage1 age2 eststage3 eststage2 eststage1 estage2 givenrate
> 1 SD SD <NA> NA <NA> <NA> <NA> NA 0.08
> 2 P1 SD <NA> NA <NA> <NA> <NA> NA 0.10
> 3 P2 P1 <NA> NA <NA> <NA> <NA> NA 0.10
> 4 P3 P2 <NA> NA <NA> <NA> <NA> NA 0.10
> 5 SL P3 <NA> NA <NA> <NA> <NA> NA 0.05
> 6 SL SL <NA> NA <NA> <NA> <NA> NA 0.05
> 7 D SL <NA> NA D D <NA> NA NA
> 8 V1 SL <NA> NA V1 D <NA> NA NA
> 9 V2 SL <NA> NA V2 D <NA> NA NA
> 10 V3 SL <NA> NA V3 D <NA> NA NA
> 11 SD rep <NA> NA <NA> <NA> <NA> NA NA
> 12 P1 rep <NA> NA <NA> <NA> <NA> NA NA
> multiplier convtype convtype_t12
> 1 NA 1 1
> 2 NA 1 1
> 3 NA 1 1
> 4 NA 1 1
> 5 NA 1 1
> 6 NA 1 1
> 7 0.7 1 1
> 8 0.7 1 1
> 9 0.7 1 1
> 10 0.7 1 1
> 11 1125.0 3 1
> 12 1125.0 3 1Now the treatment supplement.
germ_t <- 0.48
cypsupp2_trt <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL",
"D", "V1", "V2", "V3", "SD", "P1"),
stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "SL", "rep",
"rep"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA),
givenrate = c(0.08, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, sl_mult, sl_mult, sl_mult, sl_mult,
0.5 * seeds_per_fruit * germ_t, 0.5 * seeds_per_fruit * germ_t),
type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "R", "R"),
stageframe = cypframe_fb, historical = FALSE)
cypsupp2_trt
> stage3 stage2 stage1 age2 eststage3 eststage2 eststage1 estage2 givenrate
> 1 SD SD <NA> NA <NA> <NA> <NA> NA 0.08
> 2 P1 SD <NA> NA <NA> <NA> <NA> NA 0.10
> 3 P2 P1 <NA> NA <NA> <NA> <NA> NA 0.10
> 4 P3 P2 <NA> NA <NA> <NA> <NA> NA 0.10
> 5 SL P3 <NA> NA <NA> <NA> <NA> NA 0.05
> 6 SL SL <NA> NA <NA> <NA> <NA> NA 0.05
> 7 D SL <NA> NA D D <NA> NA NA
> 8 V1 SL <NA> NA V1 D <NA> NA NA
> 9 V2 SL <NA> NA V2 D <NA> NA NA
> 10 V3 SL <NA> NA V3 D <NA> NA NA
> 11 SD rep <NA> NA <NA> <NA> <NA> NA NA
> 12 P1 rep <NA> NA <NA> <NA> <NA> NA NA
> multiplier convtype convtype_t12
> 1 NA 1 1
> 2 NA 1 1
> 3 NA 1 1
> 4 NA 1 1
> 5 NA 1 1
> 6 NA 1 1
> 7 0.7 1 1
> 8 0.7 1 1
> 9 0.7 1 1
> 10 0.7 1 1
> 11 3600.0 3 1
> 12 3600.0 3 1Finally, let’s run two simulations. We are interested in developing a management plant that will allow persistence for at least 10 years. So, let’s try a replicated stochastic simulation projecting 10 time steps forward.
First the control.
set.seed(42)
trial2f_control <- f_projection3(data = cypfb, format = 3, stochastic = TRUE,
substoch = 2, stageframe = cypframe_fb, supplement = cypsupp2_control,
modelsuite = cypmodels2, times = 10, nreps = 1000, integeronly = TRUE)
> Warning: Option patch not set, so will set to first patch/population.
summary(trial2f_control, ext_time = TRUE)
>
> The input lefkoProj object covers 1 population-patches.
> It is a single projection including 10 projected steps per replicate, and 1000 replicates, respectively.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> $milepost_sums
> 1 1
> 0 1000
> 0.25 1000
> 0.5 1000
> 0.75 1000
> 1 167
>
> $extinction_times
> ext_reps ext_time
> 1 833 10.2473Next the treatment simulation.
set.seed(42)
trial2f_trt <- f_projection3(data = cypfb, format = 3, stochastic = TRUE,
substoch = 2, stageframe = cypframe_fb, supplement = cypsupp2_trt,
modelsuite = cypmodels2, times = 10, nreps = 1000, integeronly = TRUE)
> Warning: Option patch not set, so will set to first patch/population.
summary(trial2f_trt, ext_time = TRUE)
>
> The input lefkoProj object covers 1 population-patches.
> It is a single projection including 10 projected steps per replicate, and 1000 replicates, respectively.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> $milepost_sums
> 1 1
> 0 1000
> 0.25 1000
> 0.5 1000
> 0.75 1000
> 1 983
>
> $extinction_times
> ext_reps ext_time
> 1 17 11The summaries all suggest that the hand pollination treatment leads to increased persistence. Let’s plot the results.
This certainly looks better, but also perhaps a bit difficult to interpret. Our quasi-extinction analysis needs a bit more than this - a numeric metric showing us the difference in extinction probability. Let’s estimate the quasi-extinction probabilities.
control_sum <- summary(trial2f_control)
>
> The input lefkoProj object covers 1 population-patches.
> It is a single projection including 10 projected steps per replicate, and 1000 replicates, respectively.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
trt_sum <- summary(trial2f_trt)
>
> The input lefkoProj object covers 1 population-patches.
> It is a single projection including 10 projected steps per replicate, and 1000 replicates, respectively.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
control_ext <- 1 - (control_sum$milepost_sums[5] / control_sum$milepost_sums[1])
trt_ext <- 1 - (trt_sum$milepost_sums[5] / trt_sum$milepost_sums[1])
writeLines(paste("Control quasi-extinction prob: ", control_ext))
> Control quasi-extinction prob: 0.833
writeLines(paste("Treatment quasi-extinction prob: ", trt_ext))
> Treatment quasi-extinction prob: 0.017Finally, let’s estimate the quasi-extinction ratio. This is the extinction probability of the treatment relative to the control.
And so we see that our hand pollination treatment reduces the chance of extinction to only 2% of the control over the next 10 years. This would suggest that our management plan is probably worth putting into action.
Naturally, quasi-extinction analysis is not exact and strongly dependent on assumptions. It is also dependent on conditions staying as they are during the monitoring period. We can try other approaches, including adding density dependence, shifting climate, etc. We will leave that to the user to try.
13.2 Points to remember
- Population viability analyses (PVA) use projection analyses to assess the persistence of populations.
- Quasi-extinction analysis is used to assess the impact of a treatment on the chance of populatoon extinction, given a comparison to a control.
- Functions
projection3()andf_projection3()allow users the ability to run paired projections to estimate quasi-extinction rates.
References
Beissinger, S.R. & Westphal, M.I. (1998). On the use of demographic models of population viability in endangered species management. Journal of Wildlife Management, 62, 821–841.