vignettes/westBrook/fitWbLengthModels.R
In evanchildress/perform:
library(perform)
reconnect()
rivers<-c("wb jimmy","wb mitchell","wb obear","west brook")
#r<-"wb jimmy"
for(r in c("west brook")){
for(sp in c("bnt")){
# for(sp in c("ats")){
if(sp=="bnt"&r=="wb obear") next
priors<-list(bkt=list(tOptMean=0,#14.2,
tOptPrecision=0.0001,
ctMaxMean=25,#23.4,
ctMaxPrecision=0.0001,
ctUltimate=40),#25.3),
bnt=list(tOptMean=0,#15.95,
tOptPrecision=0.0001,
ctMaxMean=25,#22.5,
ctMaxPrecision=0.0001,
ctUltimate=40)#28)
)
endDate<-as.POSIXct("2016-10-01")
core<-createCoreData("electrofishing",columnsToAdd="observedWeight") %>%
data.table() %>%
.[,n:=.N,by=tag] %>%
.[n>1] %>%
.[,n:=NULL] %>%
data.frame() %>%
addTagProperties() %>%
createCmrData(dateEnd = endDate) %>%
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==118,observedLength:=218]
# core[tag=="00088cfb9e"&observedLength==223,observedLength:=NA]
# core[tag=="00088d215d"&observedLength==67,observedLength:=NA]
# core[tag=="0009f6f4d5"&observedLength==115,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]
core[,tagIndex:=match(tag,unique(tag))] %>%
setkey(river,year,season)
#
bktBiomass<-readRDS("vignettes/westBrook/bktBiomass.rds")
bntBiomass<-readRDS("vignettes/westBrook/bntBiomass.rds")
core<-bktBiomass[core] %>%
bntBiomass[.] %>%
setkey(tag,detectionDate)
# if(r=="wb obear"){
# core[,bntBiomass:=0]
# }
if(r=="wb mitchell"){
core[is.na(bntBiomass),bntBiomass:=0]
}
# if(sp=="ats"){
# core[,":="(meanBktBiomass=mean(bktBiomass,na.rm=T),
# meanBntBiomass=mean(bntBiomass,na.rm=T)),by=season] %>%
# .[is.na(bktBiomass),bktBiomass:=meanBktBiomass] %>%
# .[is.na(bntBiomass),bntBiomass:=meanBntBiomass] %>%
# .[,":="(meanBktBiomass=NULL,
# meanBntBiomass=NULL)]
# }
#core<-core[!tag %in% c("00088d1ac1","00088d2f6f")]
# .[,tagIndex:=match(tag,unique(tag))] %>%
# setkey(tagIndex,detectionDate) %>%
# .[,firstObs:=detectionDate==min(detectionDate),by=tag]
jagsData<-createJagsData(data.frame(core) %>%
addEnvironmental()) %>%
.[c("firstObsRows",
"nFirstObsRows",
"evalRows",
"nEvalRows",
"nAllRows",
"lengthDATA",
"ind",
"season")]
core<-core[,.(tag,tagIndex,detectionDate,observedLength,
bktBiomass,bntBiomass,medianFlow)]
t<-temp %>%
.[river==r] %>%
.[,date:=as.Date(datetime)]%>%
.[datetime>=min(core$detectionDate)&datetime<=max(core$detectionDate)] %>%
setkey(datetime)
#setting an arbitrary temperature for a period of NAs,
#but growth over this period is excluded from the likelihood in the model
if(sp=="ats"&r=="west brook"){
t[is.na(temperature),temperature:=23]
}
core[,':='(time=which.min(
abs(
t$datetime-as.POSIXct(paste0(detectionDate," 12:00:00"))
)
)),
by=.(tag,detectionDate)]
# t[,month:=ceiling((month(date))/3)]
# hoursPerMonth<-t[date>=as.Date("2003-01-01")&date<=as.Date("2014-12-31"),
# .(hours=.N/length(unique(year(date)))),
# by=month] %>%
# setkey(month)
#
# propMonth<-core[,.(tag,time)] %>%
# .[,startTime:=as.numeric(shift(time)),by=tag] %>%
# .[,obs:=1:nrow(.)] %>%
# .[is.na(startTime),startTime:=time-1]
#
# propMonth<-propMonth[!is.na(startTime),t[startTime:time,.N,by=month] %>%
# setkey(month) %>%
# .[hoursPerMonth] %>%
# .[is.na(N),N:=0] %>%
# .[,propMonth:=N/hours] %>%
# .[,.(propMonth,month)],
# by=.(obs)] %>%
# melt(id.vars=c("month","obs")) %>%
# acast(obs~month)
# #
# jagsData$time<-core$time<-list(lengthDATA=core$observedLength,
# firstObsRows=which(core$firstObs==1),
# nFirstObsRows=length(which(core$firstObs==1)),
# evalRows=which(core$firstObs==0),
# nEvalRows=length(which(core$firstObs==0)),
# propMonth=propMonth,
# tempDATA=t$temperature,
# nTimes=nrow(t),
# time=core$time,
# nMonths=max(hoursPerMonth$month)
# )
# jagsData$propMonth<-propMonth
# jagsData$nMonths<-max(hoursPerMonth$month)
jagsData$tempDATA<-t$temperature
jagsData$nTimes<-nrow(t)
jagsData$time<-core$time
jagsData$nInd<-max(core$tagIndex)
jagsData$isSpring<-as.numeric(jagsData$season==2)
jagsData$bktBiomassDATA<-scale(core$bktBiomass)[,1]
jagsData$bntBiomassDATA<-scale(core$bntBiomass)[,1]
jagsData$biomassDATA<-scale(core[[paste0(sp,"Biomass")]])[,1]
jagsData$flowDATA<-scale(core$medianFlow)[,1]
jagsData$tOptMean<-priors[[sp]][["tOptMean"]]
jagsData$tOptPrecision<-priors[[sp]][["tOptPrecision"]]
jagsData$ctMaxMean<-priors[[sp]][["ctMaxMean"]]
jagsData$ctMaxPrecision<-priors[[sp]][["ctMaxPrecision"]]
jagsData$ctUltimate<-priors[[sp]][["ctUltimate"]]
jagsData$meanFlow<-mean(core$medianFlow,na.rm=T)
jagsData$sdFlow<-sd(core$medianFlow,na.rm=T)
jagsData$meanBktBiomass<-mean(core$bktBiomass,na.rm=T)
jagsData$sdBktBiomass<-sd(core$bktBiomass,na.rm=T)
jagsData$meanBntBiomass<-mean(core$bntBiomass,na.rm=T)
jagsData$sdBntBiomass<-sd(core$bntBiomass,na.rm=T)
if(r=="wb obear"){
jagsData$bntBiomassDATA<-rep(0,nrow(core))
}
if(r=="wb mitchell"){
badStarts<-as.POSIXct("2007-08-02 19:00:00","2013-05-01 00:00:00")
badEnds<-as.POSIXct("2007-08-04 00:00:00","2013-10-01 00:00:00")
for(b in 1:length(badStarts)){
st<-which(t$datetime==badStarts[b])
en<-which(t$datetime==badEnds[b])
jagsData$evalRows<-jagsData$evalRows[
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]<st)|
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]>en)]
}
jagsData$nEvalRows<-length(jagsData$evalRows)
}
if(r=="west brook" & sp=="ats"){
badStarts<-as.POSIXct("1998-06-24 20:00:00")
badEnds<-as.POSIXct("1998-08-21 16:00:00")
for(b in 1:length(badStarts)){
st<-which(t$datetime==badStarts[b])
en<-which(t$datetime==badEnds[b])
jagsData$evalRows<-jagsData$evalRows[
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]<st)|
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]>en)]
}
jagsData$nEvalRows<-length(jagsData$evalRows)
}
core[,lengthInit:=approx(observedLength,n=length(observedLength))$y,by=tag]
core[!is.na(observedLength),lengthInit:=NA]
inits<-function(){list(lengthData=core$lengthInit,
beta2= -5e-05,
eps=0.003)
}
parsToMonitor<-c("beta1","beta2","beta3","beta4","beta5",
"tOpt",'ctMax',"sigma","ranInd","ranSlope",
'eps',"sigmaInd","lengthExp","b","gamma","psi")
out<-fitModel(jagsData=jagsData,inits=inits,parallel=T,params=parsToMonitor,
nb=8000,ni=18000,nt=10,modelFile="modelLengthField.R",codaOnly=c("lengthExp","ranInd","ranSlope","beta1"))
# saveRDS(out,file=paste0("vignettes/westBrook/results/out",sp,toupper(substr(r,1,1)),substr(r,2,nchar(r)),".rds"))
# saveRDS(core,file=paste0("vignettes/westBrook/results/core",sp,toupper(substr(r,1,1)),substr(r,2,nchar(r)),".rds"))
print(out)
assign(paste0("out",sp,which(r==rivers)),out)
assign(paste0("core",sp,which(r==rivers)),core)
# core[jagsData$evalRows,predictedLength:=apply(out$sims.list$lengthExp,2,mean)]
# core[,residual:=observedLength-predictedLength]
}
}
#
# for(r in c("wb jimmy","wb mitchell","wb obear","west brook")){
# for(sp in c("bkt")){
# if(sp=="bnt"&r=="wb obear") next
#
# out<-get(paste0("out",sp,which(r==rivers)))
# plot(NA,xlim=c(0,22),ylim=c(-1,1),main=paste(r,sp),xlab="temp",ylab="performance")
# for(i in sample(1:length(out$sims.list$tOpt),300,replace=T)){
# points(predictPerformance(0:22,tOpt=out$sims.list$tOpt[i],
# sigma=out$sims.list$sigma[i],
# ctMax=out$sims.list$ctMax[i])~I(0:22),
# type='l',col='gray')
# }
#
# core<-get(paste0("core",sp,which(r==rivers)))
# core$dummy<-1
# core[,firstObs:=detectionDate==min(detectionDate),by=tag] %>%
# .[firstObs==F,predictedLength:=apply(out$sims.list$lengthExp,2,mean)]
#
# assign(paste0("gr",sp,which(r==rivers)),
# core[,.(obsGrowth=diff(observedLength),
# predGrowth=predictedLength[2:length(observedLength)]-
# observedLength[1:(length(observedLength)-1)],
# date=detectionDate[2:length(detectionDate)]),
# by=tag] %>%
# .[,residual:=obsGrowth-predGrowth])
# assign(paste0("core",sp,which(r==rivers)),core)
#
# plot(predictedLength~observedLength,data=get(paste0("core",sp,which(r==rivers))))
# a<-lm(predictedLength~observedLength,get(paste0("core",sp,which(r==rivers))))
# text(75,150,bquote(R^2==.(round(summary(a)$r.squared,2))))
#
# plot(obsGrowth~predGrowth,data=get(paste0("gr",sp,which(r==rivers))),pch=19,col=gray(0.45,0.5))
# a<-lm(obsGrowth~predGrowth,get(paste0("gr",sp,which(r==rivers))))
# abline(a)
# abline(0,1,lty=2)
# text(5,15,bquote(R^2==.(round(summary(a)$r.squared,2))),col='black')
#
#
# # plot(mmPerDay~startLength,data=get(paste0("gr",which(r==rivers))),
# # col=gray(0.5,0.5),pch=19,main=r,xlab="startLength",ylab="growth (mm/day)")
# # x<-seq(60,300)
# # for(i in sample(1:length(out$sims.list$tOpt),300,replace=T)){
# # points(predictVonBert(x,
# # out$sims.list$beta[i,1],
# # out$sims.list$beta[i,2],
# # derivative=T)*24~x,
# # type='l',col='blue')
# # }
# #
# #
# #
# # predicted<-apply(out$sims.list$length,2,mean)
# # plot(predicted~get(paste0("gr",which(r==rivers)))$growth,main=r,
# # ylab="predicted growth",xlab="observed growth")
# # abline(0,1)
# }
# }
#
# ###code to estimate bias and mse
# bktSummary<-corebkt4[!is.na(predictedLength)&!is.na(observedLength),
# .(rmse=sqrt(sum(((observedLength-predictedLength))^2)/.N),
# relativeBias=sum((observedLength-predictedLength)/observedLength)/.N)]
# bntSummary<-corebnt4[!is.na(predictedLength)&!is.na(observedLength),
# .(rmse=sqrt(sum(((observedLength-predictedLength))^2)/.N),
# relativeBias=sum((observedLength-predictedLength)/observedLength)/.N)]
# bktGrowthSummary<-grbkt4[!is.na(predGrowth)&!is.na(obsGrowth),
# .(rmse=sqrt(sum(((obsGrowth-predGrowth))^2)/.N),
# relativeBias=sum((obsGrowth-predGrowth))/.N/mean(obsGrowth))]
# bntGrowthSummary<-grbnt4[!is.na(predGrowth)&!is.na(obsGrowth),
# .(rmse=sqrt(sum(((obsGrowth-predGrowth))^2)/.N),
# relativeBias=sum((obsGrowth-predGrowth))/.N/mean(obsGrowth))]
#
#
# plot(NA,xlim=c(0,25),ylim=c(-1.5,1))
# for(i in sample(1:length(out$sims.list$tOpt),300,replace=T)){
# points(predictPerformance(seq(0,23.3,length.out=100),tOpt=outbkt1$sims.list$tOpt[i],
# sigma=outbkt1$sims.list$sigma[i],
# ctMax=outbkt1$sims.list$ctMax[i])~I(seq(0,23.3,length.out=100)),
# type='l',col=palette()[1])
#
# points(predictPerformance(seq(0,25.4,length.out=100),tOpt=outbkt2$sims.list$tOpt[i],
# sigma=outbkt2$sims.list$sigma[i],
# ctMax=outbkt2$sims.list$ctMax[i])~I(seq(0,24.4,length.out=100)),
# type='l',col=palette()[2])
#
# points(predictPerformance(seq(0,21.2,length.out=100),tOpt=outbkt3$sims.list$tOpt[i],
# sigma=outbkt3$sims.list$sigma[i],
# ctMax=outbkt3$sims.list$ctMax[i])~I(seq(0,21.2,length.out=100)),
# type='l',col=palette()[3])
# }
#
evanchildress/perform documentation built on May 16, 2019, 9:35 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(perform)
reconnect()
rivers<-c("wb jimmy","wb mitchell","wb obear","west brook")
#r<-"wb jimmy"
for(r in c("west brook")){
for(sp in c("bnt")){
# for(sp in c("ats")){
if(sp=="bnt"&r=="wb obear") next
priors<-list(bkt=list(tOptMean=0,#14.2,
tOptPrecision=0.0001,
ctMaxMean=25,#23.4,
ctMaxPrecision=0.0001,
ctUltimate=40),#25.3),
bnt=list(tOptMean=0,#15.95,
tOptPrecision=0.0001,
ctMaxMean=25,#22.5,
ctMaxPrecision=0.0001,
ctUltimate=40)#28)
)
endDate<-as.POSIXct("2016-10-01")
core<-createCoreData("electrofishing",columnsToAdd="observedWeight") %>%
data.table() %>%
.[,n:=.N,by=tag] %>%
.[n>1] %>%
.[,n:=NULL] %>%
data.frame() %>%
addTagProperties() %>%
createCmrData(dateEnd = endDate) %>%
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==118,observedLength:=218]
# core[tag=="00088cfb9e"&observedLength==223,observedLength:=NA]
# core[tag=="00088d215d"&observedLength==67,observedLength:=NA]
# core[tag=="0009f6f4d5"&observedLength==115,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]
core[,tagIndex:=match(tag,unique(tag))] %>%
setkey(river,year,season)
#
bktBiomass<-readRDS("vignettes/westBrook/bktBiomass.rds")
bntBiomass<-readRDS("vignettes/westBrook/bntBiomass.rds")
core<-bktBiomass[core] %>%
bntBiomass[.] %>%
setkey(tag,detectionDate)
# if(r=="wb obear"){
# core[,bntBiomass:=0]
# }
if(r=="wb mitchell"){
core[is.na(bntBiomass),bntBiomass:=0]
}
# if(sp=="ats"){
# core[,":="(meanBktBiomass=mean(bktBiomass,na.rm=T),
# meanBntBiomass=mean(bntBiomass,na.rm=T)),by=season] %>%
# .[is.na(bktBiomass),bktBiomass:=meanBktBiomass] %>%
# .[is.na(bntBiomass),bntBiomass:=meanBntBiomass] %>%
# .[,":="(meanBktBiomass=NULL,
# meanBntBiomass=NULL)]
# }
#core<-core[!tag %in% c("00088d1ac1","00088d2f6f")]
# .[,tagIndex:=match(tag,unique(tag))] %>%
# setkey(tagIndex,detectionDate) %>%
# .[,firstObs:=detectionDate==min(detectionDate),by=tag]
jagsData<-createJagsData(data.frame(core) %>%
addEnvironmental()) %>%
.[c("firstObsRows",
"nFirstObsRows",
"evalRows",
"nEvalRows",
"nAllRows",
"lengthDATA",
"ind",
"season")]
core<-core[,.(tag,tagIndex,detectionDate,observedLength,
bktBiomass,bntBiomass,medianFlow)]
t<-temp %>%
.[river==r] %>%
.[,date:=as.Date(datetime)]%>%
.[datetime>=min(core$detectionDate)&datetime<=max(core$detectionDate)] %>%
setkey(datetime)
#setting an arbitrary temperature for a period of NAs,
#but growth over this period is excluded from the likelihood in the model
if(sp=="ats"&r=="west brook"){
t[is.na(temperature),temperature:=23]
}
core[,':='(time=which.min(
abs(
t$datetime-as.POSIXct(paste0(detectionDate," 12:00:00"))
)
)),
by=.(tag,detectionDate)]
# t[,month:=ceiling((month(date))/3)]
# hoursPerMonth<-t[date>=as.Date("2003-01-01")&date<=as.Date("2014-12-31"),
# .(hours=.N/length(unique(year(date)))),
# by=month] %>%
# setkey(month)
#
# propMonth<-core[,.(tag,time)] %>%
# .[,startTime:=as.numeric(shift(time)),by=tag] %>%
# .[,obs:=1:nrow(.)] %>%
# .[is.na(startTime),startTime:=time-1]
#
# propMonth<-propMonth[!is.na(startTime),t[startTime:time,.N,by=month] %>%
# setkey(month) %>%
# .[hoursPerMonth] %>%
# .[is.na(N),N:=0] %>%
# .[,propMonth:=N/hours] %>%
# .[,.(propMonth,month)],
# by=.(obs)] %>%
# melt(id.vars=c("month","obs")) %>%
# acast(obs~month)
# #
# jagsData$time<-core$time<-list(lengthDATA=core$observedLength,
# firstObsRows=which(core$firstObs==1),
# nFirstObsRows=length(which(core$firstObs==1)),
# evalRows=which(core$firstObs==0),
# nEvalRows=length(which(core$firstObs==0)),
# propMonth=propMonth,
# tempDATA=t$temperature,
# nTimes=nrow(t),
# time=core$time,
# nMonths=max(hoursPerMonth$month)
# )
# jagsData$propMonth<-propMonth
# jagsData$nMonths<-max(hoursPerMonth$month)
jagsData$tempDATA<-t$temperature
jagsData$nTimes<-nrow(t)
jagsData$time<-core$time
jagsData$nInd<-max(core$tagIndex)
jagsData$isSpring<-as.numeric(jagsData$season==2)
jagsData$bktBiomassDATA<-scale(core$bktBiomass)[,1]
jagsData$bntBiomassDATA<-scale(core$bntBiomass)[,1]
jagsData$biomassDATA<-scale(core[[paste0(sp,"Biomass")]])[,1]
jagsData$flowDATA<-scale(core$medianFlow)[,1]
jagsData$tOptMean<-priors[[sp]][["tOptMean"]]
jagsData$tOptPrecision<-priors[[sp]][["tOptPrecision"]]
jagsData$ctMaxMean<-priors[[sp]][["ctMaxMean"]]
jagsData$ctMaxPrecision<-priors[[sp]][["ctMaxPrecision"]]
jagsData$ctUltimate<-priors[[sp]][["ctUltimate"]]
jagsData$meanFlow<-mean(core$medianFlow,na.rm=T)
jagsData$sdFlow<-sd(core$medianFlow,na.rm=T)
jagsData$meanBktBiomass<-mean(core$bktBiomass,na.rm=T)
jagsData$sdBktBiomass<-sd(core$bktBiomass,na.rm=T)
jagsData$meanBntBiomass<-mean(core$bntBiomass,na.rm=T)
jagsData$sdBntBiomass<-sd(core$bntBiomass,na.rm=T)
if(r=="wb obear"){
jagsData$bntBiomassDATA<-rep(0,nrow(core))
}
if(r=="wb mitchell"){
badStarts<-as.POSIXct("2007-08-02 19:00:00","2013-05-01 00:00:00")
badEnds<-as.POSIXct("2007-08-04 00:00:00","2013-10-01 00:00:00")
for(b in 1:length(badStarts)){
st<-which(t$datetime==badStarts[b])
en<-which(t$datetime==badEnds[b])
jagsData$evalRows<-jagsData$evalRows[
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]<st)|
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]>en)]
}
jagsData$nEvalRows<-length(jagsData$evalRows)
}
if(r=="west brook" & sp=="ats"){
badStarts<-as.POSIXct("1998-06-24 20:00:00")
badEnds<-as.POSIXct("1998-08-21 16:00:00")
for(b in 1:length(badStarts)){
st<-which(t$datetime==badStarts[b])
en<-which(t$datetime==badEnds[b])
jagsData$evalRows<-jagsData$evalRows[
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]<st)|
(jagsData$time[jagsData$evalRows]&jagsData$time[jagsData$evalRows-1]>en)]
}
jagsData$nEvalRows<-length(jagsData$evalRows)
}
core[,lengthInit:=approx(observedLength,n=length(observedLength))$y,by=tag]
core[!is.na(observedLength),lengthInit:=NA]
inits<-function(){list(lengthData=core$lengthInit,
beta2= -5e-05,
eps=0.003)
}
parsToMonitor<-c("beta1","beta2","beta3","beta4","beta5",
"tOpt",'ctMax',"sigma","ranInd","ranSlope",
'eps',"sigmaInd","lengthExp","b","gamma","psi")
out<-fitModel(jagsData=jagsData,inits=inits,parallel=T,params=parsToMonitor,
nb=8000,ni=18000,nt=10,modelFile="modelLengthField.R",codaOnly=c("lengthExp","ranInd","ranSlope","beta1"))
# saveRDS(out,file=paste0("vignettes/westBrook/results/out",sp,toupper(substr(r,1,1)),substr(r,2,nchar(r)),".rds"))
# saveRDS(core,file=paste0("vignettes/westBrook/results/core",sp,toupper(substr(r,1,1)),substr(r,2,nchar(r)),".rds"))
print(out)
assign(paste0("out",sp,which(r==rivers)),out)
assign(paste0("core",sp,which(r==rivers)),core)
# core[jagsData$evalRows,predictedLength:=apply(out$sims.list$lengthExp,2,mean)]
# core[,residual:=observedLength-predictedLength]
}
}
#
# for(r in c("wb jimmy","wb mitchell","wb obear","west brook")){
# for(sp in c("bkt")){
# if(sp=="bnt"&r=="wb obear") next
#
# out<-get(paste0("out",sp,which(r==rivers)))
# plot(NA,xlim=c(0,22),ylim=c(-1,1),main=paste(r,sp),xlab="temp",ylab="performance")
# for(i in sample(1:length(out$sims.list$tOpt),300,replace=T)){
# points(predictPerformance(0:22,tOpt=out$sims.list$tOpt[i],
# sigma=out$sims.list$sigma[i],
# ctMax=out$sims.list$ctMax[i])~I(0:22),
# type='l',col='gray')
# }
#
# core<-get(paste0("core",sp,which(r==rivers)))
# core$dummy<-1
# core[,firstObs:=detectionDate==min(detectionDate),by=tag] %>%
# .[firstObs==F,predictedLength:=apply(out$sims.list$lengthExp,2,mean)]
#
# assign(paste0("gr",sp,which(r==rivers)),
# core[,.(obsGrowth=diff(observedLength),
# predGrowth=predictedLength[2:length(observedLength)]-
# observedLength[1:(length(observedLength)-1)],
# date=detectionDate[2:length(detectionDate)]),
# by=tag] %>%
# .[,residual:=obsGrowth-predGrowth])
# assign(paste0("core",sp,which(r==rivers)),core)
#
# plot(predictedLength~observedLength,data=get(paste0("core",sp,which(r==rivers))))
# a<-lm(predictedLength~observedLength,get(paste0("core",sp,which(r==rivers))))
# text(75,150,bquote(R^2==.(round(summary(a)$r.squared,2))))
#
# plot(obsGrowth~predGrowth,data=get(paste0("gr",sp,which(r==rivers))),pch=19,col=gray(0.45,0.5))
# a<-lm(obsGrowth~predGrowth,get(paste0("gr",sp,which(r==rivers))))
# abline(a)
# abline(0,1,lty=2)
# text(5,15,bquote(R^2==.(round(summary(a)$r.squared,2))),col='black')
#
#
# # plot(mmPerDay~startLength,data=get(paste0("gr",which(r==rivers))),
# # col=gray(0.5,0.5),pch=19,main=r,xlab="startLength",ylab="growth (mm/day)")
# # x<-seq(60,300)
# # for(i in sample(1:length(out$sims.list$tOpt),300,replace=T)){
# # points(predictVonBert(x,
# # out$sims.list$beta[i,1],
# # out$sims.list$beta[i,2],
# # derivative=T)*24~x,
# # type='l',col='blue')
# # }
# #
# #
# #
# # predicted<-apply(out$sims.list$length,2,mean)
# # plot(predicted~get(paste0("gr",which(r==rivers)))$growth,main=r,
# # ylab="predicted growth",xlab="observed growth")
# # abline(0,1)
# }
# }
#
# ###code to estimate bias and mse
# bktSummary<-corebkt4[!is.na(predictedLength)&!is.na(observedLength),
# .(rmse=sqrt(sum(((observedLength-predictedLength))^2)/.N),
# relativeBias=sum((observedLength-predictedLength)/observedLength)/.N)]
# bntSummary<-corebnt4[!is.na(predictedLength)&!is.na(observedLength),
# .(rmse=sqrt(sum(((observedLength-predictedLength))^2)/.N),
# relativeBias=sum((observedLength-predictedLength)/observedLength)/.N)]
# bktGrowthSummary<-grbkt4[!is.na(predGrowth)&!is.na(obsGrowth),
# .(rmse=sqrt(sum(((obsGrowth-predGrowth))^2)/.N),
# relativeBias=sum((obsGrowth-predGrowth))/.N/mean(obsGrowth))]
# bntGrowthSummary<-grbnt4[!is.na(predGrowth)&!is.na(obsGrowth),
# .(rmse=sqrt(sum(((obsGrowth-predGrowth))^2)/.N),
# relativeBias=sum((obsGrowth-predGrowth))/.N/mean(obsGrowth))]
#
#
# plot(NA,xlim=c(0,25),ylim=c(-1.5,1))
# for(i in sample(1:length(out$sims.list$tOpt),300,replace=T)){
# points(predictPerformance(seq(0,23.3,length.out=100),tOpt=outbkt1$sims.list$tOpt[i],
# sigma=outbkt1$sims.list$sigma[i],
# ctMax=outbkt1$sims.list$ctMax[i])~I(seq(0,23.3,length.out=100)),
# type='l',col=palette()[1])
#
# points(predictPerformance(seq(0,25.4,length.out=100),tOpt=outbkt2$sims.list$tOpt[i],
# sigma=outbkt2$sims.list$sigma[i],
# ctMax=outbkt2$sims.list$ctMax[i])~I(seq(0,24.4,length.out=100)),
# type='l',col=palette()[2])
#
# points(predictPerformance(seq(0,21.2,length.out=100),tOpt=outbkt3$sims.list$tOpt[i],
# sigma=outbkt3$sims.list$sigma[i],
# ctMax=outbkt3$sims.list$ctMax[i])~I(seq(0,21.2,length.out=100)),
# type='l',col=palette()[3])
# }
#
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.