# GenEst - Command Line Walkthrough

#### 2019-07-17

This vignette walks through an example of GenEst at the command line and was constructed using GenEst version 1.3.0 on 2019-07-17.

## Data

For this vignette, we will be using a completely generic, mock dataset provided with the GenEst package, which contains Searcher Efficiency (SE), Carcass Persistence (CP), Search Schedule (SS), Density Weighted Proportion (DWP), and Carcass Observation (CO) Data.

data(mock)
names(mock)
#> [1] "SE"  "CP"  "SS"  "DWP" "CO"

## Searcher Efficiency

### Single Searcher Efficiency Model

The central function for searcher efficiency analyses is pkm, which, in its most basic form, conducts a singular searcher efficiency analysis (i.e., a singular set of $$p$$ and $$k$$ formulae and a singular size classification of carcasses). As a first example, we will ignore the size category and use intercept-only models for both $$p$$ and $$k$$:

data_SE <- mock$SE pkModel <- pkm(formula_p = p ~ 1, formula_k = k ~ 1, data = data_SE) Here, we have taken advantage of pkm’s default behavior of selecting observation columns (see ?pkm for details). head(data_SE) #> seID Visibility HabitatType Season Size Search1 Search2 Search3 Search4 #> 1 se1 L HT1 SF S 1 NA NA NA #> 2 se2 L HT1 SF S 1 NA NA NA #> 3 se3 L HT1 WS S 0 0 0 0 #> 4 se4 L HT1 WS S 1 NA NA NA #> 5 se5 L HT1 WS S 0 1 NA NA #> 6 se6 L HT1 SF S 1 NA NA NA If we wanted to explicitly control the observations, we would use the obsCol argument: pkModel <- pkm(formula_p = p ~ 1, formula_k = k ~ 1, data = data_SE, obsCol = c("Search1", "Search2", "Search3", "Search4") ) Note that the search observations must be entered in order such that no carcasses have non-detected observations (i.e., 0) after detected observations (i.e., 1). Further, no carcasses can be detected more than once. If successfully fit, a pkm model output contains a number of elements, some printed automatically: pkModel #>$call
#> pkm0(formula_p = formula_p, formula_k = formula_k, data = data,
#>     obsCol = obsCol, kFixed = kFixed, kInit = kInit, CL = CL,
#>     quiet = quiet)
#>
#> $formula_p #> p ~ 1 #> #>$formula_k
#> k ~ 1
#>
#> $predictors #> character(0) #> #>$AICc
#> [1] 1145
#>
#> $convergence #> [1] 0 #> #>$cell_pk
#>   cell   n p_median p_lower p_upper k_median k_lower k_upper
#> 1  all 480    0.568   0.532   0.604    0.599   0.543   0.653
#>
#> $CL #> [1] 0.9 and others available upon request: names(pkModel) #> [1] "call" "data" "data0" "formula_p" #> [5] "formula_k" "predictors" "predictors_p" "predictors_k" #> [9] "AIC" "AICc" "convergence" "varbeta" #> [13] "cellMM_p" "cellMM_k" "nbeta_p" "nbeta_k" #> [17] "betahat_p" "betahat_k" "cells" "ncell" #> [21] "cell_pk" "CL" "observations" "carcCells" #> [25] "loglik" "pOnly" pkModel$cells
#>   group CellNames
#> 1   all       all

The plot function has been defined for pkm objects, such that one can simply run

plot(pkModel)

to visualize the model’s output.

You can generate random draws of the $$p$$ and $$k$$ parameters for each cell grouping (in pkModel there are no predictors, so there is one cell grouping called “all”) using the rpk function which, like other r* functions in R (e.g., rnorm, runif) takes the number of random draws (n) as the first argument:

