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.