---
title: "Script 4 Environmental cues for migration"
author: "Michio Fukushima and Pete Rand"
date: "`r format(Sys.time(),'%Y/%m/%d')`"
output: html_document
---

```{r setup, include=F}
library(kableExtra)
library(data.table)
library(directlabels)
knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)

# Function to change time zone from UTC to Asia/Tokyo
Tz <- function(x) as.POSIXct(x) %>% force_tz(tzone = "Asia/Tokyo")
```

NOTE: Run Script 1 before running this script.

# Data import
```{r Env data import, eval=TRUE, include=FALSE}
env18 <- read_csv("Karibetsu_Water_2018.csv") # Water temperature and water level in 2018
env19 <- read_csv("Karibetsu_Water_2019.csv") # Water temperature and water level in 2019
env18$Time <- Tz(env18$Time)
env19$Time <- Tz(env19$Time)

env18 <- env18 %>% rename(Tmp=Temp,Lev=Level.m) %>% add_column(.before="Time",SN=seq(nrow(env18)))
env19 <- env19 %>% rename(Tmp=Temp,Lev=Level.m) %>% add_column(.before="Time",SN=seq(nrow(env19)))
# Tmp: Water temperature recorded at Site L (Celcius)
# Lev: Water level recorded at Site L (m)
```

# Sliding window analysis
```{r Sliding window of X-hr lag and Y-hr-forward duration for Water Temperature, include=TRUE, echo=TRUE}
# WARNING: the following calculation takes several minutes or more!
ma <- function(x, X, Y){stats::filter(lag(x,6*(X - Y)),rep(1/(6*Y),6*Y),sides=1) %>% round(3)} # Moving average
xhrs <- yhrs <- seq(8,240,by=4)  # With "by=4"", it takes about 10 minutes for the loop below.
env18T <- env18 %>% select(Time,Tmp)
env19T <- env19 %>% select(Time,Tmp)

for(i in seq_along(xhrs)){
  for(j in seq_along(yhrs)){
    if(j <= i){
      nam18 <- nam19 <- paste0("T",xhrs[i],".",yhrs[j])
      assign(nam18, ma(env18T$Tmp,xhrs[i],yhrs[j])) 
      env18T <- env18T %>% mutate(!!nam18 := get(nam18))

      assign(nam19, ma(env19T$Tmp,xhrs[i],yhrs[j])) 
      env19T <- env19T %>% mutate(!!nam19 := get(nam19))
    }
  }
}
# For dynamically naming vars above, refer to: 
# https://intellipaat.com/community/11413/dplyr-mutate-use-dynamic-variable-names
```

```{r Same for Water Level, include=TRUE, echo=TRUE}
env18L <- env18 %>% select(Time,Lev)
env19L <- env19 %>% select(Time,Lev)

for(i in seq_along(xhrs)){ # "xhrs" and "yhrs" as defined above
  for(j in seq_along(yhrs)){
    if(j <= i){
      nam18 <- nam19 <- paste0("L",xhrs[i],".",yhrs[j])
      assign(nam18, ma(env18L$Lev,xhrs[i],yhrs[j])) 
      env18L <- env18L %>% mutate(!!nam18 := get(nam18))

      assign(nam19, ma(env19L$Lev,xhrs[i],yhrs[j])) 
      env19L <- env19L %>% mutate(!!nam19 := get(nam19))
    }
  }
}
```

```{r Fish data, echo=TRUE}
joined <-  joined %>% mutate(U.A18=Tz(TimeU.A18),U.T18=Tz(TimeU.T18),D.T18=Tz(TimeD.T18),D.A18=Tz(TimeD.A18),
                             U.A19=Tz(TimeU.A19),U.T19=Tz(TimeU.T19),D.T19=Tz(TimeD.T19),D.A19=Tz(TimeD.A19)) %>% 
           select(-TimeU.A18,-TimeU.T18,-TimeD.T18,-TimeD.A18,-TimeU.A19,-TimeU.T19,-TimeD.T19,-TimeD.A19)

pick <- is.na(joined$U.T18)|is.na(joined$U.T19)
joined18 <- joined %>% filter(!pick) %>% select(PIT,Sex,Utime=U.T18) # Tributary Arrival Timing is extracted. 
joined19 <- joined %>% filter(!pick) %>% select(PIT,Sex,Utime=U.T19)
```

# Arrival timings and environmental cues
Below, fish arrival timings at tributary sites (Utime) are linked to data for mean water temperature and water level within corresponding sliding windows during both 2018 and 2019 spawning seasons.  