rpk(n = 10, pkModel)
#> $all #> p k #> [1,] 0.5488693 0.6278329 #> [2,] 0.5833221 0.4930200 #> [3,] 0.6028768 0.5782386 #> [4,] 0.5712496 0.5812995 #> [5,] 0.5718182 0.5464391 #> [6,] 0.5747276 0.6408535 #> [7,] 0.6265282 0.5724079 #> [8,] 0.5929173 0.6148904 #> [9,] 0.5800293 0.6137182 #> [10,] 0.5681824 0.6018989 You can complicate the $$p$$ and $$k$$ formulae independently pkm(formula_p = p ~ Visibility, formula_k = k ~ HabitatType, data = data_SE, obsCol = c("Search1", "Search2", "Search3", "Search4") ) #>$call
#> pkm0(formula_p = formula_p, formula_k = formula_k, data = data,
#>     obsCol = obsCol, kFixed = kFixed, kInit = kInit, CL = CL,
#>     quiet = quiet)
#>
#> $formula_p #> p ~ Visibility #> #>$formula_k
#> k ~ HabitatType
#>
#> $predictors #> [1] "Visibility" "HabitatType" #> #>$AICc
#> [1] 1149.67
#>
#> $convergence #> [1] 0 #> #>$cell_pk
#>    cell  n p_median p_lower p_upper k_median k_lower k_upper
#> 1 H.HT1 80    0.564   0.503   0.623    0.560   0.482   0.636
#> 2 L.HT1 80    0.578   0.516   0.637    0.560   0.482   0.636
#> 3 M.HT1 80    0.563   0.505   0.620    0.560   0.482   0.636
#> 4 H.HT2 80    0.564   0.503   0.623    0.631   0.557   0.700
#> 5 L.HT2 80    0.578   0.516   0.637    0.631   0.557   0.700
#> 6 M.HT2 80    0.563   0.505   0.620    0.631   0.557   0.700
#>
#> $CL #> [1] 0.9 And you can fix $$k$$ at a nominal value between 0 and 1 (inclusive) using the kFixed argument pkm(formula_p = p ~ Visibility, kFixed = 0.7, data = data_SE, obsCol = c("Search1", "Search2", "Search3", "Search4")) #>$call
#> pkm0(formula_p = formula_p, formula_k = formula_k, data = data,
#>     obsCol = obsCol, kFixed = kFixed, kInit = kInit, CL = CL,
#>     quiet = quiet)
#>
#> $formula_p #> p ~ Visibility #> #>$formula_k
#> fixedk
#>    0.7
#>
#> $predictors #> [1] "Visibility" #> #>$AICc
#> [1] 1155.63
#>
#> $convergence #> [1] 0 #> #>$cell_pk
#>   cell   n p_median p_lower p_upper k_median k_lower k_upper
#> 1    H 160    0.531   0.474   0.588      0.7     0.7     0.7
#> 2    L 160    0.545   0.487   0.601      0.7     0.7     0.7
#> 3    M 160    0.538   0.482   0.593      0.7     0.7     0.7
#>
#> $CL #> [1] 0.9 ### Set of Searcher Efficiency Models If the arg allCombos = TRUE is provided, pkm fits a set of pkm models defined as all allowable models simpler than, and including, the provided model for both formulae (where “allowable” means that any interaction terms have all component terms included in the model). Consider the following model set analysis, where visibility and habitat type are included in the $$p$$ formula but only habitat type is in the $$k$$ formula. This generates a set of 10 models: pkmModSet <- pkm(formula_p = p ~ Visibility*HabitatType, formula_k = k ~ HabitatType, data = data_SE, obsCol = c("Search1", "Search2", "Search3", "Search4"), allCombos = TRUE ) class(pkmModSet) #> [1] "pkmSet" "list" names(pkmModSet) #> [1] "p ~ Visibility * HabitatType; k ~ HabitatType" #> [2] "p ~ Visibility + HabitatType; k ~ HabitatType" #> [3] "p ~ HabitatType; k ~ HabitatType" #> [4] "p ~ Visibility; k ~ HabitatType" #> [5] "p ~ 1; k ~ HabitatType" #> [6] "p ~ Visibility * HabitatType; k ~ 1" #> [7] "p ~ Visibility + HabitatType; k ~ 1" #> [8] "p ~ HabitatType; k ~ 1" #> [9] "p ~ Visibility; k ~ 1" #> [10] "p ~ 1; k ~ 1" The plot function is defined for the pkmSet class, and by default, creates a new plot window on command for each sub-model. If we want to only plot a specific single (or subset) of models from the full set, we can utilize the specificModel argument: plot(pkmModSet, specificModel = "p ~ Visibility + HabitatType; k ~ 1") The resulting model outputs can be compared in an AICc table aicc(pkmModSet) #> p Formula k Formula AICc <U+0394>AICc #> 10 p ~ 1 k ~ 1 1145.00 0.00 #> 5 p ~ 1 k ~ HabitatType 1145.70 0.70 #> 3 p ~ HabitatType k ~ HabitatType 1146.57 1.57 #> 8 p ~ HabitatType k ~ 1 1146.76 1.76 #> 9 p ~ Visibility k ~ 1 1148.96 3.96 #> 4 p ~ Visibility k ~ HabitatType 1149.67 4.67 #> 2 p ~ Visibility + HabitatType k ~ HabitatType 1150.55 5.55 #> 7 p ~ Visibility + HabitatType k ~ 1 1150.73 5.73 #> 1 p ~ Visibility * HabitatType k ~ HabitatType 1153.45 8.45 #> 6 p ~ Visibility * HabitatType k ~ 1 1153.49 8.49 ### Multiple Sizes of Animals and Sets of Searcher Efficiency Models Often, carcasses are grouped in multiple size classes, and we are interested in analyzing a set of models separately for each size class. To do so, we use the sizeCol arg to tell pkm which column in data_CP gives the carcass size class. If, in addition, allCombos = TRUE, pkm will fit a pkmSet that runs for each unique size class in the column identified by the sizeCol argument: pkmModSetSize <- pkm(formula_p = p ~ Visibility*HabitatType, formula_k = k ~ HabitatType, data = data_SE, obsCol = c("Search1", "Search2", "Search3", "Search4"), sizeCol = "Size", allCombos = TRUE) class(pkmModSetSize) #> [1] "pkmSetSize" "list" The pkmSetSize object is a list where each element corresponds to a different unique size class, and contains the associated pkmSetobject, which itself is a list of pkm outputs: names(pkmModSetSize) #> [1] "L" "M" "S" "XL" names(pkmModSetSize[[1]]) #> [1] "p ~ Visibility * HabitatType; k ~ HabitatType" #> [2] "p ~ Visibility + HabitatType; k ~ HabitatType" #> [3] "p ~ HabitatType; k ~ HabitatType" #> [4] "p ~ Visibility; k ~ HabitatType" #> [5] "p ~ 1; k ~ HabitatType" #> [6] "p ~ Visibility * HabitatType; k ~ 1" #> [7] "p ~ Visibility + HabitatType; k ~ 1" #> [8] "p ~ HabitatType; k ~ 1" #> [9] "p ~ Visibility; k ~ 1" #> [10] "p ~ 1; k ~ 1" ## Carcass Persistence ### Single Carcass Persistence Model The central function for carcass persistence analyses is cpm, which, in its simplest form, conducts a singular carcass persistence analysis (i.e., a singular set of $$l$$ and $$s$$ formulae and a singular size classification of carcasses). Note that we use $$l$$ and $$s$$ to reference $$location$$ and $$scale$$ as the parameters for survival models, following survreg, however we also provide an alternative parameterization (using parameters $$a$$ and $$b$$, referred to as “ab” or “ppersist”). As a first example, we will ignore the size category, use intercept-only models for both $$l$$ and $$s$$, and use the Weibull distribution: data_CP <- mock$CP
cpModel <- cpm(formula_l = l ~ 1, formula_s = s ~ 1, data = data_CP,
left = "LastPresentDecimalDays",
right = "FirstAbsentDecimalDays", dist = "weibull"
)

