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(arm)
library(relaimpo)
library(linkedModels)
library(arrayhelpers)
library(tidyverse)
library(gridExtra)

load("D:/projects/linkedModels/data/out/ddD_bktbntats_1997.RData")
load("D:/projects/linkedModels/data/out/ddddD_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")
ddGInAll <- rbind(ddGIn2[[1]],ddGIn2[[2]],ddGIn2[[3]]) %>%
  filter(year %in% 2000:2015, river == 'west brook') %>%
  mutate( consecOcc = ifelse( is.na(grLength), 0,1 ) )

ddGInAll %>%
  group_by(river,species) %>%
  distinct(tag) %>%
  summarize( n = n() )

#   species     n
#   <chr>   <int>
# 1 ats      5032
# 2 bkt      8556
# 3 bnt      6389

ddGInAll %>%
  group_by(river,species) %>%
  summarize( n = n() )

#   species     n
#   <chr>   <int>
# 1 ats     15791
# 2 bkt     18412
# 3 bnt     16510

ddGInAll %>%
  group_by(river,species,enc) %>%
  summarize( n = n() )

#   species   enc     n
#   <chr>   <dbl> <int>
# 1 ats         0  3515
# 2 ats         1 12276
# 3 bkt         0  3464
# 4 bkt         1 14948
# 5 bnt         0  3855
# 6 bnt         1 12655

ddGInAll %>%
  group_by(river,species,consecOcc) %>%
  summarize( n = n() ) %>%
  spread(consecOcc,n, sep="_") %>%
  mutate( prop0 = consecOcc_0/(consecOcc_0+consecOcc_1),
          prop1 = 1 - prop0)

#   river       species consecOcc_0 consecOcc_1 prop0 prop1
# 5 west brook  ats           10203        5588 0.646 0.354
# 6 west brook  bkt           13631        4781 0.740 0.260
# 7 west brook  bnt           12420        4090 0.752 0.248

length(unique(ddGInAll$sampleNumber))
ddGIn2[['bkt']] %>%
  filter(river == 'west brook') %>%
  distinct(river,year,season,nAllFishBySpeciesPBKT,nAllFishBySpeciesPBNT,nAllFishBySpeciesPATS) %>%
  arrange(year, season)
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))
allFishBySpecies$seasonGG <- factor(allFishBySpecies$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = 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 %>% filter(riverOrderedGG == "West Brook"), aes(year,nAllFishBySpecies)) +
  geom_point(aes(color = speciesGG)) +
  geom_line(aes(color = speciesGG)) +
  geom_line( data = allFish %>% filter(riverOrderedGG == "West Brook"), aes(year,nAllFish) ) +
  theme_publication() +
  labs( x = "Year", y = "Count" ) +
  facet_wrap( ~ seasonGG )#, scales = "free_y" )

ggsave('figures/3Spp_WB/counts.png')

ggplot(allFishBySpeciesYOY %>% filter(riverOrderedGG == "West Brook", year >=2000), aes(year,nAllFishBySpecies)) +
  geom_point(aes(color = speciesGG), size=3) +
#  geom_line(aes(color = speciesGG)) +
  geom_smooth(aes(color = speciesGG), se=FALSE,method='lm') +
  geom_line( data = allFish %>% filter(riverOrderedGG == "West Brook", year >=2000), aes(year,nAllFish), size=1.25, color='darkgrey', linetype = 5 ) +
  scale_y_continuous(limits = c(0,NA)) +
  theme_publication() +
  labs( x = "Year", y = "Estimated abundance" ) +
  scale_color_grey(start = 0, end = 0.9) +
  facet_wrap( ~ seasonGG , scales = "free_y" )

ggsave('figures/3Spp_WB/countsWSlope.png', width = 11, height = 9, units='in')

nSlopes <- allFishBySpeciesYOY %>%
  filter(riverOrderedGG == "West Brook", year >=2000) %>%
  group_by( species,  season ) %>%
  nest()

nModel <- function(df){
  lm( nAllFishBySpecies ~ year, data = df )
}

nSlopes <- nSlopes %>%
  mutate( model = map(nSlopes$data, nModel) )

nSlopesTidy <- nSlopes %>%
  mutate( glance = map(model, broom::tidy)) %>% 
  unnest(glance, .drop = TRUE) %>%
  filter( term == 'year')

save(nSlopesTidy,file = 'D:/projects/linkedModels/data/out/nSlopesTidy.RData')
########################
# 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" )
# sample summary table
summaryTable <- allFishBySpecies %>% 
          select(-massAllFishBySpecies) %>% 
          filter(river == 'west brook', year %in% 2000:2015) %>% 
          spread(species,nAllFishBySpecies) %>% 
          arrange(year,seasonGG)

write.csv(summaryTable, file = 'D:/projects/linkedModels/tables/rawCounts.csv')
# 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
# Has mean size changed over time?
means <- 
  cd %>% mutate( age = year - cohort ) %>%
    filter( riverOrdered == 'west brook', ageInSamples %in% 1:10, age < 3, year < 2016 ) %>%
    group_by(species,ageInSamples,year,season,age) %>%
    summarize(meanLen = mean(observedLength,na.rm = TRUE),
              sdLen = sd(observedLength,na.rm = TRUE),
              meanGR = mean(grLength,na.rm = TRUE),
              sdGR = sd(grLength,na.rm = TRUE)
              )

lenSlopes <- means %>%
  group_by( species, ageInSamples, season, age ) %>%
  nest()

lengthModel <- function(df){
  lm( meanLen ~ year, data = df )
}

grModel <- function(df){
  lm( meanGR ~ year, data = df )
}

lenSlopes <- lenSlopes %>%
  mutate( model = map(lenSlopes$data, lengthModel),
          modelGR = map(lenSlopes$data, grModel))

###################################
######### len here, gr below ######
lenSlopesTidy <- lenSlopes %>%
  mutate(glance = map(model, broom::tidy)) %>% 
  unnest(glance, .drop = TRUE) %>%
  filter(term == 'year') %>%
  mutate( sig = ifelse(p.value < 0.05,"*","")) %>%
  mutate( sig = ifelse(p.value < 0.001,"**",sig)) %>%
  mutate( sig = ifelse(p.value < 0.0001,"***",sig))

save(lenSlopesTidy,file = 'D:/projects/linkedModels/data/out/lenSlopesTidy.RData')
##
means$seasonGG <- factor(means$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
means$speciesGG <- factor(means$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)
lenSlopesTidy$seasonGG <- factor(lenSlopesTidy$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
lenSlopesTidy$speciesGG <- factor(lenSlopesTidy$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)

ggplot( means %>% filter(year >= 2000, species != 'ats'), aes(year,meanLen, group= speciesGG)) +
  geom_point( aes(color = speciesGG), size = 4) +
  geom_smooth( aes(color = speciesGG), se = FALSE, method = 'lm') +
  labs( x = "Year", y = "Mean body length (mm)" ) +
  # geom_text(data = lenSlopesTidy %>% filter(species == 'bkt'), aes(x=1999,y=220, label=paste(round(estimate,2),round(p.value,4)), color = speciesGG)) +
  # geom_text(data = lenSlopesTidy %>% filter(species == 'bnt'), aes(x=1999,y=200, label=paste(round(estimate,2),round(p.value,4)), color = speciesGG)) +
  # geom_text(data = lenSlopesTidy %>% filter(species == 'ats'), aes(x=1999,y=180, label=paste(round(estimate,2),round(p.value,4)), color = speciesGG)) +
  geom_text(data = lenSlopesTidy %>% filter(species == 'bkt', age<2), aes(x=2000,y=225, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +
  geom_text(data = lenSlopesTidy %>% filter(species == 'bnt', age<2), aes(x=2000,y=203, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +
#  geom_text(data = lenSlopesTidy %>% filter(species == 'ats', age<2), aes(x=2000,y=181, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +

    geom_text(data = lenSlopesTidy %>% filter(species == 'bkt', age==2), aes(x=2012.5,y=104, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +
  geom_text(data = lenSlopesTidy %>% filter(species == 'bnt', age==2), aes(x=2012.5,y=82, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +
#  geom_text(data = lenSlopesTidy %>% filter(species == 'ats', age==2), aes(x=2012.5,y=80, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +
  theme_publication() +
  scale_color_grey(start = 0.15, end = 0.75) +
  facet_grid(seasonGG~age)

ggsave(paste0('figures/3Spp_WB/meanLengths.png'), width = 11, height = 11, units='in')

# mean increases 
lenSlopesTidy %>%
  filter(species != 'ats') %>%
  summarise( means = mean(estimate),
             sd = sd(estimate))

# mean increases by species
lenSlopesTidy %>%
  group_by(species) %>%
  summarise( means = mean(estimate),
             sd = sd(estimate))

# do increases (slopes) increase over age? Yes
ggplot(lenSlopesTidy %>% filter(species != 'ats'), aes((ageInSamples),estimate, color=speciesGG)) +
  geom_point(aes( shape = factor(seasonGG)), size = 3) +
  geom_smooth(method='lm', se=FALSE) +
  scale_x_continuous(breaks = 1:10, labels = c(0,0,1,1,1,1,2,2,2,2)) +
  labs( x = "Age", y = "Slope of mean body length over years" ) +
  theme_publication() +
  scale_color_grey(start = 0.15, end = 0.75) 

ggsave(paste0('figures/3Spp_WB/meanLengthsOverAge.png'), width = 6, height = 6, units='in')

summary(lm(estimate~ageInSamples, data = lenSlopesTidy %>% filter(species=='bkt')))
summary(lm(estimate~ageInSamples, data = lenSlopesTidy %>% filter(species=='bnt')))


#########################
#### gr here, len above

grSlopesTidy <- lenSlopes %>%
  mutate(glance = map(modelGR, broom::tidy)) %>% 
  unnest(glance, .drop = TRUE) %>%
  filter(term == 'year') %>%
  mutate( sig = ifelse(p.value < 0.05,"*","")) %>%
  mutate( sig = ifelse(p.value < 0.001,"**",sig)) %>%
  mutate( sig = ifelse(p.value < 0.0001,"***",sig))

save(grSlopesTidy,file = 'D:/projects/linkedModels/data/out/grSlopesTidy.RData')
##
means$seasonGG <- factor(means$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
means$speciesGG <- factor(means$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)
grSlopesTidy$seasonGG <- factor(grSlopesTidy$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
grSlopesTidy$speciesGG <- factor(grSlopesTidy$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)

ggplot( means %>% filter(year >= 2000, species != 'ats'), aes(year,meanGR, group= speciesGG)) +
  geom_point( aes(color = speciesGG), size = 4) +
  geom_smooth( aes(color = speciesGG), se = FALSE, method = 'lm') +
  labs( x = "Year", y = "Mean growth rate (mm/d)" ) +
  # geom_text(data = grSlopesTidy %>% filter(species == 'bkt'), aes(x=1999,y=220, label=paste(round(estimate,2),round(p.value,4)), color = speciesGG)) +
  # geom_text(data = grSlopesTidy %>% filter(species == 'bnt'), aes(x=1999,y=200, label=paste(round(estimate,2),round(p.value,4)), color = speciesGG)) +
  # geom_text(data = grSlopesTidy %>% filter(species == 'ats'), aes(x=1999,y=180, label=paste(round(estimate,2),round(p.value,4)), color = speciesGG)) +
  geom_text(data = grSlopesTidy %>% filter(species == 'bkt', age<2), aes(x=2000,y=0.7, label=paste(round(estimate,4),sig)),hjust = 0)+#, color = speciesGG)) +
  geom_text(data = grSlopesTidy %>% filter(species == 'bnt', age<2), aes(x=2000,y=0.6, label=paste(round(estimate,4),sig)),hjust = 0)+#, color = speciesGG)) +
#  geom_text(data = grSlopesTidy %>% filter(species == 'ats', age<2), aes(x=2000,y=181, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +

    geom_text(data = grSlopesTidy %>% filter(species == 'bkt', age==2), aes(x=2011.5,y=0.66, label=paste(round(estimate,4),sig)),hjust = 0)+#, color = speciesGG)) +
  geom_text(data = grSlopesTidy %>% filter(species == 'bnt', age==2), aes(x=2011.5,y=0.56, label=paste(round(estimate,4),sig)),hjust = 0)+#, color = speciesGG)) +
#  geom_text(data = grSlopesTidy %>% filter(species == 'ats', age==2), aes(x=2012.5,y=80, label=paste(round(estimate,2),sig)),hjust = 0)+#, color = speciesGG)) +
  theme_publication() +
  scale_color_grey(start = 0.15, end = 0.75) +
  facet_grid(seasonGG~age)

ggsave(paste0('figures/3Spp_WB/meanGRs.png'), width = 11, height = 11, units='in')

# mean increases 
grSlopesTidy %>%
  filter(species != 'ats') %>%
  summarise( means = mean(estimate),
             sd = sd(estimate))

# mean increases by species
grSlopesTidy %>%
  group_by(species) %>%
  summarise( means = mean(estimate),
             sd = sd(estimate))

# do increases (slopes) increase over age? Yes
ggplot(grSlopesTidy %>% filter(species != 'ats'), aes((ageInSamples),estimate, color=speciesGG)) +
  geom_point(aes( shape = factor(seasonGG)), size = 3) +
  geom_smooth(method='lm', se=FALSE) +
  scale_x_continuous(breaks = 1:10, labels = c(0,0,1,1,1,1,2,2,2,2)) +
  labs( x = "Age", y = "Slope of mean growth rate over years" ) +
  theme_publication() +
  scale_color_grey(start = 0.15, end = 0.75) 

ggsave(paste0('figures/3Spp_WB/meanGRsOverAge.png'), width = 6, height = 6, units='in')

summary(lm(estimate~ageInSamples, data = grSlopesTidy %>% filter(species=='bkt')))
summary(lm(estimate~ageInSamples, data = grSlopesTidy %>% filter(species=='bnt')))
d2$riverOrderedGG <- factor(d2$river,
                         levels=c('west brook', 'wb jimmy', 'wb mitchell',"wb obear"),
                         labels = c("West Brook","Open Large","Open Small","Isolated Small"), ordered = T)
d2$seasonGG <- factor(d2$season, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)

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/3Spp_WB/flow.png')

lmFunctionF <- function(df){ lm( meanF ~ year, data=df ) }

flowMods <- d2 %>%
  group_by(riverOrderedGG,seasonGG,year) %>%
  summarize(meanF = mean(meanFlow, na.rm=T),
            sdF = sd(meanFlow, na.rm=T)) %>%
  group_by( riverOrderedGG, seasonGG ) %>%
  nest() %>%
  mutate( model = map(data,lmFunctionF) )

flowSlopesTidy <- flowMods %>%
  mutate(tidy = map(model,broom::tidy)) %>%
  dplyr::select(-data,-model) %>%
  unnest() %>%
  filter( term == 'year') %>%
  mutate( sig = ifelse(p.value < 0.05,"*",""))

save(flowSlopesTidy,file = 'D:/projects/linkedModels/data/out/flowSlopesTidy.RData')

#################


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) )

tempSlopesTidy <- 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)) %>%
  filter(riverOrderedGG == "West Brook") %>%
  ggplot( aes(year,meanT) ) +
  geom_smooth(method = 'lm', se=F, color = 'darkgrey') +
  geom_point( size = 3 ) +
  geom_line() +
  geom_errorbar( aes(ymin = meanT - sdT, ymax = meanT + sdT, width = 0.33 ) ) +
  labs( x = "Year", y = "Stream temperature (C)" ) +
  geom_text(data=tempSlopesTidy %>% filter(riverOrderedGG == "West Brook"),aes(x=2000.5,y=17.5, label=paste(round(estimate,2),sig))) +
  theme_publication() +
  facet_wrap(~seasonGG)

ggsave('figures/3Spp_WB/temp.png')

save(tempSlopesTidy,file = 'D:/projects/linkedModels/data/out/tempSlopesTidy.RData')
######
# for Carys Jones analysis
dJonesT <- 
  d2 %>% group_by(riverOrderedGG,seasonGG,year) %>%
    summarize(meanT = mean(meanTemperature, na.rm=T),
              sdT = sd(meanTemperature, na.rm=T)) %>%
    filter(riverOrderedGG == "West Brook", seasonGG == "Summer") %>%
    ungroup() %>%
    select(year,meanLen,sdLen)

#####
# Lengths
# run 'r mean lengths' chunk above

dJonesL <- 
  means %>%
    filter(season == 3, age == 1) %>%
    ungroup() %>%
    select(species,year,meanLen,sdLen)

######
# abundance
# run 'r adjusted counts' above

dJonesA <- 
  allFishBySpeciesYOY %>% 
    filter(riverOrderedGG == "West Brook", year >=2000, season == 3, isYOY == 2) %>%
    ungroup() %>%
    select(species,year,nAllFishBySpecies)

######
# phi
# run 'r p(detection)' below
# only has summer survival, may want annual survival

dJonesP <- phi %>%
  filter(season == 2, isYOYGG == 'Adult') %>%
  select(speciesGG,yearGG,phi) %>%
  rename(species = speciesGG, year = yearGG) %>%
  mutate(species = as.character(species))


dJonesAL <- left_join(dJonesA,dJonesL) 
dJonesALT <-  left_join(dJonesAL,dJonesT)
dJonesALTP <-  left_join(dJonesALT,dJonesP)

##
write.csv(dJonesALTP, file = 'D:/projects/linkedModels/dataForCarysJones/dJonesALTP.csv')
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)

growthYOY <- 
cd %>% group_by(speciesGG,species,riverOrderedGG,seasonGG,isYOY) %>%
  filter( riverOrdered == 'west brook', year < 2016, species != 'ats' )  %>%
  summarize(meanGR = mean(grLength, na.rm=T),
            sdGR = sd(grLength, na.rm=T)) 

ggplot(growthYOY%>%filter(isYOY == 1), aes(seasonGG,meanGR, color=speciesGG, shape=speciesGG) ) +
  geom_point( position = pd, size = 3 ) +
  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 ) +
  scale_color_grey(start = 0.0, end = 0.5) +
  theme_publication() +
  theme( panel.grid.minor = element_blank() ) 
#  facet_wrap( ~ riverOrderedGG)

ggsave('figures/3Spp_WB/meanGRYOY1.png', width = 4, height = 4, units='in')


ggplot(growthYOY%>%filter(isYOY == 2), aes(seasonGG,meanGR, color=speciesGG, shape=speciesGG) ) +
  geom_point( position = pd, size = 3 ) +
  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 ) +
  scale_color_grey(start = 0.0, end = 0.5) +
  theme_publication() +
  theme( panel.grid.minor = element_blank() ) 

ggsave('figures/3Spp_WB/meanGRYOY2.png', width = 4, height = 4, units='in')


# get ratios of growth rates for the 2 species
growthYOYWide <- growthYOY %>% ungroup() %>% select(-sdGR,-speciesGG) %>% spread(key=species,value=meanGR) %>%
  mutate(ratio = bnt/bkt)

mean(growthYOYWide$ratio)
# 1.107745
cd %>% group_by(species) %>% filter(riverOrderedGG == "West Brook",species != 'ats') %>%
  summarize(meanGR = mean(grLength, na.rm=T),
            sdGR = sd(grLength, na.rm=T))
.197/.182#bnt/bkt

cd %>% group_by(isYOY) %>% filter(riverOrderedGG == "West Brook",species != 'ats') %>%
  summarize(meanGR = mean(grLength, na.rm=T),
            sdGR = sd(grLength, na.rm=T))
.145/.244 #adult/yoy
.244/.145


cd %>% group_by(season) %>% filter(riverOrderedGG == "West Brook",species != 'ats') %>%
  summarize(meanGR = mean(grLength, na.rm=T),
            sdGR = sd(grLength, na.rm=T))
.434/((.114+.124+.067)/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) %>% filter(river == 1, year > 3) 

  p$speciesGG <- factor(p$species, levels = 1:3, labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)
  p$seasonGG <- factor(p$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
  p$isYOYGG <- factor(p$isYOY, levels = 1:2, labels = c("Juvenile","Adult"), ordered = T)
  p$yearGG <- p$year + 1996

  ggplot(p, aes(yearGG,invlogit(est),color=factor(speciesGG))) + 
    geom_point(size = 2) + 
    labs( x = 'Year', y = 'Probability of detection') +
    scale_color_grey(start = 0.0, end = 0.85) +
    theme_publication() +
    facet_grid( seasonGG ~ isYOYGG)

  # seem to need to run the ggsave outside of markdown by pasting into console
   ggsave('figures/3Spp_WB/p.png', height = 11, width = 9, units = 'in')

  p %>% group_by(species) %>%
    summarize(meanP = mean(invlogit(est)),
              sdP = sd(invlogit(est)))

  p %>% group_by(season) %>%
    summarize(meanP = mean(invlogit(est)))

  ##################### 
  phi <- array2df(ddD$q50$phiBetaInt, label.x = "est") %>%
    rename(est=est,isYOY=d1,species=d2,season=d3,river=d4,year=d5) %>% filter(river == 1, year > 3)  

  phi$speciesGG <- factor(phi$species, levels = 1:3, labels = c("bkt", "bnt", "ats"), ordered = T)
  phi$seasonGG <- factor(phi$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
  phi$isYOYGG <- factor(phi$isYOY, levels = 1:2, labels = c("Juvenile","Adult"), ordered = T)
  phi$yearGG <- phi$year + 1996
  phi$phi <- invlogit(phi$est)

  ggplot(phi, aes(yearGG,phi,color=factor(speciesGG))) +geom_point() + facet_grid(seasonGG ~ isYOYGG)

  ggsave('figures/3Spp_WB/phi.png')

  phi %>% group_by(species) %>%
    summarize(meanPhi = mean(invlogit(est)),
              sdPhi = sd(invlogit(est)))

  phi %>% group_by(season) %>%
    summarize(meanPhi = mean(invlogit(est)))

# # 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)
rBkt <- unlist(mcmc[['bkt']]$Rhat)
sort(rBkt[rBkt > 1.1])

rBnt <- unlist(mcmc[['bnt']]$Rhat)
sort(rBnt[rBnt > 1.1])
mcmc[['bnt']]$Rhat$grBeta[,,,1]
mcmc[['bnt']]$Rhat$grInt
mcmc[['bnt']]$Rhat$grIntMu

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")

hist(rHat[["bkt"]]$rHat)
hist(rHat[["bnt"]]$rHat)
rHat[["bnt"]] %>% filter( rHat >1.1 )

ggplot(rHat[["bkt"]], 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") %>% mutate(species = 'bkt')
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") %>%   mutate(species = 'bnt')
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") %>% mutate(species = 'ats')


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)

#grIntDF <- rbind(grInt[['bkt']],grInt[['bnt']],grInt[['ats']])
#tmp <- grIntDF %>% filter(river == 'west brook')
#save(tmp , file= "data/out/grIntDFAllIters.RData") #for Ron

##################################################
# as range plot

plotHPD <- function(sppToPlot){

  nRivers <- list()
  nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1

  probs <- c(0.025,0.1,0.5,0.9,0.975)
  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$seasonGG <- factor(p$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
  p$parameter <- p$seasonGG

  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)

  p$seasonGG <- factor(p$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
  #print(c(min(p$intByDayDown,na.rm=T),max(p$intByDayDown,na.rm=T)))

  #save(p, file= paste0("data/out/allInts_",sppToPlot,".RData"))  #for projections and tables

    gg <- ggplot(p %>% filter(river=='west brook'), aes(x = intByDay, y = reorder(parameter,intByDay,na.rm=TRUE), color = factor(isYOY), shape = factor(isYOY))) + 
        theme_publication() +  
        geom_point(size = 4) +
        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.00,0.55)) +
        xlab("Mean daily growth (mm/d)") + ylab("Season") +
        scale_color_grey(start = 0, end = 0.65) 
       # ggtitle(sppToPlot) 

    return(gg)
}

plotHPD('bkt'); ggsave('figures/3Spp_WB/HDP_bkt.png', height = 5, width = 5, scale = 1)
plotHPD('bnt'); ggsave('figures/3Spp_WB/HDP_bnt.png', height = 5, width = 5, scale = 1)
plotHPD('ats'); ggsave('figures/3Spp_WB/HDP_ats.png', height = 5, width = 5, scale = 1)
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.025,0.1,0.5,0.9,0.975)
  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$seasonGG <- factor(p$season, levels = 1:4, labels = c("Spring","Summer","Autumn","Winter"), ordered = T)
  p$parameter <- p$seasonGG


  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 %>% filter(river == 'west brook'), aes(x = exp(int), y = reorder(parameter,int,na.rm=TRUE), color = factor(isYOY), 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("Season") +
        scale_color_grey(start = 0, end = 0.65) 
       # ggtitle(sppToPlot) 

    return(gg)
}

plotHPDSigma('bkt'); ggsave('figures/3Spp_WB/HDPSigma_bkt.png', height = 5, width = 5, scale = 1)
plotHPDSigma('bnt'); ggsave('figures/3Spp_WB/HDPSigma_bnt.png', height = 5, width = 5, scale = 1)
plotHPDSigma('ats'); ggsave('figures/3Spp_WB/HDPsigma_ats.png', height = 5, width = 5, scale = 1)
#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)

#tmp <- grBetaDF %>% filter(river == 'west brook')
#save(tmp , file= "data/out/grBetaDFAllIters.RData") #for Ron
##################################################
# as range plot

  nRivers <- list()
  nRivers[['bkt']] <- 4; nRivers[['bnt']] <- 2; nRivers[['ats']] <- 1
  speciesList <- list('bkt','bnt','ats')

  probs <- c(0.025,0.1,0.5,0.9,0.975)
  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.025, hi = X0.975 )

  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 = q$speciesGG
q$parameterOrd <- factor(q$parameter, 
                         levels = c( "ATS", "BNT", "BKT"), 
                         labels = c( "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")) )

  save(q, file= "data/out/allBetas.RData")  #for projections and tables

  #reorder(parameter,median,na.rm=TRUE)
(min(q$median))
  (max(q$median))
#  q[q$median < -0.13,]

    ggplot(q %>% filter(river == 'west brook', species != 'ats'), aes(x = median, y = parameterOrd, color = factor(species), shape = factor(isYOY))) + 
        geom_vline(xintercept = 0,color='grey') +
        theme_publication() +  
        geom_point(size = 4) +
        geom_segment(aes(x = median, xend = hi, yend = parameterOrd), 
        size = 0.5, arrow = arrow(angle = 90, length = unit(0.1, "inches"))) + 
        geom_segment(aes(x = median, xend = lo, yend = parameterOrd), 
        size = 0.5,arrow = arrow(angle = 90, length = unit(0.1, "inches"))) +
     #   xlim(c(-0.15,0.11)) +
        scale_x_continuous( breaks=c(-0.1,0,0.1)) +
        xlab("Parameter estimate") + theme(axis.text.y = element_blank(), axis.title.y = element_blank()) +
        scale_color_grey(start = 0, end = 0.5)  +
        facet_grid(seasonGG~betaGG) 

ggsave('figures/3Spp_WB/grBetas.png', height = 8, width = 10)

save(betas,file = 'D:/projects/linkedModels/data/out/nSlopesTidy.RData')
load(file = "data/out/allInts_bkt.RData") #loads p
pBKT <- p %>% mutate(species = 'bkt')

load(file = "data/out/allInts_bnt.RData") #loads p
pBNT <- p %>% mutate(species = 'bnt')

pp <- rbind(pBKT,pBNT)

intsForTable <- pp %>% 
  filter(river == 'west brook', species %in% c('bkt','bnt')) %>%
  arrange(species,isYOY,seasonGG) %>%
  mutate(lo = round(intDown,3),mean = round(int,3),hi = round(intUp,3)) %>%
  select(species,isYOY,seasonGG,intDown,int,intUp)

write.csv(intsForTable,file = 'D:/projects/linkedModels/data/out/intsForTableAgeGroup.csv')
load(file= "data/out/allBetas.RData") #loads q

betasForTable <- q %>% 
  filter(river == 'west brook', species %in% c('bkt','bnt')) %>%
  arrange(species,betaGG,isYOY,seasonGG,betaGG) %>%
  mutate(lo = round(lo,3),median = round(median,3),hi = round(hi,3)) %>%
  select(species,betaGG,isYOY,seasonGG,betaGG,lo,median,hi)

write.csv(betasForTable,file = 'D:/projects/linkedModels/data/out/betasForTableAgeGroup.csv')
# 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')

  probs <- c(0.025,0.1,0.5,0.9,0.975)
  nSpecies <- 2


  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.025, hi = X0.975 )

  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 = q$speciesGG
q$parameterOrd <- factor(q$parameter, 
                         levels = c( "ATS", "BNT", "BKT"), 
                         labels = c( "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 %>% filter(river == 'west brook', species != 'ats'), 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") + theme(axis.text.y = element_blank(), axis.title.y = element_blank()) +
        scale_color_grey(start = 0, end = 0.5)  +
        facet_grid(seasonGG~betaGG) 

ggsave('figures/3Spp_WB/grBetaSigmas.png', height = 8, width = 8)
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 %>% filter(river == 'west brook'), "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/3Spp_WBtmp.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 %>% filter(species != 'ats'), 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]) +
      scale_color_grey(start = 0, end = 0.5)  +
      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)")
# save from rStudio as export, width=900, height=500
# ggsave(gridPred,'figures/3Spp_WB/predGRYOY1.png', width = 9, height = 7, units='in')
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 <- 1
riverGGToPlot <- "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, species != 'ats')
  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]) +
      scale_color_grey(start = 0, end = 0.5)  +
      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)")

# 0 = 875,600
# 1 = 875,500

#ggsave(ggCV, paste0('figures/3Spp_WBpredGRCV_YOY',isYOYGG,'_',riverGGToPlot,'.png'), width = 9, height = 7, units='in')
###############################
# flow*temp interaction plots
isYOYGG <- 0
riverGGToPlot <- "West Brook" #"Open Large" #
outputToPlot <- "meanPredGRSigma" #"meanPredGR" #

interact <- pred %>% filter(isYOY == isYOYGG,
                            riverGG == riverGGToPlot,
                            len == 0, cBKT == 0, cBNT == 0, cATS == 0 ) %>%
                     mutate( temp_spp = paste0(temp,speciesGG))

    ggplot( interact %>% filter(species != 'ats'), aes( flow, eval(as.name(outputToPlot)), group=temp_spp, linetype=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('Stream flow') +
      scale_y_continuous("Standard deviation of growth rate (mm/d)") +
      #scale_color_manual(values=c('blue','darkgreen'))  +
      facet_grid(~seasonGG)


#yoy0 = 875,350
#yoy1 = 875,280
###############################
# flow*temp interaction plots
isYOYGG <- 1
riverGGToPlot <- "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 %>% filter(species != 'ats'), aes( cBKT, eval(as.name(outputToPlot)), group=bnt_spp, linetype=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') +
      scale_y_continuous("Growth rate (mm/d)") +
      facet_grid(~seasonGG)

# for saving from export
#yoy0 = 875,350
#yoy1 = 875,275
# run 'Load data' chunk at top to get the 'mcmc' object

#data for running the simulation
simData <- list(
  nReps = 100, #num of mcmc iterations to grab
  numInd = 500 , # number of individual to simulate
  nOccYOY = 3,
  nOccAdult = 7,
  riverIndex = 1, #Levels: west brook < wb jimmy < wb mitchell < wb obear
  minSD = -1.5, #range and step for independent variable values in the sim
  maxSD = 1.5, 
  step = 1.5
)
simData$nOcc <- simData$nOccYOY + simData$nOccAdult
simData$indVarSeq = seq(simData$minSD,simData$maxSD,simData$step)

stats3Spp <- grStats3Spp <- anov <- elapsed <- lm <- relImp <- list()
spp <- c('bkt','bnt') #,'ats')

for( speciesForProj in spp ){

  # get mean on mean length by cohort for the fall yoy sample
  obsLenStats <-  
    ddGIn2[[speciesForProj]] %>%
      group_by(cohort,ageInSamples) %>%
      summarise(lenMean = mean(observedLength, na.rm=TRUE),
                lenSD = sd(observedLength, na.rm=TRUE)) %>%
      filter(ageInSamples == 1) %>%
      ungroup() %>%
      summarise(lenMeanMean = mean(lenMean, na.rm=TRUE),
                lenSDMean = mean(lenSD, na.rm=TRUE)) 


  simData$species = speciesForProj
  simData$meanInit = obsLenStats$lenMeanMean
  simData$sdInit = obsLenStats$lenSDMean
  simData$nSamplesAvail = mcmcInfo$nChains * mcmcInfo$nSamples

  simData <- updateSimData(simData) # add #season, itersForSim etc data to simData
  nB <- getNBetas() # get number of betas from mcmc[]
  m <- createMatrices()

  indVars <- getIndVars() # get independent variables table, -1.5, 0, 1.5

  startTime <- Sys.time()
  sizeStats <- sizeStats2 <- grStats <- grStats2 <- NULL
  #loop over conditions
  for(testRow in indVars$rowNum){
    project <- projectSizes(m,simData,testRow)

    sizeMat <- project$sizeMat
    grMat <- project$grMat

    sizeStats1 <- summarizeLens(sizeMat,testRow)
    sizeStats <- rbind(sizeStats,sizeStats1[[1]]) #means over individuals (grouped by iter, occ)
    sizeStats2 <- rbind(sizeStats2,sizeStats1[[2]]) #means by occ of means over individuals (grouped by occ)

    grStats1 <- summarizeGRs(grMat,testRow)
    grStats <- rbind(grStats,grStats1[[1]]) #means over individuals (grouped by iter, occ)
    grStats2 <- rbind(grStats2,grStats1[[2]]) #means by occ of means over individuals (grouped by occ)
  }
  endTime <- Sys.time()
  (elapsed[[speciesForProj]] <- endTime - startTime)

  stats3Spp[[speciesForProj]] <- list( sizeStats,sizeStats2 )
  grStats3Spp[[speciesForProj]] <- list( grStats,grStats2 )

  # lm[[speciesForProj]] <- lm( lenMean ~ temp + flow + bkt + bnt + ats + temp * flow + temp^2 + flow^2 + bnt*bkt, data = sizeStats %>% filter(occ == (simData$nOcc)))
  # 
  # relImp[[speciesForProj]] <- calc.relimp(lm[[speciesForProj]], type = c("lmg"), rela = TRUE)
  #   
  # anov[[speciesForProj]] <- anova( lm[[speciesForProj]] )

}

#projOut <- list(stats3Spp, lm, relImp, anov, elapsed, simData)
save( stats3Spp, grStats3Spp, elapsed, simData, file = paste0('D:/projects/linkedModels/data/out/projectionOut_N',simData$nReps,'.RData'))

# if not re-running the sim
# load( file = 'D:/projects/linkedModels/data/out/projectionOut.RData')
# load( file = paste0('D:/projects/linkedModels/data/out/projectionOut_N',simData$nReps,'.RData'))
# spp <- c('bkt','bnt','ats')

# raw sim data figs - size
ggplot(stats3Spp[['bkt']][[1]] %>% filter(ats == 0, bkt == 0, bnt == 0), aes(as.numeric(occ),lenMean)) +
  geom_jitter(aes(color = factor(ats))) +
  #geom_smooth(aes(color = factor(temp), linetype = factor(flow)))+#(bnt + 2)/4 )) +
  geom_smooth()+#(bnt + 2)/4 )) +
  geom_smooth(color = 'green', data = stats3Spp[['bnt']][[1]] %>% filter(ats == 0, bkt == 0, bnt == 0)) +
  facet_grid(temp ~ flow)


########################
# raw sim data figs - gr
ggplot(grStats3Spp[['bkt']][[1]] %>% filter(ats == 0, bkt == 0, bnt == 0), aes(as.numeric(occ),grMean)) +
  geom_jitter(aes(color = factor(ats)), width = 0.2, size = 0.5) +
  #geom_smooth(aes(color = factor(temp), linetype = factor(flow)))+#(bnt + 2)/4 )) +
  #geom_smooth()+#(bnt + 2)/4 )) +
  geom_jitter(color = 'darkblue', width = 0.1, size = 0.5, data = grStats3Spp[['bnt']][[1]] %>% filter(ats == 0, bkt == 0, bnt == 0)) +
  facet_grid(temp ~ flow)
# relative proportion stats with occ in the lm 
grStats3Spp[['bkt']][[1]]$season <- as.factor(c(3,4,1,2,3,4,1,2,3,4))
grStats3Spp[['bkt']][[1]]$temp2 <- grStats3Spp[['bkt']][[1]]$temp^2
grStats3Spp[['bkt']][[1]]$flow2 <- grStats3Spp[['bkt']][[1]]$flow^2


indepVarVals <- c(-1.5,0,1.5) # make sure symetrical around 0
indepVars <- c('temp', 'flow', 'temp:flow', 'bkt', 'bnt', 'ats', 'bkt:bnt')
rIBKTBNTAll <- NULL
# loop over occs
for(rep in 1:simData$nReps){

  print(rep)

  datOccBKTAll <- grStats3Spp[['bkt']][[1]] %>% filter(iter == rep, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals)
  lmAllBKT <- lm( grMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt + occ, data = datOccBKTAll) 
  relImpAllBKT <- calc.relimp(lmAllBKT, type = c("lmg"), rela = TRUE)

  datOccBNTAll <- grStats3Spp[['bnt']][[1]] %>% filter(iter == rep, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals)
  lmAllBNT <- lm( grMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt + occ, data = datOccBNTAll)
  relImpAllBNT <- calc.relimp(lmAllBNT, type = c("lmg"), rela = TRUE)

  lmgBKT <- relImpAllBKT$lmg[(names(relImpAllBKT$lmg) %in% indepVars)]
  rBKT <- lmgBKT / sum(lmgBKT) 

  lmgBNT <- relImpAllBNT$lmg[(names(relImpAllBNT$lmg) %in% indepVars)]
  rBNT <- lmgBNT / sum(lmgBNT) 

  rIBKTBNTAll <- rbind(rIBKTBNTAll, data.frame(cbind(rBKT,rBNT)) %>% 
                      rownames_to_column() %>% 
                      gather(key = rowname) %>% 
                      mutate(var = rep(1:length(indepVars),2),
                             varText = rep(names(rBKT),2),
                             iter = rep)
                    )
}

rIBKTBNTAllMeans <- rIBKTBNTAll %>%
  group_by(rowname, var) %>%
  summarise(rIMean = mean(value),
            rISd = sd(value),
            rIUpperCI = quantile(value, 0.95))

dodge <- position_dodge(width=0.9)
ggplot(rIBKTBNTAllMeans, aes(factor(var),rIMean, fill = factor(rowname, levels = c('rBKT','rBNT'), ordered = TRUE))) +
  geom_errorbar(aes(ymin = rIMean - 0.001, ymax = rIUpperCI), width = 0.25, position=dodge) +
  geom_col(position=dodge) +
  scale_x_discrete('Parameter',
                   labels = c('1' = 'Temperature','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")) +
  scale_fill_grey() +
  scale_y_continuous('Relative proportion of variance explained') +
  theme_publication() +
  geom_hline(yintercept = 0) +
  theme( legend.position = 'none')

ggsave(paste0('figures/3Spp_WB/relImpOccMeanModAll_GR_N',simData$nReps,'.png'), width = 9, height = 7, units='in')

# for MS
rIBKTBNTAllMeans
#    rowname   var  rIMean    rISd rIUpperCI
#    <chr>   <int>   <dbl>   <dbl>     <dbl>
#  1 rBKT        1 0.132   0.0251     0.170 
#  2 rBKT        2 0.670   0.0369     0.731 
#  3 rBKT        3 0.0116  0.0118     0.0330
#  4 rBKT        4 0.0346  0.0256     0.0703
#  5 rBKT        5 0.0951  0.0284     0.140 
#  6 rBKT        6 0.0441  0.0263     0.0919
#  7 rBKT        7 0.0121  0.0155     0.0410
#  8 rBNT        1 0.106   0.0212     0.141 
#  9 rBNT        2 0.568   0.0364     0.623 
# 10 rBNT        3 0.0885  0.0223     0.124 
# 11 rBNT        4 0.00518 0.00556    0.0163
# 12 rBNT        5 0.122   0.0285     0.169 
# 13 rBNT        6 0.0314  0.0176     0.0624
# 14 rBNT        7 0.0782  0.0270     0.122 

rIBKTBNTAll %>%
  mutate(isEnvVar = ifelse(varText %in% c('temp','flow','temp:flow'), 1,0)) %>%
  group_by(rowname,isEnvVar) %>%
  summarize(meanMean = mean(value),
            sumMean = sum(value),
            sdMean = sd(value))

#   rowname isEnvVar meanMean sumMean sdMean
#   <chr>      <dbl>    <dbl>   <dbl>  <dbl>
# 1 rBKT           0   0.0383    15.3 0.0402
# 2 rBKT           1   0.282     84.7 0.279 
# 3 rBNT           0   0.0736    29.4 0.0484
# 4 rBNT           1   0.235     70.6 0.239 

#bkt
.282/(.282+.038) # 0.88125
#bnt
.235/(.235+.074) # 0.7605178
.881/.76 # 1.159211
# relative proportion stats with occ in the lm 
grStats3Spp[['bkt']][[1]]$season <- as.factor(c(3,4,1,2,3,4,1,2,3,4))
grStats3Spp[['bkt']][[1]]$temp2 <- grStats3Spp[['bkt']][[1]]$temp^2
grStats3Spp[['bkt']][[1]]$flow2 <- grStats3Spp[['bkt']][[1]]$flow^2


ageForRI <-  'yoy' #'yoy' or 'adult'
occs <- if (ageForRI == 'yoy') 1:3 else 4:10
indepVarVals <- c(-1.5,0,1.5) # make sure symetrical around 0
indepVars <- c('temp', 'flow', 'temp:flow', 'bkt', 'bnt', 'ats', 'bkt:bnt')
rIBKTBNT <- NULL
# loop over occs
for(rep in 1:simData$nReps){

  print(rep)

  datOccBKT <- grStats3Spp[['bkt']][[1]] %>% filter(iter == rep, occ %in% occs, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals)
  lmAllBKT <- lm( grMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt + occ, data = datOccBKT) 
  relImpAllBKT <- calc.relimp(lmAllBKT, type = c("lmg"), rela = TRUE)

  datOccBNT <- grStats3Spp[['bnt']][[1]] %>% filter(iter == rep, occ %in% occs, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals)
  lmAllBNT <- lm( grMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt + occ, data = datOccBNT)
  relImpAllBNT <- calc.relimp(lmAllBNT, type = c("lmg"), rela = TRUE)

  lmgBKT <- relImpAllBKT$lmg[(names(relImpAllBKT$lmg) %in% indepVars)]
  rBKT <- lmgBKT / sum(lmgBKT) 

  lmgBNT <- relImpAllBNT$lmg[(names(relImpAllBNT$lmg) %in% indepVars)]
  rBNT <- lmgBNT / sum(lmgBNT) 

  rIBKTBNT <- rbind(rIBKTBNT, data.frame(cbind(rBKT,rBNT)) %>% 
                      rownames_to_column() %>% 
                      gather(key = rowname) %>% 
                      mutate(var = rep(1:length(indepVars),2),
                             varText = rep(names(rBKT),2),
                             iter = rep)
                    )
}

rIBKTBNTMeans <- rIBKTBNT %>%
  group_by(rowname, var) %>%
  summarise(rIMean = mean(value),
            rISd = sd(value),
            rIUpperCI = quantile(value, 0.95))

dodge <- position_dodge(width=0.9)
ggplot(rIBKTBNTMeans, aes(factor(var),rIMean, fill = factor(rowname, levels = c('rBKT','rBNT'), ordered = TRUE))) +
  geom_errorbar(aes(ymin = rIMean - 0.001, ymax = rIUpperCI), width = 0.25, position=dodge) +
  geom_col(position=dodge) +
  scale_x_discrete('Parameter',
                   labels = c('1' = 'Temperature','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")) +
  scale_fill_grey() +
  scale_y_continuous('Relative proportion of variance explained') +
  theme_publication() +
  geom_hline(yintercept = 0) +
  theme( legend.position = 'none')

ggsave(paste0('figures/3Spp_WB/relImpOccMeanModWOcc_GR',simData$nReps,'_',ageForRI,'.png'), width = 9, height = 7, units='in')

#    rowname   var  rIMean    rISd rIUpperCI
#    <chr>   <int>   <dbl>   <dbl>     <dbl>
#  1 rBKT        1 0.516   0.0639     0.619 
#  2 rBKT        2 0.123   0.0332     0.175 
#  3 rBKT        3 0.0209  0.0222     0.0568
#  4 rBKT        4 0.0579  0.0307     0.100 
#  5 rBKT        5 0.193   0.0487     0.273 
#  6 rBKT        6 0.0539  0.0374     0.126 
#  7 rBKT        7 0.0362  0.0303     0.0889
#  8 rBNT        1 0.154   0.0566     0.253 
#  9 rBNT        2 0.0106  0.00999    0.0281
# 10 rBNT        3 0.140   0.0479     0.208 
# 11 rBNT        4 0.00421 0.00467    0.0135
# 12 rBNT        5 0.463   0.0780     0.598 
# 13 rBNT        6 0.173   0.0676     0.294 
# 14 rBNT        7 0.0547  0.0340     0.121


# get isEnvVar sens for yoy and adults - need to rerun the above code for each age
rIBKTBNT %>%
  mutate(isEnvVar = ifelse(varText %in% c('temp','flow','temp:flow'), 1,0)) %>%
  group_by(rowname,isEnvVar) %>%
  summarize(meanMean = mean(value),
            sumMean = sum(value),
            sdMean = sd(value)) %>%
  mutate(age = ageForRI)

#   rowname isEnvVar meanMean sumMean sdMean age  
#   <chr>      <dbl>    <dbl>   <dbl>  <dbl> <chr>
# 1 rBKT           0   0.0769    30.8 0.0763 yoy  
# 2 rBKT           1   0.231     69.2 0.209  yoy  
# 3 rBNT           0   0.165     66.2 0.185  yoy  
# 4 rBNT           1   0.113     33.8 0.0888 yoy

#   rowname isEnvVar meanMean sumMean sdMean age  
#   <chr>      <dbl>    <dbl>   <dbl>  <dbl> <chr>
# 1 rBKT           0   0.0287    11.5 0.0272 adult
# 2 rBKT           1   0.295     88.5 0.377  adult
# 3 rBNT           0   0.0411    16.4 0.0285 adult
# 4 rBNT           1   0.279     83.6 0.343  adult




# ratio of env/(env+spp) means - yoy
.231/(.231 + .0769) # bkt 0.7502436
.165/(.165 + .113) # bnt 0.5935252

# ratio of env/(env+spp) means - adult
.295/(.295 + .0287) # bkt 0.9113377
.279/(.279 + .0411) # bnt 0.8716026

.75/.593 # 1.26
.911/.872 # 1.04
##############################################################
############################################################

# projOut <- list(stats3Spp, lm, relImp, anov, elapsed)
load(paste0('D:/projects/linkedModels/data/out/projectionOut_N',simData$nReps,'.RData'))

#spp <- c('bkt','bnt','ats')
# scenarios
# Plot subsets of the simulations
#means for each iter
stats1 <- rbind(data.frame(grStats3Spp[['bkt']][[1]]), 
                data.frame(grStats3Spp[['bnt']][[1]]), 
                data.frame(grStats3Spp[['ats']][[1]])) %>%
          mutate(spp = rep(spp, each = nrow(grStats3Spp[['bkt']][[1]])))

#means for each occ
stats2 <- rbind(data.frame(grStats3Spp[['bkt']][[2]]), 
                data.frame(grStats3Spp[['bnt']][[2]]), 
                data.frame(grStats3Spp[['ats']][[2]])) %>%
          mutate(spp = rep(spp, each = nrow(grStats3Spp[['bkt']][[2]])))

ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, ats == 0), aes(occ,grMeanMean, group = spp)) +
  geom_point(aes(color = spp) ) +
  geom_line( aes(color = spp) ) +
  geom_point(aes(color = spp), data = stats2 %>% filter(temp == 1.5, flow == 0, bkt == 0, bnt == 0, ats == 0) ) +
  geom_line( aes(color = spp), data = stats2 %>% filter(temp == 1.5, flow == 0, bkt == 0, bnt == 0, ats == 0), linetype = 2 ) 


ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, ats == 0), aes(occ,grMeanMean, group = spp)) +
  geom_point(aes(color = spp) ) +
  geom_line( aes(color = spp) ) +
  geom_point(aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 1.5, bkt == 0, bnt == 0, ats == 0) ) +
  geom_line( aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 1.5, bkt == 0, bnt == 0, ats == 0), linetype = 2 ) 

ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, ats == 0), aes(occ,grMeanMean, group = spp)) +
  geom_point(aes(color = spp) ) +
  geom_line( aes(color = spp) ) +
  geom_point(aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 0, bkt == 1.5, bnt == 0, ats == 0) ) +
  geom_line( aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 0, bkt == 1.5, bnt == 0, ats == 0), linetype = 2 ) 


