Mike

code/biology era regression hier.R

@@ -36,49 +36,75 @@ npgo$win.yr <- ifelse(npgo$month %in% 11:12, npgo$Year+1, npgo$Year)
 win.npgo <- tapply(npgo$value, npgo$win.yr, mean)
 
 # and smoothed (2yr) values of each
-npgo2 <- rollapply(win.npgo, 2, mean, align="right", fill=NA)
-pdo2 <- rollapply(win.pdo, 2, mean, align="right", fill=NA)
-names(pdo2) <- 1900:2019
+win.npgo <- rollapply(win.npgo, 2, mean, align="right", fill=NA)
+names(win.npgo) <- 1950:2019
+win.pdo <- rollapply(win.pdo, 2, mean, align="right", fill=NA)
+names(win.pdo) <- 1900:2019
 
-# load three "non-salmon" data sets: GOA crustaceans/fish, Farallon seabirds, CalCOFI ichthyo
+# load four "non-salmon" data sets: EBS groundfish recruitment, GOA crustaceans/fish, Farallon seabirds, CalCOFI ichthyo
 dat <- read.csv("data/farallon.sbrd.biol.csv", row.names = 1)
+
+# examine distributions
+look <- dat %>%
+  gather(key, value, -year)
+
+ggplot(look, aes(value)) +
+  geom_histogram() +
+  facet_wrap(~key, scales="free")
+
 # add era term
 dat$era <- as.factor(ifelse(dat$year <= 1988, 1, 2))
 
 # and pdo/npgo
-dat$pdo <- pdo2[match(dat$year, names(pdo2))]
-dat$npgo <- npgo2[match(dat$year, names(npgo2))]
+dat$pdo <- win.pdo[match(dat$year, names(win.pdo))]
+dat$npgo <- win.npgo[match(dat$year, names(win.npgo))]
 
 # reshape with year, era, and pdo and npgo as the grouping variables
 melted <- melt(dat, id.vars = c("year","pdo","era","npgo"))
 melted$variable_era = paste0(melted$era,melted$variable)
 # standardize all the time series by variable -- so slopes are on same scale
 m1 = dplyr::group_by(melted, variable) %>%
   mutate(scale_x = scale(value))
-m1$system <- "Farallon seabirds"
+m1$system <- "Central California Current"
 
 ####
 # CalCOFI
 dat <- read.csv("data/calcofi.biol.csv", row.names=1)
 
+# examine distributions
+look <- dat %>%
+  gather(key, value, -year)
+
+ggplot(look, aes(value)) +
+  geom_histogram() +
+  facet_wrap(~key, scales="free")
+
 dat$era <- as.factor(ifelse(dat$year <= 1988, 1, 2))
 
 # and pdo/npgo
-dat$pdo <- pdo2[match(dat$year, names(pdo2))]
-dat$npgo <- npgo2[match(dat$year, names(npgo2))]
+dat$pdo <- win.pdo[match(dat$year, names(win.pdo))]
+dat$npgo <- win.npgo[match(dat$year, names(win.npgo))]
 
 # reshape with year, era, and pdo and npgo as the grouping variables
 melted <- melt(dat, id.vars = c("year","pdo","era","npgo"))
 melted$variable_era = paste0(melted$era,melted$variable)
 # standardize all the time series by variable -- so slopes are on same scale
 m2 = dplyr::group_by(melted, variable) %>%
   mutate(scale_x = scale(value))
-m2$system <- "CalCOFI ichthyoplankton"
+m2$system <- "Southern California Current"
 
 #########
-
 dat <- read.csv("data/goa.biol.csv")
 colnames(dat)[1] <- "year"
+
+# examine distributions
+look <- dat %>%
+  gather(key, value, -year)
+
+ggplot(look, aes(value)) +
+  geom_histogram() +
+  facet_wrap(~key, scales="free")
+
 dat$era <- as.factor(ifelse(dat$year <= 1988, 1, 2))
 
 # and pdo/npgo
@@ -91,12 +117,41 @@ melted$variable_era = paste0(melted$era,melted$variable)
 # standardize all the time series by variable -- so slopes are on same scale
 m3 = dplyr::group_by(melted, variable) %>%
   mutate(scale_x = scale(value))