If successfully fit, a cpm model output contains a number of elements, some printed automatically:

cpModel
#> $call #> cpm0(formula_l = formula_l, formula_s = formula_s, data = data, #> left = left, right = right, dist = dist, CL = CL, quiet = quiet) #> #>$formula_l
#> l ~ 1
#>
#> $formula_s #> s ~ 1 #> #>$distribution
#> [1] "weibull"
#>
#> $predictors #> character(0) #> #>$AICc
#> [1] 2102.11
#>
#> $convergence #> [1] 0 #> #>$cell_ls
#>   cell   n l_median l_lower l_upper s_median s_lower s_upper
#> 1  all 480    2.671   2.592   2.749    0.966   0.901   1.037
#>
#> $cell_ab #> cell n pda_median pda_lower pda_upper pdb_median pdb_lower pdb_upper #> 1 all 480 1.035 1.11 0.964 14.454 13.356 15.627 #> #>$CL
#> [1] 0.9
#>
#> $cell_desc #> cell medianCP r1 r3 r7 r14 r28 #> 1 all 10.1437 0.9696732 0.9094574 0.8003887 0.6463498 0.4427757 and others available upon request: names(cpModel) #> [1] "call" "data" "formula_l" "formula_s" #> [5] "distribution" "predictors" "predictors_l" "predictors_s" #> [9] "AIC" "AICc" "convergence" "varbeta" #> [13] "cellMM_l" "cellMM_s" "nbeta_l" "nbeta_s" #> [17] "betahat_l" "betahat_s" "cells" "ncell" #> [21] "cell_ls" "cell_ab" "CL" "observations" #> [25] "carcCells" "loglik" "cell_desc" cpModel$cells
#>   group CellNames
#> 1   all       all

