## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----data---------------------------------------------------------------------

library(MRgrowth)

data(MRdata)
head(MRdata)
str(MRdata)

## ----format-------------------------------------------------------------------
formatted.data <- format.data(
  x        = MRdata,
  id.col   = "id",
  date.col = "date",
  size.col = "size",
  units    = "years",
  retain   = "max"
)

head(formatted.data)

## ----fit_fast-----------------------------------------------------------------
# Using low runs for vignette build speed
results <- calc.growth(
  x           = formatted.data,
  size0       = 50,
  age.vect    = 1:50,
  models      = c("vb", "gompertz", "logistic"),
  type        = "ML",
  runs        = 100,
  return.plot = TRUE
)

## ----params-------------------------------------------------------------------
results$parameters

## ----sizes--------------------------------------------------------------------
head(results$sizes.at.ages)

## -----------------------------------------------------------------------------
#look at model diagnostics
print(results$parameters)


# Extract parameters from the best-fitting model
best.model <- results$parameters[which(results$parameters$model=="vb"), ]


## ----size_at_age--------------------------------------------------------------
#calculate expected sizes for a 10, 20 and 30 year old turtle
#note that the model argument specifies that we are using a vb model

size.at.age(
  A     = best.model$A,
  k     = best.model$k,
  age   = c(10, 20, 30),
  size0 = 50,
  model = "vb")


## ----age_at_size--------------------------------------------------------------

#Estimate the age (years) of a 150, 250, and 350 mm turtle
#note that the model argument specifies that we are using a vb model

age.at.size(
  A     = best.model$A,
  k     = best.model$k,
  size  = c(150, 250, 350),
  size0 = 50,
  model = "vb")


## ----growth_funcs-------------------------------------------------------------
# Predict size of a 200 mm turtle after 3 years using the vb model
#note that unlike most functions in 'MRgrowth', each model as a separate function, rather than specifying the model in a single function

vb(
  A     = best.model$A,
  k     = best.model$k,
  size1 = 200,
  td    = 3
)

## -----------------------------------------------------------------------------

predicted.size2 <- vb(
  A     = best.model$A,
  k     = best.model$k,
  size1 = formatted.data$size1, #use size 1 data
  td    = formatted.data$td #use time difference in the data
)

plot(x = formatted.data$size2,
     y = predicted.size2,
     xlab = "Observed size2",
     ylab = "Estimated size2",
     pch = 16,
     col = "blue",
     bg  = "white",
     panel.first = rect(par("usr")[1], par("usr")[3],
                        par("usr")[2], par("usr")[4],
                        col = "white"))
abline(lm(predicted.size2 ~ formatted.data$size2),
       col = "black",
       lwd = 1)
r2 <- summary(lm(predicted.size2 ~ formatted.data$size2))$r.squared
legend("topleft",
       legend = bquote(R^2 == .(round(r2, 4))),
       bty    = "n")


## -----------------------------------------------------------------------------
#remove any individuals where size1 > A
observed.size2.b <- formatted.data$size2[which(formatted.data$size1 < best.model$A)]
predicted.size2.b <- predicted.size2[which(formatted.data$size1 < best.model$A)]

plot(x = observed.size2.b,
     y = predicted.size2.b,
     xlab = "Observed size2",
     ylab = "Estimated size2",
     pch = 16,
     col = "blue",
     bg  = "white",
     panel.first = rect(par("usr")[1], par("usr")[3],
                        par("usr")[2], par("usr")[4],
                        col = "white"))
abline(lm(predicted.size2.b ~ observed.size2.b),
       col = "black",
       lwd = 1)
r2 <- summary(lm(predicted.size2.b ~ observed.size2.b))$r.squared
legend("topleft",
       legend = bquote(R^2 == .(round(r2, 4))),
       bty    = "n")


## -----------------------------------------------------------------------------
predicted.size <- size.at.age(A=best.model$A,
                              k=best.model$k,
                              age=c(5,10,12,20),
                              size0=50,
                              model="vb")

#expected sizes given their ages
print(predicted.size)

#percent difference between observed and predicted sizes
print((c(200,320,330,400)-predicted.size)/c(200,320,330,400)*100)

## -----------------------------------------------------------------------------
logistic.model <- results$parameters[which(results$parameters$model=="logistic"), ]

predicted.size.logistic <- size.at.age(A=logistic.model$A,
                              k=logistic.model$k,
                              age=c(5,10,12,20),
                              size0=50,
                              model="logistic")

#predicted sizes based on the logistic model
print(predicted.size.logistic)

#percent difference between observed and predicted sizes for the logistic model
print((c(200,320,330,400)-predicted.size.logistic)/c(200,320,330,400)*100)


## -----------------------------------------------------------------------------
predicted.age <- age.at.size(A=best.model$A,
                              k=best.model$k,
                              size=c(200,320,330,400),
                              size0=50,
                              model="vb")

print(predicted.age)