-m3$system <- "GOA fish & crustaceans"
+m3$system <- "Gulf of Alaska"
+
+#######
+dat <- read.csv("data/ebs.biol.data.csv")
+dat[,2:5] <- sqrt(dat[,2:5])
+
+# examine distributions
+look <- dat %>%
+  gather(key, value, -year)
+
+ggplot(look, aes(value)) +
+  geom_histogram() +
+  facet_wrap(~key, scales="free")
+
+dat$era <- as.factor(ifelse(dat$year <= 1988, 1, 2))
+
+# and pdo/npgo
+dat$pdo <- win.pdo[match(dat$year, names(win.pdo))]
+dat$npgo <- win.npgo[match(dat$year, names(win.npgo))]
+
+# reshape with year, era, and pdo and npgo as the grouping variables
+melted <- melt(dat, id.vars = c("year","pdo","era","npgo"))
+melted$variable_era = paste0(melted$era,melted$variable)
+melted$value <- as.numeric(melted$value)
+
+# standardize all the time series by variable -- so slopes are on same scale
+m4 = dplyr::group_by(melted, variable) %>%
+  mutate(scale_x = scale(value))
+m4$system <- "Bering Sea"
+
 
 #######
 # combine
 
-melted <- rbind(m1, m2, m3)
+melted <- rbind(m1, m2, m3, m4)
 melted$year <- as.numeric(melted$year)
 melted$variable <- as.factor(melted$variable)
 melted$variable_era <- as.factor(melted$variable_era)
