Read in growth and release datasets.

Carlin_full_growth<-readRDS(paste(shared.path, "Salmon/Data/Growth/Data/Carlin_growth_fulldata.rds", sep = ""))

releasedata<- readRDS(paste(shared.path, "Salmon/Data/Mark-recap/Carlin_all_releases.rds", sep = ""))

Create “RecaptureDate” formatted as a date.

Carlin_full_growth$RecaptureDate<-do.call(paste,list(Carlin_full_growth$RecaptureYear, Carlin_full_growth$RecaptureMonthNum, Carlin_full_growth$RecaptureDay,sep="/"))

Carlin_full_growth$RecaptureDate<-as.Date(Carlin_full_growth$RecaptureDate,format="%Y/%m/%d") 

Create and format range of release dates (Date 1 and Date 2). Calculate # days between Release Date 1 and Release Date 2.

calcreleasedatesdf$ReleaseDate1<-do.call(paste,list(calcreleasedatesdf$ReleaseYear.y, calcreleasedatesdf$ReleaseMonth1, calcreleasedatesdf$ReleaseDay1,sep="/"))

calcreleasedatesdf$ReleaseDate2<-do.call(paste,list(calcreleasedatesdf$ReleaseYear.y, calcreleasedatesdf$ReleaseMonth2, calcreleasedatesdf$ReleaseDay2,sep="/"))

calcreleasedatesdf$ReleaseDate1<-as.Date(calcreleasedatesdf$ReleaseDate1,format="%Y/%m/%d") 
calcreleasedatesdf$ReleaseDate2<-as.Date(calcreleasedatesdf$ReleaseDate2,format="%Y/%m/%d") 

calcreleasedatesdf$ReleaseDate2[is.na(calcreleasedatesdf$ReleaseDay2)] = calcreleasedatesdf$ReleaseDate1[is.na(calcreleasedatesdf$ReleaseDay2)]

calcreleasedatesdf<-calcreleasedatesdf %>%
  dplyr::mutate(DateInt = lubridate::interval(lubridate::ymd(ReleaseDate1),lubridate::ymd(ReleaseDate2),tzone = "America/New_York")) %>%
  dplyr::mutate(RdayInt = lubridate::as.period(DateInt,unit = "days")) %>%
  dplyr::mutate(RdayInt = lubridate::day(RdayInt)) 

#Number of days between Release Day 1 and Release Date 2
#mean(calcreleasedatesdf$JdayInt,na.rm = TRUE) mean = 7.9
#min(calcreleasedatesdf$JdayInt,na.rm = TRUE) min = 0
#max(calcreleasedatesdf$JdayInt,na.rm = TRUE) max = 213

Write function to calculate midpoint of interval. Calculate midpoint of release interval using “int_midpoint” function. Add it to the dataframe as “EmigrationDate.”

int_midpoint <- function(interval) {
    int_start(interval) + (int_end(interval) - int_start(interval))/2
}

calcreleasedatesdf<-calcreleasedatesdf %>%
  mutate(EmigrationDate = do.call("c",lapply(calcreleasedatesdf$DateInt, int_midpoint))) %>%
  mutate(EmigrationDate = as.Date(EmigrationDate))

Add release dates dataframe to full growth dataframe. Remove fish with unknown recapture or emigration dates.

Carlin_simpledata<-Carlin_full_growth %>%
  left_join(calcreleasedatesdf,by="Batch_TagID") %>%
  dplyr::select(JoinID,Batch_TagID,RecaptureYear,RecaptureDate,SeaAge,RecaptureType,ReleaseYear,
                EmigrationDate,ends_with(".circ")) %>%
  filter(!is.na(RecaptureDate)) %>%
  filter(!is.na(EmigrationDate)) 

Add “SolsticeDateY1”, then convert to julian date (December 21). Convert emigration and recapture dates to julian date. Calculate the the number of days between emigration date and recapture date. Convert # days at sea to julian days and name JdayInt.

monthlygrowth<-Carlin_simpledata %>%
  dplyr::mutate(SolsticeDateY1=do.call(paste,list(ReleaseYear,"/12/21",sep=""))) %>%
  dplyr::mutate(SolsticeDateY1=as.Date(SolsticeDateY1,format="%Y/%m/%d")) %>%
  dplyr::mutate(SolsticeDOY=lubridate::yday(SolsticeDateY1)) %>%
  dplyr::mutate(EmigrationDOY=lubridate::yday(EmigrationDate)) %>%
  dplyr::mutate(RecaptureDOY=lubridate::yday(RecaptureDate)) %>%
  dplyr::mutate(DateInt = lubridate::interval(lubridate::ymd(EmigrationDate),lubridate::ymd(RecaptureDate),tzone = "America/New_York")) %>%
  dplyr::mutate(JdayInt = lubridate::as.period(DateInt,unit = "days")) %>%
  dplyr::mutate(JdayInt = lubridate::day(JdayInt)) %>%
  filter(JdayInt>365 & JdayInt < 1000)

