Calculate monthly growth increments
LGCarlson
5/14/2019
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") Select only important columns and remove date NAs.
monthlygrowth<-Carlin_full_growth %>%
dplyr::select(JoinID,ReleaseYear,RecaptureYear,RecaptureDate,SeaAge,RecaptureType,
FMC,M1,M2,ends_with(".circ")) %>%
filter(!is.na(RecaptureDate))
monthlygrowth<-monthlygrowth %>%
dplyr::mutate(SolsticeDateY1=do.call(paste,list(ReleaseYear,"/12/21",sep=""))) %>%
dplyr::mutate(EmigrationDate=do.call(paste,list(ReleaseYear,"/5/7",sep=""))) %>%
dplyr::mutate(SolsticeDateY1=as.Date(SolsticeDateY1,format="%Y/%m/%d")) %>%
dplyr::mutate(EmigrationDate=as.Date(EmigrationDate,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)emigration<-monthlygrowth %>%
mutate(PeriodID="emigration") %>%
mutate(circ.num = 0) %>%
select(JoinID,PeriodID,EmigrationDOY,circ.num) %>%
rename("DOY" = "EmigrationDOY")
solstice<-monthlygrowth %>%
mutate(PeriodID="solstice") %>%
mutate(circ.num=(FS.circ + FW.circ)) %>%
select(JoinID,PeriodID,SolsticeDOY,circ.num) %>%
rename("DOY" = "SolsticeDOY")
recapture<-monthlygrowth %>%
mutate(PeriodID="recapture") %>%
select(JoinID,PeriodID,JdayInt,marine.circ) %>%
rename("DOY" = "JdayInt","circ.num" = "marine.circ")
longmonthlygrowth<-rbind(emigration,solstice,recapture)testdf<-longmonthlygrowth %>%
filter(JoinID == "1963_377" | JoinID == "1965_990" | JoinID == "1966_1444" | JoinID == "1966_1842")library(devtools)
source_gist("524eade46135f6348140")## Sourcing https://gist.githubusercontent.com/kdauria/524eade46135f6348140/raw/676acaca9a0a144ef320ae2ef00a31c3daa7179d/ggplot_smooth_func.R## SHA-1 hash of file is c0b163b9fd2d7fe7bd5541e3266d8d36ff3b895dggplot(testdf,aes(circ.num,DOY)) + geom_point(aes(color = PeriodID),size=3) +
stat_smooth(method = "lm", formula = y ~ x + I(x^2), size = 1) + facet_wrap(~JoinID) +
geom_point(aes(color = PeriodID),size=3) +
stat_smooth_func(geom="text",method="lm",hjust=0,parse=TRUE)## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
Create a model for each fish (as shown in the example dataset above).
longmonthlygrowth$circ.num2<-(longmonthlygrowth$circ.num)^2
JoinIDs<-as.vector(monthlygrowth$JoinID)
modellist <- c()
for(i in JoinIDs){
modellist[[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.
listoftibs<-list()
for(i in JoinIDs){
listoftibs[[i]]<-dplyr::as_tibble(matrix(NA,nrow = 1,ncol = 1)) %>%
mutate(interceptterm = as.numeric(unlist(modellist[[i]][["coefficients"]][1]))) %>%
mutate(slopeterm = as.numeric(unlist(modellist[[i]][["coefficients"]][2]))) %>%
mutate(quadterm = as.numeric(unlist(modellist[[i]][["coefficients"]][3]))) %>%
mutate(V1 = i)
}## Warning: `as_tibble.matrix()` requires a matrix with column names or a `.name_repair` argument. Using compatibility `.name_repair`.
## This warning is displayed once per session.modelcoefficients = do.call(rbind, listoftibs)
colnames(modelcoefficients) <- c("JoinID","interceptterm","slopeterm","quadterm") Reformat to long form for plotting.
longmodelling<-longmonthlygrowth %>%
left_join(modelcoefficients,by="JoinID")Create and check lm to use marine circulus number to predict marine increment for DOY.
marine.lm<-lm(Carlin_full_growth$marine.incr~Carlin_full_growth$marine.circ); summary(marine.lm)##
## Call:
## lm(formula = Carlin_full_growth$marine.incr ~ Carlin_full_growth$marine.circ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.97102 -0.15928 0.00031 0.15927 0.73753
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.0047422 0.0283816 -0.167
## Carlin_full_growth$marine.circ 0.0534032 0.0006249 85.461
## Pr(>|t|)
## (Intercept) 0.867
## Carlin_full_growth$marine.circ <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2476 on 1659 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8149, Adjusted R-squared: 0.8148
## F-statistic: 7304 on 1 and 1659 DF, p-value: < 0.00000000000000022ggplot(Carlin_full_growth,aes(marine.circ,marine.incr)) + geom_point(color='deepskyblue4') +
geom_abline(intercept = marine.lm$coefficients[1], slope = marine.lm$coefficients[2],size=1.2,
color = 'darkblue') + theme_bw() + labs(title = "",x = "Marine circuli",
y = "Marine increment (mm)") +
annotate("text",x=20,y=3.5,label = "y = 0.534x - 0.004")## Warning: Removed 1 rows containing missing values (geom_point).
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, therefore calculating monthly increments.
marine.lm<-lm(Carlin_full_growth$marine.incr~Carlin_full_growth$marine.circ)
fittedlist<-list()
for(i in JoinIDs){
q <- seq(from=0, to=max(monthlygrowth$marine.circ[monthlygrowth$JoinID == i]), by=1)
fittedlist[[i]]<-dplyr::as_tibble(matrix(NA,nrow = length(q),ncol = 2)) %>%
mutate(fittedDOY=((as.numeric(unlist(modellist[[i]][["coefficients"]][2]))*q) +
(as.numeric(unlist(modellist[[i]][["coefficients"]][3]))*(q^2)) +
as.numeric(unlist(modellist[[i]][["coefficients"]][1])))) %>%
mutate(V1 = i) %>%
mutate(V2 = q) %>%
mutate(pred.marine.inc = (marine.lm$coefficients[2]*q + marine.lm$coefficients[1])) %>%
rename("JoinID" = "V1","circulivalues" = "V2")
}
modelpredictions = do.call(rbind, fittedlist)hist(modelpredictions$fittedDOY)
hist(modelpredictions$circulivalues)
write.csv(modelpredictions,"preds.csv")testdf<-longmonthlygrowth %>%
filter(JoinID=="1968_2201" | JoinID=="1969_8327" | JoinID=="1970_38846" | JoinID=="1971_87680")library(devtools)
source_gist("524eade46135f6348140")## Sourcing https://gist.githubusercontent.com/kdauria/524eade46135f6348140/raw/676acaca9a0a144ef320ae2ef00a31c3daa7179d/ggplot_smooth_func.R## SHA-1 hash of file is c0b163b9fd2d7fe7bd5541e3266d8d36ff3b895dggplot(testdf,aes(circ.num,DOY)) + geom_point(aes(color = PeriodID),size=3) +
stat_smooth(method = "lm", formula = y ~ x + I(x^2), size = 1) + facet_wrap(~JoinID) +
geom_point(aes(color = PeriodID),size=3) +
stat_smooth_func(geom="text",method="lm",hjust=0,parse=TRUE)## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning in predict.lm(model, newdata = data.frame(x = xseq), se.fit = se, :
## prediction from a rank-deficient fit may be misleading## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).