# trajectories across indep var values for different combinations. Trying to relate to relImp results
#temp
p1 <- ggplot( stats2 %>% filter(flow == 0, bkt == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(temp))) +
  geom_line(aes(color = factor(temp))) +
  geom_point(aes(color = factor(temp))) +
  facet_wrap(~spp) 

#flow
p2 <- ggplot( stats2 %>% filter(temp == 0, bkt == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(flow))) +
  geom_line(aes(color = factor(flow))) +
  geom_point(aes(color = factor(flow))) +
  facet_wrap(~spp)  

#bkt
p3 <- ggplot( stats2 %>% filter(temp == 0, flow == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(bkt))) +
  geom_line(aes(color = factor(bkt))) +
  geom_point(aes(color = factor(bkt))) +
  facet_wrap(~spp) 

#bnt
p4 <- ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(bnt))) +
  geom_line(aes(color = factor(bnt))) +
  geom_point(aes(color = factor(bnt))) +
  facet_wrap(~spp)

#ats
p5 <- ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(ats))) +
  geom_line(aes(color = factor(ats))) +
  geom_point(aes(color = factor(ats))) +
  facet_wrap(~spp)

multiplot(p1, p2, p3, p4, p5, cols=2) # multiplot function is below

#temp*flow
ggplot( stats2 %>% filter(bkt == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(flow))) +
  geom_line(aes(color = factor(flow))) +
  geom_point(aes(color = factor(flow))) +
  facet_grid(temp~spp)

