```

This tutorial provides an introduction to the array of command line tools **GenEst** provides for estimating carcass arrivals and detection probabilities. The approach is to walk through analyses of realistic but fictitious data sets representing studies of bird and bat mortality at a wind power facility.

The general steps in the analysis are:

- Construct a model for Searcher Efficiency
- Construct a model for Carcass Persistance
- Estimate mortality
- Specify the type of summary desired (e.g., mortality by species and season)

Data required for a full analysis include results of searcher efficiency trials, results of carcass persistence trials, search schedules for all units searched, the search coverage within each unit, and results of periodic carcass surveys. More information about the kinds of data required can be found in the User Guide. For convenience, data required in this tutorial are available in R as packaged lists, which can readily be loaded as described below. Alternatively, the data may be downloaded into zipped .csv files from the GUI.

To perform the analyses illustrated in the tutorial, begin by starting R and loading GenEst from the command line: `library(GenEst)`

Searcher efficiency and carcass persistence would be expected to vary with carcass size (sparrow, eagle, bat), ground characteristics (road & pad, cleared field, vegetation type), season, etc. In this first example, we limit the analysis to one carcass size (`bat`

) and one ground visibility class (`RP`

= road and pad). A more complicated scenario is analyzed in example 2.

The required data is stored in `wind_RPbat`

, a list of data frames with results for searcher efficiency (`SE`

) and carcass persistence trials (`CP`

), search schedules for all turbines (`SS`

), the search coverage or density weighted proportion (`DWP`

) of area searched at each turbine (i.e., the fraction of carcasses expected to fall in the search plots), and carcass observation (`CO`

) data.

Load the full data set into R:

To streamline the notation, extract the data from the `wind_RPbat`

list into its components:

```
data_SE <- wind_RPbat$SE
data_CP <- wind_RPbat$CP
data_SS <- wind_RPbat$SS
data_DWP <- wind_RPbat$DWP
data_CO <- wind_RPbat$CO
```

`SE`

)Searcher efficiency trials were conducted on roads and pads, with a total of 60 fresh carcasses placed in the field over the course of the entire monitoring period, evenly divided among seasons (spring, summer, fall). Carcasses that were later discovered by search teams during the course of normal carcass surveys were removed from the field. Carcasses were left in the field for up to 5 searches after carcass placement.

Results of the SE field trials are stored in the `data_SE`

data frame:

```
head(data_SE)
#> pkID Season s1 s2 s3 s4 s5
#> 1 pk1 spring 0 NA NA NA NA
#> 2 pk2 spring 0 0 NA NA NA
#> 3 pk3 spring 1 NA NA NA NA
#> 4 pk4 spring 1 NA NA NA NA
#> 5 pk5 spring 0 NA NA NA NA
#> 6 pk6 spring 1 NA NA NA NA
```

Columns `s1, s2, ..., s5`

show the fate of carcass `pkID`

on the 1st, 2nd, … 5th searches after the carcass was placed. A 1 indicates that the carcass was discovered, a 0 indicates that the carcass was present but not discovered, and NA indicates that the carcass was no longer present for discovery or no search was conducted.

`CP`

)Carcass persistence trials were conducted on roads and pads, with a total of 60 fresh carcasses placed in the field over the course of the entire monitoring period, evenly divided among seasons (spring, summer, fall). Carcasses were checked approximately 1, 2, 3, 4, 7, 10, 14, 21, and 28 days after placement in the field. Exact times were entered as decimal fractions of days after placement.

Results of the SE field trials are stored in the `data_CP`

data frame:

```
head(data_CP)
#> cpID Season LastPresent FirstAbsent
#> 1 cp1 spring 0.00 0.93
#> 2 cp2 spring 0.98 1.97
#> 3 cp3 spring 0.00 1.01
#> 4 cp4 spring 13.99 21.13
#> 5 cp5 spring 0.00 1.17
#> 6 cp6 spring 20.95 27.97
```

Exact persistence times are not known, but a carcass that was present at one check and absent at the next check is assumed to have been removed at some point in the interval. The left endpoint of the interval was entered as `LastPresent`

and the right endpoint as `FirstAbsent`

. For carcasses that had not been scavenged by the end of the study, `LastPresent`

is the time of the last check and `FirstAbsent`

is `Inf`

. For carcasses whose removal time is known exactly (e.g., scavenging was recorded by camera), `LastPresent = FirstAbsent`

. The `Season`

column gives the season at the time the carcass was placed in the field.

`SS`

)Carcass searches were conducted on roads and pads within a 120 m radius from all 100 turbines at a fictitious wind power facility. Monitoring began on 1955-04-15 and continued through 1955-11-01. Searches spanned 3 seasons: spring, summer, fall. Search intervals varied by turbine and by time of year, ranging from daily searches at some turbines in the fall and searches once every 12 days in the spring at some other turbines. Search schedules for all turbines are stored in `data_SS`

, which is a data frame with a column for search dates (including all dates that any turbine was searched); a column of 0s and 1s for each turbine, indication whether it was searched on the given date; and zero or more optional columns giving additional information about the date (e.g., season).

```
head(data_SS[, 1:10])
#> SearchDate Season t1 t2 t3 t4 t5 t6 t7 t8
#> 1 1955-04-15 spring 0 0 0 1 0 0 0 1
#> 2 1955-04-18 spring 1 0 0 0 1 0 0 0
#> 3 1955-04-21 spring 0 1 0 0 0 1 0 0
#> 4 1955-04-24 spring 0 0 1 0 0 0 1 0
#> 5 1955-04-27 spring 0 0 0 1 0 0 0 1
#> 6 1955-04-30 spring 1 0 0 0 1 0 0 0
```

Note that we have only displayed a few of the turbine columns - there are 100 turbine columns altogether (t1, …, t100).

`DWP`

)The density-weighted proportion (`DWP`

) is the expected fraction of carcasses that fell in the searched area. Carcass density is not the same at all distances from a turbine, but typically rises over a short distance then decreases eventually to 0. Searches were conducted on roads and pads within a 120 m radius from all 100 turbines to provide sufficient data with which to model the change in density with distance and from this, accurately calculate the fraction of all carcasses that are expected to land on road and pad surrounding each turbine (density-weighted proportion or DWP). The configuration of the roads and pads differs among turbines, hence the DWP must be calculated for each turbine. DWPs for bats at each turbine are stored in `data_DWP`

, which is a data frame with a column for turbine name (note that turbine names also must be included among column names in `data_SS`

, which gives the search schedule at each turbine) and a column of DWP labeled `bat`

. In other studies where, for example, mortality of birds might be of interest, `DWP`

would be expected vary with carcass size or species (e.g., the spatial distributions of bats and large birds around a turbine would likely differ from one another). In that case, each carcass size class (`large`

, `medium`

, `small`

, `bat`

) would have its own column.

`CO`

)Information about each (non-trial) carcass observed during searches is stored in `CO_data`

which is a data frame with at least 4 columns: carcass ID, the turbine (or unit) at which it was found, the date it was found and its distance from the turbine center. In `data_CO`

we also have turbine type, species, and species group variables by which we will later summarize mortality estimates.

```
head(data_CO)
#> carcID Turbine TurbineType DateFound Species SpeciesGroup Distance
#> 1 x79 t86 Z 1955-05-03 BA bat1 7.1
#> 2 x37 t19 X 1955-05-06 BB bat1 5.3
#> 3 x175 t37 Y 1955-05-12 BB bat1 9.8
#> 4 x181 t91 Z 1955-05-12 BA bat1 2.7
#> 5 x203 t38 Y 1955-05-15 BA bat1 3.8
#> 6 x291 t7 X 1955-05-18 BB bat1 112.3
```

Searcher efficiency and carcass persistence parameters are estimated by fitting models using functions `pkm`

and `cpm`

(carcass persistence model) which are patterned after familar R functions such as `lm`

and `glm`

. The “pk” in `pkm`

refers to GenEst’s model of searcher efficiency which includes two parameters: `p`

, which is the initial searcher efficiency for carcasses on the first search after they have arrived, and `k`

, which is a parameter governing the decrease in searcher efficiency in later searches. In this relatively simple example, our SE and CP field trials were conducted for one carcass size (bat) on one type of terrain (roads and pads) in three seasons (spring, summer, fall). The only potential predictor variable we have is `Season`

, which is entered as a column in both `data_SE`

and `data_CP`

.

Searcher efficiency is the probability of detection of a carcass that is present in the searched area at the time of search. Searcher efficiency typically decreases with carcass age because older carcasses tend to become harder to find as they accumulate dust or debris, fall deeper into vegetation, get blown against objects or into holes, decay, or get partially scavenged. In addition, carcasses missed in one search tend to be more likely to be missed in subsequent searches because the relatively easy-to-find carcasses are preferentially removed in the first searches after carcass arrival, leaving mostly the harder-to-find carcasses available in subsequences searches. GenEst accounts for a non-constant searcher efficiency using two parameters, `p`

(searcher efficiency on the first search after carcass arrivals) and `k`

(proportional change in searcher efficiency with each successive search). The `k`

parameter can be estimated from field trials if carcasses that are not discovered in the first search after arrival are left in the field for possible discovery in later searches.

`p`

and `k`

may both depend on covariates such as season, visibility class, or carcass size, and GenEst allows for them to be modeled as functions of different covariate combinations.

```
model_SE <- pkm(p ~ Season, k ~ 1, data = data_SE)
class(model_SE)
#> [1] "pkm" "list"
model_SE
#> $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 ~ Season
#>
#> $formula_k
#> k ~ 1
#>
#> $predictors
#> [1] "Season"
#>
#> $AICc
#> [1] 104.42
#>
#> $convergence
#> [1] 0
#>
#> $cell_pk
#> cell n p_median p_lower p_upper k_median k_lower k_upper
#> 1 fall 20 0.618 0.431 0.776 0.921 0.177 0.998
#> 2 spring 20 0.440 0.280 0.613 0.921 0.177 0.998
#> 3 summer 20 0.619 0.437 0.773 0.921 0.177 0.998
#>
#> $CL
#> [1] 0.9
```

NOTE: The `pkm`

family of functions by default interprets columns with names that begin with and “s” or “S” and end with a number contain search results data (carcass found = 0, not found = 1). A user can override the auto-parsing by explicitly listing the names of the search data columns in a vector of character strings in the `obsCol`

argument.

The probability of a carcass persisting a given length of time without being removed by scavengers (or other factors) is modeled as a Weibull, lognormal, loglogistic, or exponential distribution. Like the `p`

and `k`

parameters for searcher efficiency, the location and scale parameters (Therneau 2015) of the persistence distribution may depend on covariates. GenEst allows for them to be modeled as separate functions of predictor combinations.

```
model_CP <- cpm(l ~ Season, s ~ Season, data = data_CP, dist = "weibull",
left = "LastPresent", right = "FirstAbsent")
class(model_CP)
#> [1] "cpm" "list"
model_CP
#> $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 ~ Season
#>
#> $formula_s
#> s ~ Season
#>
#> $distribution
#> [1] "weibull"
#>
#> $predictors
#> [1] "Season"
#>
#> $AICc
#> [1] 235.51
#>
#> $convergence
#> [1] 0
#>
#> $cell_ls
#> cell n l_median l_lower l_upper s_median s_lower s_upper
#> 1 fall 20 1.163 0.612 1.715 1.360 0.932 1.986
#> 2 spring 20 0.857 0.169 1.545 1.679 1.198 2.354
#> 3 summer 20 1.246 0.561 1.930 1.688 1.155 2.467
#>
#> $cell_ab
#> cell n pda_median pda_lower pda_upper pdb_median pdb_lower pdb_upper
#> 1 fall 20 0.735 1.073 0.504 3.200 1.844 5.557
#> 2 spring 20 0.596 0.835 0.425 2.356 1.184 4.688
#> 3 summer 20 0.592 0.866 0.405 3.476 1.752 6.890
#>
#> $CL
#> [1] 0.9
#>
#> $cell_desc
#> cell medianCP r1 r3 r7 r14 r28
#> 1 fall 1.943506 0.7877998 0.5968500 0.4035271 0.2504831 0.1362684
#> 2 spring 1.273808 0.6946923 0.5067622 0.3406993 0.2163474 0.1221627
#> 3 summer 1.871560 0.7459730 0.5779176 0.4168571 0.2838980 0.1714068
```

The model summary shows descriptive statistics for the cellwise estimates of the `l`

and `s`

parameters. The `location`

and `scale`

parameterization is common in survival analysis, but the `pda`

and `pdb`

parameterization (Dalthorp and Huso, 2014) is also shown. These parameterizations are convenient to work with in statistical calculations but are not as convenient for giving users quick insight into the distributions, so a third set of summary statistics about the fitted distributions is also given. Namely, the median persistence time and the `r`

statistic, which is the probability that a carcass persists until the first search after arrival (assuming uniformly distributed arrival times within the interval). Clearly, `r`

depends on the length of the search interval, and the table shows `r`

for intervals of 1, 3, 7, 14, and 28 days. A rough, back-of-the-envelope calculation for the probability of observing a carcass that arrives at a site during the monitored period would be `DWP * r * p * f`

, which is `DWP`

= fraction of carcasses that arrive in the area searched at a unit, `r`

= fraction of carcasses that persist until a search, `p`

= fraction of carcasses found on the first search after arrival (given that they persisted), and `f`

= fraction of carcasses that arrive at the units searched.

In other scenarios we might consider other predictors, like the visibility of the ground searched or the search team. We might also be interested in carcasses of different sizes (e.g., large, medium, and small birds instead of or in addition to bats). We are not restricted to using the same predictors of both SE and CP. The modeling complexity increases with each additional predictor, but, in theory, any number of predictors can be used. The only rule is that sufficient numbers of trial carcasses must be placed in each cell combination of factor levels among the selected predictors. For example, if we were to place 15 carcasses for each cell for predictors that include season (spring, summer, fall, winter), size (S, M, L, B), visibility (RP, M, D), search team (dogs, humans), and turbine type (small, medium, large), we’d need 15 x 4 x 4 x 3 x 2 x 2 = 2880 carcasses. Typically, the number of predictors is limited to a few key variables.

Each carcass’s contribution to the total mortality in each search interval is estimated using the `estM`

function.

```
Mhat <- estM(nsim = 1000, data_CO = data_CO, data_SS = data_SS,
data_DWP = data_DWP, model_SE = model_SE, model_CP = model_CP,
unitCol = "Turbine", COdate = "DateFound")
summary(Mhat)
#> median 5% 95%
#> 1447.02 1021.36 2149.11
plot(Mhat)
```

Mortality estimates may be partitioned or split into desired categories, such as species, season, or turbine type. Splits may be performed according to characteristics of the carcasses or where they were found (e.g., species, turbine or other variable found in `data_CO`

) or when they were found (e.g., season or other variable associated with search schedule and found in `data_SS`

, or a vector of specific times).

Mortality by `Species`

(a CO split because it is a column in the CO file):

```
M_species <- calcSplits(M = Mhat, split_CO = "Species", data_CO = data_CO)
summary(M_species)
#> X 5% 25% 50% 75% 95%
#> BA 42 564.53025 706.85569 817.56556 982.35781 1277.15937
#> BB 14 180.30993 272.46342 344.49886 437.00362 639.70012
#> BC 9 75.82602 119.99603 156.58053 203.71929 289.95029
#> BD 5 24.99493 57.17510 83.61704 120.39107 177.47340
#> BE 1 1.00000 3.36694 11.95894 23.34428 44.74911
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "Species"
#> attr(,"type")
#> [1] "CO"
plot(M_species)
```

Mortality estimates may also be split by temporal variables that are represented as columns in `data_SS`

or as numeric vectors spanning the monitoring season (from day 0 to length of monitoring season). If several temporal splits are to be calculated, creating a specially formatted `prepSS`

object for the search schedule can streamline the calculations.

Mortality by `Season`

(an SS split because it is a column in the SS file):

```
M_season <- calcSplits(M = Mhat, split_SS = "Season", data_SS = SSdat,
split_CO = NULL, data_CO = data_CO)
summary(M_season)
#> X 5% 25% 50% 75% 95%
#> spring 18.10738 315.36044 457.6462 585.9600 768.8649 1145.3525
#> summer 10.18962 97.79268 153.2999 201.7058 268.9006 395.0101
#> fall 42.70300 414.51466 522.9091 615.3305 726.2513 962.1429
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "Season"
#> attr(,"type")
#> [1] "SS"
#> attr(,"times")
#> [1] 0 60 130 200
plot(M_season)
```

Mortality by month (a temporal split that spans the monitoring period):

```
M_month <- calcSplits(M = Mhat, split_time = seq(0, max(SSdat$days), by = 28),
data_SS = SSdat, data_CO = data_CO)
summary(M_month)
#> X 5% 25% 50% 75% 95%
#> 28 6.273750 89.096298 157.36243 221.45797 313.3490 518.5814
#> 56 11.706212 175.050123 271.19429 348.31998 464.5671 709.0288
#> 84 2.119896 5.634015 37.18257 65.80332 107.5328 193.2661
#> 112 4.568000 23.807485 48.91660 74.86874 106.8389 171.0078
#> 140 10.020857 75.257779 112.23405 145.17049 182.3748 254.8024
#> 168 27.861476 265.748155 341.39260 409.31692 496.5294 671.4778
#> 196 8.449810 53.104713 84.25840 112.53674 144.5494 205.9294
#> 200 0.000000 0.000000 0.00000 0.00000 0.0000 0.0000
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "time"
#> attr(,"type")
#> [1] "time"
#> attr(,"times")
#> [1] 0 28 56 84 112 140 168 196 200
plot(M_month)
```

Temporal splits that divide the monitoring season into separate time intervals (like season or month) can be plotted as the number per interval (`rate = FALSE`

, which is the default arg in `calcSplits`

) or the number per unit time (`rate = TRUE`

).

```
M_various_times <- calcSplits(M = Mhat,
split_time = c(seq(0, 90, by = 15), 120, 150, seq(155, 200, by = 5)),
data_SS = SSdat, data_CO = data_CO)
plot(M_various_times)
```

Finally, splits can be calculated for combinations of splitting covariates, like species by season or species group by turbine type. No more than two splitting covariates may be used in one call to `calcSplits`

and at most one temporal split may be used (whether it is an SS split or a vector of times).

Thorough searches out to a radius of 60 m from each turbine were conducted at 23 out of 100 turbines. The searched area was divided into three visibility classes (`RP`

, `M`

, `D`

) according the difficulty of finding carcasses.

Searcher efficiency and carcass persistence would be expected to vary with carcass size (sparrow, eagle, bat), ground characteristics (road & pad, cleared field, vegetation type), season, etc. In this example, we perform a full analysis of scenario with four classes of carcass (`lrg`

, `med`

, `sml`

, and `bat`

), three visibility classes (difficult = `D`

, moderate = `M`

, and road & pad = `RP`

), and three seasons (`spring`

, `summer`

, and `fall`

).

The required data is stored in `wind_cleared`

, a list of data frames with results for searcher efficiency (`SE`

) and carcass persistence trials (`CP`

), search schedules for all turbines (`SS`

), the search coverage or density weighted proportion (`DWP`

) of area searched at each turbine (i.e., the fraction of carcasses expected to fall in the search plots), and carcass observation (`CO`

) data.

Load the full data set into R:

To streamline the notation, extract the data from the `wind_cleared`

list into its components:

```
data_SE <- wind_cleared$SE
data_CP <- wind_cleared$CP
data_SS <- wind_cleared$SS
data_DWP <- wind_cleared$DWP
data_CO <- wind_cleared$CO
```

In searcher efficiency and carcass persistence trials, 15 trial carcasses were placed in each combination of visibility class (`D`

, `M`

, `RP`

), season (`spring`

, `summer`

, `fall`

), and size class (`lrg`

, `med`

, `sml`

, `bat`

). Data formats are like those of example 1:

```
head(data_SE)
#> pkID Size Season Visibility s1 s2 s3 s4 s5
#> 1 pk1 bat spring RP 0 NA NA NA NA
#> 2 pk2 bat spring RP 0 0 NA NA NA
#> 3 pk3 bat spring RP 1 NA NA NA NA
#> 4 pk4 bat spring RP 1 NA NA NA NA
#> 5 pk5 bat spring RP 0 NA NA NA NA
#> 6 pk6 bat spring RP 1 NA NA NA NA
head(data_CP)
#> cpID Size Season Visibility LastPresent FirstAbsent
#> 1 cp1 bat spring RP 0.00 0.93
#> 2 cp2 bat spring RP 0.98 1.97
#> 3 cp3 bat spring RP 0.00 1.01
#> 4 cp4 bat spring RP 13.99 21.13
#> 5 cp5 bat spring RP 0.00 1.17
#> 6 cp6 bat spring RP 20.95 27.97
```

With 36 combinations of covariate levels (3 visibilities x 3 seasons x 4 sizes) and two parameters (`p`

and `k`

), the number of possible models to consider for searcher efficiency is unwieldy using simple calls to `pkm`

, but `pkm`

has powerful model building and model selection capabilities that can be accessed via the arg list: `allCombos`

and `sizeCol`

, which are discussed below.

When `allCombos = TRUE`

, `pkm`

fits the set of submodels of the given covariate combinations, including the full model, the null model, and everything in between. For example, if the parameter models are `p ~ Visibility * Season`

and `k ~ Visibility`

, `pkm`

with `allCombos = TRUE`

would fit all combinations of possible `p`

models (`p ~ Visibility * Season`

, `p ~ Visibility + Season`

, `p ~ Visibility`

, `p ~ Season`

, and `p ~ 1`

) and possible `k`

models (`k ~ Visibility`

and `k ~ 1`

), or 10 models in all.

Carcasses in different size classes would be expected to have different searcher efficiency and carcass persistence parameters and would even be likely to be affected by covariates in different ways. When analyzing data with carcasses in different size classes, it is recommended that separate CP and SE models be fit for size classes separately. This can be accomplished using the `sizeCol`

argument in `pkm`

and `cpm`

, where `sizeCol`

gives the name of the column that gives the carcass size classes in `data_SE`

.

```
pkModels <- pkm(p ~ Visibility * Season, k ~ Visibility * Season, data = data_SE,
allCombos = TRUE, sizeCol = "Size")
class(pkModels)
#> [1] "pkmSetSize" "list"
names(pkModels)
#> [1] "bat" "lrg" "med" "sml"
```

When `allCombos = TRUE`

and `sizeCol`

is defined, `pkm`

returns a list of sets of models for each size class. The sets of models for each size class include the full spectrum of models that can be constructed using simple combinations of the covariates.

```
names(pkModels[["sml"]])
#> [1] "p ~ Visibility * Season; k ~ Visibility * Season" "p ~ Visibility + Season; k ~ Visibility * Season"
#> [3] "p ~ Season; k ~ Visibility * Season" "p ~ Visibility; k ~ Visibility * Season"
#> [5] "p ~ 1; k ~ Visibility * Season" "p ~ Visibility * Season; k ~ Visibility + Season"
#> [7] "p ~ Visibility + Season; k ~ Visibility + Season" "p ~ Season; k ~ Visibility + Season"
#> [9] "p ~ Visibility; k ~ Visibility + Season" "p ~ 1; k ~ Visibility + Season"
#> [11] "p ~ Visibility * Season; k ~ Season" "p ~ Visibility + Season; k ~ Season"
#> [13] "p ~ Season; k ~ Season" "p ~ Visibility; k ~ Season"
#> [15] "p ~ 1; k ~ Season" "p ~ Visibility * Season; k ~ Visibility"
#> [17] "p ~ Visibility + Season; k ~ Visibility" "p ~ Season; k ~ Visibility"
#> [19] "p ~ Visibility; k ~ Visibility" "p ~ 1; k ~ Visibility"
#> [21] "p ~ Visibility * Season; k ~ 1" "p ~ Visibility + Season; k ~ 1"
#> [23] "p ~ Season; k ~ 1" "p ~ Visibility; k ~ 1"
#> [25] "p ~ 1; k ~ 1"
class(pkModels[["sml"]])
#> [1] "pkmSet" "list"
```

To estimate mortality, one model for each size class must be selected from the long list of models fit. GenEst provides several tools for guiding the selection. First, the models can be listed by AICc, which gives a score for the quality of the model for the given data. Complicated models that use many parameters may fit the data more closely than a simpler model but are penalized because of their complexity and relative instability. The scores have meaning only in comparison with other models’. AICc provides a rough but useful guide for model selection, but should in no way be relied upon as definitive. Its utility is in identifying relatively poor models and in narrowing the choice of plausible models to a manageable number.

The `aicc`

function lists the fitted models in order of \({\small \Delta}\)AICc for each size class (if applicable). For this discussion, we will focus on `sml`

only, but for a full analysis, all size classes would need to be similarly analyzed.

```
aicc(pkModels[["sml"]])
#> p Formula k Formula AICc <U+0394>AICc
#> 24 p ~ Visibility k ~ 1 318.55 0.00
#> 19 p ~ Visibility k ~ Visibility 319.69 1.14
#> 14 p ~ Visibility k ~ Season 321.42 2.87
#> 21 p ~ Visibility * Season k ~ 1 321.53 2.98
#> 22 p ~ Visibility + Season k ~ 1 322.17 3.62
#> 16 p ~ Visibility * Season k ~ Visibility 322.98 4.43
#> 17 p ~ Visibility + Season k ~ Visibility 323.14 4.59
#> 9 p ~ Visibility k ~ Visibility + Season 323.44 4.89
#> 11 p ~ Visibility * Season k ~ Season 323.97 5.42
#> 12 p ~ Visibility + Season k ~ Season 324.33 5.78
#> 7 p ~ Visibility + Season k ~ Visibility + Season 326.33 7.78
#> 6 p ~ Visibility * Season k ~ Visibility + Season 326.51 7.96
#> 4 p ~ Visibility k ~ Visibility * Season 327.99 9.44
#> 2 p ~ Visibility + Season k ~ Visibility * Season 331.07 12.52
#> 1 p ~ Visibility * Season k ~ Visibility * Season 333.47 14.92
#> 20 p ~ 1 k ~ Visibility 339.51 20.96
#> 10 p ~ 1 k ~ Visibility + Season 342.41 23.86
#> 18 p ~ Season k ~ Visibility 342.74 24.19
#> 8 p ~ Season k ~ Visibility + Season 344.10 25.55
#> 25 p ~ 1 k ~ 1 346.45 27.90
#> 5 p ~ 1 k ~ Visibility * Season 347.47 28.92
#> 15 p ~ 1 k ~ Season 347.66 29.11
#> 23 p ~ Season k ~ 1 350.35 31.80
#> 13 p ~ Season k ~ Season 350.44 31.89
#> 3 p ~ Season k ~ Visibility * Season 350.49 31.94
```

Preference should normally be given to models with \({\small \Delta}\)AICc less than 6 or 7. Models with differences of less than 3 or 4 are generally considered indistinguishable by this measure. Choices among such models should be based on other criteria.

Diagnostic plots can be used to identify potential problems with model fits and to help distinguish between models with similar AICc scores. The `plot`

function is defined for `pkm`

(one model) and `pkmSet`

(set of models for a given size class) objects. For example, `plot(pkModels[["sml"]][[1]])`

would produce the single figure for the first model for the `sml`

size class. `plot(pkModels[["sml"]])`

would create plots for each model fit for the `sml`

size class. To plot a specific single model from the full set, use the `specificModel`

argument. For example, diagnostic plots for the model with the lowest AICc score are shown below:

The top row shows box plots of estimated `p`

and `k`

parameters for all cells (i.e., combinations of covariate levels, like `D.fall`

for difficult visibility in the fall) for both the selected model (black) and the full model (gray). With the full model, the fits for each cell are based solely on data from that specific cell. The advantage is that each cell’s estimates are untainted by data from other cells. The disadvantage is that the sample size for each estimate is relatively small and the error bars large. In the reduced models, estimates for one cell borrow strength from estimates in related cells. This gives smaller error bars but can lead to errors if the model structure does not properly reflect the dependence of searcher efficiency parameters on cell characteristics.

In the figure, the estimates of `p`

from the selected model are markedly less variable than the estimates from the full model, while the locations of the boxes are very similar for the two models. Thus, this selected model (the one with the lowest AICc) appears to be an improvement over the full cell model for estimating `p`

.

The \({\small \Delta}\)AICc for the full model is 12.52, which indicates a serious deficiency in comparison to the model with the best fit, `"p ~ Visibility; k ~ 1"`

. The boxplots for `k`

highlight one particular problem with the fit of the reference model. In some of the cells, the boxes extend from 0 to 1, which suggests that the reference model is unable to estimate `k`

for the given cell. Selecting a simpler model for `k`

often remedies this problem. If all models display this 0-1 phenomenon, a fixed `k`

of 1 is appropriate if a smaller proportion of carcasses was found on the first search occasion than on later searches.

By comparison, the model with the highest \({\small \Delta}\)AICc (= `{r max(aicc(pkModels[["sml"]]))}`

routinely estimates `p`

and `k`

either well above or well below the reference model but has fairly tight error bars–bad estimates but quite confident about them!

A similar model selection exercise gives the same form of model (`p ~ Visibility; k ~ 1`

) for large birds, medium birds, and bats. These can all be collated into a list for later analysis of detection probabilities and mortality rates.

The work flow for carcass persistence modeling is similar to that for searcher efficiency except that in addition to selecting covariates for two different parameters (location = `l`

and scale = `s`

), there are four model forms to choose from: Weibull, lognormal, loglogistic, and exponential.

```
cpModels <- cpm(
l ~ Visibility * Season, s ~ Visibility * Season,
data = data_CP, left = "LastPresent", right = "FirstAbsent",
dist = c( "weibull", "lognormal", "loglogistic", "exponential"),
allCombos = TRUE, sizeCol = "Size"
)
```

The list of models is long:

```
lapply(aicc(cpModels), nrow)
#> $bat
#> [1] 80
#>
#> $lrg
#> [1] 80
#>
#> $med
#> [1] 80
#>
#> $sml
#> [1] 80
aicc(cpModels[["sml"]])
#> Distribution Location Formula Scale Formula AICc <U+0394>AICc
#> 25 weibull l ~ 1 s ~ 1 780.00 0.00
#> 24 weibull l ~ Visibility s ~ 1 780.16 0.16
#> 15 weibull l ~ 1 s ~ Season 780.38 0.38
#> 19 weibull l ~ Visibility s ~ Visibility 780.97 0.97
#> 20 weibull l ~ 1 s ~ Visibility 781.32 1.32
#> 10 weibull l ~ 1 s ~ Visibility + Season 782.00 2.00
#> 14 weibull l ~ Visibility s ~ Season 782.04 2.04
#> 22 weibull l ~ Visibility + Season s ~ 1 782.55 2.55
#> 74 loglogistic l ~ Visibility s ~ 1 782.57 2.57
#> 21 weibull l ~ Visibility * Season s ~ 1 782.77 2.77
#> 49 lognormal l ~ Visibility s ~ 1 782.84 2.84
#> 23 weibull l ~ Season s ~ 1 782.88 2.88
#> 71 loglogistic l ~ Visibility * Season s ~ 1 783.06 3.06
#> 13 weibull l ~ Season s ~ Season 783.16 3.16
#> 50 lognormal l ~ 1 s ~ 1 783.18 3.18
#> 75 loglogistic l ~ 1 s ~ 1 783.30 3.30
#> 9 weibull l ~ Visibility s ~ Visibility + Season 783.33 3.33
#> 16 weibull l ~ Visibility * Season s ~ Visibility 783.42 3.42
#> 46 lognormal l ~ Visibility * Season s ~ 1 783.86 3.86
#> 47 lognormal l ~ Visibility + Season s ~ 1 783.86 3.86
#> 40 lognormal l ~ 1 s ~ Season 783.88 3.88
#> 11 weibull l ~ Visibility * Season s ~ Season 783.95 3.95
#> 72 loglogistic l ~ Visibility + Season s ~ 1 784.00 4.00
#> 17 weibull l ~ Visibility + Season s ~ Visibility 784.09 4.09
#> 48 lognormal l ~ Season s ~ 1 784.13 4.13
#> 69 loglogistic l ~ Visibility s ~ Visibility 784.28 4.28
#> 65 loglogistic l ~ 1 s ~ Season 784.30 4.30
#> 12 weibull l ~ Visibility + Season s ~ Season 784.48 4.48
#> 73 loglogistic l ~ Season s ~ 1 784.54 4.54
#> 44 lognormal l ~ Visibility s ~ Visibility 784.58 4.58
#> 18 weibull l ~ Season s ~ Visibility 784.85 4.85
#> 66 loglogistic l ~ Visibility * Season s ~ Visibility 784.85 4.85
#> 64 loglogistic l ~ Visibility s ~ Season 784.95 4.95
#> 39 lognormal l ~ Visibility s ~ Season 785.01 5.01
#> 38 lognormal l ~ Season s ~ Season 785.09 5.09
#> 45 lognormal l ~ 1 s ~ Visibility 785.14 5.14
#> 6 weibull l ~ Visibility * Season s ~ Visibility + Season 785.29 5.29
#> 70 loglogistic l ~ 1 s ~ Visibility 785.31 5.31
#> 8 weibull l ~ Season s ~ Visibility + Season 785.43 5.43
#> 61 loglogistic l ~ Visibility * Season s ~ Season 785.47 5.47
#> 41 lognormal l ~ Visibility * Season s ~ Visibility 785.68 5.68
#> 63 loglogistic l ~ Season s ~ Season 785.71 5.71
#> 36 lognormal l ~ Visibility * Season s ~ Season 785.89 5.89
#> 35 lognormal l ~ 1 s ~ Visibility + Season 786.07 6.07
#> 42 lognormal l ~ Visibility + Season s ~ Visibility 786.15 6.15
#> 37 lognormal l ~ Visibility + Season s ~ Season 786.29 6.29
#> 67 loglogistic l ~ Visibility + Season s ~ Visibility 786.30 6.30
#> 7 weibull l ~ Visibility + Season s ~ Visibility + Season 786.48 6.48
#> 43 lognormal l ~ Season s ~ Visibility 786.48 6.48
#> 60 loglogistic l ~ 1 s ~ Visibility + Season 786.56 6.56
#> 62 loglogistic l ~ Visibility + Season s ~ Season 786.56 6.56
#> 68 loglogistic l ~ Season s ~ Visibility 786.89 6.89
#> 59 loglogistic l ~ Visibility s ~ Visibility + Season 786.94 6.94
#> 34 lognormal l ~ Visibility s ~ Visibility + Season 787.01 7.01
#> 56 loglogistic l ~ Visibility * Season s ~ Visibility + Season 787.41 7.41
#> 33 lognormal l ~ Season s ~ Visibility + Season 787.71 7.71
#> 31 lognormal l ~ Visibility * Season s ~ Visibility + Season 787.87 7.87
#> 58 loglogistic l ~ Season s ~ Visibility + Season 788.35 8.35
#> 32 lognormal l ~ Visibility + Season s ~ Visibility + Season 788.83 8.83
#> 57 loglogistic l ~ Visibility + Season s ~ Visibility + Season 789.12 9.12
#> 5 weibull l ~ 1 s ~ Visibility * Season 789.98 9.98
#> 4 weibull l ~ Visibility s ~ Visibility * Season 791.45 11.45
#> 3 weibull l ~ Season s ~ Visibility * Season 793.01 13.01
#> 1 weibull l ~ Visibility * Season s ~ Visibility * Season 794.00 14.00
#> 30 lognormal l ~ 1 s ~ Visibility * Season 794.47 14.47
#> 2 weibull l ~ Visibility + Season s ~ Visibility * Season 794.50 14.50
#> 55 loglogistic l ~ 1 s ~ Visibility * Season 794.95 14.95
#> 29 lognormal l ~ Visibility s ~ Visibility * Season 795.13 15.13
#> 54 loglogistic l ~ Visibility s ~ Visibility * Season 795.14 15.14
#> 28 lognormal l ~ Season s ~ Visibility * Season 795.65 15.65
#> 53 loglogistic l ~ Season s ~ Visibility * Season 796.37 16.37
#> 51 loglogistic l ~ Visibility * Season s ~ Visibility * Season 796.80 16.80
#> 26 lognormal l ~ Visibility * Season s ~ Visibility * Season 797.20 17.20
#> 27 lognormal l ~ Visibility + Season s ~ Visibility * Season 797.36 17.36
#> 52 loglogistic l ~ Visibility + Season s ~ Visibility * Season 797.70 17.70
#> 76 exponential l ~ Visibility * Season NULL 835.23 55.23
#> 79 exponential l ~ Visibility NULL 839.78 59.78
#> 77 exponential l ~ Visibility + Season NULL 840.96 60.96
#> 80 exponential l ~ 1 NULL 842.24 62.24
#> 78 exponential l ~ Season NULL 844.75 64.75
```

It is not uncommon to see the fits for the exponential distribution at the bottom of the AIC list. The exponential distribution has only one parameter and does not have nearly as much flexibility as the others, which each have two-parameters. An implicit assumption of the exponential model is that the scavenging rate is constant, regardless of carcass age. When that assumption is not met, the exponential provides an inferior fit.

To compare among a set of cp models with plausible AICs, use, for example, `plot(cpModels[["sml"]])`

to browse through all the models or `plot(cpModels[["sml"]][["p ~ Visibility; k ~ 1"]])`

to view a single model specified by name. It can be seen from the AICc table (`aicc(cpModels[["sml"]])`

) that the top 10 models according to AICc are:

```
cp_smlCandidates <-
names(cpModels[["sml"]])[c(25, 24, 15, 19, 20, 10, 14, 22, 74, 21)]
cp_smlCandidates
#> [1] "dist: weibull; l ~ 1; s ~ 1" "dist: weibull; l ~ Visibility; s ~ 1"
#> [3] "dist: weibull; l ~ 1; s ~ Season" "dist: weibull; l ~ Visibility; s ~ Visibility"
#> [5] "dist: weibull; l ~ 1; s ~ Visibility" "dist: weibull; l ~ 1; s ~ Visibility + Season"
#> [7] "dist: weibull; l ~ Visibility; s ~ Season" "dist: weibull; l ~ Visibility + Season; s ~ 1"
#> [9] "dist: loglogistic; l ~ Visibility; s ~ 1" "dist: weibull; l ~ Visibility * Season; s ~ 1"
```

These can be compared in graphs as follows:

The figure shows the raw persistence data (fraction of carcasses remaining after the given time) for each cell, as a black stair case with Kaplan-Meier confidence intervals as dashed lines. In addition, the fitted curves for each of the distributions are shown in color, with the `specificModel`

distribution having a thicker than the others. The two-parameter models (Weibull, lognormal, and loglogistic) tend to be very similar and a relatively close fit to the data. The exponential model tends to be somewhat removed from the others. Clicking on the graphing window brings up the next set of figures.

We’re looking for a good fit between the selected model and the data in as many cells as possible. The first several models among the `cp_smlCandidates`

seems to provide a reasonably good fit in all cells, although there seems to be a trade-off between fitting the RP.spring cell well or fitting the RP.summer cell well. `Visibility`

occurs frequently in the top models, while `Season`

appears more frequently in the bottom models. This indicates that `Season`

is probably not a strong predictor of carcass persistence, while `Visibility`

is.

Selecting the top AICc model for use in mortality estimation:

Following a similar model selection process for the other size classes, we select the following:

```
cp_med <- cpModels[["med"]][["dist: weibull; l ~ Visibility; s ~ Season"]]
cp_lrg <- cpModels[["lrg"]][["dist: exponential; l ~ Visibility + Season; NULL"]]
cp_bat <- cpModels[["bat"]][["dist: weibull; l ~ Visibility + Season; s ~ 1"]]
```

NOTE: For large carcasses, the exponential distribution was at the top of the AICc list. The exponential requires only one parameter (`l`

), so no scale parameter is provided in the model.

Again, we collate the models into a list for later analysis of detection probabilities and mortality rates.

Each carcass’s contribution to the total mortality in each search interval is estimated using the `estM`

function. The function call is largely similar that used in the simple scenario discussed in example 1. However, there are some important differences. First, `model_SE`

and `model_CP`

are lists of models, one element for each size class. In addition, the name of the size class variable is provided as `sizeCol = "Size"`

. Finally, the `frac`

argument represents the sampling fraction or the fraction of carcasses expected to fall at the units that were searched. In this example, 23 out of 100 turbines were searched, and `frac`

is set equal to 0.23 under the assumption that the mortality rates at the unsearched turbines did not differ substantially from the rates at the searched turbines.

```
Mhat <- estM(nsim = 1000, data_CO = data_CO, data_SS = data_SS, frac = 0.23,
data_DWP = data_DWP, model_SE = pkMods, model_CP = cpMods,
sizeCol = "Size", unitCol = "Turbine", COdate = "DateFound")
summary(Mhat)
#> median 5% 95%
#> 1926.32 1505.99 2463.93
plot(Mhat)
```

This estimate is for the total number of fatalities among all size classes combined, from hummingbirds and bats to eagles and may be too vague to be very useful. Fortunately, mortality estimates may be partitioned or split into desired categories, such as species, size, or season. Splits may be performed according to characteristics of the carcasses or where they were found (e.g., species, turbine or other variable found in `data_CO`

) or when they were found (e.g., season or other variable associated with search schedule and found in `data_SS`

, or a vector of specific times).

Carcasses were categorized not only by size but also by species, species group, the type of turbine they were found at, the visibility class of the ground where they were found, and distance from nearest turbine.

Although species groups may extend across different size classes (e.g., `raptors`

could include kestrels, red-tailed hawks, and golden eagles; `passerines`

could include sparrows and ravens), splits according to species group can easily be accomplished:

```
M_speciesGroup <- calcSplits(M = Mhat,
split_CO = "SpeciesGroup", data_CO = data_CO)
summary(M_speciesGroup)
#> X 5% 25% 50% 75% 95%
#> bat0 4 18.33396 41.96731 62.61233 89.03721 133.69702
#> bat1 49 809.73453 973.78664 1113.19755 1249.40451 1535.40347
#> brd1 19 309.12235 425.37161 524.22698 659.57902 888.08233
#> brd2 10 70.26812 110.08250 139.17522 175.88369 233.38252
#> brd3 5 18.82971 39.50241 54.78262 71.39421 99.09671
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "SpeciesGroup"
#> attr(,"type")
#> [1] "CO"
plot(M_speciesGroup)
```

Split by species and season:

```
M_speciesseason <- calcSplits(M = Mhat,
split_CO = "Species", data_CO = data_CO, split_SS = "Season", data_SS = data_SS)
summary(M_speciesseason)
#> $BA
#> X 5% 25% 50% 75% 95%
#> spring 7.000 78.03770 139.53775 188.74411 248.1228 353.9261
#> summer 3.298 13.43985 46.95295 76.29841 112.7901 175.9622
#> fall 22.702 276.78301 359.31352 426.39835 507.4478 630.3254
#>
#> $BB
#> X 5% 25% 50% 75% 95%
#> spring 3.000 19.70609 70.38817 116.35144 176.97157 265.68982
#> summer 1.149 1.14900 8.10605 23.96345 45.31941 88.96974
#> fall 6.851 56.19216 94.96797 127.16241 164.02256 221.68390
#>
#> $BC
#> X 5% 25% 50% 75% 95%
#> spring 1.0281538 1.0281538 9.992437 30.49797 62.54395 123.53807
#> summer 0.9738462 0.9738462 5.967803 19.65440 37.60093 72.23311
#> fall 2.9980000 7.4416836 27.907840 49.48437 73.58833 117.18915
#>
#> $BD
#> X 5% 25% 50% 75% 95%
#> spring 0.000 0.000 0.00000 0.00000 0.00000 0.00000
#> summer 0.002 0.002 0.00200 0.00200 0.00200 0.00200
#> fall 1.998 1.998 14.98986 29.37361 47.05278 84.49484
#>
#> $BE
#> X 5% 25% 50% 75% 95%
#> spring 0 0 0.00000 0.00000 0.00000 0.00000
#> summer 0 0 0.00000 0.00000 0.00000 0.00000
#> fall 2 2 16.05041 29.68113 47.63888 78.11209
#>
#> $LA
#> X 5% 25% 50% 75% 95%
#> spring 2 2 9.848526 18.76283 27.74579 42.98474
#> summer 0 0 0.000000 0.00000 0.00000 0.00000
#> fall 2 2 11.146153 19.75408 30.57606 47.18836
#>
#> $LB
#> X 5% 25% 50% 75% 95%
#> spring 1 1 3.105801 8.709779 16.31249 29.54316
#> summer 0 0 0.000000 0.000000 0.00000 0.00000
#> fall 0 0 0.000000 0.000000 0.00000 0.00000
#>
#> $LC
#> X 5% 25% 50% 75% 95%
#> spring 1.0141 1.0141 2.518306 8.597191 15.13264 28.75196
#> summer 0.9859 0.9859 2.739716 8.068422 14.30869 26.06010
#> fall 0.0000 0.0000 0.000000 0.000000 0.00000 0.00000
#>
#> $LE
#> X 5% 25% 50% 75% 95%
#> spring 0.066 0.066 0.066000 0.066000 0.06600 2.303007
#> summer 0.934 0.934 1.495151 7.689636 14.60354 26.788470
#> fall 0.000 0.000 0.000000 0.000000 0.00000 0.000000
#>
#> $MA
#> X 5% 25% 50% 75% 95%
#> spring 1.001 1.001 4.819965 14.97960 27.86713 52.06967
#> summer 0.999 0.999 3.720563 11.89201 21.81530 37.81837
#> fall 1.000 1.000 3.518491 11.52242 20.62347 36.45960
#>
#> $MB
#> X 5% 25% 50% 75% 95%
#> spring 0.000 0.000 0.000000 0.00000 0.00000 0.00000
#> summer 0.021 0.021 0.021000 0.02100 0.02100 0.02100
#> fall 0.979 0.979 2.709526 10.08919 18.58794 33.23729
#>
#> $MC
#> X 5% 25% 50% 75% 95%
#> spring 1.284 1.284 7.414389 18.877355 34.83163 61.73823
#> summer 0.716 0.716 0.716000 8.800436 24.90795 50.61471
#> fall 0.000 0.000 0.000000 0.000000 0.00000 0.00000
#>
#> $MD
#> X 5% 25% 50% 75% 95%
#> spring 1 1 3.956838 12.17275 22.46417 41.2197
#> summer 0 0 0.000000 0.00000 0.00000 0.0000
#> fall 0 0 0.000000 0.00000 0.00000 0.0000
#>
#> $ME
#> X 5% 25% 50% 75% 95%
#> spring 0.000 0.000 0.00000 0.00000 0.00000 0.00000
#> summer 0.004 0.004 0.00400 0.00400 0.00400 0.00400
#> fall 0.996 0.996 3.36441 10.02639 19.41026 32.97884
#>
#> $MF
#> X 5% 25% 50% 75% 95%
#> spring 0 0.000000 0.00000 0.000 0.00000 0.0000
#> summer 0 0.000000 0.00000 0.000 0.00000 0.0000
#> fall 3 7.132314 25.67748 42.902 62.21671 100.0018
#>
#> $MH
#> X 5% 25% 50% 75% 95%
#> spring 1 1 4.150007 12.699201 22.70208 40.67432
#> summer 1 1 3.225510 9.947665 18.07888 33.39145
#> fall 0 0 0.000000 0.000000 0.00000 0.00000
#>
#> $SA
#> X 5% 25% 50% 75% 95%
#> spring 1 1 10.82479 34.24883 66.24702 121.3315
#> summer 0 0 0.00000 0.00000 0.00000 0.0000
#> fall 0 0 0.00000 0.00000 0.00000 0.0000
#>
#> $SB
#> X 5% 25% 50% 75% 95%
#> spring 0.456 0.456 0.456 0.456 29.49794 95.32605
#> summer 0.544 0.544 0.544 0.544 29.03662 80.06962
#> fall 0.000 0.000 0.000 0.000 0.00000 0.00000
#>
#> $SC
#> X 5% 25% 50% 75% 95%
#> spring 3 17.80458 63.35057 97.27212 147.2447 251.5965
#> summer 0 0.00000 0.00000 0.00000 0.0000 0.0000
#> fall 0 0.00000 0.00000 0.00000 0.0000 0.0000
#>
#> $SD
#> X 5% 25% 50% 75% 95%
#> spring 1.1359231 1.1359231 9.282457 25.70477 48.83796 103.72860
#> summer 0.8640769 0.8640769 4.558270 23.10571 47.74858 96.25288
#> fall 0.0000000 0.0000000 0.000000 0.00000 0.00000 0.00000
#>
#> $SE
#> X 5% 25% 50% 75% 95%
#> spring 0.003 0.003 0.003000 0.00300 0.00300 0.00300
#> summer 0.997 0.997 8.239688 27.83813 53.56048 95.42225
#> fall 0.000 0.000 0.000000 0.00000 0.00000 0.00000
#>
#> $SG
#> X 5% 25% 50% 75% 95%
#> spring 3.0415 19.76785 68.02091 108.67316 165.10378 259.1624
#> summer 1.9585 1.95850 29.08550 56.83186 92.25347 169.1607
#> fall 0.0000 0.00000 0.00000 0.00000 0.00000 0.0000
#>
#> attr(,"class")
#> [1] "splitSummary"
#> attr(,"CL")
#> [1] 0.9
#> attr(,"vars")
#> [1] "Season" "Species"
#> attr(,"type")
#> [1] "SS" "CO"
#> attr(,"times")
#> [1] 0 60 130 200
plot(M_speciesseason)
```

There are so many vertical panels that it is difficult to glean any useful information out of the graph. However, the panels may be transposed and graphed for better interpretability: