R/FatalFcns.R

In josephstatwick/EagleModelCalc: Eagle fatality calculator for wind farms

#' @title Fatality simulator (deprecated)
#' @author Mark Otto (FWS)
#'
#' @description This function calculates the annual predicted fatalities at a location. It has been deprecated. Use simFatalCPr() instead
#' @param Emin Total number of eagle minutes observed in pre-construction monitoring
#' @param SmpHrKM2 Total surveyed hours * area surveyed
#' @param ExpFac Expansion factor
#' @param CPrPriors Collision probability priors, either USFWS or Bay. If neither, custom priors must be provided.
#' @param aPriExp Alpha parameter for the exposure prior. Defaults to 0.9776543.
#' @param bPriExp Beta parameter for the exposure prior. Defaults to 2.777427.
#' @param aPriCpr Alpha parameter for the exposure prior. Must be supplied if neither USFWS nor Bay priors are selected.
#' @param bPriCpr Beta parameter for the exposure prior. Must be supplied if neither USFWS nor Bay priors are selected.
#' @export
#' @import rv


### Fatality Functions - 23 Apr 2013 ###
simFatal<-function(EMin,SmpHrKM2,ExpFac,CPrPriors = "USFWS",
  aPriExp=0.9776543,bPriExp=2.777427,aPriCPr=NULL, bPriCPr=NULL){

  if(CPrPriors == "USFWS"){
    aPriCPr=2.31
    bPriCPr=396.69
  } else if(CPrPriors == "Bay"){
    aPriCPr<-9.28
    bPriCPr<-3224.51
  }


# Update the exposure prior
  aPostExp<-aPriExp+EMin
  bPostExp<-bPriExp+SmpHrKM2
  print(c(aPostExp=aPostExp,bPostExp=bPostExp,aPriCPr=aPriCPr,bPriCPr=bPriCPr))

  Exp<-rvgamma(n=1,aPostExp,bPostExp)
  CPr<-if(bPriCPr==0){
   aPriCPr
  } else {
   rvbeta(n=1,aPriCPr,bPriCPr)
  }
  Fatalities<-ExpFac*Exp*CPr
  attr(Fatalities,"Exp")<-cbind(Mean=rvmean(Exp),SD=rvsd(Exp))
  return(Fatalities)
}

#' @title Fatality simulator with collison probaility
#' @author Joe Statwick (Arcadis)
#' @author Mark Otto (FWS)
#'
#' @description This function calculates the annual predicted fatalities at a location, when post construction mortality monitoring is completed.
#' @param Emin Total number of eagle minutes observed in pre-construction monitoring. Entering a negative value indicates no surveys performed.
#' @param SmpHrKM2 Total surveyed hours * area surveyed in pre-construction monitoring
#' @param EOutmin Total number of eagle carcasses observed in post-construction mortality monitoring. Entering a negative value indicates no surveys performed.
#' @param CprHrKM2 Total surveyed hours * hazardous area surveyed in post-construction mortality monitoring
#' @param ExpFac Expansion factor
#' @param CPrPriors Collision probability priors, either USFWS or Bay. If neither, custom priors must be provided.
#' @param aPriExp Alpha parameter for the exposure prior. Defaults to 0.9776543.
#' @param bPriExp Beta parameter for the exposure prior. Defaults to 2.777427.
#' @param aPriCpr Alpha parameter for the exposure prior. Must be supplied if neither USFWS nor Bay priors are selected.
#' @param bPriCpr Beta parameter for the exposure prior. Must be supplied if neither USFWS nor Bay priors are selected.
#' @export
#' @import rv

simFatalCPr <- function(EMin, EOutMin, SmpHrKM2, CprHrKM2, ExpFac, CPrPriors = "USFWS", aPriExp=0.9776543,bPriExp=2.777427, aPriCPr=NULL, bPriCPr=NULL){

  # EMin: observed number of eagle minutes
  # EOutMin: annual eagle fatalities on an operational wind facility
  # SmpHrKM2:total hr km2 surveyed for eagle minutes
  # ExpFac:expansion factor
  # aPriExp: alpha parameter for the prior on lambda
  # bPriExp: beta parameter for the prior on lambda
  # aPriCPr: alpha parameter for the prior on C
  # bPriCPr: beta parameter for the prior on C

  # Entering a negative value for EMin or EOutMin indicates that no data
  # were collected for those model inputs

  if(CPrPriors == "USFWS"){
    aPriCPr=2.31
    bPriCPr=396.69
  } else if(CPrPriors == "Bay"){
    aPriCPr<-9.28
    bPriCPr<-3224.51
  }

  # Update the exposure prior
  if(EMin>=0){
    aPostExp <- aPriExp + EMin
    bPostExp <- bPriExp + SmpHrKM2
  }else{
    aPostExp <- aPriExp
    bPostExp <- bPriExp}

  Exp <- rvgamma(n=1, aPostExp, bPostExp)

  nsurveys<-CprHrKM2*Exp #hazardous space-time surveyed times exposure

  # Update the collisions prior
  if(EOutMin>=0){
    aPostCPr <- aPriCPr + EOutMin
    bPostCPr <- bPriCPr + nsurveys
  }else{
    aPostCPr <- aPriCPr
    bPostCPr <- bPriCPr}

  CPr <- rvbeta(n=1, aPostCPr, bPostCPr)

  Fatalities             <- ExpFac * Exp * CPr
  attr(Fatalities,"Exp") <- c(Mean=rvmean(Exp), SD=rvsd(Exp))
  attr(Fatalities,"CPr") <- c(Mean=rvmean(CPr), SD=rvsd(CPr))

  return(Fatalities)}

#' @title Fatality plotter
#' @author Mark Otto (FWS)
#'
#' @description This function plots the annual predicted fatalities at a location.
#' @param Fatalities Predicted distribution of fatalities, as generated by simFatalCpr
#' @param probs Desired upper confidence limit(s) on the mean estimate of fatalities. Defaults to 0.8
#' @param ... further arguments passed to the function hist
#' @export
#' @import rv

plotFatal<-function(Fatalities,probs=0.8,PlotHist=TRUE,
                    xlim=NULL,xlab="Collisions",ylab="Density",# bty="o",
                    col="red",add=FALSE,...){

  Names<-if(is.null(names(Fatalities))) 1:length(Fatalities) else
    names(Fatalities)
  Smry<-RVSmry(Names,Fatalities,probs=probs)
  ColIdx<-grepl("CI",colnames(Smry))
  CIs<-Smry[,ColIdx]

  if(!add){
    if(is.null(xlim)) xlim<-c(0,1.1*rvquantile(Fatalities,probs=0.99))
    rvhist(Fatalities,xlab=xlab,ylab=ylab,
           xlim=xlim,freq=FALSE,# bty=bty,
           ...
    )
  }

  lines(density(as.numeric(Fatalities[[1]],bw="sj")),col=if(add) col else "blue")
  abline(v=Smry$Mean,col=if(add) col else "black")
  abline(v=CIs,col=col)
  invisible(NULL)
}