#bkt*bnt
ggplot( stats2 %>% filter(temp == 0, flow == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,grMeanMean, group = factor(bnt))) +
  geom_line(aes(color = factor(bnt))) +
  geom_point(aes(color = factor(bnt))) +
  facet_grid(bkt~spp)




# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}



# size=100
# gr=10
# sd=1
# mean(rnorm(10000,size+gr,sd))
# mean(size+rnorm(10000,gr,sd))
# make figures of relative importance for each sampling occasion
rIAll <- cAll <- NULL

for(ageGroupLoop in 0:1){
  for(seasonLoop in 1:4){
#for( occSelected in 1:(simData$nOcc - 1)){
  if(ageGroupLoop == 0 & seasonLoop == 2) next  

  lmAll <- ggRIAll <- relImpAll <- coefsAll <- list()

  for(spp1 in spp){
    print(c(ageGroupLoop, seasonLoop, spp1))  

    datAgeGroupSeason <- grStats3Spp[[spp1]][[1]] %>%
         mutate(ageGroup = ifelse(occ %in% 1:3, 0, 1),
         season = rep(c(3,4,1,2,3,4,1,2,3,4), nrow(grStats3Spp[[spp1]][[1]])/simData$nOcc))

      lmAll[[spp1]] <- lm( grMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt, data = datAgeGroupSeason %>% filter(ageGroup == ageGroupLoop, season == seasonLoop, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals))

      relImpAll[[spp1]] <- calc.relimp(lmAll[[spp1]], type = c("lmg"), rela = TRUE)

      coefsAll[[spp1]] <- lmAll[[spp1]]$coefficients
  }

  # make bar chart of relative importance of each variable for each spp
  rI <- cbind(data.frame(bktRI = relImpAll[['bkt']]$lmg), 
              data.frame(bntRI = relImpAll[['bnt']]$lmg)) %>%
        rownames_to_column() %>% 
        gather(key=rowname) %>% 
        mutate(var = rep(1:length(relImpAll[['bkt']]$lmg),length(spp)),
               ageGroup = ageGroupLoop,
               season = seasonLoop
               )
  rI$varText <- rep(names(relImpAll[['bkt']]$lmg),2)

  rIAll <- rbind(rIAll,rI)

  #rIAll$age <- ifelse(rIAll$occ %in% 1:2, 0, ifelse(rIAll$occ %in% 3:6, 1, 2))

  }  
}