The plot function has been defined for cpm objects, such that one can simply run

plot(cpModel)

to visualize the model’s output.

You can generate random draws of the $$l$$ and $$s$$ (or $$a$$ and $$b$$) parameters for each cell grouping (in cpModel there are no predictors, so there is one cell grouping called “all”) using the rcp function which, like other r* functions in R (e.g., rnorm) takes the number of random draws (n) as the first argument:

rcp(n = 10, cpModel)
#> $all #> l s #> [1,] 2.614048 1.0117197 #> [2,] 2.657523 0.9668281 #> [3,] 2.617342 0.9494785 #> [4,] 2.622189 0.9330116 #> [5,] 2.658662 1.0266295 #> [6,] 2.623847 0.9310500 #> [7,] 2.758454 0.9853990 #> [8,] 2.727480 0.9377525 #> [9,] 2.644048 0.9857435 #> [10,] 2.651328 0.9221541 rcp(n = 10, cpModel, type = "ppersist") #>$all
#>             pda      pdb
#>  [1,] 0.9864405 12.70573
#>  [2,] 1.0349598 14.97980
#>  [3,] 1.0226685 15.44520
#>  [4,] 1.0402250 14.03074
#>  [5,] 1.0296836 14.24607
#>  [6,] 1.0748220 14.67369
#>  [7,] 0.9962308 13.78008
#>  [8,] 1.0394273 15.13397
#>  [9,] 1.0982423 14.86640
#> [10,] 0.9900430 14.55487

You can complicate the $$l$$ and $$s$$ formulae independently