## Data preparation
```{r Merging migration timings with mean environmental signals for each sliding window, echo=TRUE}
setDT(joined18)
setDT(joined19)
setDT(env18T)
setDT(env19T)
setDT(env18L)
setDT(env19L)

joined18[,time := Utime]
joined19[,time := Utime]
env18T[,  time := Time]
env19T[,  time := Time]
env18L[,  time := Time]
env19L[,  time := Time]

T8  <- env18T[joined18, on = "time", roll = "nearest"] %>% arrange(PIT) %>% select(-Tmp,-time,-Time,-Utime) 
T9  <- env19T[joined19, on = "time", roll = "nearest"] %>% arrange(PIT) %>% select(-Tmp,-time,-Time,-Utime)
L8  <- env18L[joined18, on = "time", roll = "nearest"] %>% arrange(PIT) %>% select(-Lev,-time,-Time,-Utime) 
L9  <- env19L[joined19, on = "time", roll = "nearest"] %>% arrange(PIT) %>% select(-Lev,-time,-Time,-Utime)
```

## Correlations between pairs of all sliding windows  
Function "cor.coef()" is applied to calculate the time window-specific Pearson correlation coefficients of both water temperature and water level experienced by females and males separately.

```{r Sliding window implementation, echo=TRUE}
cor.coef <- function(df1,df2,sex) {# calculates correlations between all pairs of sliding windows
  df1. <- df1 %>% filter(Sex == `sex`) %>% select(-PIT,-Sex)
  df2. <- df2 %>% filter(Sex == `sex`) %>% select(-PIT,-Sex)
  sapply(1:ncol(df1.), function(x) cor(df1.[,x,with=FALSE], df2.[,x,with=FALSE], use="compl")) %>% set_names(names(df1.))
  } # See: http://datatable.r-forge.r-project.org/datatable-faq.pdf with regard to "with=FALSE" above.

R.T.F <- cor.coef(df1=T8,df2=T9,sex="Female") # water temp
R.T.M <- cor.coef(df1=T8,df2=T9,sex="Male")
R.L.F <- cor.coef(df1=L8,df2=L9,sex="Female") # water level
R.L.M <- cor.coef(df1=L8,df2=L9,sex="Male")
```

The resulting correlations are displayed using heatmaping against time lag and duration.  

```{r Heatmapping correlation coefs for water temerature, echo=TRUE}
trans <- function(vec) str_split(str_sub(names(vec),2,-1),"\\.",simplify = TRUE) %>% as.character() %>% as.numeric() %>% 
  matrix(ncol=2) %>% as.data.frame() %>% add_column(vec) %>% set_names(c("X","Y","Coef")) # See p.174 of Wickham's book for "\\."

R.T <- bind_rows(add_column(trans(R.T.F),Sex="Female"), add_column(trans(R.T.M),Sex="Male"))
brks <- c(seq(0.1,0.6,0.1), 0.65)

p <- ggplot(R.T, aes(X,Y,z=Coef)) +
  geom_raster(aes(fill=Coef)) +
  scale_x_reverse() + 
  scale_fill_distiller(palette = "Spectral", direction = -1) +
  geom_contour(aes(colour = ..level..), breaks=brks) +
  labs(x = "Time lag before PIT detection (hours)", y = "Duration (hours)") +
  coord_fixed() +
  facet_wrap(~Sex)

direct.label(p, list("far.from.others.borders", colour='black'))
```

```{r Heatmapping correlation coefs for water level, echo=TRUE}
R.L <- bind_rows(add_column(trans(R.L.F),Sex="Female"), add_column(trans(R.L.M),Sex="Male")) 
brks <- c(0.1, 0.2, 0.3, 0.4, 0.45, 0.5, 0.55, 0.60)

p <- ggplot(R.L, aes(X,Y,z=Coef)) + 
  geom_raster(aes(fill=Coef)) +
  scale_x_reverse() + 
  scale_fill_distiller(palette = "Spectral", direction = -1) +
  geom_contour(aes(colour = ..level..), breaks=brks) +
  labs(x = "Time lag before PIT detection (hours)", y = "Duration (hours)") +
  coord_fixed() +
  facet_wrap(~Sex)

direct.label(p, list("far.from.others.borders", colour='black'))
```

## Maximum correlations for the enviromental cues and their significance   
```{r Maximum Correlations}
nmTF <- names(R.T.F)[R.T.F == max(R.T.F)] # (T200.40)  for females
nmTM <- names(R.T.M)[R.T.M == max(R.T.M)] # (T124.12)  for males
nmLF <- names(R.L.F)[R.L.F == max(R.L.F)] # (L224.144) for females
nmLM <- names(R.L.M)[R.L.M == max(R.L.M)] # (L216.32)  for males
```

The significance of these correlations were tested using randomization (van de Pol 2016). 