######################################################
# relImp over occasions for each spp and parameter
labels <- c('1' = 'Temp','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")

# ggplot(rIAll %>% filter(rowname %in% c('bktRI','bntRI')), aes((occ), value, color = factor(rowname))) +
#   geom_point() +
#   geom_line() +
#   geom_vline(xintercept = 2.5) +
#   scale_x_continuous('Sampling occasion',
#                     breaks = 1:9,
#                     labels = c('Autumn','Winter','Spring','Summer', 'Autumn','Winter','Spring','Summer', 'Autumn')
#                     ) +
#   scale_y_continuous('Relative proportion of variance explained') +
#   theme_publication() +
#   theme(axis.text.x  = element_text(angle=90, vjust=0.25, size=10)) +
#   facet_wrap(~ var, labeller=labeller(var = labels))
# 
# ggsave(paste0('figures/3Spp_WB/relImpOccAll_GR.png'), width = 11, height = 9, units='in')


rIAllMinus <- rIAll %>% mutate(value2 = ifelse(rowname == 'bktRI',value * -1, value))

labelsSeason <- c('1'='Autumn','2'='Winter','3'='Spring','4'='Summer','5'='Autumn','6'='Winter','7'='Spring','8'='Summer','9'='Autumn')
labelsSeason2 <- c('1'='Spring','2'='Summer','3'='Autumn','4'='Winter')
labelAge <- c('0' = 'Growth year 0/1', '1' = 'Growth Year 1+')

### get prop env
summ1 <- rIAll %>%
  mutate(isEnvVar = ifelse(varText %in% c('temp','flow','temp:flow'), 1,0)) %>%
  group_by(rowname,isEnvVar,ageGroup, season) %>%
  summarize(mean=mean(value)) %>%
  arrange(ageGroup,season,rowname,isEnvVar)

summ2 <- summ1 %>%
  group_by(rowname,ageGroup,season) %>%
  summarize(sum=sum(mean))

summ1a <- left_join(summ1,summ2) %>%
  mutate(prop = mean/sum,
         xAxisLoc = ifelse(rowname == 'bktRI',-0.85, 0.85)) %>%
  filter(isEnvVar == 1) 
###

####
# add in betas to rIAllMinus
load("D:/projects/linkedModels/data/out/allBetas.RData") #loads df 'q'
qFig1 <- q %>% filter(riverGG == 'WB', species %in% c('bkt','bnt'))
qFig1$season <- as.numeric(as.character(qFig1$season))
qFig1$isYOY <- as.numeric(as.character(qFig1$isYOY))
qFig <- qFig1 %>%
  mutate( varText = tolower(betaGG)) %>%
  mutate( varText = ifelse(varText == 'temp*flow','temp:flow', varText)) %>%
  mutate( varText = ifelse(varText == 'bnt*bkt','bkt:bnt', varText)) %>%
  select(isYOY,season,species,median,varText)

# prepare rIAllMinus
rIAllMinusQ <- rIAllMinus %>%
  mutate(species = str_sub(rowname,1,3),
         isYOY = ageGroup) %>%
  left_join(qFig) %>%
  mutate(median01 = ifelse(median > 0, 1, 0),
         medianBreaks = cut(median, breaks = c(-1, -0.025, 0 , 0.025, 1)))

#

ggplot(rIAllMinus, aes(var, y = value2, fill = rowname)) +
  geom_col() +
  coord_flip() +
  scale_x_reverse('Variable', 
                     breaks = 1:7,
                     labels = c('Temp', 'Flow', 'Bkt', 'Bnt','Ats','Temp*Flow',  "Bkt*Bnt")) +
  scale_y_continuous('Relative proportion of variance explained',
                     breaks = c(-1.0, -0.5, 0.0, 0.5),
                     labels = c('1.0', '0.5', '0.0', '0.5')) +
  geom_text(data = summ1a, aes(y=xAxisLoc, x=7, label=paste(round(prop,2)), color=rowname),size = 5) +
  scale_color_grey(start = 0.15, end = 0.75) +
  theme_publication() +
  scale_fill_grey(start = 0.15, end = 0.75) +
  facet_grid(season ~ ageGroup, labeller=labeller(season = labelsSeason2, ageGroup = labelAge))

ggsave(paste0('figures/3Spp_WB/relImpOccAll_GR_HorizBars.png'), width = 11, height = 9, units='in')

# bars colered by beta
ggplot(rIAllMinusQ, aes(var, y = value2)) +
  geom_col(aes(fill = median)) +
  coord_flip() +
  scale_x_reverse('Variable', 
                     breaks = 1:7,
                     labels = c('Temp', 'Flow', 'Bkt', 'Bnt','Ats','Temp*Flow',  "Bkt*Bnt")) +
  scale_y_continuous('Relative proportion of variance explained',
                     breaks = c(-0.9, -0.25, -0.5, 0.0, 0.25, 0.5, 0.9),
                     labels = c('0.9', '0.25', '0.5\n\nBrook Trout\n', '0.0', '0.25', '0.5\n\nBrown Trout\n', '0.9')) +
  geom_text(data = summ1a, aes(y=xAxisLoc, x=7, label=paste(round(prop,2))),size = 5) +
  geom_hline(yintercept = 0) +
  #scale_color_grey(start = 0.15, end = 0.75) +
  theme_publication() +
 #  scale_fill_continuous() +
  scale_fill_distiller(name="Parameter\nestimate", palette = "RdBu", direction = 1) +
  theme(legend.position="right",             
        #legend.key = element_rect(colour = NA),
        legend.direction = "vertical",
        legend.key.size= unit(0.25, "in"),
        #    legend.margin = unit(0, "cm"),
        legend.title = element_text(face="plain")#element_blank()
        ) +
  facet_grid(season ~ ageGroup, labeller=labeller(season = labelsSeason2, ageGroup = labelAge))

ggsave(paste0('figures/3Spp_WB/relImpOccAll_GR_HorizBarsColor.png'), width = 11, height = 9, units='in')


#######################
# coefs
labelsC <- c('1' = 'intercept', '2' = 'Temp','3' = 'Flow','4' = 'Bkt','5' = 'Bnt','6' = 'Ats','7' = 'Temp*Flow', '8' = "Bkt*Bnt")


ggplot(cAll %>% filter(var %in% 2:8, rowname %in% c('bktC','bntC')), aes((occ), value, color = factor(rowname))) +
  geom_point(size=3) +
  geom_line() +
  geom_vline(xintercept = 3.5) +
  scale_x_continuous('Sampling occasion',
                    breaks = 2:10,
                    labels = c('Winter','Spring','Summer', 'Autumn','Winter','Spring','Summer', 'Autumn','Winter')
                    ) +
  scale_y_continuous('Coefficients') +
  theme_publication() +
  theme(axis.text.x  = element_text(angle=90, vjust=0.25, size=10)) +
  facet_wrap(~ var, labeller=labeller(var = labelsC))#, scales = 'free')

ggsave(paste0('figures/3Spp_WB/coefsOccAll.png'), width = 11, height = 9, units='in')

# intercept only
cAll2 <- rbind(cAll,data.frame(rowname=c('bktC','bntC','atsC'),value=simData$meanInit,var=1,occ=1,season=3))
ggplot(cAll2 %>% filter(var == 1), aes((occ), value, color = factor(rowname))) +
  geom_point(size=4) +
  geom_line() +
  geom_vline(xintercept = 3.5) +
  scale_x_continuous('Sampling occasion',
                    breaks = 1:10,
                    labels = c('Autumn','Winter','Spring','Summer', 'Autumn','Winter','Spring','Summer', 'Autumn','Winter')
                    ) +
  scale_y_continuous('Intercept') +
  theme_publication() +
  theme(axis.text.x  = element_text(angle=90, vjust=0.25, size=10))

ggsave(paste0('figures/3Spp_WB/coefsInterceptOccAll.png'), width = 11, height = 9, units='in')




#####################################
# no need to start sims with diff sizes by species
yoyMeans <- 
cd %>% filter( riverOrdered == 'west brook', ageInSamples == 1 ) %>%
  group_by(species) %>%
  summarize(meanLen = mean(observedLength,na.rm = TRUE),
           sdLen = sd(observedLength,na.rm = TRUE))

#   species meanLen sdLen
#   <chr>     <dbl> <dbl>
# 1 ats        69.0  5.86
# 2 bkt        72.7  9.94
# 3 bnt        72.9  8.73
#####################################
# make figures of relative importance for each sampling occasion
rIAll <- cAll <- NULL
for( occSelected in 2:simData$nOcc){
  print(occSelected)
  lmAll <- ggRIAll <- relImpAll <- coefsAll <- list()

  for(spp1 in spp){

      lmAll[[spp1]] <- lm( lenMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt, data = stats3Spp[[spp1]][[1]] %>% filter(occ == occSelected))

      relImpAll[[spp1]] <- calc.relimp(lmAll[[spp1]], type = c("lmg"), rela = TRUE)

      coefsAll[[spp1]] <- lmAll[[spp1]]$coefficients
  }

  # make bar chart of relative importance of each variable for each spp
  rI <- cbind(data.frame(bktRI = relImpAll[['bkt']]$lmg), 
              data.frame(bntRI = relImpAll[['bnt']]$lmg), 
              data.frame(atsRI = relImpAll[['ats']]$lmg)) %>%
        rownames_to_column() %>% 
        gather(key=rowname) %>% 
        mutate(var = rep(1:length(relImpAll[['bkt']]$lmg),length(spp)),
               occ = occSelected,
               season = simData$seasons[occSelected])

  rIAll <- rbind(rIAll,rI)

  ggplot(rI, aes(factor(var),value, fill = factor(rowname, levels = c('bktRI','bntRI','atsRI'), ordered = TRUE))) +
    geom_col(position='dodge') +
    scale_x_discrete('Parameter',
                     labels = c('1' = 'Temp','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")) +
    scale_fill_grey() +
    scale_y_continuous('Relative proportion of variance explained') +
    ggtitle(paste0("occ = ",occSelected,"  Season = ",simData$seasons[occSelected])) +
    theme_publication() +
    theme( legend.position = 'none')

##############################  ggsave(paste0('figures/3Spp_WB/relImpOcc',occSelected,'.png'), width = 9, height = 7, units='in')

  # coefficients
  coefs <- cbind(data.frame(bktC = coefsAll[['bkt']]), 
                 data.frame(bntC = coefsAll[['bnt']]), 
                 data.frame(atsC = coefsAll[['ats']])) %>%
        rownames_to_column() %>% 
        gather(key=rowname) %>% 
        mutate(var = rep(1:length(coefsAll[['bkt']]),length(spp)),
               occ = occSelected,
               season = simData$seasons[occSelected])

  cAll <- rbind(cAll,coefs)
}  


######################################################
# relImp over occasions for each spp and parameter
labels <- c('1' = 'Temp','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")

# leaving out bkt*bnt, too small
ggplot(rIAll %>% filter(rowname %in% c('bktRI','bntRI')), aes((occ), value, color = factor(rowname))) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = 3.5) +
  scale_x_continuous('Sampling occasion',
                    breaks = 2:10,
                    labels = c('Winter','Spring','Summer', 'Autumn','Winter','Spring','Summer', 'Autumn','Winter')
                    ) +
  scale_y_continuous('Relative proportion of variance explained') +
  theme_publication() +
  theme(axis.text.x  = element_text(angle=90, vjust=0.25, size=10)) +
  facet_wrap(~ var, labeller=labeller(var = labels))

ggsave(paste0('figures/3Spp_WB/relImpOccAll.png'), width = 11, height = 9, units='in')

rIMeans <- rIAll %>%
  group_by(rowname,var) %>%
  summarize( mean = mean(value),
             sd = sd(value),
             stderr = sd(value)/sqrt(length(value))) %>%
  mutate( up = mean + stderr )

dodge <- position_dodge(width=0.9)
ggplot(rIMeans %>% filter(rowname %in% c('bktRI','bntRI')), aes(factor(var),mean, fill = factor(rowname, levels = c('bktRI','bntRI','atsRI'), ordered = TRUE))) +
  geom_errorbar(aes(ymin = mean-0.01, ymax = up), width = 0.25, position=dodge) +
  geom_col(position=dodge) +
  scale_x_discrete('Parameter',
                   labels = c('1' = 'Temperature','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")) +
  scale_fill_grey() +
  scale_y_continuous('Relative proportion of variance explained') +
  theme_publication() +
  theme( legend.position = 'none')

ggsave(paste0('figures/3Spp_WB/relImpOccMean.png'), width = 9, height = 7, units='in')

# mean rel Imp for env vs spp variables
rIMeans %>%
  mutate(isEnvVar = ifelse(var %in% c(1,2,6), 1,0)) %>%
  group_by(rowname,isEnvVar) %>%
  summarize(meanMean = mean(mean),
            sumMean = sum(mean),
            sdMean = sd(mean))

# ratio of env/(env+spp) means
.247/(.246+.0648) # bkt
.122/(.122+.159) # bnt

# ratio of env/(env+spp) sums - donT need to do ratio, just use sum for env
.635/.259 # den-dep bnt/bkt
.741/.365 # env bkt/bnt

#######################
# coefs
labelsC <- c('1' = 'intercept', '2' = 'Temp','3' = 'Flow','4' = 'Bkt','5' = 'Bnt','6' = 'Ats','7' = 'Temp*Flow', '8' = "Bkt*Bnt")


ggplot(cAll %>% filter(var %in% 2:8, rowname %in% c('bktC','bntC')), aes((occ), value, color = factor(rowname))) +
  geom_point(size=3) +
  geom_line() +
  geom_vline(xintercept = 3.5) +
  scale_x_continuous('Sampling occasion',
                    breaks = 2:10,
                    labels = c('Winter','Spring','Summer', 'Autumn','Winter','Spring','Summer', 'Autumn','Winter')
                    ) +
  scale_y_continuous('Coefficients') +
  theme_publication() +
  theme(axis.text.x  = element_text(angle=90, vjust=0.25, size=10)) +
  facet_wrap(~ var, labeller=labeller(var = labelsC))#, scales = 'free')

ggsave(paste0('figures/3Spp_WB/coefsOccAll.png'), width = 11, height = 9, units='in')

# intercept only
cAll2 <- rbind(cAll,data.frame(rowname=c('bktC','bntC','atsC'),value=simData$meanInit,var=1,occ=1,season=3))
ggplot(cAll2 %>% filter(var == 1), aes((occ), value, color = factor(rowname))) +
  geom_point(size=4) +
  geom_line() +
  geom_vline(xintercept = 3.5) +
  scale_x_continuous('Sampling occasion',
                    breaks = 1:10,
                    labels = c('Autumn','Winter','Spring','Summer', 'Autumn','Winter','Spring','Summer', 'Autumn','Winter')
                    ) +
  scale_y_continuous('Intercept') +
  theme_publication() +
  theme(axis.text.x  = element_text(angle=90, vjust=0.25, size=10))

ggsave(paste0('figures/3Spp_WB/coefsInterceptOccAll.png'), width = 11, height = 9, units='in')




#####################################
# no need to start sims with diff sizes by species
yoyMeans <- 
cd %>% filter( riverOrdered == 'west brook', ageInSamples == 1 ) %>%
  group_by(species) %>%
  summarize(meanLen = mean(observedLength,na.rm = TRUE),
           sdLen = sd(observedLength,na.rm = TRUE))

#   species meanLen sdLen
#   <chr>     <dbl> <dbl>
# 1 ats        69.0  5.86
# 2 bkt        72.7  9.94
# 3 bnt        72.9  8.73
#####################################
round(mcmc$summary$all.chains[startsWith(attributes(mcmc$summary$all.chains)$dimnames[[1]], "grBeta["),],2)
# relative proportion stats with occ in the lm 
stats3Spp[['bkt']][[1]]$season <- as.factor(c(3,4,1,2,3,4,1,2,3,4))
stats3Spp[['bkt']][[1]]$temp2 <- stats3Spp[['bkt']][[1]]$temp^2
stats3Spp[['bkt']][[1]]$flow2 <- stats3Spp[['bkt']][[1]]$flow^2


ageForRI <-  'adult' #'yoy' or 'adult'
occs <- if (ageForRI == 'yoy') 1:3 else 4:10
indepVarVals <- c(-1.5,0,1.5) # make sure symetrical around 0
rIBKTBNT <- NULL
# loop over occs
for(rep in 1:simData$nReps){

  print(rep)

  datOccBKT <- stats3Spp[['bkt']][[1]] %>% filter(iter == rep, occ %in% occs, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals)
  lmAllBKT <- lm( lenMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt + occ, data = datOccBKT) 
  relImpAllBKT <- calc.relimp(lmAllBKT, type = c("lmg"), rela = TRUE)

  datOccBNT <- stats3Spp[['bnt']][[1]] %>% filter(iter == rep, occ %in% occs, temp %in% indepVarVals, flow %in% indepVarVals, bkt %in% indepVarVals, bnt %in% indepVarVals, ats %in% indepVarVals)
  lmAllBNT <- lm( lenMean ~ temp + flow + temp * flow + temp^2 + flow^2 + bkt + bnt + ats + bnt * bkt + occ, data = datOccBNT)
  relImpAllBNT <- calc.relimp(lmAllBNT, type = c("lmg"), rela = TRUE)

  rBKT <- relImpAllBKT$lmg[2:8] / sum(relImpAllBKT$lmg[2:8]) 
  rBNT <- relImpAllBNT$lmg[2:8] / sum(relImpAllBNT$lmg[2:8])

  rIBKTBNT <- rbind(rIBKTBNT, data.frame(cbind(rBKT,rBNT)) %>% 
                      rownames_to_column() %>% 
                      gather(key = rowname) %>% 
                      mutate(var = rep(1:7,2),
                             varText = rep(names(rBKT),2),
                             iter = rep)
                    )

}


rIBKTBNTMeans <- rIBKTBNT %>%
  group_by(rowname, var) %>%
  summarise(rIMean = mean(value),
            rISd = sd(value),
            rIUpperCI = quantile(value, 0.95))

dodge <- position_dodge(width=0.9)
ggplot(rIBKTBNTMeans, aes(factor(var),rIMean, fill = factor(rowname, levels = c('rBKT','rBNT'), ordered = TRUE))) +
  geom_errorbar(aes(ymin = rIMean - 0.001, ymax = rIUpperCI), width = 0.25, position=dodge) +
  geom_col(position=dodge) +
  scale_x_discrete('Parameter',
                   labels = c('1' = 'Temperature','2' = 'Flow','3' = 'Bkt','4' = 'Bnt','5' = 'Ats','6' = 'Temp*Flow', '7' = "Bkt*Bnt")) +
  scale_fill_grey() +
  scale_y_continuous('Relative proportion of variance explained') +
  theme_publication() +
  geom_hline(yintercept = 0) +
  theme( legend.position = 'none')

ggsave(paste0('figures/3Spp_WB/relImpOccMeanModWOcc_Size',simData$nReps,'_',ageForRI,'.png'), width = 9, height = 7, units='in')

# uses 42 G - does not run
#lmAll <- lm( lenMean ~ (temp + flow +  temp^2 + flow^2 + bkt + bnt + ats )*occ + temp * flow + bnt * bkt, data = stats3Spp[['bkt']][[1]] %>%               filter(iter %in% 1) )

rIBKTBNT %>%
  mutate(isEnvVar = ifelse(var %in% c(1,2,6), 1,0)) %>%
  group_by(rowname,isEnvVar) %>%
  summarize(meanMean = mean(value),
            sumMean = sum(value),
            sdMean = sd(value))
# 25 reps
#   rowname isEnvVar meanMean sumMean sdMean
#   <chr>      <dbl>    <dbl>   <dbl>  <dbl>
# 1 rBKT           0   0.0664   0.265  0.105
# 2 rBKT           1   0.245    0.735  0.249
# 3 rBNT           0   0.147    0.588  0.104
# 4 rBNT           1   0.137    0.412  0.158

# 100 reps
#   rowname isEnvVar meanMean sumMean sdMean
#   <chr>      <dbl>    <dbl>   <dbl>  <dbl>
# 1 rBKT           0   0.0676    27.0 0.0914
# 2 rBKT           1   0.243     73.0 0.200 
# 3 rBNT           0   0.146     58.2 0.0922
# 4 rBNT           1   0.139     41.8 0.131

# ratio of env/(env+spp) means
.243/(.243 + .0676) # bkt 0.7823567
.139/(.139 + .146) # bnt 0.4877193

.782/.488 # 1.60
##############################################################
############################################################

# projOut <- list(stats3Spp, lm, relImp, anov, elapsed)
load(paste0('D:/projects/linkedModels/data/out/projectionOut_N',simData$nReps,'.RData'))

spp <- c('bkt','bnt','ats')
# scenarios
# Plot subsets of the simulations
#means for each iter
stats1 <- rbind(data.frame(stats3Spp[['bkt']][[1]]), 
                data.frame(stats3Spp[['bnt']][[1]]), 
                data.frame(stats3Spp[['ats']][[1]])) %>%
          mutate(spp = rep(spp, each = nrow(stats3Spp[['bkt']][[1]])))

#means for each occ
stats2 <- rbind(data.frame(stats3Spp[['bkt']][[2]]), 
                data.frame(stats3Spp[['bnt']][[2]]), 
                data.frame(stats3Spp[['ats']][[2]])) %>%
          mutate(spp = rep(spp, each = nrow(stats3Spp[['bkt']][[2]])))

ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, ats == 0), aes(occ,lenMeanMean, group = spp)) +
  geom_point(aes(color = spp) ) +
  geom_line( aes(color = spp) ) +
  geom_point(aes(color = spp), data = stats2 %>% filter(temp == 1.5, flow == 0, bkt == 0, bnt == 0, ats == 0) ) +
  geom_line( aes(color = spp), data = stats2 %>% filter(temp == 1.5, flow == 0, bkt == 0, bnt == 0, ats == 0), linetype = 2 ) 


ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, ats == 0), aes(occ,lenMeanMean, group = spp)) +
  geom_point(aes(color = spp) ) +
  geom_line( aes(color = spp) ) +
  geom_point(aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 1.5, bkt == 0, bnt == 0, ats == 0) ) +
  geom_line( aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 1.5, bkt == 0, bnt == 0, ats == 0), linetype = 2 ) 

ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, ats == 0), aes(occ,lenMeanMean, group = spp)) +
  geom_point(aes(color = spp) ) +
  geom_line( aes(color = spp) ) +
  geom_point(aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 0, bkt == 1.5, bnt == 0, ats == 0) ) +
  geom_line( aes(color = spp), data = stats2 %>% filter(temp == 0, flow == 0, bkt == 1.5, bnt == 0, ats == 0), linetype = 2 ) 