Load function from github that write figure equations. Select a few fish to view curve between emigration, solstice, and recapture date.

Create a model for each fish (as shown in the example dataset above).

#square circuli for quadratic 
longmonthlygrowth$circ.num2<-(longmonthlygrowth$circ.num)^2

#remove fish with unknown circuli at time of recapture. 
longmonthlygrowth<-longmonthlygrowth %>%
  filter(!is.na(circ.num))

#loop through JoinIDs to create unique models for each individual 
JoinIDs<-as.vector(longmonthlygrowth$JoinID)
quadmodellist <- list()
for(i in JoinIDs){
  quadmodellist[[i]] <- lm(longmonthlygrowth$DOY[longmonthlygrowth$JoinID == i] ~ longmonthlygrowth$circ.num[longmonthlygrowth$JoinID == i] + longmonthlygrowth$circ.num2[longmonthlygrowth$JoinID == i])
}

Pull coefficients from model for each fish. Store in a list, then flatten the list to a df.

Use quadratic model coefficients to predict DOY estimates for each marine circulus. Use lm coefficients to predict marine circulus increment for each marine circulus. Therefore, DOY can be used to find marine circulus increment. Save to each individuals df to a list, then flatten to a df.

#create empty list
quadfittedlist<-list()

#loop through JoinIDs
JoinIDs<-as.vector(monthlygrowth$JoinID)

for(i in JoinIDs){
  #q is a sequence from 0 to the number of marine circuli each fish has 
         #(does not project beyond what is known)
q <- seq(from=0, to=max(monthlygrowth$marine.circ[monthlygrowth$JoinID == i]), by=1)

quadfittedlist[[i]]<-dplyr::as_tibble(matrix(NA,nrow = length(q),ncol = 2)) %>%
  mutate(fittedDOY=((as.numeric(unlist(quadmodellist[[i]][["coefficients"]][2]))*q) + 
        (as.numeric(unlist(quadmodellist[[i]][["coefficients"]][3]))*(q^2)) + 
        as.numeric(unlist(quadmodellist[[i]][["coefficients"]][1])))) %>%
  mutate(V1 = i) %>%
  mutate(V2 = q) %>% 
  mutate(pred.marine.inc = (marine.lm$coefficients[2]*q + marine.lm$coefficients[1]) + 
           Carlin_full_growth$smolt.incr[Carlin_full_growth$JoinID == i]) %>%
  rename("JoinID" = "V1","circulivalues" = "V2")
}
quad_modelpredictions = do.call(rbind, quadfittedlist)