cpm(formula_l = l ~ Visibility * GroundCover, formula_s = s ~ 1, data = data_CP,
left = "LastPresentDecimalDays", right = "FirstAbsentDecimalDays",
dist = "weibull"
)
#> $call #> cpm0(formula_l = formula_l, formula_s = formula_s, data = data, #> left = left, right = right, dist = dist, CL = CL, quiet = quiet) #> #>$formula_l
#> l ~ Visibility * GroundCover
#>
#> $formula_s #> s ~ 1 #> #>$distribution
#> [1] "weibull"
#>
#> $predictors #> [1] "Visibility" "GroundCover" #> #>$AICc
#> [1] 2109.36
#>
#> $convergence #> [1] 0 #> #>$cell_ls
#>   cell  n l_median l_lower l_upper s_median s_lower s_upper
#> 1  H.A 80    2.647   2.457   2.836    0.964   0.898   1.034
#> 2  L.A 80    2.618   2.427   2.809    0.964   0.898   1.034
#> 3  M.A 80    2.536   2.350   2.722    0.964   0.898   1.034
#> 4  H.B 80    2.789   2.597   2.981    0.964   0.898   1.034
#> 5  L.B 80    2.710   2.521   2.900    0.964   0.898   1.034
#> 6  M.B 80    2.716   2.527   2.905    0.964   0.898   1.034
#>
#> $cell_ab #> cell n pda_median pda_lower pda_upper pdb_median pdb_lower pdb_upper #> 1 H.A 80 1.037 1.114 0.967 14.112 11.670 17.047 #> 2 L.A 80 1.037 1.114 0.967 13.708 11.325 16.593 #> 3 M.A 80 1.037 1.114 0.967 12.629 10.486 15.211 #> 4 H.B 80 1.037 1.114 0.967 16.265 13.423 19.708 #> 5 L.B 80 1.037 1.114 0.967 15.029 12.441 18.174 #> 6 M.B 80 1.037 1.114 0.967 15.120 12.516 18.265 #> #>$CL
#> [1] 0.9
#>
#> $cell_desc #> cell medianCP r1 r3 r7 r14 r28 #> 1 H.A 9.910449 0.9691191 0.9076888 0.7965531 0.6402617 0.4355206 #> 2 L.A 9.626732 0.9681953 0.9050529 0.7912752 0.6323916 0.4266794 #> 3 M.A 8.868981 0.9654396 0.8972323 0.7757825 0.6096922 0.4019134 #> 4 H.B 11.422439 0.9732702 0.9196215 0.8208024 0.6773283 0.4789728 #> 5 L.B 10.554432 0.9710322 0.9131703 0.8076196 0.6569919 0.4547593 #> 6 M.B 10.618339 0.9712095 0.9136797 0.8086542 0.6585719 0.4566077 Given that the exponential only has one parameter ($$l$$, location), a model for scale (formula_s) is not required: cpModExp <- cpm(formula_l = l ~ Visibility * GroundCover, data = data_CP, left = "LastPresentDecimalDays", right = "FirstAbsentDecimalDays", dist = "exponential" ) ### Set of Carcass Persistence Models If the arg allCombos = TRUE is provided, cpm fits a set of cpm models defined as all allowable models simpler than, and including, the provided model formulae (where “allowable” means that any interaction terms have all component terms included in the model). In addition, cpm with allCombos can include any subset of the four base distributions (exponential, weibull, lognormal, loglogistic) and crosses them with the predictor models. Consider the following model set analysis, where Visibility and Season are included in the $$l$$ formula but only Visibility is in the $$s$$ formula, and only the exponential and lognormal distributions are included. This generates a set of 15 models: cpmModSet <- cpm(formula_l = l ~ Visibility * Season, formula_s = s ~ Visibility, data = data_CP, left = "LastPresentDecimalDays", right = "FirstAbsentDecimalDays", dist = c("exponential", "lognormal"), allCombos = TRUE ) class(cpmModSet) #> [1] "cpmSet" "list" names(cpmModSet) #> [1] "dist: exponential; l ~ Visibility * Season; NULL" #> [2] "dist: exponential; l ~ Visibility + Season; NULL" #> [3] "dist: exponential; l ~ Season; NULL" #> [4] "dist: exponential; l ~ Visibility; NULL" #> [5] "dist: exponential; l ~ 1; NULL" #> [6] "dist: lognormal; l ~ Visibility * Season; s ~ Visibility" #> [7] "dist: lognormal; l ~ Visibility + Season; s ~ Visibility" #> [8] "dist: lognormal; l ~ Season; s ~ Visibility" #> [9] "dist: lognormal; l ~ Visibility; s ~ Visibility" #> [10] "dist: lognormal; l ~ 1; s ~ Visibility" #> [11] "dist: lognormal; l ~ Visibility * Season; s ~ 1" #> [12] "dist: lognormal; l ~ Visibility + Season; s ~ 1" #> [13] "dist: lognormal; l ~ Season; s ~ 1" #> [14] "dist: lognormal; l ~ Visibility; s ~ 1" #> [15] "dist: lognormal; l ~ 1; s ~ 1" The resulting model outputs can be compared in an AICc table aicc(cpmModSet) #> Distribution Location Formula Scale Formula AICc <U+0394>AICc #> 5 exponential l ~ 1 NULL 2100.72 0.00 #> 3 exponential l ~ Season NULL 2101.08 0.36 #> 4 exponential l ~ Visibility NULL 2104.16 3.44 #> 2 exponential l ~ Visibility + Season NULL 2104.48 3.76 #> 1 exponential l ~ Visibility * Season NULL 2108.17 7.45 #> 15 lognormal l ~ 1 s ~ 1 2159.24 58.52 #> 13 lognormal l ~ Season s ~ 1 2160.78 60.06 #> 10 lognormal l ~ 1 s ~ Visibility 2162.15 61.43 #> 14 lognormal l ~ Visibility s ~ 1 2163.19 62.47 #> 8 lognormal l ~ Season s ~ Visibility 2163.67 62.95 #> 12 lognormal l ~ Visibility + Season s ~ 1 2164.74 64.02 #> 9 lognormal l ~ Visibility s ~ Visibility 2166.13 65.41 #> 11 lognormal l ~ Visibility * Season s ~ 1 2166.91 66.19 #> 7 lognormal l ~ Visibility + Season s ~ Visibility 2167.66 66.94 #> 6 lognormal l ~ Visibility * Season s ~ Visibility 2169.79 69.07 The plot function is defined for the cpmSet class, and by default, creates a new plot window on command for each sub-model. If we want to only plot a specific single (or subset) of models from the full set, we can utilize the specificModel argument: plot(cpmModSet, specificModel = "dist: lognormal; l ~ Visibility * Season; s ~ Visibility" ) ### Multiple Sizes of Animals and Sets of Carcass Persistence Models Often, carcasses are grouped in multiple size classes, and we are interested in analyzing a set of models separately for each size class. To do so, we furnish cpm with sizeCol, which is the name of the column in data_CP that gives the size classes of the carcasses. If, in addition, allCombos = TRUE, then cpm returns a cpmSet for each unique size class in the column identified by the sizeCol argument: cpmModSetSize <- cpm(formula_l = l ~ Visibility * Season, formula_s = s ~ Visibility, data = data_CP, left = "LastPresentDecimalDays", right = "FirstAbsentDecimalDays", dist = c("exponential", "lognormal"), sizeCol = "Size", allCombos = TRUE) class(cpmModSetSize) #> [1] "cpmSetSize" "list" The cpmSetSize object is a list where each element corresponds to a different unique size class, and contains the associated cpmSeto bject, which itself is a list of cpm outputs: names(cpmModSetSize) #> [1] "L" "M" "S" "XL" names(cpmModSetSize[[1]]) #> [1] "dist: exponential; l ~ Visibility * Season; NULL" #> [2] "dist: exponential; l ~ Visibility + Season; NULL" #> [3] "dist: exponential; l ~ Season; NULL" #> [4] "dist: exponential; l ~ Visibility; NULL" #> [5] "dist: exponential; l ~ 1; NULL" #> [6] "dist: lognormal; l ~ Visibility * Season; s ~ Visibility" #> [7] "dist: lognormal; l ~ Visibility + Season; s ~ Visibility" #> [8] "dist: lognormal; l ~ Season; s ~ Visibility" #> [9] "dist: lognormal; l ~ Visibility; s ~ Visibility" #> [10] "dist: lognormal; l ~ 1; s ~ Visibility" #> [11] "dist: lognormal; l ~ Visibility * Season; s ~ 1" #> [12] "dist: lognormal; l ~ Visibility + Season; s ~ 1" #> [13] "dist: lognormal; l ~ Season; s ~ 1" #> [14] "dist: lognormal; l ~ Visibility; s ~ 1" #> [15] "dist: lognormal; l ~ 1; s ~ 1" class(cpmModSetSize[[1]]) #> [1] "cpmSet" "list" ## Generic Detection Probability For the purposes of mortality estimation, we calculate carcass-specific detection probabilities (see below), which may be difficult to generalize, given the specific history of each observed carcass. Thus, we also provide a simple means to calculate generic detection probabilities that are cell-specific, rather than carcass-specific. For any estimation of detection probability ($$\hat{g}$$), we need to have singular SE and CP models to use for each of the size classes. Here, we use the best-fit of the models for each size class: pkMods <- c("S" = "p ~ 1; k ~ 1", "L" = "p ~ 1; k ~ 1", "M" = "p ~ 1; k ~ 1", "XL" = "p ~ 1; k ~ HabitatType" ) cpMods <- c("S" = "dist: exponential; l ~ Season; NULL", "L" = "dist: exponential; l ~ 1; NULL", "M" = "dist: exponential; l ~ 1; NULL", "XL" = "dist: exponential; l ~ 1; NULL" ) The estgGenericSize function produces n random draws of generic (i.e., cell-specific, not carcass-sepecific) detection probabilities for each of the possible carcass cell combinations across the selected SE and CP models across the size classes. estgGeneric is a single-size-class version of function and estgGenericSize actually loops over estgGeneric. The generic $$\hat{g}$$ is estimated according to a particular search schedule. When we pass averageSS a full data_SS table like we have here, it will assume that columns filled exclusively with 0s and 1s represent search schedules for units and will create the average search schedule across the units. data_SS <- mock$SS
avgSS <- averageSS(data_SS)

