R/fit.R

In ggcostoya/limon: Evolutionary Constraints among Thermal Performance Traits

### FUNCTIONS TO FIT DATA

### * Fit TPD

#' Fit TPD
#'
#' Fits a thermal-performance dataset (TPD) using a the best (lowest AIC scored) thermal-performance model (TPM) through a Non-linear least squares (nls) approach
#'
#' @param tpd A thermal-performance dataset (TPD)
#'
#' @return Model results in tibble format divided in 3 columns: AIC score, TPM name and model results (nested)
#'
#' @examples
#' tpd <- gen_base(topt = 30, tb = 5, skw = -2, ctmin = 15, ctmax = 35, pmax = 10, pmin = 0.1)
#' fit_tpd(tpd)
#' fit_tpd(tpd) %>% select(results) %>% unnest()
#'
#' @export

fit_tpd <- function(tpd){

  # run model fits
  fit_flinn <- fit_flinn(tpd)
  fit_gaussian <- fit_gaussian(tpd)
  fit_spain <- fit_spain(tpd)
  fit_weibull <- fit_weibull(tpd)

  # overcome errors caused by not being able to calculate AIC
  aic_flinn <- ifelse(is.null(fit_flinn), NA, AIC(fit_flinn))
  aic_gaussian <- ifelse(is.null(fit_gaussian), NA, AIC(fit_gaussian))
  aic_spain <- ifelse(is.null(fit_spain), NA, AIC(fit_spain))
  aic_weibull <- ifelse(is.null(fit_weibull),NA, AIC(fit_weibull))

  # determine the best model
  best_tpm <- suppressWarnings(tibble(aic = c(aic_flinn,
                                              aic_gaussian,
                                              aic_spain,
                                              aic_weibull), # get aic scores
                                      tpm = c("flinn",
                                              "gaussian",
                                              "spain",
                                              "weibull"), # get tpm names
                                      results = c(nest(tidy(fit_flinn)),
                                                  nest(tidy(fit_gaussian)),
                                                  nest(tidy(fit_spain)),
                                                  nest(tidy(fit_weibull))))) %>% # # get model results
    filter(aic == min(aic,na.rm = TRUE)) %>% unnest(results) # select best model result

  return(best_tpm)
}

### * Fit TPD using Gaussian TPM

#' Fit a TPD using a Gaussian TPM
#'
#' Fits a thermal-performance dataset (TPD) using a Gaussian thermal-performance model (TPM) through a Non-linear least squares (nls) approach.
#'
#' @param tpd A thermal-performance dataset (TPD)
#'
#' @return An nls model object
#'
#' @examples
#' tpd <- gen_base_tpd(topt = 30, tb = 5, skw = -2, ctmin = 15, ctmax = 35, pmax = 10, pmin = 0.1)
#' fit_gaussian(tpd)
#'
#' @export

fit_gaussian <- function(tpd){

  # set starting values
  starting_values <- c(s = max(tpd$p, na.rm = TRUE), # ~ pmax
                       a = mean(tpd$t[tpd$p == max(tpd$p, na.rm = TRUE)]), # ~ topt
                       b = (max(tpd$t,na.rm = TRUE) - min(tpd$t, na.rm = TRUE))/2) # ~ ctmax - ctmin

  # set lower limit values
  lower_limit_values <- c(s = min(tpd$p, na.rm = TRUE), # ~ pmin
                          a = min(tpd$t, na.rm = TRUE), # ~ ctmin
                          b = 0) # following rTPC instructions

  # set upper limit values
  upper_limit_values <- c(s = max(tpd$p, na.rm = TRUE) * 10, # ~ pmax * 10
                          a = max(tpd$t, na.rm = TRUE), # ~ ctmax
                          b = (max(tpd$t, na.rm = TRUE) - min(tpd$t, na.rm = TRUE))*10) # ~ (ctmax - ctmin)*10

  # fitting model using nlsLM
  fit <- nls.multstart::nls_multstart(formula = p ~ gaussian(t = t, s, a, b),
                                      data = tpd,
                                      iter = 500,
                                      start_lower = starting_values -1,
                                      start_upper = starting_values + 1,
                                      lower = lower_limit_values,
                                      upper = upper_limit_values,
                                      supp_errors = "Y")
  return(fit)
}

### * Fit TPD using a Weibull TPM

