SST x growth EDA

Read in SST and growth marker data.

sst<-read_csv(paste(sst.path,"/migrationareasERSST.csv",sep = ""))
Carlin_full_growth<-readRDS(paste(growth.path,"/Carlin_growth_fulldata.rds",sep = ""))

Clip SST data to applicable years
Calculate mean annual SST for years in question

sst_summarized_long<-sst %>% 
  filter(Dates >= "1963-01-01", 
         Dates < "1992-01-01") %>% 
  mutate(SmoltYear = year(Dates)) %>% 
  group_by(SmoltYear, Area) %>% 
  dplyr::summarize(AnnualSST = mean(ERSST))

Select growth markers of interest
Calculate annual mean growth marker values by smolt year

growth_summarized_long<- Carlin_full_growth %>% 
  filter(!is.na(ReleaseYear)) %>% 
  dplyr::select(JoinID, ReleaseYear, RecaptureYear, FMC, M1, M2, marine.incr, FS.incr, FW.incr, PS.incr, 
         smolt.circ, FS.circ, FW.circ, PS.circ) %>% 
  gather(key = "marker", value = "value", 4:14) %>% 
  group_by(ReleaseYear, marker) %>% 
  dplyr::summarize(meanmarker = mean(value, na.rm = TRUE)) %>% 
  dplyr::rename("SmoltYear" = "ReleaseYear")

Pairwise-combination of all growth marker data and all SST data

modelinglong<- growth_summarized_long %>% 
  left_join(sst_summarized_long, by = "SmoltYear") %>% 
  dplyr::rename("responsevar" = "marker", "predictvar" = "Area", "meanSST" = "AnnualSST")

Loop through all possible model combinations (singular pred var and no interactions)

responses<-as.vector(unique(as.character(modelinglong$responsevar)))
predictors<-as.vector(unique(as.character(modelinglong$predictvar)))

AIClist<-list()
r2list<-list()
pvlist<-list()

for(i in responses){
  for(j in predictors){
    
temp<-modelinglong %>% 
  filter(responsevar == i & predictvar == j)

glancesummary<-glance(lm(meanmarker ~ meanSST, data = temp))

AIClist[[i]][j]<-glancesummary$AIC
r2list[[i]][j]<-glancesummary$r.squared
pvlist[[i]][j]<-glancesummary$p.value
  }
}

AIClist<-do.call(rbind,AIClist)
r2list<-do.call(rbind,r2list)
pvlist<-do.call(rbind,pvlist)

#Reformat model summary outputs

pvdf<-pvlist %>% 
  as_tibble() %>% 
  mutate(Marker = unique(as.character(modelinglong$responsevar))) %>% 
  gather(key = "Area", value = "p.value", 1:6)

AICdf<-AIClist %>% 
  as_tibble() %>% 
  mutate(Marker = unique(as.character(modelinglong$responsevar))) %>% 
  gather(key = "Area", value = "AIC", 1:6)

r2df<-r2list %>% 
  as_tibble() %>% 
  mutate(Marker = unique(as.character(modelinglong$responsevar))) %>% 
  gather(key = "Area", value = "Rsquared", 1:6)

model_summary_df<-AICdf %>% 
  left_join(r2df, by = c("Marker", "Area")) %>% 
  left_join(pvdf, by = c("Marker", "Area")) %>% 
  filter(p.value < 0.05 | AIC == min(AIC)) %>% 
  dplyr::rename("responsevar" = "Marker", "predictvar" = "Area") %>% 
  left_join(modelinglong, by = c("responsevar","predictvar"))

Combine all AIC, r2, and p-value into a single model summary data frame. Keep only models with p-values less than 0.05 or models which have the lowest AIC for each marker (redundant for my data, but may not be for yours).
Add back in the actual values associated with each model (join modelinglong).

df<-data.frame(
  responsevar = c(rep("marine.incr",5),rep("PS.incr",5),rep("FS.incr",5),rep("FW.incr",5)), 
  predictvar = c(rep(c("n.labsea","nf","northatlantic","s.labsea","wgreen"),4)),
  x = c(rep(c(5.1, 5.3, 5, 7.8, 2.6),4)),
  y = c(rep(4,5),rep(2,5),rep(0.51,5), rep(1.5,5)), 
  pval = c(0.10,0.07,0.38,0.84,0.05,0.02,0.08,0.18,0.64,0.05,0.92,0.16,0.56,0.04,0.22,0.03,0.02,0.30,0.73,0.02)
) %>% 
  mutate(pval2 = paste("p = ",pval, sep = ""))

