Using acmeR
acme.summary
This function provides all the information that most users will require. acme.summary() has only one required parameter - fname, the data described in the previous section. fname can be either a character string file path that points to a CSV file, or a data frame. If fname is not specified, this function will throw an error and then create an example dataset titled ‘altamont.csv’ in the working directory.
acme.summary()## Error: Missing csv file 'fname'. Example 'altamont.csv' created in
## working directory. Try acme.summary('altamont.csv').Running acme.summary() with no file name specified also allows a user to examine an example dataset to get a better understand of the data required. (Note, the same data is attached with the library and can be accessed as a data frame by data(altamont), but some users might find the CSV format more convenient).
The following parameters are optional:
- spec allows the user to subset the analysis to certain species, specified by a character vector of all-captial 4-letter American Ornithologists’ Union (AOU) codes. For example, spec = c(“BHCO”, “HOWR”) will subset to Brown-Headed Cowbirds and House Wrens. A full list of codes can be found at The Institute for Bird Populations
- blind is a logical that is by default set to TRUE. This assumes that Field Technicians (FT) are unware of carcasses. If FALSE, this sets all subsequent FT discoveries (after the initial discovery) of a carcass as PFM checks.
- ps allows the user to generate a postscript file - if left as default empty, no file will be generated
- plot_scav and plot_srch are logicals which allow the user to include persistence and search efficiency plots, respectively. For each plot, the estimated distributions with MLE parameters are overlayed on empirical distributions.
The output from acme.summary() is two-fold - the details of the model printed to the console, and the model estimates returned as a list. For most users, the information printed to the cosole will be sufficient. However, those intent on finding posterior mortality probability through acme.post() or acme.table() will find use for the returned values.
For example, suppose we are using the Altamont dataset. To use the csv that we generated above (which should be in the user’s working directory), we run the following line:
acme.val <- acme.summary(fname = 'altamont.csv',spec=c("BHCO", "HOWR"))acme.summary() can also be run with a data frame rather than a csv. The following line of code will run the function on a data frame version of the csv we created earlier, which (again) can be accessed by data(altamont). (Note, CSV files given by a character string and data frames given by R objects are currently the only two classes for fname that are accepted.)
acme.val <- acme.summary(fname=altamont,spec=c("BHCO","HOWR"))## Data set 'altamont' includes records of 104 placed carcasses
## of 1 species:
## Brown-Headed Cowbird (BHCO:104)
##
## Field technicians discovered 32 of these 104 placed carcasses
## (30.77%) on at least one occasion. Carcass placements ranged
## from Fri, Jan 07, 2011 to Fri, Apr 01, 2011.
##
## FT Search intervals averaged 7.44 days, with sd = 2.05 days.
## PFM Check intervals averaged 2.51 days, with sd = 1.15 days.
##
## The estimated scavenger removal survival function is:
##
## P[ remain > t days ] = exp(-(0.0764*t)^0.4986),
##
## (Weibull distribution) with mean removal time of 26.32 days.
## Weibull shape parameter alpha=0.4986 is 8.8 standard errors away
## from 1, making the Exponential Distribution completely implausible
## (P-value < 2.1e-18).
##
## Field Technicians tried 222 times to find carcasses known to be present,
## and succeeded 38 times (17.12% of the time). Of the total number 104 of
## carcasses placed, 32 were removed before the first FT search. Of the
## remaining 72 carcasses ever present during a search, 32 were ever
## discovered (44.44%), 23 of those (10.36%) on the first try.
## The estimated FT search proficiency is
##
## P[ discovery | carcass age = t days ] = exp(-1.0067-0.0771*t),
##
## starting at 36.54% for fresh carcasses, falling to 20.58% after one
## 7.4-day FT search interval, 11.59% after two, and 6.52% after three.
##
## Carcasses were later discovered after an initial miss 9 times, and
## after some earlier miss 12 times, suggesting bleed-through occurred.
## Estimated bleed-through rate is 99.89%.
##
## ACME mortality estimates for 7-day search intervals that discover C
## carcasses are:
## M* = 3.87 * C = C/R*, for R* = 0.2584.acme.val## $params
## alp.alp rho.rho a b bt
## 0.49864804 0.07637650 1.00667650 0.07713975 0.99886031
##
## $Rstar
## [1] 0.2583747
##
## $T
## [1] 0.1807887
##
## $I
## mu
## 7As is shown the output to the console is quite extensive. The top gives a summary of the data, with our subset species listed (here we see Altamont includes no House Wren carcasses) as well as useful information such as search and check intervals. The middle paragraphs describe the model and MLE parameters, while the bottom of the console output displays the ACME inverse-inflation factor R*. In this case, R* is 0.2584, so the ACME estimate for mortality given the number of carcasses C would be C/0.2584.
The other output is the list that the function returns. The components of the list are a 5-element numeric vector titled params that contains the MLE parameters calculated to run the model, a numeric titled Rstar containing the same as R* printed to the console, a numeric titled T containing the first component of R* as calculated in equation 7b in Wolpert (2015), and a numeric titled I containing the mean mean interval time between searches. These values will be useful when running acme.post() or acme.table().
acme.post
acme.post() allows the user to plot the posterior distribution of true mortality M given observed counts C. The model specifies the mortality as being drawn from a Poisson distribution with mean \(mI\),where \(I\) is the interval length and m is a parameter with prior distribution of Gamma(\(\xi\), \(\lambda\)). Appendix A2 of Wolpert (2015) gives the full model description, and calculates the posterior of M given C to be a scaled Gauss’ hypergeometric function.
acme.post() has 10 parameters, with defaults relating to the example from the paper (and thus this Vignette). The parameters are as follows:
- C is the observed mortality counts in the current search.
- Rstar, T, and I are values returned from acme.summary()
- gam is a 2-element vector of coverage probabilities for highest posterior density (hpd) credible intervals.
- xlim is 2-element vector of plotting ranges (if the maximum value of the larger hpd credible interval is xmax, then the default value for xlim is c(0,xmax + 2)).
- Mmax is the largest value of M for which the posterior probability is calculated, for purposes of finding the posterior mean.
- xi and lam are the prior parameters \(\xi\) and \(\lambda\), respectively, with default values relating to the Objective Bayes approach.
- ps is the postscript message, with default empty string that supresses output.
- plotit is a boolean which allows the user to suppress plotting.
#Carcass count of 5, with output from acme.summary as parameters
acme.post(C=5, Rstar=acme.val$Rstar, T=acme.val$T, I = acme.val$I,
xi = 1/2, lam = 0)
## gam cvg lo hi
## [1,] 0.5 0.5114542 13 23
## [2,] 0.9 0.9057251 8 34Here we see up to 38 posterior probability values for the mortality count M given the observed count C. In this case the observed fatality count is 5, the ACME estimator (\(\hat{M}\)) is 19.35, and the posterior mean (\(\overline{M}\)) is 20.94. This plot makes explicit that while the model may estimate close to 20 mortalities, there are a wide range of true mortalities with non-negligable posterior probabilities. Even if the count is 5, there is a small, nonzero chance that the true number of mortalities in this search is fewer than 5 (i.e. that bleed-through occurred).
The 90% hpd credible interval is covered by M=8 to M=34 (as shown by blue circles), and we see this range actually covers 90.57% of the density (because the distribution is discrete, exact coverage probabilities cannot be obtained, and hpd credible intervals will always cover more than the exact specified amount). The red squres show the 50% hpd credible interval from 13 to 24 in this case.
#Carcass count at 3, different hpd credible intervals values
acme.post(C=3,Rstar=acme.val$Rstar, T=acme.val$T, I = acme.val$I,
xi = 1/2, lam = 0, gam = c(.9,.95))
## gam cvg lo hi
## [1,] 0.90 0.9094126 3 23
## [2,] 0.95 0.9526068 3 27This graph depicts a different observed count (here C = 3), and a different set of values for hpd credible intervals (here we have 0.9 and 0.95 rather than the default 0.5 and 0.9).
acme.table
If the user is interested in summary statistics for multiple observed counts, but does not need plots for each one, acme.table() is valuable. This function serves as a wrapper for acme.post() with plotit=FALSE, while concatinating different values of C (acme.post() can only be run with a single value for C). The function outputs a matrix with a row for each observed count, along with ACME estimate, posterior mean, and specified highest posterior density credible intervals. With the exception of the counts (and lack of plotit parameter), the parameters of acme.table() are identical to those of acme.post().
acme.table(C=0:5, Rstar=acme.val$Rstar, T=acme.val$T, I = acme.val$I,
xi = 1/2, lam = 0, gam=c(.5,.9))## C M_hat Post_Mean 0.5_low 0.5_hi 0.9_low 0.9_hi
## [1,] 0 0.00 1.59 0 1 0 5
## [2,] 1 3.87 5.46 1 5 0 11
## [3,] 2 7.74 9.33 3 9 2 18
## [4,] 3 11.61 13.20 6 14 3 23
## [5,] 4 15.48 17.07 10 19 5 28
## [6,] 5 19.35 20.94 13 23 8 34References
Wolpert, Robert. 2015. “ACME A Partially Periodic Estimator of Avian & Chiropteran Mortality at Wind Turbines.” ArXiv E-Prints. http://arxiv.org/abs/1507.00749.