gsGeneric <- estgGenericSize(nsim = 1000, days = avgSS,
modelSetSize_SE = pkmModSetSize,
modelSetSize_CP = cpmModSetSize,
modelSizeSelections_SE = pkMods,
modelSizeSelections_CP = cpMods
)

The output from estgGeneric can be simply summarized

summary(gsGeneric)
#> $L #> Group 5% 25% 50% 75% 95% #> 1 all 0.379 0.41 0.433 0.453 0.482 #> #>$M
#>   Group   5%   25%   50%   75%  95%
#> 1   all 0.35 0.383 0.405 0.429 0.46
#>
#> $S #> Season 5% 25% 50% 75% 95% #> 1 SF 0.355 0.397 0.425 0.452 0.491 #> 2 WS 0.297 0.333 0.359 0.382 0.421 #> #>$XL
#>   HabitatType    5%   25%   50%   75%   95%
#> 1         HT1 0.316 0.349 0.370 0.392 0.424
#> 2         HT2 0.328 0.362 0.384 0.407 0.443

or plotted.

plot(gsGeneric)

## Mortality Estimation

When estimating mortality, detection probability is determined for individual carcasses based on the dates when they are observed, size class values, associated covariates, the searcher efficiency and carcass persistence models, and the search schedule. The carcass-specific detection probabilities (as opposed to the generic/cell-specific detection probabilities above) are therefore calculated before estimating the total mortality. Although it is possible to estimate these detection probabilities separately, they are best interpreted in the context of a full mortality estimation.