#' Fit a TPD using a Weibull TPM
#'
#' Fits a thermal-performance dataset (TPD) using a Gaussian thermal-performance model (TPM) through a Non-linear least squares (NLS) approach.
#'
#' @param tpd A thermal-performance dataset (TPD)
#'
#' @return An nls model object
#'
#' @examples
#' tpd <- gen_base_tpd(topt = 30, tb = 5, skw = -2, ctmin = 15, ctmax = 35, pmax = 10, pmin = 0.1)
#' fit_weibull(tpd)
#'
#' @export

fit_weibull <- function(tpd){

  # set starting values
  starting_values <- c(s = max(tpd$p, na.rm = T), # ~ pmax
                       a = mean(tpd$t[tpd$p == max(tpd$p, na.rm = T)]), # ~ topt
                       b = max(tpd$t,na.rm = T) - min(tpd$t, na.rm = T), # ~ ctmax - ctmin
                       c = 4 ) # following by rTPC

  # set lower limit values
  lower_limit_values <- c(s = 1, # recommended by rTPC
                          a = min(tpd$t, na.rm = T), # ~ ctmin
                          b = 0.1, c = 0.1) # following rTPC instructions

  # set upper limit values
  upper_limit_values <- c(s = 100, # rare to find something higher than 100.
                          a = max(tpd$t, na.rm = TRUE), # ~ ctmax
                          b = Inf, c = Inf) # following rTPC instructions

  # fitting model using nlsLM
  fit <- nls.multstart::nls_multstart(formula = p ~ weibull(t = t, s, a, b, c),
                                      data = tpd,
                                      iter = 500,
                                      start_lower = starting_values -1,
                                      start_upper = starting_values + 1,
                                      lower = lower_limit_values,
                                      upper = upper_limit_values,
                                      supp_errors = "Y")

  return(fit)

}

### * Fit TPD using a Flinn TPM

#' Fit a TPD using a Flinn TPM
#'
#' Fits a thermal-performance dataset (TPD) using a Flinn thermal-performance model (TPM) through a Non-linear least squares (NLS) approach.
#'
#' @param tpd A thermal-performance dataset (TPD)
#'
#' @return An nls model object
#'
#' @examples
#' tpd <- gen_base_tpd(topt = 30, tb = 5, skw = -2, ctmin = 15, ctmax = 35, pmax = 10, pmin = 0.1)
#' fit_flinn(tpd)
#'
#' @export

fit_flinn <- function(tpd){

  # set starting values
  b = (-2*0.005*tpd$t[tpd$p == max(tpd$p)])[1]
  a = -min(b*tpd$t +0.005*(tpd$t^2))
  starting_values <- c(b = b, a = a, c = 1) # following rTPC instructions

  # set lower limit values
  lower_limit_values <- c(b = -100, a = -100, c = 0) # following rTPC instructions

  # set upper limit values
  upper_limit_values <- c(b = 100, a = 100, c = 10) # following rTPC instructions

  # fitting model using nlsLM
  fit <- nls.multstart::nls_multstart(formula = p ~ flinn(t = t, a, b, c),
                                      data = tpd,
                                      iter = 500,
                                      start_lower = starting_values -1,
                                      start_upper = starting_values + 1,
                                      lower = lower_limit_values,
                                      upper = upper_limit_values,
                                      supp_errors = "Y")
  return(fit)
}

### * Fit TPD using a Spain TPM

#' Fit a TPD using a Spain TPM
#'
#' Fits a thermal-performance dataset (TPD) using a Spain thermal-performance model (TPM) through a Non-linear least squares (NLS) approach.
#'
#' @param tpd A thermal-performance dataset (TPD)
#'
#' @return An nls model object
#'
#' @examples
#' tpd <- gen_base(topt = 30, tb = 5, skw = -2, ctmin = 15, ctmax = 35, pmax = 10, pmin = 0.1)
#' fit_spain(tpd)
#'
#' @export

fit_spain <- function(tpd){

  #set starting values
  starting_values <- c(s = min(tpd$p,na.rm = T), a = 0, b = 0, c = 0)

  # set lower limit values
  lower_limit_values <- c(s = abs(min(tpd$p,na.rm = T))*-10, a = -20, b = -20, c = -20)

  # set upper limit values
  upper_limit_values <- c(s = max(tpd$p,na.rm = T),a = 20, b = 20, c = 20)

  # fitting the model using nls_multstart
  fit <- nls.multstart::nls_multstart(formula = p ~ spain(t = t, s, a, b, c),
                                      data = tpd,
                                      iter = 500,
                                      start_lower = starting_values -1,
                                      start_upper = starting_values + 1,
                                      lower = lower_limit_values,
                                      upper = upper_limit_values,
                                      supp_errors = "Y")
  return(fit)


}