In natureanalytics-ca/BycatchEstimator: Generic Model-Based Bycatch Estimation Procedure

#outDir<-"C:/Users/Elizabeth A Babcock/Box/Shrimp trawl bycatch/bycatchEstimatorRuns/Output SawfishSetNoPool"
setupObj<-readRDS(file=paste0(outDir,"/",Sys.Date(),"_BycatchModelSpecification.rds"))
library(tidyverse)
library(kableExtra)
theme_set(theme_bw())

#Unpack setupObj
  modelTry<-obsdat<-logdat<-yearVar<-obsEffort<-logEffort<-logUnsampledEffort<-
    includeObsCatch<-matchColumn<-factorNames<-randomEffects<-randomEffects2<-
    EstimateIndex<-EstimateBycatch<-logNum<-sampleUnit<-complexModel<-simpleModel<-indexModel<-
    designMethods<-designVars<-designPooling<-poolTypes<-pooledVar<-adjacentNum<-
    minStrataUnit<-
    baseDir<-runName<-runDescription<-
    common<-sp<-obsCatch<-catchUnit<-catchType<-NULL

  numSp<-modelTable<-modelSelectTable<-modFits<-modPredVals<-modIndexVals<-residualTab<-bestmod<-predbestmod<-indexbestmod<-allmods<-allindex<-
    modelFail<-rmsetab<-metab<-dat<-yearSum<-requiredVarNames<-allVarNames<-indexDat<-strataSum<-NumCores<-NULL

  for(r in 1:NROW(setupObj$bycatchInputs)) assign(names(setupObj$bycatchInputs)[r], setupObj$bycatchInputs[[r]])
  for(r in 1:NROW(setupObj$bycatchOutputs)) assign(names(setupObj$bycatchOutputs)[r],setupObj$bycatchOutputs[[r]])

numericalVars<-c("Effort",allVarNames[!allVarNames %in% factorNames])
obsdat<-mutate(obsdat,SampleUnits=1)
##Function to plot by species
makeObserverPlots<-function(spNum,figurenum) {
## Show where species was caught in observer data
obsdatTemp <-obsdat %>% rename(Catch=!!obsCatch[spNum]) %>%
  mutate(Present=factor(ifelse(Catch>0,1,0)),
         CPUE=Catch/Effort)
# Presence/absence across years
cat("\n Figure ", figurenum,". Presence and absence of",common[spNum], "in each year. \n",sep="")
print(ggplot(obsdatTemp,aes(x=Year,fill=Present))+
  geom_bar()+
  scale_fill_manual(values=c("grey","red"))+
    ylab(paste0("Count of observed ",sampleUnit)))
figurenum<-figurenum+1
# Presence/absence across all factors
cat("\n")
cat("\n Figure ", figurenum,". Presence and absence of ",common[spNum], " across factor variables. \n",sep="")
temp<-obsdatTemp %>% pivot_longer(cols=all_of(factorNames),names_to="Variable",values_to = "Level") 
print(ggplot(temp,aes(x=Level,fill=Present))+
  geom_bar()+
  facet_wrap(Variable~.,scales="free")+
  scale_fill_manual(values=c("grey","red"))+
    ylab(paste0("Count of observed ",sampleUnit)))
figurenum<-figurenum+1

# Presence/absence across numerical variables and effort
cat("\n")
cat("\n Figure ", figurenum,". Presence and absence of ",common[spNum], " across numerical variables. \n",sep="")
temp<-obsdatTemp %>% pivot_longer(cols=all_of(c("Effort",numericalVars)),names_to="Variable",values_to = "Value") 
print(ggplot(temp,aes(x=Value,fill=Present))+
  geom_histogram()+
  facet_wrap(Variable~.,scales="free")+
  scale_fill_manual(values=c("grey","red"))+
    ylab(paste0("Count of observed ",sampleUnit)))
figurenum<-figurenum+1

#CPUE
cat("\n")
cat("\n Figure ",figurenum,". Observed CPUE of ",common[spNum], " by factor levels.\n",sep="")
temp<-obsdatTemp %>% pivot_longer(cols=all_of(factorNames),names_to="Variable",values_to = "Level") 
print(ggplot(temp,aes(x=Level,y=CPUE))+
  geom_violin(fill="grey")+
  facet_wrap(Variable~.,scales="free")+
  stat_summary())
figurenum<-figurenum+1
cat("\n")
cat("\n Figure ",figurenum,". Observed CPUE of",common[spNum], " across numerical variables. \n",sep="")
temp<-obsdatTemp %>% pivot_longer(cols=all_of(numericalVars),names_to="Variable",values_to = "Value") 
print(ggplot(temp,aes(x=Value,y=CPUE))+
  geom_point(alpha=0.3)+
  stat_smooth()+
  facet_wrap(Variable~.,scales="free"))
  figurenum
}