# trajectories across indep var values for different combinations. Trying to relate to relImp results
#temp
p1 <- ggplot( stats2 %>% filter(flow == 0, bkt == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(temp))) +
  geom_line(aes(color = factor(temp))) +
  geom_point(aes(color = factor(temp))) +
  facet_wrap(~spp) 

#flow
p2 <- ggplot( stats2 %>% filter(temp == 0, bkt == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(flow))) +
  geom_line(aes(color = factor(flow))) +
  geom_point(aes(color = factor(flow))) +
  facet_wrap(~spp)  

#bkt
p3 <- ggplot( stats2 %>% filter(temp == 0, flow == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(bkt))) +
  geom_line(aes(color = factor(bkt))) +
  geom_point(aes(color = factor(bkt))) +
  facet_wrap(~spp) 

#bnt
p4 <- ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(bnt))) +
  geom_line(aes(color = factor(bnt))) +
  geom_point(aes(color = factor(bnt))) +
  facet_wrap(~spp)

#ats
p5 <- ggplot( stats2 %>% filter(temp == 0, flow == 0, bkt == 0, bnt == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(ats))) +
  geom_line(aes(color = factor(ats))) +
  geom_point(aes(color = factor(ats))) +
  facet_wrap(~spp)

multiplot(p1, p2, p3, p4, p5, cols=2) # multiplot function is below