Fig 1. Correlation between growth markers and migration areas

Temporal trends of growth patterns are only weakly related to ERSST trends in migration areas of importance to salmon.

axislabs<-as_labeller(c("marine.incr" = "Marine increment","PS.incr" = "Post-smolt increment", 
                        "FS.incr" = "First summer increment", "FW.incr" = "First winter increment",
                        "n.labsea" = "North Labrador Sea", "nf" = "Newfoundland", 
                        "northatlantic" = "North Atlantic", "s.labsea" = "South Labrador Sea",
                        "wgreen" = "West Greenland"))

AICdf %>% 
  left_join(r2df, by = c("Marker", "Area")) %>% 
  left_join(pvdf, by = c("Marker", "Area")) %>% 
  dplyr::rename("responsevar" = "Marker", "predictvar" = "Area") %>% 
  left_join(modelinglong, by = c("responsevar","predictvar")) %>% 
  filter(responsevar == "PS.incr" | responsevar == "marine.incr" | responsevar == "FS.incr" | responsevar == "FW.incr") %>% 
  filter(predictvar != "gom") %>% 
  mutate(responsevar = fct_reorder(responsevar, meanmarker)) %>% 
  mutate(predictvar = fct_reorder(predictvar, desc(meanSST))) %>% 
  ggplot(aes(meanSST,meanmarker)) + geom_point() +
  geom_text(data = df,aes(x = x, y = y, label = pval2)) + 
  facet_grid(rows=vars(responsevar), cols = vars(predictvar), scales = "free", labeller = axislabs) + 
  geom_smooth(method = "lm") +  
  theme(panel.grid = element_blank()) + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 14)) + labs(x = "Annual mean SST C", y = "Growth marker annual mean", title = "Growth marker x SST") + 
  theme(strip.background = element_blank(), strip.text = element_text(size = 12, color = "black"))

#ggsave(filename = "growthxsst_full.pdf",plot = last_plot(), width = 11, height = 9.5)

sst_summarized_month<-sst %>% 
  filter(Dates >= "1963-01-01", 
         Dates < "1992-01-01") %>% 
  mutate(SmoltYear = year(Dates)) %>% 
  mutate(Month = month(Dates)) %>% 
  group_by(Month, SmoltYear,Area) %>% 
  dplyr::summarize(MonthlySST = mean(ERSST))

Select growth markers of interest
Calculate annual mean growth marker values by smolt year
Pairwise-combine all growth marker data and all SST data

growth_summarized_long<- Carlin_full_growth %>% 
  dplyr::select(JoinID, ReleaseYear, RecaptureYear,FS.incr, FW.incr, PS.incr, SS.incr, SW.incr, marine.incr) %>% 
  gather(key = "marker", value = "value", 4:9) %>% 
  group_by(ReleaseYear, marker) %>% 
  dplyr::summarize(meanmarker = mean(value, na.rm = TRUE)) %>% 
  dplyr::rename("SmoltYear" = "ReleaseYear")


modelinglong<- growth_summarized_long %>% 
  left_join(sst_summarized_month, by = "SmoltYear") %>% 
  dplyr::rename("responsevar" = "marker", "predictvar" = "Area", "meanSST" = "MonthlySST") %>% 
  filter(!is.na(predictvar))

Loop through all possible model combinations (singular pred var and no interactions)

responses<-as.vector(unique(as.character(modelinglong$responsevar)))
predictors<-as.vector(unique(as.character(modelinglong$predictvar)))
month<-as.character(seq(1,12,1))

AIClist<-list()
r2list<-list()
pvlist<-list()
bigAIC<-list()
bigr2<-list()
bigpv<-list()

for(m in month){
  
month_df<-modelinglong %>% 
    filter(Month == m)

for(i in responses){
  for(j in predictors){
    
temp<-month_df %>% 
  filter(responsevar == i & predictvar == j)

glancesummary<-glance(lm(meanmarker ~ meanSST, data = temp))

AIClist[[i]][j] <-glancesummary$AIC
r2list[[i]][j]<-glancesummary$r.squared
pvlist[[i]][j]<-glancesummary$p.value

}
}

bigAIC[[m]]<-do.call(rbind,AIClist)
bigr2[[m]]<-do.call(rbind,r2list)
bigpv[[m]]<-do.call(rbind,pvlist)

  }

