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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.