#Compare observed and logbook  
makeEffortPlots<-function(figurenum) {
allData<-bind_rows(list(observer=select(obsdat,all_of(c(allVarNames,"Effort","SampleUnits"))),
                  effort=select(logdat,all_of(c(allVarNames,"Effort","SampleUnits")))),
                  .id = "Source")%>%
  pivot_longer(cols=all_of(factorNames),names_to="Variable",values_to = "Level") 

cat("\n")
cat("\n Figure ",figurenum,". Total Effort by factor levels, observed and total. \n",sep="")
print(ggplot(allData,aes(x=Level,weight=Effort,fill=Source))+
  geom_bar(position=position_dodge())+
  facet_wrap(Variable~.,ncol=2,scales="free")+
  scale_fill_manual(values=c("grey","black"))+
    ylab(paste0("Sum of effort")))
figurenum<-figurenum+1
cat("\n")
cat("\n Figure ",figurenum,". Count of ",sampleUnit," by factor levels, observed and total. \n",sep="")
print(ggplot(allData,aes(x=Level,weight=SampleUnits,fill=Source))+
  geom_bar(position=position_dodge())+
  facet_wrap(Variable~.,ncol=2,scales="free")+
  scale_fill_manual(values=c("grey","black"))+
    ylab(paste0("Count of ",sampleUnit)))
figurenum<-figurenum+1

logdatTemp<-logdat %>%
  mutate(Effort=Effort/SampleUnits) %>%
  uncount(SampleUnits)

allData<-bind_rows(list(observer=select(obsdat,all_of(c(allVarNames,"Effort"))),
                  effort=select(logdatTemp,all_of(c(allVarNames,"Effort")))),
                  .id = "Source")%>%
  pivot_longer(cols=all_of(numericalVars),names_to="Variable",values_to = "Value")
cat("\n")
cat("\n Figure ",figurenum,". Effort across numerical variables. \n",sep="")
print(ggplot(allData,aes(x=Value,fill=Source))+
  geom_histogram(position=position_dodge())+
  facet_wrap(Variable~.,ncol=2,scales="free")+
  scale_fill_manual(values=c("grey","black"))+
    ylab(paste0("Count of ",sampleUnit)))
figurenum<-figurenum+1
figurenum
}

Summary of data for r runDescription , r Sys.Date()

if(!"Year" %in% factorNames) {
  obsdat$Year<-obsdat$Year+startYear
  logdat$Year<-logdat$Year+startYear
}
fignum<-makeEffortPlots(1)

Table 1. Observed and unobserved effort by r sampleUnit.

temp<-bind_rows(list(observer=select(obsdat,Effort),
                  effort=select(logdat,Effort)),
                  .id = "Source") %>% group_by(Source) %>%
  summarize(`Sample units` = n(),
            `Total effort`=sum(Effort))%>%
  mutate(`Proportion units`=round(`Sample units`/sum(`Sample units`),3),
         `Proportion effort`=round(`Total effort`/sum(`Total effort`),3)) 
kbl(temp,format="simple")

for(i in 1:numSp) {
 fignum<-makeObserverPlots(i,fignum)
}

natureanalytics-ca/BycatchEstimator documentation built on Oct. 15, 2024, 10:28 p.m.