bigAIC_df<-do.call(rbind,bigAIC)

monthno<-tibble(rep(seq(1,12,1),length(unique(rownames(bigAIC_df))))) %>% 
  dplyr::rename("month" = "rep(seq(1, 12, 1), length(unique(rownames(bigAIC_df))))") %>% 
  arrange(month)

finalAIC <- bigAIC_df %>% 
  as_tibble() %>% 
  mutate(Marker = rownames(bigAIC_df)) %>% 
  cbind(monthno) %>% 
  gather(key = "Area", value = "AIC", 1:6)


bigr2_df<-do.call(rbind,bigr2)

finalr2 <- bigr2_df %>% 
  as_tibble() %>% 
  mutate(Marker = rownames(bigr2_df)) %>% 
  cbind(monthno) %>% 
  gather(key = "Area", value = "Rsquared", 1:6)


bigpv_df<-do.call(rbind,bigpv)

finalpv <- bigpv_df %>% 
  as_tibble() %>% 
  mutate(Marker = rownames(bigpv_df)) %>% 
  cbind(monthno) %>% 
  gather(key = "Area", value = "p.value", 1:6)

model_summary_df<-finalAIC %>% 
  left_join(finalr2, by = c("Marker", "Area", "month")) %>% 
  left_join(finalpv, by = c("Marker", "Area", "month")) %>% 
  dplyr::rename("responsevar" = "Marker", "predictvar" = "Area") %>% 
  left_join(modelinglong, by = c("responsevar","predictvar"))

notinarea<-tribble(
  ~predictvar, ~month, ~presence,
  "gom", 1,"absent",
  "gom", 2,"absent",
  "gom", 7,"absent",
  "gom", 8,"absent",
  "gom", 9,"absent",
  "wgreen", 1,"absent",
  "wgreen", 2,"absent",
  "wgreen", 3,"absent",
  "wgreen", 4,"absent",
  "wgreen", 11,"absent",
  "wgreen", 12,"absent",
  "s.labsea", 7,"absent",
  "s.labsea", 8,"absent"
)


monthname<-tribble(
  ~Moname, ~month,
  "January",1,
  "February",2,
  "March",3,
  "April",4,
  "May",5,
  "June",6,
  "July",7,
  "August",8,
  "September",9,
  "October",10,
  "November",11,
  "December",12
)

Fig 2. R2 and p-value for correlation between migration areas and growth markers

model_summary_df %>% 
  mutate_at(6,round,2) %>% 
  filter(responsevar != "FS.circ", 
         responsevar != "FW.circ", 
         responsevar != "PS.circ", 
         responsevar != "SS.incr", 
         responsevar != "SW.incr", 
         predictvar != "gom") %>% 
  left_join(notinarea, by = c("predictvar","month")) %>% 
  mutate(Rsquared = ifelse(is.na(presence) == TRUE, Rsquared, NA)) %>% 
    mutate(responsevar = forcats::fct_relevel(responsevar, c("FS.incr", "FW.incr","PS.incr", "marine.incr"))) %>% 
    mutate(predictvar = forcats::fct_relevel(predictvar,"s.labsea","nf","n.labsea","northatlantic","wgreen")) %>% 
  mutate(addtext = ifelse(p.value <= 0.05, p.value,NA)) %>% 
  left_join(monthname, by = "month") %>% 
  mutate(Moname = fct_reorder(Moname,month)) %>% 
  ggplot(aes(responsevar,predictvar, label = addtext)) + geom_tile(aes(fill = Rsquared)) + viridis::scale_fill_viridis(na.value = "white") + facet_wrap(~Moname) + geom_text(color = "white", size = 3.4) + 
  theme(axis.text.y = element_text(size = 15), axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
        axis.title = element_text(size = 16), legend.text = element_text(size = 14), 
        legend.title = element_text(size = 15)) + labs(x = "Annual mean growth marker", y = "Migration area mean monthly SST") + theme(panel.grid = element_blank(), strip.background = element_blank(), strip.text = element_text(size = 14, color = "black"))

#ggsave(filename = "monthlycorrelations.pdf", plot = last_plot(), height = 9, width = 10)

SST x growth EDA

LGCarlson

6/27/2019

Fig 1. Correlation between growth markers and migration areas

Fig 2. R2 and p-value for correlation between migration areas and growth markers