vignettes/westBrook/elliottGrowthEstimates.R

In evanchildress/perform:

library(perform)
library(plotHacks)
reconnect()

rivers<-c("wb jimmy","wb mitchell","wb obear","west brook")
#r<-"wb jimmy"
r<-"west brook"
sp<-"bnt"

    core<-createCoreData("electrofishing",columnsToAdd="observedWeight") %>%
      data.table() %>%
      .[,n:=.N,by=tag] %>%
      .[n>1] %>%
      .[,n:=NULL] %>%
      data.frame() %>%
      addTagProperties() %>%
      createCmrData() %>%
      addKnownZ() %>%
      filter(knownZ==1) %>%
      fillSizeLocation(size=F) %>%
      addSampleProperties() %>%
      filter(river==r&species==sp) %>%
      addEnvironmental(funName="mean") %>%
      rename(medianFlow=meanFlow)%>%
      select(-meanTemperature,-lagDetectionDate) %>%
      data.table() %>%
      .[,diffSample:=c(NA,diff(sampleIndex)),by=tag]

    if(sp=="bnt"){
      core[tag=="257c59a7cc"&observedLength==258,observedLength:=NA]
      core[tag=="1c2c582218"&observedLength==240,observedLength:=NA]
      core[tag=="00088d22ae"&observedLength==70,observedLength:=NA]
      core[tag=="257c59b698"&observedLength==117,observedLength:=NA]
      core[tag=="1bf0fec1d7"&observedLength==66,observedLength:=NA]
    }
    if(sp=="bkt"){
      core[tag=="257c67ca48"&observedLength==224,observedLength:=NA]
    }

    movers<-unique(core[diffSample>1,tag])
    for(t in movers){
      core[tag==t&
             sampleIndex %in%
             sampleIndex[min(which(diffSample>1)):length(sampleIndex)],
           observedLength:=NA]
    }
    core[,diffSample:=NULL]
    core<-core[,nMeasures:=sum(!is.na(observedLength)),by=tag] %>%
      .[nMeasures>1] %>%
      .[,nMeasures:=NULL] %>%
      .[,lastObs:=detectionDate==max(detectionDate),by=tag] %>%
      .[!lastObs|!is.na(observedLength)] %>%
      .[,lastObs:=NULL] %>%
      .[,firstDate:=min(detectionDate[!is.na(observedLength)]),tag] %>%
      .[detectionDate>=firstDate] %>%
      .[,firstDate:=NULL]

    t<-temp %>%
      .[river==r] %>%
      .[,date:=as.Date(datetime)]%>%
      .[datetime>=min(core$detectionDate)&datetime<=max(core$detectionDate)] %>%
      setkey(datetime)

    core[,':='(time=which.min(
      abs(
        t$datetime-as.POSIXct(paste0(detectionDate," 12:00:00"))
      )
    )),
    by=.(tag,detectionDate)]

    lengthWeight<-lm(log(observedLength)~log(observedWeight),data=core)
    core<-core[!is.na(observedWeight)]
    core[,timeLagged:=shift(time),by=tag]
    core[,wExp:=getElliottGrowth(timeLagged,time,observedWeight,tempData=t),by=tag]
    core[,lExp:=exp(predict(lengthWeight,data.frame(observedWeight=wExp)))]
    setkey(core,river,year,season)


    bktBiomass<-readRDS("vignettes/westBrook/bktBiomass.rds")
    bntBiomass<-readRDS("vignettes/westBrook/bntBiomass.rds")

    core<-bktBiomass[core] %>%
      bntBiomass[.] %>%
      setkey(tag,detectionDate)
    core[,residual:=observedWeight-wExp]

    gr<-core[,.(gObsWeight=diff(observedWeight),
                 gExpWeight=wExp[2:length(wExp)]-observedWeight[1:(length(wExp)-1)],
                 gObs=diff(observedLength),
                 gExp=lExp[2:length(lExp)]-observedLength[1:(length(lExp)-1)],
                 startLength=observedLength[1:(length(wExp)-1)],
                 startDate=detectionDate[1:(length(wExp)-1)],
                 endDate=detectionDate[2:length(wExp)]),
              by=tag]
    lengthSumStats<-core[!is.na(observedLength)&!is.na(lExp),
               .(relativeBias=sum(lExp-observedLength)/.N/mean(observedLength),
                 rmse=sqrt(sum((lExp-observedLength)^2)/.N))]
    weightSumStats<-core[!is.na(observedWeight)&!is.na(wExp),
                 .(relativeBiasWeight=sum(wExp-observedWeight)/.N/mean(observedWeight),
                 rmseWeight=sqrt(sum((wExp-observedWeight)^2)/.N))]
    growthSumStats<-gr[!is.na(gObs)&!is.na(gExp)&!is.na(gObsWeight)&!is.na(gExpWeight),
                         .(relativeBias=sum(gExp-gObs)/.N/mean(gObs),
                           rmse=sqrt(sum((gExp-gObs)^2)/.N),
                           relativeBiasWeight=sum(gExpWeight-gObsWeight)/.N/mean(gObsWeight),
                           rmseWeight=sqrt(sum((gExpWeight-gObsWeight)^2)/.N))]

    tiff.par("vignettes/westBrook/results/figures/elliotGrowth.tif")
    plot(gExp~gObs,data=gr,pch=16,col=gray(0.5,0.5),
         ylab="Elliott Model Predicted Growth (g)",
         xlab="Observed Growth (g)")
    abline(0,1)
dev.off()