The estM function is the general wrapper function for estimating M, whether for a single size class or multiple size classes. Prior to estimation, we need to reduce the model-set-size complexed to just a single chosen model per size class, corresponding to the pkMods and cpMods vectors given above. To reduce the model set complexity, we can use the trimSetSize function:

pkmModSize <- trimSetSize(pkmModSetSize, pkMods)
cpmModSize <- trimSetSize(cpmModSetSize, cpMods)

In addition to the models and search schedule data, estM requires density-weighted proportion (DWP) and carcass observation (CO) data. If more than one size class is represented in the data, a required input is also the column names associated with the DWP value for each size class (argument DWPCol in estM):

data_CO <- mock$CO data_DWP <- mock$DWP
#>    Unit    S    M    L   XL
#> 1 Unit1 0.70 0.70 0.60 0.60
#> 2 Unit2 0.70 0.70 0.60 0.60
#> 3 Unit3 0.56 0.56 0.48 0.48
#> 4 Unit4 0.56 0.56 0.48 0.48
#> 5 Unit5 0.70 0.70 0.60 0.60
DWPcolnames <- names(pkmModSize)

eM <- estM(data_CO = data_CO, data_SS = data_SS, data_DWP = data_DWP,
frac = 1, model_SE = pkmModSize, model_CP = cpmModSize,
seed_SE = NULL, seed_CP = NULL, seed_g = NULL, seed_M = NULL,
unitCol = "Unit", COdate = "DateFound",
SSdate = "DateSearched", sizeCol = "Size", nsim = 1000)