# Randomization test of correlations in water temp and level between 2018 and 2019
## Initialization 
```{r data and function, echo=TRUE}
T8f <- T8 %>% filter(Sex == "Female") %>% select(-PIT,-Sex) %>% as.data.frame()
T9f <- T9 %>% filter(Sex == "Female") %>% select(-PIT,-Sex) %>% as.data.frame()
T8m <- T8 %>% filter(Sex == "Male")   %>% select(-PIT,-Sex) %>% as.data.frame()
T9m <- T9 %>% filter(Sex == "Male")   %>% select(-PIT,-Sex) %>% as.data.frame()

L8f <- L8 %>% filter(Sex == "Female") %>% select(-PIT,-Sex) %>% as.data.frame()
L9f <- L9 %>% filter(Sex == "Female") %>% select(-PIT,-Sex) %>% as.data.frame()
L8m <- L8 %>% filter(Sex == "Male")   %>% select(-PIT,-Sex) %>% as.data.frame()
L9m <- L9 %>% filter(Sex == "Male")   %>% select(-PIT,-Sex) %>% as.data.frame()

# Randomization function
cor.coef.rand <- function(df1,df2,n) {# this func works best with data.frame (tibble is slow)
  Summary.cor <- data.frame(
    matrix(nrow=10,ncol=6,
           dimnames=list(c(),c("Min","1st Qu.","Median","Mean","3rd Qu.","Max"))))
  for(i in seq(n)){
    cors <- sapply(1:ncol(df1), function(x) cor(df1[,x], df2[sample(nrow(df2)),x], use="compl"))
    Summary.cor[i,] <- summary(cors)
    cat(i,"\n")
  }
  Summary.cor
} 
```

## Execute randomization
```{r Randomizations of Temperature and Water Level for Females and Males, echo=TRUE}
set.seed(100)

RTF.rand <- cor.coef.rand(df1=T8f,df2=T9f,n=999);   cat("End Temperature for females","\n")
RTM.rand <- cor.coef.rand(df1=T8m,df2=T9m,n=999);   cat("End Temperature for males","\n")
RLF.rand <- cor.coef.rand(df1=L8f,df2=L9f,n=999);   cat("End Water Level for females","\n")
RLM.rand <- cor.coef.rand(df1=L8m,df2=L9m,n=999);   cat("End Water Level for males","\n")
```

## P values for the test
```{r P values from randomization tests, echo=TRUE}
P.Tf <- sum(RTF.rand$Max >= max(R.T.F))/(1 + nrow(RTF.rand)) # P.Tf = 0.001 !
P.Tm <- sum(RTM.rand$Max >= max(R.T.M))/(1 + nrow(RTM.rand)) # P.Tm = 0.999
P.Lf <- sum(RLF.rand$Max >= max(R.L.F))/(1 + nrow(RLF.rand)) # P.Lf = 0.014 !
P.Lm <- sum(RLM.rand$Max >= max(R.L.M))/(1 + nrow(RLM.rand)) # P.Lm = 0.999
```

```{r Summary table of the critical time windows}
Tab1 <- data.frame(
  Var = c("Water temp", "Water temp", "Water level", "Water level"),
  Sex = c("Female", "Male", "Female", "Male"),
  Window=c(nmTF,nmTM,nmLF,nmLM),
  'Cor coef'=c(max(R.T.F),max(R.T.M),max(R.L.F),max(R.L.M)) %>% round(.,3),
  P = c(P.Tf,P.Tm,P.Lf,P.Lm)
)
  
knitr::kable(Tab1,align="r",caption = "Table X. Critical windows for temperature and water level, their between-year correlations, and statistical significance.") %>% 
  kable_styling(full_width = F, position="left") %>% 
  column_spec(1:5,width="2.5cm")
```


```{r Critical time windows and approximate time of departure from the river mouth}
Windows <- data.frame(window=c("Water level","Water temperature"),start=c(-9.3,-8.3),end=c(-3.3,-6.7))
# Migration distance = 24.2 km, Female up-migration speed = 1.98 km/d (95% CI = 0.68 - 5.79 km/d)
# (see 'Migration Speed' Section of "Script 2")

ggplot(Windows,aes(x=start,y=c(1,2)))+
  annotate("rect",xmin=-24.2/0.68,xmax=-24.2/5.79,ymin=.2,ymax=2.8,fill="grey90")+ #
  geom_vline(xintercept = -24.2/1.98,col="black",linetype="dashed")+  # 12.2 days of travel time
  geom_hline(yintercept = c(1,2),col="grey60",linetype="dotted")+
  geom_segment(aes(xend=end,yend=c(1,2)),size=4,col="black")+
  labs(x="Days before arrival",y="")+
  scale_x_continuous(breaks=seq(-40,0,5))+
  scale_y_continuous(breaks=c(1,2),label=c("Water level","Water temperature"))+
  coord_cartesian(xlim=c(-40,0),ylim=c(.4,2.6))+
  theme(aspect.ratio = .5,
        axis.title.x = element_text(margin = margin(t = 10)))
```