#temp*flow
ggplot( stats2 %>% filter(bkt == 0, bnt == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(flow))) +
  geom_line(aes(color = factor(flow))) +
  geom_point(aes(color = factor(flow))) +
  facet_grid(temp~spp)

#bkt*bnt
ggplot( stats2 %>% filter(temp == 0, flow == 0, ats == 0, spp %in% c('bkt','bnt')), aes(occ,lenMeanMean, group = factor(bnt))) +
  geom_line(aes(color = factor(bnt))) +
  geom_point(aes(color = factor(bnt))) +
  facet_grid(bkt~spp)




# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}



# size=100
# gr=10
# sd=1
# mean(rnorm(10000,size+gr,sd))
# mean(size+rnorm(10000,gr,sd))
pd <- position_dodge(0.65)
d2$speciesGG <- factor(d2$species, levels = c('bkt','bnt','ats'), labels = c("Brook Trout", "Brown Trout", "Atlantic Salmon"), ordered = T)


# isYOY1
d2 %>% group_by(speciesGG,riverOrderedGG,seasonGG,year) %>%
  filter(isYOYDATA == 1, riverOrderedGG == 'West Brook') %>%
  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_wrap(~seasonGG )

ggsave('figures/3Spp_WB/grByYearYOY1.png')

# isYOY2
d2 %>% group_by(speciesGG,riverOrderedGG,seasonGG,year) %>%
  filter(isYOYDATA == 2, riverOrderedGG == 'West Brook') %>%
  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_wrap(~seasonGG )


ggsave('figures/3Spp_WB/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.