estM returns an object that contains the random draws of pkm and cpm parameters (named pk and ab, respectively) and the estimated carcass-level detection parameters (g), arrival intervals (Aj), and associated total mortality (Mhat) values for each simulation. These Mhat values should be considered in combination, and can be summarized and plotted simply:

summary(eM)
#>  median      5%     95%
#> 1804.35 1635.32 1982.09
plot(eM)

### Splitting Mortality Estimations

It is possible to split the resulting mortality estimation into components that are denoted according to covariates in either the search schedule or carcass observation data sets.

First, a temporal split:

M_season <- calcSplits(M = eM, split_SS = "Construction",
split_CO = NULL, data_SS = data_SS, data_CO = data_CO
)
summary(M_season)
#>               X        5%       25%      50%       75%      95%
#> Before 161.5319  590.7384  638.9584  672.168  708.0378  774.293
#> After  264.4681 1012.8606 1078.1314 1129.559 1178.0742 1261.744
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "Construction"
#> attr(,"type")
#> [1] "SS"
#> attr(,"times")
#> [1]   0 127 349
plot(M_season)

Next, a carcass split:

M_class <- calcSplits(M = eM, split_SS = NULL,
split_CO = "Split", data_SS = data_SS, data_CO = data_CO
)
summary(M_class)
#>      X       5%      25%      50%       75%       95%
#> C1 196 724.6863 781.9022 823.8686  868.8898  925.8283
#> C2 230 874.2922 929.9648 977.6811 1024.1010 1095.6859
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "Split"
#> attr(,"type")
#> [1] "CO"
plot(M_class)

And finally, if two splits are included, the mortality estimation is expanded fully factorially:

M_SbyC <- calcSplits(M = eM, split_SS = "Construction",
split_CO = "Split", data_SS = data_SS, data_CO = data_CO
)
summary(M_SbyC)
#> $C1 #> X 5% 25% 50% 75% 95% #> Before 80.24343 274.4449 308.7490 331.9646 354.9848 394.7901 #> After 115.75657 418.3971 463.8586 493.5380 520.2003 571.0686 #> #>$C2
#>               X       5%      25%      50%      75%      95%
#> Before  81.2885 288.2849 319.3385 341.2645 364.3408 400.0452
#> After  148.7115 551.5801 601.2314 635.1906 671.1242 726.5226
#>
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "Construction" "Split"
#> attr(,"type")
#> [1] "SS" "CO"
#> attr(,"times")
#> [1]   0 127 349
plot(M_SbyC)