@@ -132,7 +187,7 @@ for(s in levels.syst) {
   pars = rstan::extract(mod,permuted=TRUE)
 
   model.data = rbind(model.data,
-    data.frame(system=unique(melted$system)[s],
+    data.frame(system=s,
       ratio=100*exp(pars$mu_ratio)))
 
   temp$pred = apply(pars$pred,2,mean)
@@ -171,22 +226,22 @@ for(s in levels.syst) {
 
 }
 
+# add a placeholder 
+placeholder <- data.frame(system="Northern California Current", ratio=NA)
+model.data <- rbind(model.data, placeholder)
+
 # order the systems north-south
-model.data$order <- ifelse(model.data$system=="GOA fish & crustaceans", 1,
-                           ifelse(model.data$system=="Farallon seabirds", 2, 3))
-model.data$system <- reorder(model.data$system, model.data$order)
+model.data$order <- ifelse(model.data$system=="Bering Sea", 1,
+                           ifelse(model.data$system=="Gulf of Alaska", 2, 
+                                  ifelse(model.data$system=="Northern California Current", 3, 
+                                         ifelse(model.data$system=="Central California Current", 4, 5))))
 
-cb <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
+model.data$system <- reorder(model.data$system, model.data$order)
 
 pdo.data <- model.data
 
-pdo.plot <- ggplot(pdo.data, aes(ratio/100)) + # just removing % for now
-  theme_linedraw() +
-  geom_density(fill=cb[3]) + xlab("Avg ratio: Era 1 slope / Era 2 slope") +
-  facet_wrap(~system, ncol=1) +
-  xlim(c(0,2)) +
-  geom_vline(xintercept = 1)+
-  ggtitle("a) PDO")
+# save for future reference
+write.csv(pdo.data, "models/pdo_biology_model_data.csv")
 
 #################
 ## and the same thing for npgo
@@ -254,25 +309,22 @@ for(s in levels.syst) {
   dev.off()
 }
 
+# add a placeholder 
+placeholder <- data.frame(system="Northern California Current", ratio=NA)
+model.data <- rbind(model.data, placeholder)
+
 # order the systems north-south
-model.data$order <- ifelse(model.data$system=="GOA fish & crustaceans", 1,
-                           ifelse(model.data$system=="Farallon seabirds", 2, 3))
+model.data$order <- ifelse(model.data$system=="Bering Sea", 1,
+                           ifelse(model.data$system=="Gulf of Alaska", 2, 
+                                  ifelse(model.data$system=="Northern California Current", 3, 
+                                         ifelse(model.data$system=="Central California Current", 4, 5))))
+
 model.data$system <- reorder(model.data$system, model.data$order)
 
 npgo.data <- model.data
 
-npgo.plot <- ggplot(npgo.data, aes(ratio/100)) +
-  theme_linedraw() +
-  geom_density(fill=cb[3]) + xlab("Avg ratio: Era 1 slope / Era 2 slope") +
-  facet_wrap(~system, ncol=1) +
-  geom_vline(xintercept=1) +
-  xlim(c(0,4.5)) +
-  ggtitle("b) NPGO")
-
-png("biol regression change pdo-npgo slope.png", 7, 7, units="in", res=300)
-ggarrange(pdo.plot, npgo.plot, ncol=2)
-dev.off()
-
+# save for future reference
+write.csv(npgo.data, "models/npgo_biology_model_data.csv")
 
 # Caterpillar Plot ===============================
 # Helper Functions
@@ -286,25 +338,35 @@ head(pdo.data)
 npgo.data$var <- "NPGO"
 pdo.data$var <- "PDO"
 all.data <- rbind(pdo.data, npgo.data)
+all.data$var.order <- ifelse(all.data$var=="PDO", 1, 2)
+all.data$var <- reorder(all.data$var, all.data$var.order)
+all.data$log.ratio <- log(all.data$ratio/100, 10)
+
+# colorblind...
+cb <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
 
-cat.plt <- ggplot(all.data, aes(x=system, y=ratio/100, fill=system)) +
+biol.plt <- ggplot(all.data, aes(x=reorder(system, desc(system)), y=log.ratio, fill=system)) +
              theme_linedraw() +
+            scale_fill_manual(values=cb[c(6,3,4,2,8)], 
+                              labels=c("Bering Sea", "Gulf of Alaska", "Central Cal. Curr.", "Southern Cal. Curr.")) +
              # scale_fill_colorblind() +
-             scale_fill_tableau() +
+             # scale_fill_tableau() +
              # scale_fill_brewer(c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")) +
              # geom_eye() +
-
+  
              geom_violin(alpha = 0.75, lwd=0.1, scale='width') +
              stat_summary(fun.y="q.95", colour="black", geom="line", lwd=0.75) +
              stat_summary(fun.y="q.50", colour="black", geom="line", lwd=1.5) +
              stat_summary(fun.y="median", colour="black", size=2, geom="point", pch=21) +
              facet_wrap(~var, ncol=1) +
-             ylab("Avg ratio: Era 1 slope / Era 2 slope") +
-             theme(axis.text.y = element_blank()) +
-             geom_hline(aes(yintercept=1), color="red", linetype="dotted", size=1) +
-             coord_flip(ylim=c(0,7))
-
-cat.plt
+             ylab("Log ratio: Era 1 slope / Era 2 slope") +
+             theme(axis.text.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_line(size=0),
+                   legend.title = element_blank(), legend.position = c(0.15, 0.15)) +
+             geom_hline(aes(yintercept=0), color="red", linetype="dotted", size=1) +
+             coord_flip(ylim=c()) 
+
+
+biol.plt
 ggsave("biol regression change pdo-npgo slope_cater.png", plot=cat.plt,
          height=7, width=7, units="in", dpi=300)
 

data/ebs.biol.data.csv

@@ -0,0 +1,57 @@
+year,yellowfin,turbot,flathead.sole,pollock
+1960,,46073,,
+1961,,82675,,
+1962,,185807,,
+1963,,437867,,
+1964,733,484027,455615,6294
+1965,749,204888,435985,20951
+1966,1568,85070,422282,15011
+1967,2395,44596,418335,25590
+1968,2451,29672,427371,22171
+1969,3227,24010,448528,26244
+1970,3581,23147,470738,23562
+1971,2817,26270,476659,14401
+1972,1951,35638,476513,11755
+1973,2620,59529,427138,27009
+1974,3078,119004,501793,19772
+1975,2027,209981,609171,16758
+1976,2551,147543,621030,12831
+1977,1670,93480,558901,13276
+1978,1067,49205,804832,24458
+1979,2065,17637,1065140,58151
+1980,1541,6022,1223000,25871
+1981,4464,1025,2262170,29475
+1982,827,1999,976982,16662
+1983,3697,3289,702171,50955
+1984,1278,6185,1352310,14281
+1985,982,20233,2121820,34580
+1986,1341,5317,919482,14019
+1987,1837,5756,2409590,7524
+1988,1836,5897,579668,5573
+1989,917,15803,1063110,10770
+1990,1030,3848,1202660,47657
+1991,2285,1131,1032960,25052
+1992,1361,743,1003370,22026
+1993,1148,606,795886,45352
+1994,1156,954,771220,15169
+1995,2850,3771,693230,10371
+1996,1231,1612,1207050,22524
+1997,1020,1641,1154200,30739
+1998,1249,2145,1049860,15088
+1999,1767,8332,751622,16347
+2000,1124,9648,788530,25565
+2001,1547,11358,1095690,34834
+2002,1503,1677,1362720,23412
+2003,2388,614,1275090,14270
+2004,1042,547,333892,6451
+2005,1258,801,819649,4629
+2006,1607,7292,737612,11769
+2007,1446,22137,555717,25225
+2008,1550,48122,766015,13623
+2009,1947,33061,721013,51287
+2010,1482,6103,489730,21394
+2011,,4003,1828830,12826
+2012,,1904,714452,11188
+2013,,2430,1429300,56864
+2014,,2063,1613410,42237
+2015,,2458,,12693