Create a factor called “RecordID” that is a unique ID for each circulus (JoinID_circuli#) Create a df called “loopdf” that contains the emigration date for each fish. Join it to the model predictions dataframe so that each row now contains a “EmigrationDOY” column. Rename to “modelpredictions.”

quad_modelpredictions$recordID<-as.factor(paste(quad_modelpredictions$JoinID,quad_modelpredictions$circulivalues,sep = "."))

loopdf<-as.data.frame(matrix(NA,nrow = length(unique(quad_modelpredictions$JoinID)),ncol = 2))
loopdf$V1<-unique(quad_modelpredictions$JoinID)
loopdf$V2<-quad_modelpredictions$fittedDOY[quad_modelpredictions$circulivalues == 0]
colnames(loopdf)<-c("JoinID","EmigrationDOY")

modelpredictions<-quad_modelpredictions %>%
  full_join(loopdf,by="JoinID") %>%
  filter(!is.na(fittedDOY)) %>%
  filter(!is.na(EmigrationDOY))

Model predictions should align perfectly with the previous curve. Blue = original model, black = predicted values from model.

Group adjusted DOY into “growth periods” by julian date as in Izzo et al. 2017

#this loop also takes forever
#create blank column called "season"
modelpredictions$season = NA

#loop over record IDs
recordIDs<-as.vector(modelpredictions$recordID)
for(i in recordID){

#days before 213 occur in the Y1 summer
  if(modelpredictions$adjDOY[modelpredictions$recordID == i] < 213){
    modelpredictions$season[modelpredictions$recordID == i] = "Y1S"
  }
#days between 213 and 304 occur in the Y1 fall
  else if(modelpredictions$adjDOY[modelpredictions$recordID == i] > 212 & 
          modelpredictions$adjDOY[modelpredictions$recordID == i] <= 304){
    modelpredictions$season[modelpredictions$recordID == i] = "Y1F"
  }
#days between 305 and 356 occur in the Y1 winter
  else if(modelpredictions$adjDOY[modelpredictions$recordID == i] > 304 & 
          modelpredictions$adjDOY[modelpredictions$recordID == i] <= 356){
    modelpredictions$season[modelpredictions$recordID == i] = "Y1W"
  }
#days between 357 and 447 occur in the Y2 winter
  else if(modelpredictions$adjDOY[modelpredictions$recordID == i] > 356 & 
          modelpredictions$adjDOY[modelpredictions$recordID == i] <= 447){
    modelpredictions$season[modelpredictions$recordID == i] = "Y2W"
  }
#days between 448 and 638 occur in the Y2 summer
  else if(modelpredictions$adjDOY[modelpredictions$recordID == i] > 447 & 
          modelpredictions$adjDOY[modelpredictions$recordID == i] <= 638){
    modelpredictions$season[modelpredictions$recordID == i] = "Y2S"
  }
#anything larger than 638 will be Y2 fall
  else{
    modelpredictions$season[modelpredictions$recordID == i] = "Y2F"
  }
}

Load in circuli increment data.

Carlin_circuli<-read.csv(paste(shared.path,"Salmon/Data/Circuli/Data/Carlin_circuli.csv", 
                               sep = ""))

circuli_long<-Carlin_circuli %>%
  filter(!duplicated(JoinID)) %>%
  select(JoinID,Circ..,starts_with("c")) %>%
  gather(key = "circnum",value = "increment",3:126)

Loop to calculate increment between each circulus and use smoothing function.

spacinglist<-list()
JoinIDs<-as.vector(unique(circuli_long$JoinID))

for(i in JoinIDs){
ionly<-circuli_long %>%
  filter(JoinID == i) %>% 
  mutate(circuli = seq(1,124, by = 1)) %>%
  filter(!is.na(increment)) %>%
  filter(circnum != "c1") %>%
  select(JoinID, circuli,increment)

tempdf<-circuli_long %>% 
  filter(!is.na(increment)) %>%
  filter(JoinID == i)

incrementdiff <- as.data.frame(diff(tempdf$increment))
colnames(incrementdiff) <-c("incrementdiff")
smootheddiff <- as.data.frame(c(NA,NA,rollmean(incrementdiff$incrementdiff,5),NA,NA))
colnames(smootheddiff) <-c("smootheddiff")

spacinglist[[i]]<-cbind(ionly,incrementdiff,smootheddiff)
}

circuli_spacing = do.call(rbind, spacinglist)

Plot smoothed increment by circulus and total marine increment.

Loop through “JoinIDs” and “seasons” to calculate growth increment during each growth period. Store result in list and flatten.

full_spacing<-modelpredictions %>%
  left_join(circuli_spacing,by = "recordID") %>%
  rename("JoinID" = "JoinID.x")

JoinIDs<-as.vector(unique(full_spacing$JoinID))
seasons<-c("Y1S","Y1F","Y1W","Y2W","Y2S","Y2F")

incrementlist<-list()

for(i in JoinIDs){
  for(j in seasons){
temp<-full_spacing %>%
  filter(JoinID == i) %>%
  filter(season == j)

incrementlist[[i]][j] = (max(temp$increment,na.rm = TRUE) - min(temp$increment,na.rm = TRUE))
  }
}

increments_period = do.call(rbind, incrementlist)

Size of growth increment during each growth period.

Plot mean circulus increment by smolt year for each growth period. The increments represent: Y1S = before JD 213 - before Aug 1 Y1F = between 213-304 - until the end of Oct Y1W = between 305-356 - through winter solstice (Dec 21) Y2W = between 357-447 - between winter solstice and March 31 Y2S = between 448-638 - through end of September (Y2) Y2F = between 639-712 - through winter solstice

Calculate the number of days between each circuli using the “diff” function and looping through each fish, then saving the #days between to a dataframe that begins at circ 1 rather than circ 0 (diff is full - 1). Save as list and flatten to df.

dayscirclist<-list()
for(i in JoinIDs){
ionly<-modelpredictions %>%
  filter(JoinID == i) %>%
  filter(circulivalues != 0) 

tempdf<-modelpredictions %>% 
  filter(JoinID == i)

dayspercirc <- as.data.frame(abs(diff(tempdf$adjDOY)))
colnames(dayspercirc) <-c("dayspercirc")

dayscirclist[[i]]<-cbind(ionly,dayspercirc)
}

dayspercircfull = do.call(rbind, dayscirclist)

Plot intercirculi spacing by day at sea for “Early growth.” Defined “early growth” as circuli accumulated within 60 days of emigration.

## Warning: Removed 24 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing missing values (geom_path).

Write lm for early growth of each fish. Defined “early growth” as circuli accumulated within 60 days of emigration.

earlygrowthdf<-full_spacing %>%
  mutate(earlygrowth = (EmigrationDOY + 60)) %>%
  filter(adjDOY <= earlygrowth)

#loop through JoinIDs to create unique models for each individual 
JoinIDs<-as.vector(unique(earlygrowthdf$JoinID))
earlymodellist <- list()

for(i in JoinIDs){
  earlymodellist[[i]] <- lm(earlygrowthdf$incrementdiff[earlygrowthdf$JoinID == i] ~ earlygrowthdf$adjDOY[earlygrowthdf$JoinID == i])
}

Pull coefficients from model for each fish. Store in a list, then flatten the list to a df.

earlylistoftibs<-list()
for(i in JoinIDs){
earlylistoftibs[[i]]<-dplyr::as_tibble(matrix(NA,nrow = 1,ncol = 1)) %>%
  mutate(interceptterm = as.numeric(unlist(earlymodellist[[i]][["coefficients"]][1]))) %>%
  mutate(slopeterm = as.numeric(unlist(earlymodellist[[i]][["coefficients"]][2]))) %>%
  mutate(V1 = i)
}
early_modelcoefficients = do.call(rbind, earlylistoftibs)
colnames(early_modelcoefficients) <- c("JoinID","interceptterm","slopeterm") 

Create subset data by each of first 12 months. Sum the growth increment within each month.

#loop through JoinIDs to create unique models for each individual 
JoinIDs<-as.vector(unique(full_spacing$JoinID))

#DOY df
DOY<-as.data.frame(matrix(NA,nrow = 12,ncol = 1))

#create list to store in 
monthlylist <- list()
list2 <- list()


for(i in JoinIDs){
  for(j in 1:12){
tempdf<-full_spacing %>%
  filter(JoinID == i)
DOY[j,1]<-(tempdf$EmigrationDOY[tempdf$circulivalues == 0] + 30*j)
colnames(DOY)<-c("months")

monthlylist[[i]][j]<-tempdf %>%
  filter(adjDOY < DOY$months[j]) %>%
  select(incrementdiff)

list2[[i]][j]<-sum(unlist(monthlylist[[i]][j]),na.rm = TRUE)
  }
}

monthlydf = do.call(rbind, list2)
monthlydf<-as_tibble(monthlydf)
monthlydf$JoinID<-JoinIDs
colnames(monthlydf) <- c("mo01","mo02","mo03","mo04","mo05","mo06","mo07","mo08",
                                       "mo09","mo10","mo11","mo12","JoinID") 

Calculate monthly growth increment.

Visualize monthly growth increment for a random sample of 12 fish.

randomsample<-sample_n(monthlydflong,12) %>%
  select(JoinID)
plotdf<-monthlydflong %>%
  right_join(randomsample,by = "JoinID")

ggplot(plotdf,aes(month,increment,group = 1)) + geom_path() + facet_wrap(~JoinID) +
  scale_x_discrete(breaks=c("month01","month02","month03","month04","month05","month06","month07","month08",
                            "month09","month10","month11","month12"),
        labels=c("1","2","3","4","5","6","7","8","9","10","11","12")) + labs(x = "Months at Sea",
                            y = "Growth increment (mm)")

Calculate mean monthly growth increment by smolt year.

Reproduction of Fig 7 in “Monthly indices of the post-smolt growth of Atlantic salmon from the Drammen River, Norway” by McCarthy et al 2008 which found that mean circuli spacing for proportional month segments decreased over time (especially for early growth: months 1-6). The same does not appear to be true for these data.

ggplot(monthlygrowth_full,aes(ReleaseYear,mean)) + geom_point(aes(color = nbins)) + 
  facet_wrap(~ month, labeller = molabs) + 
  geom_smooth(method = "lm",se = FALSE, color = "black") +
  theme(axis.text.x = element_text(angle = 35,hjust = 1)) + 
  labs(x = "Smolt year", y = "Circuli spacing, normalized", title = "Months at sea") + 
  scale_color_manual(values = bluegradient) 

#+   geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width = 0.3) + ylim(-0.05,0.25)