In bletcher/linkedModels: Link movement, growth, and survival models
riverOrderedIn <- factor(c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),labels = c("WB","OL","OS","IL"), ordered = T)
Plots for the manuscript
library(linkedModels)
library(arrayhelpers)
library(tidyverse)
library(gridExtra)
load("D:/projects/linkedModels/data/out/ddD_bktbntats_1997.RData")
load("D:/projects/linkedModels/data/out/rawCountsAndMasses.RData")
mcmc <- list()
dGIn <- list()
ddGIn2 <- list()
load("D:/projects/linkedModels/data/out/dG_1530971738_crossValFALSE_bkt1997_Nimble.RData")
# input data
d1 <- ddG[[1]][[1]]
d2 <- ddG[[1]][[2]]
meanSampleIntervalMean <- data.frame( meanSampleInterval = colMeans(d1$sampleIntervalMean), season = 1:d1$nSeasons )
#RECALCULATE gr for dataset. Filtering by maxSampleInterval in the creat ddddG step leaves in trailing lengths from the long samp intervals
d2 <- d2 %>%
group_by(tag) %>%
# arrange(tag,sampleNumber) %>%
mutate( lagObservedWeight = lead(observedWeight),
lagObservedLength = lead(observedLength),
grWeight = exp(lagObservedWeight - observedWeight)/as.numeric((lagDetectionDate - detectionDate)),
grLength = (lagObservedLength - observedLength)/as.numeric((lagDetectionDate - detectionDate))
) %>%
ungroup()
#d2 %>% filter(tag=='00088cc059') %>% select(tag,riverOrdered,riverN,lagRiverN,sampleNumber,observedLength,lagObservedLength)
#input data frame
ddGIn2[['bkt']] <- ddG[[1]][[2]]
dGIn[['bkt']] <- dG
# output data
mcmc[['bkt']] <- mcmcProcessed
load("D:/projects/linkedModels/data/out/dG_1530971744_crossValFALSE_bnt1997_Nimble.RData")
ddGIn2[['bnt']] <- ddG[[1]][[2]]
dGIn[['bnt']] <- dG
mcmc[['bnt']] <- mcmcProcessed
load("D:/projects/linkedModels/data/out/dG_1530993703_crossValFALSE_ats1997_Nimble.RData")
ddGIn2[['ats']] <- ddG[[1]][[2]]
dGIn[['ats']] <- dG
mcmc[['ats']] <- mcmcProcessed
rm(ddG,dG,mcmcProcessed)
load("D:/projects/linkedModels/data/cd_west_1997.RData")
ggplot(d2, aes(cBNTStd)) +
geom_histogram() +
facet_grid(season~riverOrdered)
ggplot(d2, aes(cBNTStd, grLength)) +
geom_point(alpha=0.1) +
geom_smooth(method='lm') +
facet_grid(isYOY+season~riverOrdered)
ggplot(d2, aes(cBKTStd, grLength)) +
geom_point(alpha=0.1) +
geom_smooth(method='lm') +
facet_grid(isYOY+season~riverOrdered)
len <- list()
for( spp in c("bkt","bnt","ats") ){
# len2 <- array2df(mcmc[[spp]]$sims.list$lengthExp, label.x = "est") %>%
# rename(iter=d1,obs=d2)
offset <- 1
if(spp == 'ats') offset <- 0 # not sure why this is needed. REs seem meesed up in the graphs below
inds <- data.frame(ind = dGIn[[spp]][[1]]$constants$ind[1:(length(dGIn[[spp]][[1]]$constants$ind) - offset)],
observedLength = dGIn[[spp]][[1]]$data$lengthDATA,
sampleNumber = ddGIn2[[spp]]$sampleNumber[1:(length(ddGIn2[[spp]]$sampleNumber) - offset)],
cohort = ddGIn2[[spp]]$cohort[1:(length(ddGIn2[[spp]]$cohort) - offset)],
ageInSamples = ddGIn2[[spp]]$ageInSamples[1:(length(ddGIn2[[spp]]$ageInSamples) - offset)],
riverOrdered = ddGIn2[[spp]]$riverOrdered[1:(length(ddGIn2[[spp]]$riverOrdered) - offset)],
lagObservedLength = ddGIn2[[spp]]$lagObservedLength[1:(length(ddGIn2[[spp]]$lagObservedLength) - offset)],
detectionDate = ddGIn2[[spp]]$detectionDate[1:(length(ddGIn2[[spp]]$detectionDate) - offset)],
lagDetectionDate = ddGIn2[[spp]]$lagDetectionDate[1:(length(ddGIn2[[spp]]$lagDetectionDate) - offset)])
inds$obs <- 1:(length(dGIn[[spp]][[1]]$constants$ind) - offset)
inds <- inds %>%
group_by(ind) %>%
mutate(first = min(obs),
firstObs = obs == first,
observed = !is.na(observedLength),
sampleNumber = sampleNumber,
nObs = n(),
grLength = (lagObservedLength - observedLength)/as.numeric((lagDetectionDate - detectionDate))
) %>%
ungroup()
# len1 <- left_join(len2,inds)
# rm(len2,inds)
# means
lenMean1 <- array2df(mcmc[[spp]]$mean$lengthExp, label.x = "estMean") %>%
rename(obs=d1)
len[[spp]] <- left_join(inds,lenMean1) %>%
mutate( estMean = ifelse( firstObs, observedLength, estMean ) )
grIndREMean1 <- array2df(mcmc[[spp]]$mean$grIndRE, label.x = "grIndRE") %>%
rename(ind=d1)
len[[spp]] <- left_join(len[[spp]],grIndREMean1)
# len <- left_join(len1,lenMean1)
# rm(len1,lenMean1)
rm(inds,lenMean1)
gc()
}
len[['bkt']] %>% filter( nObs == 7 ) %>%
ggplot(aes(ageInSamples,observedLength)) +
geom_point(aes(color = factor(ind)), size=3) +
geom_point(aes(ageInSamples, estMean, color=factor(ind))) +
geom_line(aes(ageInSamples, estMean, color = factor(ind))) +
geom_point(aes(ageInSamples, estMean, color=factor(ind)), size=2) +
geom_point(aes(ageInSamples, estMean, color=factor(ind)), shape = 1, size = 2, color = 'black') +
theme(legend.position="none") +
facet_wrap(~cohort)
focalFish <- 60
len[['bkt']] %>% filter( cohort == 2001 ) %>%
ggplot(aes(ageInSamples,observedLength)) +
geom_point(aes(color = riverOrdered), size = 3) +
geom_point(color = 'black', shape = 1, size = 3) +
geom_point(aes(ageInSamples, estMean, color = riverOrdered)) +
geom_line( aes(ageInSamples, estMean, group = factor(ind), color = riverOrdered)) +
geom_point(aes(ageInSamples, estMean, color = riverOrdered), size = 2) +
geom_point(aes(ageInSamples, estMean, color = riverOrdered), shape = 1, size = 2, color = 'black') +
# theme(legend.position="none") +
facet_wrap(~nObs)
len[['bkt']] %>% filter( cohort == 2001 ) %>%
ggplot(aes(ageInSamples,observedLength)) +
geom_point(aes(color = grIndRE), size = 3) +
# geom_point(color = 'black', shape = 1, size = 3) +
geom_point(aes(ageInSamples, estMean, color = grIndRE)) +
geom_line( aes(ageInSamples, estMean, group = factor(ind), color = grIndRE)) +
geom_point(aes(ageInSamples, estMean, color = grIndRE), size = 2) +
# geom_point(aes(ageInSamples, estMean, color = grIndRE), shape = 1, size = 2, color = 'black') +
# scale_colour_distiller(palette = "Spectral") +
scale_colour_gradient(low = "white", high = "black") +
# theme(legend.position="none") +
facet_wrap(~nObs)
len[['bkt']] %>% #filter( cohort == 2001 ) %>%
ggplot(aes(ageInSamples,observedLength)) +
geom_point(aes(color = grIndRE), size = 3) +
# geom_point(color = 'black', shape = 1, size = 3) +
geom_point(aes(ageInSamples, estMean, color = grIndRE)) +
geom_line( aes(ageInSamples, estMean, group = factor(ind), color = grIndRE)) +
geom_point(aes(ageInSamples, estMean, color = grIndRE), size = 2) +
# geom_point(aes(ageInSamples, estMean, color = grIndRE), shape = 1, size = 2, color = 'black') +
# scale_colour_distiller(palette = "Spectral") +
scale_colour_gradient(low = "white", high = "black") +
# theme(legend.position="none") +
facet_grid(riverOrdered~cohort)
len[['ats']] %>% filter( nObs > 0 ) %>%
group_by(ind) %>%
summarize( meanGRLength = mean(grLength,na.rm=TRUE),
meanGRIndRE = mean(grIndRE,na.rm=TRUE),
meanNObs = mean(nObs,na.rm=TRUE)) %>%
ggplot(aes(meanGRLength,meanGRIndRE)) +
geom_point() +
geom_smooth(method='lm') +
# scale_colour_gradient(low = "yellow", high = "red") +
facet_wrap(~meanNObs)
ggplot(ddG[[1]][[2]], aes(observedLength,grLength)) +
geom_point() +
facet_grid(river + isYOY ~ season+species)
allFishBySpeciesYOY$riverOrderedGG <- factor(allFishBySpeciesYOY$river,
levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T)
allFishBySpeciesYOY$seasonGG <- factor(allFishBySpeciesYOY$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
allFishBySpeciesYOY$speciesGG <- factor(allFishBySpeciesYOY$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)
#sum across YOY
allFishBySpecies <- allFishBySpeciesYOY %>%
group_by(river,season,year,species) %>%
summarize( nAllFishBySpecies = sum(nAllFishBySpeciesYOY, na.rm=T),
massAllFishBySpecies = sum(massAllFishBySpeciesYOY, na.rm=T))
allFishBySpeciesYOY <- left_join(allFishBySpeciesYOY,allFishBySpecies)
# sum across YOY and species
allFish <- allFishBySpeciesYOY %>%
group_by(river,season,year) %>%
summarize( nAllFish = sum(nAllFishBySpeciesYOY, na.rm=T),
massAllFish = sum(massAllFishBySpeciesYOY, na.rm=T))
cd <- left_join(cd,allFish)
allFish$riverOrderedGG <- factor(allFish$river,
levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T)
allFish$seasonGG <- factor(allFish$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
#########
ggplot(allFishBySpeciesYOY, aes(year,nAllFishBySpecies)) +
geom_point(aes(color = speciesGG)) +
geom_line(aes(color = speciesGG)) +
geom_line( data = allFish, aes(year,nAllFish) ) +
theme_publication() +
labs( x = "Year", y = "Count" ) +
facet_grid( seasonGG ~ riverOrderedGG)#, scales = "free_y" )
ggsave('figures/rawCounts.png')
# add in sampleNumber
tmp <- d2 %>% distinct(season,year,sampleNumber)
allFishBySpeciesYOY <- left_join(allFishBySpeciesYOY,tmp)
allFish <- left_join(allFish,tmp)
ggplot(allFishBySpeciesYOY, aes(sampleNumber,nAllFishBySpeciesYOY)) +
geom_point(aes(color = speciesGG)) +
geom_line(aes(color = speciesGG)) +
geom_line( data = allFish, aes(sampleNumber,nAllFish) ) +
theme_publication() +
labs( x = "SampleNumber", y = "Count" ) +
facet_wrap( ~riverOrderedGG )#, scales = "free_y" )
# Correlation between age classes
# justification for not using separate age class abundances
# interesting that it is less for ATS
# From the bottom of adjustCounts()
#
# nAllFishBySpeciesPStdBKT_yoy1 nAllFishBySpeciesPStdBKT_yoy2 nAllFishBySpeciesPStdBNT_yoy1 nAllFishBySpeciesPStdBNT_yoy2
# nAllFishBySpeciesPStdBKT_yoy1 1.00 0.96 -0.46 -0.43
# nAllFishBySpeciesPStdBKT_yoy2 **0.96** 1.00 -0.44 -0.38
# nAllFishBySpeciesPStdBNT_yoy1 -0.46 -0.44 1.00 0.94
# nAllFishBySpeciesPStdBNT_yoy2 -0.43 -0.38 * *0.94** 1.00
# nAllFishBySpeciesPStdATS_yoy1 -0.13 -0.17 -0.21 0.00
# nAllFishBySpeciesPStdATS_yoy2 0.06 -0.03 -0.05 0.15
# nAllFishBySpeciesPStdATS_yoy1 nAllFishBySpeciesPStdATS_yoy2
# nAllFishBySpeciesPStdBKT_yoy1 -0.13 0.06
# nAllFishBySpeciesPStdBKT_yoy2 -0.17 -0.03
# nAllFishBySpeciesPStdBNT_yoy1 -0.21 -0.05
# nAllFishBySpeciesPStdBNT_yoy2 0.00 0.15
# nAllFishBySpeciesPStdATS_yoy1 1.00 0.47
# nAllFishBySpeciesPStdATS_yoy2 *0.47* 1.00
d2 %>% group_by(riverOrderedGG,seasonGG,year) %>%
filter( riverOrderedGG == "West Brook" ) %>%
summarize(meanF = mean(meanFlow, na.rm=T),
sdF = sd(meanFlow, na.rm=T)) %>%
ggplot( aes(year,meanF) ) +
geom_point( size = 3 ) +
geom_line() +
# geom_smooth(method = 'lm', se=F) +
geom_errorbar( aes(ymin = meanF - sdF, ymax = meanF + sdF, width = 0.33 ) ) +
labs( x = "Year", y = "Stream flow (x/x)" ) +
theme_publication() +
facet_wrap(~seasonGG)
ggsave('figures/flow.png')
lmFunction <- function(df){ lm( meanT ~ year, data=df ) }
tempMods <- d2 %>%
group_by(riverOrderedGG,seasonGG,year) %>%
summarize(meanT = mean(meanTemperature, na.rm=T),
sdT = sd(meanTemperature, na.rm=T)) %>%
group_by( riverOrderedGG, seasonGG ) %>%
nest() %>%
mutate( model = map(data,lmFunction) )
tempModsTidy <- tempMods %>%
mutate(tidy = map(model,broom::tidy)) %>%
dplyr::select(-data,-model) %>%
unnest() %>%
filter( term == 'year') %>%
mutate( sig = ifelse(p.value < 0.05,"*",""))
d2 %>% group_by(riverOrderedGG,seasonGG,year) %>%
summarize(meanT = mean(meanTemperature, na.rm=T),
sdT = sd(meanTemperature, na.rm=T)) %>%
ggplot( aes(year,meanT) ) +
geom_smooth(method = 'lm', se=F, color = 'darkgrey') +
geom_point( ) +
geom_line() +
geom_errorbar( aes(ymin = meanT - sdT, ymax = meanT + sdT, width = 0.33 ) ) +
labs( x = "Year", y = "Stream temperature (C)" ) +
geom_text(data=tempModsTidy,aes(x=1999,y=18, label=paste(round(estimate,2),sig))) +
theme_publication() +
facet_grid(seasonGG ~ riverOrderedGG)
ggsave('figures/temp.png')
cd$riverOrderedGG <- factor(cd$river,
levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T)
cd$seasonGG <- factor(cd$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
cd$speciesGG <- factor(cd$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)
pd <- position_dodge(0.65)
cd %>% group_by(speciesGG,riverOrderedGG,seasonGG) %>%
filter(isYOY == 1) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T)) %>%
ggplot( aes(seasonGG,meanGR, color=speciesGG, shape=speciesGG) ) +
geom_point( position = pd, size = 2 ) +
geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) +
ggtitle( "isYOY = 1" ) +
labs( x = "Season", y = "Growth rate (mm/d)" ) +
#scale_x_continuous( breaks = 1997:2015 ) +
theme_publication() +
theme( panel.grid.minor = element_blank() ) +
facet_wrap( ~ riverOrderedGG)
ggsave('figures/meanGRYOY1.png', width = 7, height = 7, units='in')
cd %>% group_by(speciesGG,riverOrderedGG,seasonGG) %>%
filter(isYOY == 2) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T)) %>%
ggplot( aes(seasonGG,meanGR, color=speciesGG, shape=speciesGG) ) +
geom_point( position = pd,size = 2 ) +
geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) +
ggtitle( "isYOY = 2" ) +
labs( x = "Season", y = "Growth rate (mm/d)" ) +
#scale_x_continuous( breaks = 1997:2015 ) +
theme_publication() +
theme( panel.grid.minor = element_blank() ) +
facet_wrap( ~ riverOrderedGG)
ggsave('figures/meanGRYOY2.png', width = 7, height = 7, units='in')
# mean growth by species
cd %>% group_by(species) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T))
.162/.194 #bkt/bnt
.151/.194 #ats/bnt
cd %>% group_by(isYOY) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T))
.119/.215 #juv/yoy
.215/.119
cd %>% group_by(river) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T))
# river meanGR sdGR
# <chr> <dbl> <dbl>
# 1 wb jimmy 0.132 0.131
# 2 wb mitchell 0.179 0.160
# 3 wb obear 0.137 0.129
# 4 west brook 0.170 0.175
((.132+.137)/2)/((.17+.179)/2)
((.17+.179)/2)/((.132+.137)/2)
cd %>% group_by(season) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T))
((.0464+.0894+.0914)/3)/(.417)
.417/((.0464+.0894+.0914)/3)
cd %>% filter() %>%
ggplot(aes(observedLength,grLength)) +
geom_point() +
facet_grid(season~species)
# ggPInt <- array2df(ddD$sims.list$pBetaInt, label.x = "est")
#
# ggPInt$chain <- rep(1:ddD$mcmc.info$n.chains, each = ddD$mcmc.info$n.samples/ddD$mcmc.info$n.chains)
# ggPInt$iter <- 1:as.numeric(ddD$mcmc.info$n.samples/ddD$mcmc.info$n.chains)
#
# ggplot(filter(ggPInt), aes(iter,est)) +
# geom_point( aes(color=factor(chain)), size = 0.1 ) +
# ylim(-1.5,1.5) +
# facet_grid(d2+d4~d5+d3)
p <- array2df(ddD$q50$pBetaInt, label.x = "est") %>%
rename(est=est,isYOY=d1,species=d2,season=d3,river=d4,year=d5)
ggplot(p, aes(year,invlogit(est),color=factor(species))) +geom_point() + facet_grid(river ~ season + isYOY)
# seem to need to run the ggsave outside of markdown by pasting into console
ggsave('figures/p.png')
phi <- array2df(ddD$q50$phiBetaInt, label.x = "est") %>%
rename(est=est,isYOY=d1,species=d2,season=d3,river=d4,year=d5)
ggplot(phi, aes(year,invlogit(est),color=factor(species))) +geom_point() + facet_grid(river ~ season + isYOY)
ggsave('figures/phi.png')
# # of skunked sections
tmp <- dddD[[2]] %>%
group_by(riverOrdered,year,season,section) %>%
summarize( s = sum(enc) ) %>%
filter( s == 0 )
table(tmp$riverOrdered,tmp$year,tmp$season)
rHat <- list()
rHat[['bkt']] <- array2df(mcmc[['bkt']]$Rhat$grBeta, levels = list(beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="rHat")
rHat[['bnt']] <- array2df(mcmc[['bnt']]$Rhat$grBeta, levels = list(beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="rHat")
rHat[['ats']] <- array2df(mcmc[['ats']]$Rhat$grBeta, levels = list(beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="rHat")
# rHat[['bkt']] <- array2df(mcmc[['bkt']]$Rhat$grBeta, levels = list(beta=1:8,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="rHat")
# rHat[['bnt']] <- array2df(mcmc[['bnt']]$Rhat$grBeta, levels = list(beta=1:8,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="rHat")
# rHat[['ats']] <- array2df(mcmc[['ats']]$Rhat$grBeta, levels = list(beta=1:8,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="rHat")
ggplot(rHat[["ats"]], aes(beta,rHat)) +
geom_point() +
facet_grid(river+isYOY~season)
grInt <- list()
grInt[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int")
grInt[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int")
grInt[['ats']] <- array2df(mcmc[['ats']]$sims.list$grInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int")
ggplot( grInt[['bkt']] %>% filter(river %in% c('west brook','wb jimmy')), aes(int)) +
geom_freqpoly() +
geom_freqpoly(aes((int)),color = 'brown', data = grInt[['bnt']]) +
geom_freqpoly(aes((int)),color = 'lightblue', data = grInt[['ats']]) +
# xlim(-0.5,1.5) +
theme_publication() +
facet_grid(season ~ river + isYOY)
##################################################
# as range plot
plotHPD <- function(sppToPlot){
nRivers <- list()
nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1
probs <- c(0.05,0.1,0.5,0.9,0.95)
q <- grInt[[sppToPlot]] %>%
group_by(isYOY,season,river) %>%
summarize(
p = list(probs),
q = list(quantile(int,probs))
) %>%
unnest() %>%
spread(key = p, value = q)
p <- array2df(mcmc[[sppToPlot]]$mean$grInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="int") %>%
mutate( season = as.numeric(season) ) %>%
left_join(.,meanSampleIntervalMean, by = 'season')
p$intByDay <- p$int/p$meanSampleInterval
p$intByDay <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDay)
p$parameter <- paste(p$river,p$season,sep="_")
pUp <- array2df(mcmc[[sppToPlot]]$q97.5$grInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intUp") %>%
mutate( season = as.numeric(season) )
pDown <- array2df(mcmc[[sppToPlot]]$q2.5$grInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intDown") %>%
mutate( season = as.numeric(season) )
p <- left_join(p,pUp); p <- left_join(p,pDown)
p$intByDayUp <- p$intUp/p$meanSampleInterval; p$intByDayDown <- p$intDown/p$meanSampleInterval
p$intByDay <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDay)
p$intByDayUp <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDayUp)
p$intByDayDown <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDayDown)
#print(c(min(p$intByDayDown,na.rm=T),max(p$intByDayDown,na.rm=T)))
gg <- ggplot(p, aes(x = intByDay, y = reorder(parameter,intByDay,na.rm=TRUE), color = factor(season), shape = factor(isYOY))) +
theme_publication() +
geom_point(size = 3) +
geom_segment(aes(x = intByDay, xend = intByDayUp, yend = reorder(parameter, intByDay)),
size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
geom_segment(aes(x = intByDay, xend = intByDayDown, yend = reorder(parameter, intByDay)),
size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
xlim(c(-0.0015,0.01)) +
xlab("Mean daily growth (mm/d)") + ylab("River_Season")
# ggtitle(sppToPlot)
return(gg)
}
plotHPD('bkt'); ggsave('figures/HDP_bkt.png', height = 3, width = 3, scale = 2)
plotHPD('bnt'); ggsave('figures/HDP_bnt.png', height = 3, width = 3, scale = 2)
plotHPD('ats'); ggsave('figures/HDP_ats.png', height = 3, width = 3, scale = 2)
sigmaInt <- list()
sigmaInt[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$sigmaInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int")
sigmaInt[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$sigmaInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int")
sigmaInt[['ats']] <- array2df(mcmc[['ats']]$sims.list$sigmaInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int")
bw <- 0.05
ggplot( sigmaInt[['bkt']] %>% filter(river %in% c('west brook','wb jimmy')), aes(int)) +
geom_freqpoly(binwidth = bw) +
geom_freqpoly(aes(int),binwidth = bw,color = 'brown', data = sigmaInt[['bnt']]) +
geom_freqpoly(aes(int),binwidth = bw,color = 'lightblue', data = sigmaInt[['ats']]) +
#xlim(-4.5,-1) +
xlim(-5,2.5) +
theme_publication() +
facet_grid(season ~ river + isYOY)
##################################################
# as range plot
plotHPDSigma <- function(sppToPlot){
nRivers <- list()
nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1
probs <- c(0.05,0.1,0.5,0.9,0.95)
q <- sigmaInt[[sppToPlot]] %>%
group_by(isYOY,season,river) %>%
summarize(
p = list(probs),
q = list(quantile(int,probs))
) %>%
unnest() %>%
spread(key = p, value = q)
p <- array2df(mcmc[[sppToPlot]]$mean$sigmaInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="int") %>%
mutate( season = as.numeric(season) ) %>%
left_join(.,meanSampleIntervalMean, by = 'season')
# p$intByDay <- p$int/p$meanSampleInterval
# p$intByDay <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDay)
p$parameter <- paste(p$river,p$season,sep="_")
pUp <- array2df(mcmc[[sppToPlot]]$q97.5$sigmaInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intUp") %>%
mutate( season = as.numeric(season) )
pDown <- array2df(mcmc[[sppToPlot]]$q2.5$sigmaInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intDown") %>%
mutate( season = as.numeric(season) )
p <- left_join(p,pUp); p <- left_join(p,pDown)
# p$intByDayUp <- p$intUp/p$meanSampleInterval; p$intByDayDown <- p$intDown/p$meanSampleInterval
p$int <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$int)
p$intUp <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intUp)
p$intDown <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intDown)
#print(c(min(p$intDown,na.rm=T),max(p$intDown,na.rm=T)))
gg <- ggplot(p, aes(x = exp(int), y = reorder(parameter,int,na.rm=TRUE), color = factor(season), shape = factor(isYOY))) +
theme_publication() +
geom_point(size = 3) +
geom_segment(aes(x = exp(int), xend = exp(intUp), yend = reorder(parameter, int)),
size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
geom_segment(aes(x = exp(int), xend = exp(intDown), yend = reorder(parameter, int)),
size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
#xlim(c(-0.0015,0.01)) +
xlab("Standard deviation of daily growth (mm/d)") + ylab("River_Season")
# ggtitle(sppToPlot)
return(gg)
}
plotHPDSigma('bkt'); ggsave('figures/HDPSigma_bkt.png', height = 3, width = 3, scale = 2)
plotHPDSigma('bnt'); ggsave('figures/HDPSigma_bnt.png', height = 3, width = 3, scale = 2)
plotHPDSigma('ats'); ggsave('figures/HDPsigma_ats.png', height = 3, width = 3, scale = 2)
#riverOrderedIn <- factor(c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),labels = c("west brook","wb jimmy","wb mitchell","wb obear"), ordered = T)
grIntMu <- list()
grIntMu[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grIntMu, levels = list(iter=NA,season=1:4), label.x="int")
grIntMu[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grIntMu, levels = list(iter=NA,season=1:4), label.x="int")
grIntMu[['ats']] <- array2df(mcmc[['ats']]$sims.list$grIntMu, levels = list(iter=NA,season=1:4), label.x="int")
bw <- 0.05
ggplot( grIntMu[['bkt']] , aes(int)) +
geom_freqpoly(binwidth = bw, size=1.2) +
geom_freqpoly(aes(int),binwidth = bw,color = 'brown', data = grIntMu[['bnt']]) +
geom_freqpoly(aes(int),binwidth = bw,color = 'lightblue', data = grIntMu[['ats']]) +
#xlim(-4.5,-1) +
xlim(-0.25,0.75) +
# ylim(0,2500) +
theme_publication() +
facet_wrap(~season )
# get all betas in a DF
grBeta <- list()
grBeta[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grBeta, levels = list(iter=NA,beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>%
mutate( species = 'bkt')
grBeta[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grBeta, levels = list(iter=NA,beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>%
mutate( species = 'bnt')
grBeta[['ats']] <- array2df(mcmc[['ats']]$sims.list$grBeta, levels = list(iter=NA,beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>%
mutate( species = 'ats')
grBetaDF1 <- rbind(grBeta[['bkt']],grBeta[['bnt']],grBeta[['ats']])
grBetaBNT <- list()
grBetaBNT[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grBetaBNT, levels = list(iter=NA,beta=8:9,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>%
mutate( species = 'bkt')
grBetaBNT[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grBetaBNT, levels = list(iter=NA,beta=8:9,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>%
mutate( species = 'bnt')
grBetaBNT[['ats']] <- array2df(mcmc[['ats']]$sims.list$grBetaBNT, levels = list(iter=NA,beta=8:9,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>%
mutate( species = 'ats')
grBetaBNTDF <- rbind(grBetaBNT[['bkt']],grBetaBNT[['bnt']],grBetaBNT[['ats']])
grBetaATS <- list()
grBetaATS[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grBetaATS, levels = list(iter=NA,beta=10,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>%
mutate( species = 'bkt')
grBetaATS[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grBetaATS, levels = list(iter=NA,beta=10,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>%
mutate( species = 'bnt')
grBetaATS[['ats']] <- array2df(mcmc[['ats']]$sims.list$grBetaATS, levels = list(iter=NA,beta=10,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>%
mutate( species = 'ats')
grBetaATSDF <- rbind(grBetaATS[['bkt']],grBetaATS[['bnt']],grBetaATS[['ats']])
grBetaDF <- rbind(grBetaDF1,grBetaBNTDF,grBetaATSDF)
##################################################
# as range plot
nRivers <- list()
nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1
speciesList <- list('bkt','bnt','ats')
probs <- c(0.05,0.1,0.5,0.9,0.95)
nSpecies <- 3
q <- grBetaDF %>%
group_by(beta,isYOY,season,river,species) %>%
summarize(
p = list(probs),
q = list(quantile(int,probs))
) %>%
unnest() %>%
spread(key = p, value = q) %>%
data.frame() %>%
rename( median = X0.5, lo = X0.05, hi = X0.95 )
q <- q %>% filter( !(season == 2 & isYOY == 0))
q$riverGG <- factor(q$river,
levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
labels = c("WB","OL","OS","IL"), ordered = T)
q$seasonGG <- factor(q$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
q$speciesGG <- factor(q$species, levels = c('bkt','bnt','ats'), labels = c("BKT", "BNT", "ATS"), ordered = T)
q$betaGG <- factor(q$beta, levels = 1:10,
levels = c("Temp","Flow","Temp*Flow","Temp^2","Flow^2","BKT","Length","BNT","BNT*BKT","ATS"), # make sure this order matches the nimble code
labels = c("Temp","Flow","Temp*Flow","Temp^2","Flow^2","BKT","Length","BNT","BNT*BKT","ATS"),ordered = T)
q$parameter = paste0(q$riverGG,"_",q$speciesGG)
q$parameterOrd <- factor(q$parameter,
levels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"),
labels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), ordered = T)
q <- q %>% filter( !(betaGG == 'ATS' & (riverGG == "OS" | riverGG == "IL" | riverGG == "OL")) )
q <- q %>% filter( !(betaGG == 'BNT' & (riverGG == "OS" | riverGG == "IL")) )
q <- q %>% filter( !(betaGG == 'BNT*BKT' & (riverGG == "OS" | riverGG == "IL")) )
# q <- q %>% filter( !((riverGG == "OS" | riverGG == "IL")) )
#reorder(parameter,median,na.rm=TRUE)
(min(q$median))
(max(q$median))
# q[q$median < -0.13,]
ggplot(q, aes(x = median, y = parameterOrd, color = factor(speciesGG), shape = factor(isYOY))) +
geom_vline(xintercept = 0,color='black') +
theme_publication() +
geom_point(size = 3) +
geom_segment(aes(x = median, xend = hi, yend = parameterOrd),
size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
geom_segment(aes(x = median, xend = lo, yend = parameterOrd),
size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
xlim(c(-0.15,0.11)) +
xlab("Parameter estimate") + ylab("River_Species") +
facet_grid(seasonGG~betaGG)
ggsave('figures/grBetas.png', height = 7, width = 10, scale = 2)
# get all betas in a DF
grBetaSigma <- list()
grBetaSigma[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$sigmaBeta, levels = list(iter=NA,beta=1:3,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>%
mutate( species = 'bkt')
grBetaSigma[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$sigmaBeta, levels = list(iter=NA,beta=1:3,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>%
mutate( species = 'bnt')
grBetaSigma[['ats']] <- array2df(mcmc[['ats']]$sims.list$sigmaBeta, levels = list(iter=NA,beta=1:3,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>%
mutate( species = 'ats')
grBetaSigmaDF <- rbind(grBetaSigma[['bkt']],grBetaSigma[['bnt']],grBetaSigma[['ats']])
##################################################
# as range plot
nRivers <- list()
nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1
speciesList <- list('bkt','bnt','ats')
probs <- c(0.05,0.1,0.5,0.9,0.95)
nSpecies <- 3
q <- grBetaSigmaDF %>%
group_by(beta,isYOY,season,river,species) %>%
summarize(
p = list(probs),
q = list(quantile(int,probs))
) %>%
unnest() %>%
spread(key = p, value = q) %>%
data.frame() %>%
rename( median = X0.5, lo = X0.05, hi = X0.95 )
q <- q %>% filter( !(season == 2 & isYOY == 0))
q$riverGG <- factor(q$river,
levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
labels = c("WB","OL","OS","IL"), ordered = T)
q$seasonGG <- factor(q$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
q$speciesGG <- factor(q$species, levels = c('bkt','bnt','ats'), labels = c("BKT", "BNT", "ATS"), ordered = T)
q$betaGG <- factor(q$beta, levels = 1:10, labels = c("Temp", "Flow", "Temp*Flow","Temp^2","Flow^2","BKT","Length","BNT","BNT*BKT","ATS"), ordered = T)
q$parameter = paste0(q$riverGG,"_",q$speciesGG)
q$parameterOrd <- factor(q$parameter, levels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), labels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), ordered = T)
q <- q %>% filter( !(betaGG == 'ATS' & (riverGG == "OS" | riverGG == "IL" | riverGG == "OL")) )
q <- q %>% filter( !(betaGG == 'BNT' & (riverGG == "OS" | riverGG == "IL")) )
q <- q %>% filter( !(betaGG == 'BNT*BKT' & (riverGG == "OS" | riverGG == "IL")) )
# q <- q %>% filter( !((riverGG == "OS" | riverGG == "IL")) )
#reorder(parameter,median,na.rm=TRUE)
(min(q$median))
(max(q$median))
# q[q$median < -0.13,]
ggplot(q, aes(x = median, y = parameterOrd, color = factor(speciesGG), shape = factor(isYOY))) +
geom_vline(xintercept = 0,color='black') +
theme_publication() +
geom_point(size = 3) +
geom_segment(aes(x = median, xend = hi, yend = parameterOrd),
size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
geom_segment(aes(x = median, xend = lo, yend = parameterOrd),
size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) +
xlim(c(-0.35,0.5)) +
xlab("Parameter estimate") + ylab("River_Species") +
facet_grid(seasonGG~betaGG)
ggsave('figures/grBetaSigmas.png', height = 7, width = 7, scale = 2)
grIndRE <- list()
grIndRE[['bkt']] <- as.data.frame(mcmc[['bkt']]$q50$grIndRE); names(grIndRE[['bkt']]) <- 'value'
grIndRE[['bnt']] <- as.data.frame(mcmc[['bnt']]$q50$grIndRE); names(grIndRE[['bnt']]) <- 'value'
grIndRE[['ats']] <- as.data.frame(mcmc[['ats']]$q50$grIndRE); names(grIndRE[['ats']]) <- 'value'
table(mcmc[['bkt']]$q50$grIndRE==0)/length(mcmc[['bkt']]$q50$grIndRE)
table(mcmc[['bnt']]$q50$grIndRE==0)/length(mcmc[['bnt']]$q50$grIndRE)
table(mcmc[['ats']]$q50$grIndRE==0)/length(mcmc[['ats']]$q50$grIndRE)
ggplot( grIndRE[[1]], aes(value,..density..) ) +
geom_freqpoly(binwidth = 0.00025) +
geom_freqpoly(aes(value), color = "brown", binwidth = 0.00025, data = grIndRE[[2]]) +
geom_freqpoly(aes(value), color = "lightblue", binwidth = 0.00025, data = grIndRE[[3]]) +
ylim(c(0,100))
#
#####run if get new predictions - takes a while to run ####
# load('D:/projects/linkedModels/data/out/P_ForMike_bkt.RData')
# pBoth$species <- 'bkt'
# pAll <- pBoth; rm(pBoth)
#
# load('D:/projects/linkedModels/data/out/P_ForMike_bnt.RData')
# pBoth$species <- 'bnt'
# pAll <- rbind(pAll,pBoth); rm(pBoth)
#
# load('D:/projects/linkedModels/data/out/P_ForMike_ats.RData')
# pBoth$species <- 'ats'
# pAll <- rbind(pAll,pBoth); rm(pBoth)
#
# pred <- pAll %>% group_by(species,len,flow,temp,cBKT,cBNT,cATS,isYOY,season,river) %>%
# summarize(meanPredGR = mean(predGr, na.rm=T),
# sdPredGR = sd(predGr, na.rm=T),
# meanPredGRSigma = mean(predGrSigma, na.rm=T),
# sdPredGRSigma = sd(predGrSigma, na.rm=T),
# meanPredGRCV = mean(predCV, na.rm=T),
# sdPredGRCV = sd(predCV, na.rm=T)) %>%
# ungroup()
#
# rm(pAll)
# save(pred,file = 'D:/projects/linkedModels/data/out/P_ForMike_allMeans.RData')
load('D:/projects/linkedModels/data/out/P_ForMike_allMeans.RData')
pred <- pred %>% filter( !(isYOY == 0 & season == 2) )
#pred <- pred %>% filter( !(isYOY == 0 & river == 'wb jimmy' & season == 2) )
#graph function, in analyzeOutputFunctions.R
pred$riverGG <- factor(pred$river,
levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T)
pred$seasonGG <- factor(pred$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
pred$speciesGG <- factor(pred$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)
plotPredMeans(pred, "meanPredGR", "len", 0)
plotPredMeans(pred, "meanPredGR", "len", 1)
plotPredMeans(pred, "meanPredGR", "flow", 0)
plotPredMeans(pred, "meanPredGR", "flow", 1)
plotPredMeans(pred, "meanPredGR", "temp", 0)
plotPredMeans(pred, "meanPredGR", "temp", 1)
plotPredMeans(pred, "meanPredGR", "cBKT", 0)
plotPredMeans(pred, "meanPredGR", "cBKT", 1)
plotPredMeans(pred, "meanPredGR", "cBNT", 0)
plotPredMeans(pred, "meanPredGR", "cBNT", 1)
plotPredMeans(pred, "meanPredGR", "cATS", 0)
plotPredMeans(pred, "meanPredGR", "cATS", 1)
plotPredMeans(pred, "meanPredGR", c("temp","flow"), 0)
plotPredMeans(pred, "meanPredGR", c("temp","flow"), 1)
plotPredMeans(pred, "meanPredGR", c("flow","temp"), 0)
plotPredMeans(pred, "meanPredGR", c("flow","temp"), 1)
ggsave('figures/tmp.png', width = 9, height = 7, units='in')
all <- c('len','temp','flow','cBKT','cBNT','cATS')
allLabels <- c('Length (mm)','Temperature (C)','Flow (m3/s)','Brook Trout abundance','Brown Trout abundance','Atlantic Salmon abundance')
all <- c('len','cBKT','temp','cBNT','flow','cATS')
allLabels <- c('Length (mm)','Brook Trout abundance','Stream temperature (C)','Brown Trout abundance','Stream flow (m3/s)','Atlantic Salmon abundance')
remFacetLabel <- c(3,4,5,6)
selectSingleVar <- function(p, all, depVar, varsToPlot, isYOYGG){
notPlot <- NA
notPlot[1] <- all[!(all %in% varsToPlot)][1]
notPlot[2] <- all[!(all %in% varsToPlot)][2]
notPlot[3] <- all[!(all %in% varsToPlot)][3]
notPlot[4] <- all[!(all %in% varsToPlot)][4]
notPlot[5] <- all[!(all %in% varsToPlot)][5]
pGG <- p %>% filter(isYOY == isYOYGG,
eval(as.name(notPlot[1])) == 0,
eval(as.name(notPlot[2])) == 0,
eval(as.name(notPlot[3])) == 0,
eval(as.name(notPlot[4])) == 0,
eval(as.name(notPlot[5])) == 0
)#%>%
# distinct(species, isYOY, riverGG, seasonGG, speciesGG, .keep_all = TRUE)
return(pGG)
}
##### by species
isYOYGG <- 0
riverGGToPlot <- "West Brook" #"Open Large"
outputToPlot <- "meanPredGR"
ggSingle <- list()
singleDat <- list()
for(i in 1:length(all))
#https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r
local({ i <- i
singleDat[[i]] <- selectSingleVar(pred, all, outputToPlot, all[i], isYOYGG)
print(c(i,all[i]))
# print(str(singleDat[[i]]))
tmp <- singleDat[[i]] %>% filter(riverGG == riverGGToPlot)
ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=speciesGG ) ) +
geom_line( size = 1.225 ) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title.y = element_blank(),
legend.position="none") +
geom_hline(yintercept = 0) +
ylim(-0.15,0.8) +
scale_x_continuous(allLabels[i]) +
facet_grid(~seasonGG)
if( i %in% remFacetLabel ) {
ggSingle[[i]] <<- ggSingle[[i]] +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
}
})
grid.arrange(grobs=ggSingle, left = "Growth rate (mm/d)")
###### by river for bkt
isYOYGG <- 1
speciesGGToPlot <- "Brook Trout"
outputToPlot <- "meanPredGR"
ggSingle <- list()
singleDat <- list()
for(i in 1:length(all))
#https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r
local({ i <- i
singleDat[[i]] <- selectSingleVar(pred, all, outputToPlot, all[i], isYOYGG)
print(c(i,all[i]))
# print(str(singleDat[[i]]))
tmp <- singleDat[[i]] %>% filter(speciesGG == speciesGGToPlot)
ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=riverGG ) ) +
geom_line( size = 1.225 ) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title.y = element_blank(),
legend.position="none"
) +
geom_hline(yintercept = 0) +
ylim(-0.15,0.8) +
scale_x_continuous(allLabels[i]) +
facet_grid(~seasonGG)
if( i %in% remFacetLabel ) {
ggSingle[[i]] <<- ggSingle[[i]] +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
}
})
grid.arrange(grobs=ggSingle, left = "Growth rate (mm/d)")
selectPairVars <- function(p, all, depVar, varsToPlot, isYOYGG){
notPlot <- NA
notPlot[1] <- all[!(all %in% varsToPlot)][1]
notPlot[2] <- all[!(all %in% varsToPlot)][2]
notPlot[3] <- all[!(all %in% varsToPlot)][3]
notPlot[4] <- all[!(all %in% varsToPlot)][4]
# notPlot[5] <- all[!(all %in% varsToPlot)][5]
pGG <- p %>% filter( isYOY == isYOYGG,
eval(as.name(notPlot[1])) == 0,
eval(as.name(notPlot[2])) == 0,
eval(as.name(notPlot[3])) == 0,
eval(as.name(notPlot[4])) == 0 )
return(pGG)
}
h <- list(); h[['bkt']] <- 8.5; h[['bnt']] <- 4.5; h[['ats']] <- 2.75
for(i in c('bkt','bnt', 'ats')){
for(j in 0:1){
print(c(i,j))
p <- selectPairVars(pred, all, "meanPredGR", c("flow","temp"), j) %>%
filter(species == i)
ggplot(p, aes(x = flow, y = temp, z = meanPredGR)) +
geom_tile(aes(fill = meanPredGR)) +
stat_contour(bins = 15, color = 'darkgrey') +
ggtitle(paste(j,i)) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
#axis.title.y = element_blank(),
legend.position="none"
) +
facet_grid( riverGG ~ seasonGG )
ggsave(paste0('figures/predGRGrid',"_",i,"_",j,'.png'), width = 9, height = h[[i]], units='in')
}
}
#raw data plot - not very useful
ggplot(d2, aes(x=flowStd,y=tempStd,z=grLength)) +
# stat_density_2d( color = 'darkgrey') +
stat_density_2d(geom = "raster", aes(fill = ..density..), contour = FALSE,interpolate=TRUE,h=2.5) +
facet_grid(season+isYOY~river+species)
ggplot(d2, aes(x=tempStd,y=grLength)) +
geom_smooth() +
geom_smooth(method='lm',color='black',se=FALSE) +
ylim(c(0,0.8)) +
facet_grid(river+isYOY~season)
all <- c('temp','flow')
allLabels <- c('Stream temperature (C)','Stream flow (m3/s)')
remFacetLabel <- c(2)
selectSingleVarSigma <- function(p, all, depVar, varsToPlot, isYOYGG){
notPlot <- NA
notPlot[1] <- all[!(all %in% varsToPlot)][1]
pGG <- p %>% filter(isYOY == isYOYGG,
eval(as.name(notPlot[1])) == 0,
len == 0, cBKT == 0, cBNT == 0, cATS == 0
)#%>%
# distinct(species, isYOY, riverGG, seasonGG, speciesGG, .keep_all = TRUE)
return(pGG)
}
##### by species
isYOYGG <- 0
riverGGToPlot <- "Open Large" #"West Brook" #
outputToPlot <- "meanPredGRSigma"
ggSingle <- list()
singleDat <- list()
for(i in 1:length(all))
#https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r
local({ i <- i
singleDat[[i]] <- selectSingleVarSigma(pred, all, outputToPlot, all[i], isYOYGG)
print(c(i,all[i]))
# print(str(singleDat[[i]]))
tmp <- singleDat[[i]] %>% filter(riverGG == riverGGToPlot)
ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=speciesGG ) ) +
geom_line( size = 1.225 ) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title.y = element_blank(),
legend.position="none") +
geom_hline(yintercept = 0) +
# ylim(-0.15,0.8) +
scale_x_continuous(allLabels[i]) +
facet_grid(~seasonGG)
if( i %in% remFacetLabel ) {
ggSingle[[i]] <<- ggSingle[[i]] +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
}
})
# add thrid row to graph for interaction of tem and flow?
# run this and export
grid.arrange(grobs=ggSingle, left = "Standard deviation of growth rate (mm/d)")
#ggsave(ggCV, paste0('figures/predGRCV_YOY',isYOYGG,'_',riverGGToPlot,'.png'), width = 9, height = 7, units='in')
###### by river for bkt
isYOYGG <- 0
speciesGGToPlot <- "Brook Trout"
outputToPlot <- "meanPredGRSigma"
ggSingle <- list()
singleDat <- list()
for(i in 1:length(all))
#https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r
local({ i <- i
singleDat[[i]] <- selectSingleVarSigma(pred, all, outputToPlot, all[i], isYOYGG)
print(c(i,all[i]))
# print(str(singleDat[[i]]))
tmp <- singleDat[[i]] %>% filter(speciesGG == speciesGGToPlot)
ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=riverGG ) ) +
geom_line( size = 1.225 ) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title.y = element_blank(),
legend.position="none"
) +
geom_hline(yintercept = 0) +
#ylim(-0.15,0.8) +
scale_x_continuous(allLabels[i]) +
facet_grid(~seasonGG)
if( i %in% remFacetLabel ) {
ggSingle[[i]] <<- ggSingle[[i]] +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
}
})
grid.arrange(grobs=ggSingle, left = "Standard deviation of growth rate (mm/d)")
###############################
# flow*temp interaction plots
isYOYGG <- 1
riverGGToPlot <- "West Brook" #"Open Large" #
outputToPlot <- "meanPredGR" #"meanPredGRSigma"
interact <- pred %>% filter(isYOY == isYOYGG,
riverGG == riverGGToPlot,
len == 0, cBKT == 0, cBNT == 0, cATS == 0 ) %>%
mutate( temp_spp = paste0(temp,speciesGG))
ggplot( interact, aes( flow, eval(as.name(outputToPlot)), group=temp_spp, color=speciesGG ) ) +
geom_line( aes(alpha = temp + 3), size = 1.225 ) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title.y = element_blank()
# legend.position="none"
) +
geom_hline(yintercept = 0) +
# ylim(-0.15,0.8) +
scale_x_continuous('Flow') +
facet_grid(~seasonGG)
###############################
# flow*temp interaction plots
isYOYGG <- 1
riverGGToPlot <- "Open Large" #"West Brook" #
outputToPlot <- "meanPredGR" #"meanPredGRSigma"
interact <- pred %>% filter(isYOY == isYOYGG,
riverGG == riverGGToPlot,
len == 0, temp == 0, flow == 0, cATS == 0 ) %>%
mutate( bnt_spp = paste0(cBNT,speciesGG))
ggplot( interact, aes( cBKT, eval(as.name(outputToPlot)), group=bnt_spp, color=speciesGG ) ) +
geom_line( aes(alpha = cBNT + 3), size = 1.225 ) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title.y = element_blank()
# legend.position="none"
) +
geom_hline(yintercept = 0) +
# ylim(-0.15,0.8) +
scale_x_continuous('Brook Trout abundance') +
facet_grid(~seasonGG)
round(mcmc$summary$all.chains[startsWith(attributes(mcmc$summary$all.chains)$dimnames[[1]], "grBeta["),],2)
pd <- position_dodge(0.65)
# isYOY1
d2 %>% group_by(speciesGG,riverOrderedGG,seasonGG,year) %>%
filter(isYOYDATA == 1) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T)) %>%
ggplot( aes(year,meanGR, color=speciesGG, shape=speciesGG) ) +
geom_point( position=pd ) +
geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) +
ggtitle( "isYOY = 1" ) +
labs( x = "Year", y = "Growth rate (mm/d)" ) +
scale_x_continuous( breaks = 1997:2015 ) +
theme_publication() +
theme( panel.grid.minor = element_blank(),
axis.text.x = element_text(angle=90, vjust=0.5, size=8) ) +
facet_grid(seasonGG ~ riverOrderedGG)
ggsave('figures/grByYearYOY1.png')
# isYOY2
d %>% group_by(speciesGG,riverOrderedGG,seasonGG,year) %>%
filter(isYOYDATA == 2) %>%
summarize(meanGR = mean(grLength, na.rm=T),
sdGR = sd(grLength, na.rm=T)) %>%
ggplot( aes(year,meanGR, color=speciesGG, shape=speciesGG) ) +
geom_point( position=pd ) +
geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) +
ggtitle( "isYOY = 2" ) +
labs( x = "Year", y = "Growth rate (mm/d)" ) +
scale_x_continuous( breaks = 1997:2015 ) +
theme_publication() +
theme( panel.grid.minor = element_blank(),
axis.text.x = element_text(angle=90, vjust=0.5, size=8) ) +
facet_grid(seasonGG ~ riverOrderedGG)
ggsave('figures/grByYearYOY2.png')
d2$isYOY <- ifelse( d2$ageInSamples <= 3, 1, 2 )
d2 <- d2 %>%
group_by(tag) %>%
mutate( lagRiverN = lead(riverN),
riverNMove = paste(riverN,lagRiverN, sep="_")
#gr = (observedLength - lagObservedLength)/as.integer(detectionDate - lagDetectionDate)
)
#d2 %>% filter(tag=='00088cbf44') %>% select(tag,riverN,lagRiverN,sampleNumber,observedLength,lagObservedLength)
table(d2$riverNMove)
cdMeansByisYOYriverNMove <- d2 %>%
#mutate( gr = (observedLength - lagObservedLength)/as.integer(detectionDate - lagDetectionDate) ) %>%
group_by(species,riverN,lagRiverN,season,isYOY) %>%
summarize( meanGr = mean( grLength, na.rm = TRUE ),
sdGr = sd( grLength, na.rm = TRUE ),
n = n()) %>%
mutate( meanGrUp = meanGr + sdGr,
meanGrDn = meanGr - sdGr
)
#
ggplot( cdMeansByisYOYriverNMove %>% filter(!(season == 2 & isYOY == 1), species == "bkt", isYOY == 2) , aes(season, meanGr, color = factor(lagRiverN)) ) +
geom_errorbar(aes(ymin = meanGrDn, ymax = meanGrUp), width = 0.25, position = pd) +
geom_point( size = 2 ) +
geom_line()+
facet_wrap(~riverN)
d2 %>% filter(riverN == 1 & lagRiverN == 4) %>% select(tag,river,sampleNumber)
d2 %>% filter(riverN == 2, is.na(lagRiverN), !is.na(grLength)) %>% select(tag,riverN,sampleNumber,lagRiverN,observedLength,lagObservedLength,gr,grLength)
d2 %>% filter(tag=='00088cc059') %>% select(tag,riverOrdered,riverN,lagRiverN,sampleNumber,observedLength,lagObservedLength,gr,grLength)
bletcher/linkedModels documentation built on May 14, 2019, 5:24 p.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.
riverOrderedIn <- factor(c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),labels = c("WB","OL","OS","IL"), ordered = T)
Plots for the manuscript
library(linkedModels) library(arrayhelpers) library(tidyverse) library(gridExtra) load("D:/projects/linkedModels/data/out/ddD_bktbntats_1997.RData") load("D:/projects/linkedModels/data/out/rawCountsAndMasses.RData") mcmc <- list() dGIn <- list() ddGIn2 <- list() load("D:/projects/linkedModels/data/out/dG_1530971738_crossValFALSE_bkt1997_Nimble.RData") # input data d1 <- ddG[[1]][[1]] d2 <- ddG[[1]][[2]] meanSampleIntervalMean <- data.frame( meanSampleInterval = colMeans(d1$sampleIntervalMean), season = 1:d1$nSeasons ) #RECALCULATE gr for dataset. Filtering by maxSampleInterval in the creat ddddG step leaves in trailing lengths from the long samp intervals d2 <- d2 %>% group_by(tag) %>% # arrange(tag,sampleNumber) %>% mutate( lagObservedWeight = lead(observedWeight), lagObservedLength = lead(observedLength), grWeight = exp(lagObservedWeight - observedWeight)/as.numeric((lagDetectionDate - detectionDate)), grLength = (lagObservedLength - observedLength)/as.numeric((lagDetectionDate - detectionDate)) ) %>% ungroup() #d2 %>% filter(tag=='00088cc059') %>% select(tag,riverOrdered,riverN,lagRiverN,sampleNumber,observedLength,lagObservedLength) #input data frame ddGIn2[['bkt']] <- ddG[[1]][[2]] dGIn[['bkt']] <- dG # output data mcmc[['bkt']] <- mcmcProcessed load("D:/projects/linkedModels/data/out/dG_1530971744_crossValFALSE_bnt1997_Nimble.RData") ddGIn2[['bnt']] <- ddG[[1]][[2]] dGIn[['bnt']] <- dG mcmc[['bnt']] <- mcmcProcessed load("D:/projects/linkedModels/data/out/dG_1530993703_crossValFALSE_ats1997_Nimble.RData") ddGIn2[['ats']] <- ddG[[1]][[2]] dGIn[['ats']] <- dG mcmc[['ats']] <- mcmcProcessed rm(ddG,dG,mcmcProcessed) load("D:/projects/linkedModels/data/cd_west_1997.RData")
ggplot(d2, aes(cBNTStd)) + geom_histogram() + facet_grid(season~riverOrdered) ggplot(d2, aes(cBNTStd, grLength)) + geom_point(alpha=0.1) + geom_smooth(method='lm') + facet_grid(isYOY+season~riverOrdered) ggplot(d2, aes(cBKTStd, grLength)) + geom_point(alpha=0.1) + geom_smooth(method='lm') + facet_grid(isYOY+season~riverOrdered)
len <- list() for( spp in c("bkt","bnt","ats") ){ # len2 <- array2df(mcmc[[spp]]$sims.list$lengthExp, label.x = "est") %>% # rename(iter=d1,obs=d2) offset <- 1 if(spp == 'ats') offset <- 0 # not sure why this is needed. REs seem meesed up in the graphs below inds <- data.frame(ind = dGIn[[spp]][[1]]$constants$ind[1:(length(dGIn[[spp]][[1]]$constants$ind) - offset)], observedLength = dGIn[[spp]][[1]]$data$lengthDATA, sampleNumber = ddGIn2[[spp]]$sampleNumber[1:(length(ddGIn2[[spp]]$sampleNumber) - offset)], cohort = ddGIn2[[spp]]$cohort[1:(length(ddGIn2[[spp]]$cohort) - offset)], ageInSamples = ddGIn2[[spp]]$ageInSamples[1:(length(ddGIn2[[spp]]$ageInSamples) - offset)], riverOrdered = ddGIn2[[spp]]$riverOrdered[1:(length(ddGIn2[[spp]]$riverOrdered) - offset)], lagObservedLength = ddGIn2[[spp]]$lagObservedLength[1:(length(ddGIn2[[spp]]$lagObservedLength) - offset)], detectionDate = ddGIn2[[spp]]$detectionDate[1:(length(ddGIn2[[spp]]$detectionDate) - offset)], lagDetectionDate = ddGIn2[[spp]]$lagDetectionDate[1:(length(ddGIn2[[spp]]$lagDetectionDate) - offset)]) inds$obs <- 1:(length(dGIn[[spp]][[1]]$constants$ind) - offset) inds <- inds %>% group_by(ind) %>% mutate(first = min(obs), firstObs = obs == first, observed = !is.na(observedLength), sampleNumber = sampleNumber, nObs = n(), grLength = (lagObservedLength - observedLength)/as.numeric((lagDetectionDate - detectionDate)) ) %>% ungroup() # len1 <- left_join(len2,inds) # rm(len2,inds) # means lenMean1 <- array2df(mcmc[[spp]]$mean$lengthExp, label.x = "estMean") %>% rename(obs=d1) len[[spp]] <- left_join(inds,lenMean1) %>% mutate( estMean = ifelse( firstObs, observedLength, estMean ) ) grIndREMean1 <- array2df(mcmc[[spp]]$mean$grIndRE, label.x = "grIndRE") %>% rename(ind=d1) len[[spp]] <- left_join(len[[spp]],grIndREMean1) # len <- left_join(len1,lenMean1) # rm(len1,lenMean1) rm(inds,lenMean1) gc() } len[['bkt']] %>% filter( nObs == 7 ) %>% ggplot(aes(ageInSamples,observedLength)) + geom_point(aes(color = factor(ind)), size=3) + geom_point(aes(ageInSamples, estMean, color=factor(ind))) + geom_line(aes(ageInSamples, estMean, color = factor(ind))) + geom_point(aes(ageInSamples, estMean, color=factor(ind)), size=2) + geom_point(aes(ageInSamples, estMean, color=factor(ind)), shape = 1, size = 2, color = 'black') + theme(legend.position="none") + facet_wrap(~cohort) focalFish <- 60 len[['bkt']] %>% filter( cohort == 2001 ) %>% ggplot(aes(ageInSamples,observedLength)) + geom_point(aes(color = riverOrdered), size = 3) + geom_point(color = 'black', shape = 1, size = 3) + geom_point(aes(ageInSamples, estMean, color = riverOrdered)) + geom_line( aes(ageInSamples, estMean, group = factor(ind), color = riverOrdered)) + geom_point(aes(ageInSamples, estMean, color = riverOrdered), size = 2) + geom_point(aes(ageInSamples, estMean, color = riverOrdered), shape = 1, size = 2, color = 'black') + # theme(legend.position="none") + facet_wrap(~nObs) len[['bkt']] %>% filter( cohort == 2001 ) %>% ggplot(aes(ageInSamples,observedLength)) + geom_point(aes(color = grIndRE), size = 3) + # geom_point(color = 'black', shape = 1, size = 3) + geom_point(aes(ageInSamples, estMean, color = grIndRE)) + geom_line( aes(ageInSamples, estMean, group = factor(ind), color = grIndRE)) + geom_point(aes(ageInSamples, estMean, color = grIndRE), size = 2) + # geom_point(aes(ageInSamples, estMean, color = grIndRE), shape = 1, size = 2, color = 'black') + # scale_colour_distiller(palette = "Spectral") + scale_colour_gradient(low = "white", high = "black") + # theme(legend.position="none") + facet_wrap(~nObs) len[['bkt']] %>% #filter( cohort == 2001 ) %>% ggplot(aes(ageInSamples,observedLength)) + geom_point(aes(color = grIndRE), size = 3) + # geom_point(color = 'black', shape = 1, size = 3) + geom_point(aes(ageInSamples, estMean, color = grIndRE)) + geom_line( aes(ageInSamples, estMean, group = factor(ind), color = grIndRE)) + geom_point(aes(ageInSamples, estMean, color = grIndRE), size = 2) + # geom_point(aes(ageInSamples, estMean, color = grIndRE), shape = 1, size = 2, color = 'black') + # scale_colour_distiller(palette = "Spectral") + scale_colour_gradient(low = "white", high = "black") + # theme(legend.position="none") + facet_grid(riverOrdered~cohort) len[['ats']] %>% filter( nObs > 0 ) %>% group_by(ind) %>% summarize( meanGRLength = mean(grLength,na.rm=TRUE), meanGRIndRE = mean(grIndRE,na.rm=TRUE), meanNObs = mean(nObs,na.rm=TRUE)) %>% ggplot(aes(meanGRLength,meanGRIndRE)) + geom_point() + geom_smooth(method='lm') + # scale_colour_gradient(low = "yellow", high = "red") + facet_wrap(~meanNObs)
ggplot(ddG[[1]][[2]], aes(observedLength,grLength)) + geom_point() + facet_grid(river + isYOY ~ season+species)
allFishBySpeciesYOY$riverOrderedGG <- factor(allFishBySpeciesYOY$river, levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"), labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T) allFishBySpeciesYOY$seasonGG <- factor(allFishBySpeciesYOY$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T) allFishBySpeciesYOY$speciesGG <- factor(allFishBySpeciesYOY$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T) #sum across YOY allFishBySpecies <- allFishBySpeciesYOY %>% group_by(river,season,year,species) %>% summarize( nAllFishBySpecies = sum(nAllFishBySpeciesYOY, na.rm=T), massAllFishBySpecies = sum(massAllFishBySpeciesYOY, na.rm=T)) allFishBySpeciesYOY <- left_join(allFishBySpeciesYOY,allFishBySpecies) # sum across YOY and species allFish <- allFishBySpeciesYOY %>% group_by(river,season,year) %>% summarize( nAllFish = sum(nAllFishBySpeciesYOY, na.rm=T), massAllFish = sum(massAllFishBySpeciesYOY, na.rm=T)) cd <- left_join(cd,allFish) allFish$riverOrderedGG <- factor(allFish$river, levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"), labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T) allFish$seasonGG <- factor(allFish$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T) ######### ggplot(allFishBySpeciesYOY, aes(year,nAllFishBySpecies)) + geom_point(aes(color = speciesGG)) + geom_line(aes(color = speciesGG)) + geom_line( data = allFish, aes(year,nAllFish) ) + theme_publication() + labs( x = "Year", y = "Count" ) + facet_grid( seasonGG ~ riverOrderedGG)#, scales = "free_y" ) ggsave('figures/rawCounts.png') # add in sampleNumber tmp <- d2 %>% distinct(season,year,sampleNumber) allFishBySpeciesYOY <- left_join(allFishBySpeciesYOY,tmp) allFish <- left_join(allFish,tmp) ggplot(allFishBySpeciesYOY, aes(sampleNumber,nAllFishBySpeciesYOY)) + geom_point(aes(color = speciesGG)) + geom_line(aes(color = speciesGG)) + geom_line( data = allFish, aes(sampleNumber,nAllFish) ) + theme_publication() + labs( x = "SampleNumber", y = "Count" ) + facet_wrap( ~riverOrderedGG )#, scales = "free_y" )
# Correlation between age classes # justification for not using separate age class abundances # interesting that it is less for ATS # From the bottom of adjustCounts() # # nAllFishBySpeciesPStdBKT_yoy1 nAllFishBySpeciesPStdBKT_yoy2 nAllFishBySpeciesPStdBNT_yoy1 nAllFishBySpeciesPStdBNT_yoy2 # nAllFishBySpeciesPStdBKT_yoy1 1.00 0.96 -0.46 -0.43 # nAllFishBySpeciesPStdBKT_yoy2 **0.96** 1.00 -0.44 -0.38 # nAllFishBySpeciesPStdBNT_yoy1 -0.46 -0.44 1.00 0.94 # nAllFishBySpeciesPStdBNT_yoy2 -0.43 -0.38 * *0.94** 1.00 # nAllFishBySpeciesPStdATS_yoy1 -0.13 -0.17 -0.21 0.00 # nAllFishBySpeciesPStdATS_yoy2 0.06 -0.03 -0.05 0.15 # nAllFishBySpeciesPStdATS_yoy1 nAllFishBySpeciesPStdATS_yoy2 # nAllFishBySpeciesPStdBKT_yoy1 -0.13 0.06 # nAllFishBySpeciesPStdBKT_yoy2 -0.17 -0.03 # nAllFishBySpeciesPStdBNT_yoy1 -0.21 -0.05 # nAllFishBySpeciesPStdBNT_yoy2 0.00 0.15 # nAllFishBySpeciesPStdATS_yoy1 1.00 0.47 # nAllFishBySpeciesPStdATS_yoy2 *0.47* 1.00
d2 %>% group_by(riverOrderedGG,seasonGG,year) %>% filter( riverOrderedGG == "West Brook" ) %>% summarize(meanF = mean(meanFlow, na.rm=T), sdF = sd(meanFlow, na.rm=T)) %>% ggplot( aes(year,meanF) ) + geom_point( size = 3 ) + geom_line() + # geom_smooth(method = 'lm', se=F) + geom_errorbar( aes(ymin = meanF - sdF, ymax = meanF + sdF, width = 0.33 ) ) + labs( x = "Year", y = "Stream flow (x/x)" ) + theme_publication() + facet_wrap(~seasonGG) ggsave('figures/flow.png') lmFunction <- function(df){ lm( meanT ~ year, data=df ) } tempMods <- d2 %>% group_by(riverOrderedGG,seasonGG,year) %>% summarize(meanT = mean(meanTemperature, na.rm=T), sdT = sd(meanTemperature, na.rm=T)) %>% group_by( riverOrderedGG, seasonGG ) %>% nest() %>% mutate( model = map(data,lmFunction) ) tempModsTidy <- tempMods %>% mutate(tidy = map(model,broom::tidy)) %>% dplyr::select(-data,-model) %>% unnest() %>% filter( term == 'year') %>% mutate( sig = ifelse(p.value < 0.05,"*","")) d2 %>% group_by(riverOrderedGG,seasonGG,year) %>% summarize(meanT = mean(meanTemperature, na.rm=T), sdT = sd(meanTemperature, na.rm=T)) %>% ggplot( aes(year,meanT) ) + geom_smooth(method = 'lm', se=F, color = 'darkgrey') + geom_point( ) + geom_line() + geom_errorbar( aes(ymin = meanT - sdT, ymax = meanT + sdT, width = 0.33 ) ) + labs( x = "Year", y = "Stream temperature (C)" ) + geom_text(data=tempModsTidy,aes(x=1999,y=18, label=paste(round(estimate,2),sig))) + theme_publication() + facet_grid(seasonGG ~ riverOrderedGG) ggsave('figures/temp.png')
cd$riverOrderedGG <- factor(cd$river, levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"), labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T) cd$seasonGG <- factor(cd$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T) cd$speciesGG <- factor(cd$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T) pd <- position_dodge(0.65) cd %>% group_by(speciesGG,riverOrderedGG,seasonGG) %>% filter(isYOY == 1) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) %>% ggplot( aes(seasonGG,meanGR, color=speciesGG, shape=speciesGG) ) + geom_point( position = pd, size = 2 ) + geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) + ggtitle( "isYOY = 1" ) + labs( x = "Season", y = "Growth rate (mm/d)" ) + #scale_x_continuous( breaks = 1997:2015 ) + theme_publication() + theme( panel.grid.minor = element_blank() ) + facet_wrap( ~ riverOrderedGG) ggsave('figures/meanGRYOY1.png', width = 7, height = 7, units='in') cd %>% group_by(speciesGG,riverOrderedGG,seasonGG) %>% filter(isYOY == 2) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) %>% ggplot( aes(seasonGG,meanGR, color=speciesGG, shape=speciesGG) ) + geom_point( position = pd,size = 2 ) + geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) + ggtitle( "isYOY = 2" ) + labs( x = "Season", y = "Growth rate (mm/d)" ) + #scale_x_continuous( breaks = 1997:2015 ) + theme_publication() + theme( panel.grid.minor = element_blank() ) + facet_wrap( ~ riverOrderedGG) ggsave('figures/meanGRYOY2.png', width = 7, height = 7, units='in') # mean growth by species
cd %>% group_by(species) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) .162/.194 #bkt/bnt .151/.194 #ats/bnt cd %>% group_by(isYOY) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) .119/.215 #juv/yoy .215/.119 cd %>% group_by(river) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) # river meanGR sdGR # <chr> <dbl> <dbl> # 1 wb jimmy 0.132 0.131 # 2 wb mitchell 0.179 0.160 # 3 wb obear 0.137 0.129 # 4 west brook 0.170 0.175 ((.132+.137)/2)/((.17+.179)/2) ((.17+.179)/2)/((.132+.137)/2) cd %>% group_by(season) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) ((.0464+.0894+.0914)/3)/(.417) .417/((.0464+.0894+.0914)/3)
cd %>% filter() %>% ggplot(aes(observedLength,grLength)) + geom_point() + facet_grid(season~species)
# ggPInt <- array2df(ddD$sims.list$pBetaInt, label.x = "est") # # ggPInt$chain <- rep(1:ddD$mcmc.info$n.chains, each = ddD$mcmc.info$n.samples/ddD$mcmc.info$n.chains) # ggPInt$iter <- 1:as.numeric(ddD$mcmc.info$n.samples/ddD$mcmc.info$n.chains) # # ggplot(filter(ggPInt), aes(iter,est)) + # geom_point( aes(color=factor(chain)), size = 0.1 ) + # ylim(-1.5,1.5) + # facet_grid(d2+d4~d5+d3) p <- array2df(ddD$q50$pBetaInt, label.x = "est") %>% rename(est=est,isYOY=d1,species=d2,season=d3,river=d4,year=d5) ggplot(p, aes(year,invlogit(est),color=factor(species))) +geom_point() + facet_grid(river ~ season + isYOY) # seem to need to run the ggsave outside of markdown by pasting into console ggsave('figures/p.png') phi <- array2df(ddD$q50$phiBetaInt, label.x = "est") %>% rename(est=est,isYOY=d1,species=d2,season=d3,river=d4,year=d5) ggplot(phi, aes(year,invlogit(est),color=factor(species))) +geom_point() + facet_grid(river ~ season + isYOY) ggsave('figures/phi.png') # # of skunked sections tmp <- dddD[[2]] %>% group_by(riverOrdered,year,season,section) %>% summarize( s = sum(enc) ) %>% filter( s == 0 ) table(tmp$riverOrdered,tmp$year,tmp$season)
rHat <- list() rHat[['bkt']] <- array2df(mcmc[['bkt']]$Rhat$grBeta, levels = list(beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="rHat") rHat[['bnt']] <- array2df(mcmc[['bnt']]$Rhat$grBeta, levels = list(beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="rHat") rHat[['ats']] <- array2df(mcmc[['ats']]$Rhat$grBeta, levels = list(beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="rHat") # rHat[['bkt']] <- array2df(mcmc[['bkt']]$Rhat$grBeta, levels = list(beta=1:8,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="rHat") # rHat[['bnt']] <- array2df(mcmc[['bnt']]$Rhat$grBeta, levels = list(beta=1:8,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="rHat") # rHat[['ats']] <- array2df(mcmc[['ats']]$Rhat$grBeta, levels = list(beta=1:8,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="rHat") ggplot(rHat[["ats"]], aes(beta,rHat)) + geom_point() + facet_grid(river+isYOY~season)
grInt <- list() grInt[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") grInt[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") grInt[['ats']] <- array2df(mcmc[['ats']]$sims.list$grInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") ggplot( grInt[['bkt']] %>% filter(river %in% c('west brook','wb jimmy')), aes(int)) + geom_freqpoly() + geom_freqpoly(aes((int)),color = 'brown', data = grInt[['bnt']]) + geom_freqpoly(aes((int)),color = 'lightblue', data = grInt[['ats']]) + # xlim(-0.5,1.5) + theme_publication() + facet_grid(season ~ river + isYOY) ################################################## # as range plot plotHPD <- function(sppToPlot){ nRivers <- list() nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1 probs <- c(0.05,0.1,0.5,0.9,0.95) q <- grInt[[sppToPlot]] %>% group_by(isYOY,season,river) %>% summarize( p = list(probs), q = list(quantile(int,probs)) ) %>% unnest() %>% spread(key = p, value = q) p <- array2df(mcmc[[sppToPlot]]$mean$grInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="int") %>% mutate( season = as.numeric(season) ) %>% left_join(.,meanSampleIntervalMean, by = 'season') p$intByDay <- p$int/p$meanSampleInterval p$intByDay <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDay) p$parameter <- paste(p$river,p$season,sep="_") pUp <- array2df(mcmc[[sppToPlot]]$q97.5$grInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intUp") %>% mutate( season = as.numeric(season) ) pDown <- array2df(mcmc[[sppToPlot]]$q2.5$grInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intDown") %>% mutate( season = as.numeric(season) ) p <- left_join(p,pUp); p <- left_join(p,pDown) p$intByDayUp <- p$intUp/p$meanSampleInterval; p$intByDayDown <- p$intDown/p$meanSampleInterval p$intByDay <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDay) p$intByDayUp <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDayUp) p$intByDayDown <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDayDown) #print(c(min(p$intByDayDown,na.rm=T),max(p$intByDayDown,na.rm=T))) gg <- ggplot(p, aes(x = intByDay, y = reorder(parameter,intByDay,na.rm=TRUE), color = factor(season), shape = factor(isYOY))) + theme_publication() + geom_point(size = 3) + geom_segment(aes(x = intByDay, xend = intByDayUp, yend = reorder(parameter, intByDay)), size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + geom_segment(aes(x = intByDay, xend = intByDayDown, yend = reorder(parameter, intByDay)), size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + xlim(c(-0.0015,0.01)) + xlab("Mean daily growth (mm/d)") + ylab("River_Season") # ggtitle(sppToPlot) return(gg) } plotHPD('bkt'); ggsave('figures/HDP_bkt.png', height = 3, width = 3, scale = 2) plotHPD('bnt'); ggsave('figures/HDP_bnt.png', height = 3, width = 3, scale = 2) plotHPD('ats'); ggsave('figures/HDP_ats.png', height = 3, width = 3, scale = 2)
sigmaInt <- list() sigmaInt[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$sigmaInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") sigmaInt[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$sigmaInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") sigmaInt[['ats']] <- array2df(mcmc[['ats']]$sims.list$sigmaInt, levels = list(iter=NA,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") bw <- 0.05 ggplot( sigmaInt[['bkt']] %>% filter(river %in% c('west brook','wb jimmy')), aes(int)) + geom_freqpoly(binwidth = bw) + geom_freqpoly(aes(int),binwidth = bw,color = 'brown', data = sigmaInt[['bnt']]) + geom_freqpoly(aes(int),binwidth = bw,color = 'lightblue', data = sigmaInt[['ats']]) + #xlim(-4.5,-1) + xlim(-5,2.5) + theme_publication() + facet_grid(season ~ river + isYOY) ################################################## # as range plot plotHPDSigma <- function(sppToPlot){ nRivers <- list() nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1 probs <- c(0.05,0.1,0.5,0.9,0.95) q <- sigmaInt[[sppToPlot]] %>% group_by(isYOY,season,river) %>% summarize( p = list(probs), q = list(quantile(int,probs)) ) %>% unnest() %>% spread(key = p, value = q) p <- array2df(mcmc[[sppToPlot]]$mean$sigmaInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="int") %>% mutate( season = as.numeric(season) ) %>% left_join(.,meanSampleIntervalMean, by = 'season') # p$intByDay <- p$int/p$meanSampleInterval # p$intByDay <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intByDay) p$parameter <- paste(p$river,p$season,sep="_") pUp <- array2df(mcmc[[sppToPlot]]$q97.5$sigmaInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intUp") %>% mutate( season = as.numeric(season) ) pDown <- array2df(mcmc[[sppToPlot]]$q2.5$sigmaInt, levels = list(isYOY=c(1,2),season=1:4,river=riverOrderedIn[1:nRivers[[sppToPlot]]]), label.x="intDown") %>% mutate( season = as.numeric(season) ) p <- left_join(p,pUp); p <- left_join(p,pDown) # p$intByDayUp <- p$intUp/p$meanSampleInterval; p$intByDayDown <- p$intDown/p$meanSampleInterval p$int <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$int) p$intUp <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intUp) p$intDown <- ifelse(p$isYOY == 1 & p$season == 2, NA, p$intDown) #print(c(min(p$intDown,na.rm=T),max(p$intDown,na.rm=T))) gg <- ggplot(p, aes(x = exp(int), y = reorder(parameter,int,na.rm=TRUE), color = factor(season), shape = factor(isYOY))) + theme_publication() + geom_point(size = 3) + geom_segment(aes(x = exp(int), xend = exp(intUp), yend = reorder(parameter, int)), size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + geom_segment(aes(x = exp(int), xend = exp(intDown), yend = reorder(parameter, int)), size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + #xlim(c(-0.0015,0.01)) + xlab("Standard deviation of daily growth (mm/d)") + ylab("River_Season") # ggtitle(sppToPlot) return(gg) } plotHPDSigma('bkt'); ggsave('figures/HDPSigma_bkt.png', height = 3, width = 3, scale = 2) plotHPDSigma('bnt'); ggsave('figures/HDPSigma_bnt.png', height = 3, width = 3, scale = 2) plotHPDSigma('ats'); ggsave('figures/HDPsigma_ats.png', height = 3, width = 3, scale = 2)
#riverOrderedIn <- factor(c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),labels = c("west brook","wb jimmy","wb mitchell","wb obear"), ordered = T) grIntMu <- list() grIntMu[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grIntMu, levels = list(iter=NA,season=1:4), label.x="int") grIntMu[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grIntMu, levels = list(iter=NA,season=1:4), label.x="int") grIntMu[['ats']] <- array2df(mcmc[['ats']]$sims.list$grIntMu, levels = list(iter=NA,season=1:4), label.x="int") bw <- 0.05 ggplot( grIntMu[['bkt']] , aes(int)) + geom_freqpoly(binwidth = bw, size=1.2) + geom_freqpoly(aes(int),binwidth = bw,color = 'brown', data = grIntMu[['bnt']]) + geom_freqpoly(aes(int),binwidth = bw,color = 'lightblue', data = grIntMu[['ats']]) + #xlim(-4.5,-1) + xlim(-0.25,0.75) + # ylim(0,2500) + theme_publication() + facet_wrap(~season )
# get all betas in a DF grBeta <- list() grBeta[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grBeta, levels = list(iter=NA,beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>% mutate( species = 'bkt') grBeta[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grBeta, levels = list(iter=NA,beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>% mutate( species = 'bnt') grBeta[['ats']] <- array2df(mcmc[['ats']]$sims.list$grBeta, levels = list(iter=NA,beta=1:7,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>% mutate( species = 'ats') grBetaDF1 <- rbind(grBeta[['bkt']],grBeta[['bnt']],grBeta[['ats']]) grBetaBNT <- list() grBetaBNT[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grBetaBNT, levels = list(iter=NA,beta=8:9,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>% mutate( species = 'bkt') grBetaBNT[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grBetaBNT, levels = list(iter=NA,beta=8:9,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>% mutate( species = 'bnt') grBetaBNT[['ats']] <- array2df(mcmc[['ats']]$sims.list$grBetaBNT, levels = list(iter=NA,beta=8:9,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>% mutate( species = 'ats') grBetaBNTDF <- rbind(grBetaBNT[['bkt']],grBetaBNT[['bnt']],grBetaBNT[['ats']]) grBetaATS <- list() grBetaATS[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$grBetaATS, levels = list(iter=NA,beta=10,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>% mutate( species = 'bkt') grBetaATS[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$grBetaATS, levels = list(iter=NA,beta=10,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>% mutate( species = 'bnt') grBetaATS[['ats']] <- array2df(mcmc[['ats']]$sims.list$grBetaATS, levels = list(iter=NA,beta=10,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>% mutate( species = 'ats') grBetaATSDF <- rbind(grBetaATS[['bkt']],grBetaATS[['bnt']],grBetaATS[['ats']]) grBetaDF <- rbind(grBetaDF1,grBetaBNTDF,grBetaATSDF) ################################################## # as range plot nRivers <- list() nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1 speciesList <- list('bkt','bnt','ats') probs <- c(0.05,0.1,0.5,0.9,0.95) nSpecies <- 3 q <- grBetaDF %>% group_by(beta,isYOY,season,river,species) %>% summarize( p = list(probs), q = list(quantile(int,probs)) ) %>% unnest() %>% spread(key = p, value = q) %>% data.frame() %>% rename( median = X0.5, lo = X0.05, hi = X0.95 ) q <- q %>% filter( !(season == 2 & isYOY == 0)) q$riverGG <- factor(q$river, levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"), labels = c("WB","OL","OS","IL"), ordered = T) q$seasonGG <- factor(q$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T) q$speciesGG <- factor(q$species, levels = c('bkt','bnt','ats'), labels = c("BKT", "BNT", "ATS"), ordered = T) q$betaGG <- factor(q$beta, levels = 1:10, levels = c("Temp","Flow","Temp*Flow","Temp^2","Flow^2","BKT","Length","BNT","BNT*BKT","ATS"), # make sure this order matches the nimble code labels = c("Temp","Flow","Temp*Flow","Temp^2","Flow^2","BKT","Length","BNT","BNT*BKT","ATS"),ordered = T) q$parameter = paste0(q$riverGG,"_",q$speciesGG) q$parameterOrd <- factor(q$parameter, levels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), labels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), ordered = T) q <- q %>% filter( !(betaGG == 'ATS' & (riverGG == "OS" | riverGG == "IL" | riverGG == "OL")) ) q <- q %>% filter( !(betaGG == 'BNT' & (riverGG == "OS" | riverGG == "IL")) ) q <- q %>% filter( !(betaGG == 'BNT*BKT' & (riverGG == "OS" | riverGG == "IL")) ) # q <- q %>% filter( !((riverGG == "OS" | riverGG == "IL")) ) #reorder(parameter,median,na.rm=TRUE) (min(q$median)) (max(q$median)) # q[q$median < -0.13,] ggplot(q, aes(x = median, y = parameterOrd, color = factor(speciesGG), shape = factor(isYOY))) + geom_vline(xintercept = 0,color='black') + theme_publication() + geom_point(size = 3) + geom_segment(aes(x = median, xend = hi, yend = parameterOrd), size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + geom_segment(aes(x = median, xend = lo, yend = parameterOrd), size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + xlim(c(-0.15,0.11)) + xlab("Parameter estimate") + ylab("River_Species") + facet_grid(seasonGG~betaGG) ggsave('figures/grBetas.png', height = 7, width = 10, scale = 2)
# get all betas in a DF grBetaSigma <- list() grBetaSigma[['bkt']] <- array2df(mcmc[['bkt']]$sims.list$sigmaBeta, levels = list(iter=NA,beta=1:3,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:4]), label.x="int") %>% mutate( species = 'bkt') grBetaSigma[['bnt']] <- array2df(mcmc[['bnt']]$sims.list$sigmaBeta, levels = list(iter=NA,beta=1:3,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:2]), label.x="int") %>% mutate( species = 'bnt') grBetaSigma[['ats']] <- array2df(mcmc[['ats']]$sims.list$sigmaBeta, levels = list(iter=NA,beta=1:3,isYOY=c(0,1),season=1:4,river=riverOrderedIn[1:1]), label.x="int") %>% mutate( species = 'ats') grBetaSigmaDF <- rbind(grBetaSigma[['bkt']],grBetaSigma[['bnt']],grBetaSigma[['ats']]) ################################################## # as range plot nRivers <- list() nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1 speciesList <- list('bkt','bnt','ats') probs <- c(0.05,0.1,0.5,0.9,0.95) nSpecies <- 3 q <- grBetaSigmaDF %>% group_by(beta,isYOY,season,river,species) %>% summarize( p = list(probs), q = list(quantile(int,probs)) ) %>% unnest() %>% spread(key = p, value = q) %>% data.frame() %>% rename( median = X0.5, lo = X0.05, hi = X0.95 ) q <- q %>% filter( !(season == 2 & isYOY == 0)) q$riverGG <- factor(q$river, levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"), labels = c("WB","OL","OS","IL"), ordered = T) q$seasonGG <- factor(q$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T) q$speciesGG <- factor(q$species, levels = c('bkt','bnt','ats'), labels = c("BKT", "BNT", "ATS"), ordered = T) q$betaGG <- factor(q$beta, levels = 1:10, labels = c("Temp", "Flow", "Temp*Flow","Temp^2","Flow^2","BKT","Length","BNT","BNT*BKT","ATS"), ordered = T) q$parameter = paste0(q$riverGG,"_",q$speciesGG) q$parameterOrd <- factor(q$parameter, levels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), labels = c("IL_BKT", "OS_BKT", "OL_BNT", "OL_BKT", "WB_ATS", "WB_BNT", "WB_BKT"), ordered = T) q <- q %>% filter( !(betaGG == 'ATS' & (riverGG == "OS" | riverGG == "IL" | riverGG == "OL")) ) q <- q %>% filter( !(betaGG == 'BNT' & (riverGG == "OS" | riverGG == "IL")) ) q <- q %>% filter( !(betaGG == 'BNT*BKT' & (riverGG == "OS" | riverGG == "IL")) ) # q <- q %>% filter( !((riverGG == "OS" | riverGG == "IL")) ) #reorder(parameter,median,na.rm=TRUE) (min(q$median)) (max(q$median)) # q[q$median < -0.13,] ggplot(q, aes(x = median, y = parameterOrd, color = factor(speciesGG), shape = factor(isYOY))) + geom_vline(xintercept = 0,color='black') + theme_publication() + geom_point(size = 3) + geom_segment(aes(x = median, xend = hi, yend = parameterOrd), size = 0.5, arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + geom_segment(aes(x = median, xend = lo, yend = parameterOrd), size = 0.5,arrow = arrow(angle = 90, length = unit(0.05, "inches"))) + xlim(c(-0.35,0.5)) + xlab("Parameter estimate") + ylab("River_Species") + facet_grid(seasonGG~betaGG) ggsave('figures/grBetaSigmas.png', height = 7, width = 7, scale = 2)
grIndRE <- list() grIndRE[['bkt']] <- as.data.frame(mcmc[['bkt']]$q50$grIndRE); names(grIndRE[['bkt']]) <- 'value' grIndRE[['bnt']] <- as.data.frame(mcmc[['bnt']]$q50$grIndRE); names(grIndRE[['bnt']]) <- 'value' grIndRE[['ats']] <- as.data.frame(mcmc[['ats']]$q50$grIndRE); names(grIndRE[['ats']]) <- 'value' table(mcmc[['bkt']]$q50$grIndRE==0)/length(mcmc[['bkt']]$q50$grIndRE) table(mcmc[['bnt']]$q50$grIndRE==0)/length(mcmc[['bnt']]$q50$grIndRE) table(mcmc[['ats']]$q50$grIndRE==0)/length(mcmc[['ats']]$q50$grIndRE) ggplot( grIndRE[[1]], aes(value,..density..) ) + geom_freqpoly(binwidth = 0.00025) + geom_freqpoly(aes(value), color = "brown", binwidth = 0.00025, data = grIndRE[[2]]) + geom_freqpoly(aes(value), color = "lightblue", binwidth = 0.00025, data = grIndRE[[3]]) + ylim(c(0,100))
# #####run if get new predictions - takes a while to run #### # load('D:/projects/linkedModels/data/out/P_ForMike_bkt.RData') # pBoth$species <- 'bkt' # pAll <- pBoth; rm(pBoth) # # load('D:/projects/linkedModels/data/out/P_ForMike_bnt.RData') # pBoth$species <- 'bnt' # pAll <- rbind(pAll,pBoth); rm(pBoth) # # load('D:/projects/linkedModels/data/out/P_ForMike_ats.RData') # pBoth$species <- 'ats' # pAll <- rbind(pAll,pBoth); rm(pBoth) # # pred <- pAll %>% group_by(species,len,flow,temp,cBKT,cBNT,cATS,isYOY,season,river) %>% # summarize(meanPredGR = mean(predGr, na.rm=T), # sdPredGR = sd(predGr, na.rm=T), # meanPredGRSigma = mean(predGrSigma, na.rm=T), # sdPredGRSigma = sd(predGrSigma, na.rm=T), # meanPredGRCV = mean(predCV, na.rm=T), # sdPredGRCV = sd(predCV, na.rm=T)) %>% # ungroup() # # rm(pAll) # save(pred,file = 'D:/projects/linkedModels/data/out/P_ForMike_allMeans.RData') load('D:/projects/linkedModels/data/out/P_ForMike_allMeans.RData') pred <- pred %>% filter( !(isYOY == 0 & season == 2) ) #pred <- pred %>% filter( !(isYOY == 0 & river == 'wb jimmy' & season == 2) ) #graph function, in analyzeOutputFunctions.R pred$riverGG <- factor(pred$river, levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"), labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T) pred$seasonGG <- factor(pred$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T) pred$speciesGG <- factor(pred$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T) plotPredMeans(pred, "meanPredGR", "len", 0) plotPredMeans(pred, "meanPredGR", "len", 1) plotPredMeans(pred, "meanPredGR", "flow", 0) plotPredMeans(pred, "meanPredGR", "flow", 1) plotPredMeans(pred, "meanPredGR", "temp", 0) plotPredMeans(pred, "meanPredGR", "temp", 1) plotPredMeans(pred, "meanPredGR", "cBKT", 0) plotPredMeans(pred, "meanPredGR", "cBKT", 1) plotPredMeans(pred, "meanPredGR", "cBNT", 0) plotPredMeans(pred, "meanPredGR", "cBNT", 1) plotPredMeans(pred, "meanPredGR", "cATS", 0) plotPredMeans(pred, "meanPredGR", "cATS", 1) plotPredMeans(pred, "meanPredGR", c("temp","flow"), 0) plotPredMeans(pred, "meanPredGR", c("temp","flow"), 1) plotPredMeans(pred, "meanPredGR", c("flow","temp"), 0) plotPredMeans(pred, "meanPredGR", c("flow","temp"), 1) ggsave('figures/tmp.png', width = 9, height = 7, units='in')
all <- c('len','temp','flow','cBKT','cBNT','cATS') allLabels <- c('Length (mm)','Temperature (C)','Flow (m3/s)','Brook Trout abundance','Brown Trout abundance','Atlantic Salmon abundance') all <- c('len','cBKT','temp','cBNT','flow','cATS') allLabels <- c('Length (mm)','Brook Trout abundance','Stream temperature (C)','Brown Trout abundance','Stream flow (m3/s)','Atlantic Salmon abundance') remFacetLabel <- c(3,4,5,6) selectSingleVar <- function(p, all, depVar, varsToPlot, isYOYGG){ notPlot <- NA notPlot[1] <- all[!(all %in% varsToPlot)][1] notPlot[2] <- all[!(all %in% varsToPlot)][2] notPlot[3] <- all[!(all %in% varsToPlot)][3] notPlot[4] <- all[!(all %in% varsToPlot)][4] notPlot[5] <- all[!(all %in% varsToPlot)][5] pGG <- p %>% filter(isYOY == isYOYGG, eval(as.name(notPlot[1])) == 0, eval(as.name(notPlot[2])) == 0, eval(as.name(notPlot[3])) == 0, eval(as.name(notPlot[4])) == 0, eval(as.name(notPlot[5])) == 0 )#%>% # distinct(species, isYOY, riverGG, seasonGG, speciesGG, .keep_all = TRUE) return(pGG) } ##### by species isYOYGG <- 0 riverGGToPlot <- "West Brook" #"Open Large" outputToPlot <- "meanPredGR" ggSingle <- list() singleDat <- list() for(i in 1:length(all)) #https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r local({ i <- i singleDat[[i]] <- selectSingleVar(pred, all, outputToPlot, all[i], isYOYGG) print(c(i,all[i])) # print(str(singleDat[[i]])) tmp <- singleDat[[i]] %>% filter(riverGG == riverGGToPlot) ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=speciesGG ) ) + geom_line( size = 1.225 ) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.title.y = element_blank(), legend.position="none") + geom_hline(yintercept = 0) + ylim(-0.15,0.8) + scale_x_continuous(allLabels[i]) + facet_grid(~seasonGG) if( i %in% remFacetLabel ) { ggSingle[[i]] <<- ggSingle[[i]] + theme( strip.background = element_blank(), strip.text.x = element_blank() ) } }) grid.arrange(grobs=ggSingle, left = "Growth rate (mm/d)") ###### by river for bkt isYOYGG <- 1 speciesGGToPlot <- "Brook Trout" outputToPlot <- "meanPredGR" ggSingle <- list() singleDat <- list() for(i in 1:length(all)) #https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r local({ i <- i singleDat[[i]] <- selectSingleVar(pred, all, outputToPlot, all[i], isYOYGG) print(c(i,all[i])) # print(str(singleDat[[i]])) tmp <- singleDat[[i]] %>% filter(speciesGG == speciesGGToPlot) ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=riverGG ) ) + geom_line( size = 1.225 ) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.title.y = element_blank(), legend.position="none" ) + geom_hline(yintercept = 0) + ylim(-0.15,0.8) + scale_x_continuous(allLabels[i]) + facet_grid(~seasonGG) if( i %in% remFacetLabel ) { ggSingle[[i]] <<- ggSingle[[i]] + theme( strip.background = element_blank(), strip.text.x = element_blank() ) } }) grid.arrange(grobs=ggSingle, left = "Growth rate (mm/d)")
selectPairVars <- function(p, all, depVar, varsToPlot, isYOYGG){ notPlot <- NA notPlot[1] <- all[!(all %in% varsToPlot)][1] notPlot[2] <- all[!(all %in% varsToPlot)][2] notPlot[3] <- all[!(all %in% varsToPlot)][3] notPlot[4] <- all[!(all %in% varsToPlot)][4] # notPlot[5] <- all[!(all %in% varsToPlot)][5] pGG <- p %>% filter( isYOY == isYOYGG, eval(as.name(notPlot[1])) == 0, eval(as.name(notPlot[2])) == 0, eval(as.name(notPlot[3])) == 0, eval(as.name(notPlot[4])) == 0 ) return(pGG) } h <- list(); h[['bkt']] <- 8.5; h[['bnt']] <- 4.5; h[['ats']] <- 2.75 for(i in c('bkt','bnt', 'ats')){ for(j in 0:1){ print(c(i,j)) p <- selectPairVars(pred, all, "meanPredGR", c("flow","temp"), j) %>% filter(species == i) ggplot(p, aes(x = flow, y = temp, z = meanPredGR)) + geom_tile(aes(fill = meanPredGR)) + stat_contour(bins = 15, color = 'darkgrey') + ggtitle(paste(j,i)) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), #axis.title.y = element_blank(), legend.position="none" ) + facet_grid( riverGG ~ seasonGG ) ggsave(paste0('figures/predGRGrid',"_",i,"_",j,'.png'), width = 9, height = h[[i]], units='in') } } #raw data plot - not very useful ggplot(d2, aes(x=flowStd,y=tempStd,z=grLength)) + # stat_density_2d( color = 'darkgrey') + stat_density_2d(geom = "raster", aes(fill = ..density..), contour = FALSE,interpolate=TRUE,h=2.5) + facet_grid(season+isYOY~river+species) ggplot(d2, aes(x=tempStd,y=grLength)) + geom_smooth() + geom_smooth(method='lm',color='black',se=FALSE) + ylim(c(0,0.8)) + facet_grid(river+isYOY~season)
all <- c('temp','flow') allLabels <- c('Stream temperature (C)','Stream flow (m3/s)') remFacetLabel <- c(2) selectSingleVarSigma <- function(p, all, depVar, varsToPlot, isYOYGG){ notPlot <- NA notPlot[1] <- all[!(all %in% varsToPlot)][1] pGG <- p %>% filter(isYOY == isYOYGG, eval(as.name(notPlot[1])) == 0, len == 0, cBKT == 0, cBNT == 0, cATS == 0 )#%>% # distinct(species, isYOY, riverGG, seasonGG, speciesGG, .keep_all = TRUE) return(pGG) } ##### by species isYOYGG <- 0 riverGGToPlot <- "Open Large" #"West Brook" # outputToPlot <- "meanPredGRSigma" ggSingle <- list() singleDat <- list() for(i in 1:length(all)) #https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r local({ i <- i singleDat[[i]] <- selectSingleVarSigma(pred, all, outputToPlot, all[i], isYOYGG) print(c(i,all[i])) # print(str(singleDat[[i]])) tmp <- singleDat[[i]] %>% filter(riverGG == riverGGToPlot) ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=speciesGG ) ) + geom_line( size = 1.225 ) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.title.y = element_blank(), legend.position="none") + geom_hline(yintercept = 0) + # ylim(-0.15,0.8) + scale_x_continuous(allLabels[i]) + facet_grid(~seasonGG) if( i %in% remFacetLabel ) { ggSingle[[i]] <<- ggSingle[[i]] + theme( strip.background = element_blank(), strip.text.x = element_blank() ) } }) # add thrid row to graph for interaction of tem and flow? # run this and export grid.arrange(grobs=ggSingle, left = "Standard deviation of growth rate (mm/d)") #ggsave(ggCV, paste0('figures/predGRCV_YOY',isYOYGG,'_',riverGGToPlot,'.png'), width = 9, height = 7, units='in') ###### by river for bkt isYOYGG <- 0 speciesGGToPlot <- "Brook Trout" outputToPlot <- "meanPredGRSigma" ggSingle <- list() singleDat <- list() for(i in 1:length(all)) #https://stackoverflow.com/questions/31993704/storing-ggplot-objects-in-a-list-from-within-loop-in-r local({ i <- i singleDat[[i]] <- selectSingleVarSigma(pred, all, outputToPlot, all[i], isYOYGG) print(c(i,all[i])) # print(str(singleDat[[i]])) tmp <- singleDat[[i]] %>% filter(speciesGG == speciesGGToPlot) ggSingle[[i]] <<- ggplot( tmp, aes( eval(as.name(all[i])), eval(as.name(outputToPlot)), color=riverGG ) ) + geom_line( size = 1.225 ) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.title.y = element_blank(), legend.position="none" ) + geom_hline(yintercept = 0) + #ylim(-0.15,0.8) + scale_x_continuous(allLabels[i]) + facet_grid(~seasonGG) if( i %in% remFacetLabel ) { ggSingle[[i]] <<- ggSingle[[i]] + theme( strip.background = element_blank(), strip.text.x = element_blank() ) } }) grid.arrange(grobs=ggSingle, left = "Standard deviation of growth rate (mm/d)")
############################### # flow*temp interaction plots isYOYGG <- 1 riverGGToPlot <- "West Brook" #"Open Large" # outputToPlot <- "meanPredGR" #"meanPredGRSigma" interact <- pred %>% filter(isYOY == isYOYGG, riverGG == riverGGToPlot, len == 0, cBKT == 0, cBNT == 0, cATS == 0 ) %>% mutate( temp_spp = paste0(temp,speciesGG)) ggplot( interact, aes( flow, eval(as.name(outputToPlot)), group=temp_spp, color=speciesGG ) ) + geom_line( aes(alpha = temp + 3), size = 1.225 ) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.title.y = element_blank() # legend.position="none" ) + geom_hline(yintercept = 0) + # ylim(-0.15,0.8) + scale_x_continuous('Flow') + facet_grid(~seasonGG)
############################### # flow*temp interaction plots isYOYGG <- 1 riverGGToPlot <- "Open Large" #"West Brook" # outputToPlot <- "meanPredGR" #"meanPredGRSigma" interact <- pred %>% filter(isYOY == isYOYGG, riverGG == riverGGToPlot, len == 0, temp == 0, flow == 0, cATS == 0 ) %>% mutate( bnt_spp = paste0(cBNT,speciesGG)) ggplot( interact, aes( cBKT, eval(as.name(outputToPlot)), group=bnt_spp, color=speciesGG ) ) + geom_line( aes(alpha = cBNT + 3), size = 1.225 ) + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.title.y = element_blank() # legend.position="none" ) + geom_hline(yintercept = 0) + # ylim(-0.15,0.8) + scale_x_continuous('Brook Trout abundance') + facet_grid(~seasonGG)
round(mcmc$summary$all.chains[startsWith(attributes(mcmc$summary$all.chains)$dimnames[[1]], "grBeta["),],2)
pd <- position_dodge(0.65) # isYOY1 d2 %>% group_by(speciesGG,riverOrderedGG,seasonGG,year) %>% filter(isYOYDATA == 1) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) %>% ggplot( aes(year,meanGR, color=speciesGG, shape=speciesGG) ) + geom_point( position=pd ) + geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) + ggtitle( "isYOY = 1" ) + labs( x = "Year", y = "Growth rate (mm/d)" ) + scale_x_continuous( breaks = 1997:2015 ) + theme_publication() + theme( panel.grid.minor = element_blank(), axis.text.x = element_text(angle=90, vjust=0.5, size=8) ) + facet_grid(seasonGG ~ riverOrderedGG) ggsave('figures/grByYearYOY1.png') # isYOY2 d %>% group_by(speciesGG,riverOrderedGG,seasonGG,year) %>% filter(isYOYDATA == 2) %>% summarize(meanGR = mean(grLength, na.rm=T), sdGR = sd(grLength, na.rm=T)) %>% ggplot( aes(year,meanGR, color=speciesGG, shape=speciesGG) ) + geom_point( position=pd ) + geom_errorbar( aes(ymin = meanGR - sdGR, ymax = meanGR + sdGR), width = 0.33, position=pd ) + ggtitle( "isYOY = 2" ) + labs( x = "Year", y = "Growth rate (mm/d)" ) + scale_x_continuous( breaks = 1997:2015 ) + theme_publication() + theme( panel.grid.minor = element_blank(), axis.text.x = element_text(angle=90, vjust=0.5, size=8) ) + facet_grid(seasonGG ~ riverOrderedGG) ggsave('figures/grByYearYOY2.png')
d2$isYOY <- ifelse( d2$ageInSamples <= 3, 1, 2 ) d2 <- d2 %>% group_by(tag) %>% mutate( lagRiverN = lead(riverN), riverNMove = paste(riverN,lagRiverN, sep="_") #gr = (observedLength - lagObservedLength)/as.integer(detectionDate - lagDetectionDate) ) #d2 %>% filter(tag=='00088cbf44') %>% select(tag,riverN,lagRiverN,sampleNumber,observedLength,lagObservedLength) table(d2$riverNMove) cdMeansByisYOYriverNMove <- d2 %>% #mutate( gr = (observedLength - lagObservedLength)/as.integer(detectionDate - lagDetectionDate) ) %>% group_by(species,riverN,lagRiverN,season,isYOY) %>% summarize( meanGr = mean( grLength, na.rm = TRUE ), sdGr = sd( grLength, na.rm = TRUE ), n = n()) %>% mutate( meanGrUp = meanGr + sdGr, meanGrDn = meanGr - sdGr ) # ggplot( cdMeansByisYOYriverNMove %>% filter(!(season == 2 & isYOY == 1), species == "bkt", isYOY == 2) , aes(season, meanGr, color = factor(lagRiverN)) ) + geom_errorbar(aes(ymin = meanGrDn, ymax = meanGrUp), width = 0.25, position = pd) + geom_point( size = 2 ) + geom_line()+ facet_wrap(~riverN) d2 %>% filter(riverN == 1 & lagRiverN == 4) %>% select(tag,river,sampleNumber) d2 %>% filter(riverN == 2, is.na(lagRiverN), !is.na(grLength)) %>% select(tag,riverN,sampleNumber,lagRiverN,observedLength,lagObservedLength,gr,grLength) d2 %>% filter(tag=='00088cc059') %>% select(tag,riverOrdered,riverN,lagRiverN,sampleNumber,observedLength,lagObservedLength,gr,grLength)
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.