R/fit.R

In ggcostoya/momentum:

### * FUNCTION FIT & AIC ESTIMATOR FUNCTIONS

### * FIT TPD

#' Fit a TPD using the lowest AIC scoring model
#'
#' @param tpd A thermal performance dataset with "t" and "p" as temperature and performance.
#'
#' @return A data frame with the model used and paramter value estimations
#'
#' @examples
#'
#' @export

fit_tpd <- function(tpd){

  # AIC of the fit using all data
  aic <- c(AIC(fit_gaussian(tpd)), AIC(fit_emg(tpd)), AIC(fit_weibull(tpd)))
  model <- c("gaussian", "emg", "weibull")
  comparison <- data.frame(aic,model)

  # Select best model
  best <- comparison %>% filter(aic == min(aic))

  # Get the best fit and preparte return dataset
  if(best$model == "gaussian"){
    coefs <- data.frame(as.list(coef(fit_gaussian(tpd))))
    model <- "gaussian"
    fit <- cbind(model,coefs)}
  if(best$model == "emg"){
    coefs <- data.frame(as.list(coef(fit_emg(tpd))))
    model <- "emg"
    fit <- cbind(model,coefs)}
  if(best$model == "weibull"){
    coefs <- data.frame(as.list(coef(fit_weibull(tpd))))
    model <- "weibull"
    fit <- cbind(model,coefs)}

  return(fit)

}

### * FIT & AIC GAUSSIAN

#' Fit TPD using Gaussian function and get AIC
#'
#' @param tpd A thermal performance dataset with "t" and "p" as temperature and performance.
#'
#' @return A small table with parameter estimates and AIC score
#'
#' @examples
#'
#' @export

fit_gaussian <- function(d){

  # Set starting values for param estimation depending on the carachteristics of the data
  s_s <- max(d$p, na.rm = T) # Starting value for the scale parameter is Pmax
  s_a <- mean(d$t[d$p == max(d$p,na.rm = TRUE)]) #bStarting value for a paramter is Topt
  s_b <- max(d$t,na.rm = TRUE) - min(d$t, na.rm = TRUE) # Starting value for b paramter is Tbr
  starting_values <- c(s = s_s, a = s_a, b = s_b)

  # Set lower limit values for param estimation.
  low_s <- 1 # Lower limit is 1, there cannot be negative TPCs
  low_a <- min(d$t, na.rm = T) # Lower limit is CTmin
  low_b <- 0.1 # Should be zero but just to keep it positive
  lower_values <- c(low_s, low_a, low_b)

  # Set upper limit values for param estimation
  up_s <- 100 # Rare to find something higher than this in TPCs
  up_a <- max(d$t, na.rm = T) # Upper limit is CTmax
  up_b <- (max(d$t,na.rm = TRUE) - min(d$t, na.rm = TRUE))*10
  upper_values <- c(up_s, up_a, up_b)

  # Fitting model using non-linear least squares
  fit <- nlsLM(p ~ pred_gaussian(x=t, s, a, b),
             data = d,
             start = starting_values,
             lower =  lower_values,
             upper = upper_values,
             algorithm = "port")

  return(fit)
}

### * FIT EXPONENTIALLY MODIFIED GAUSSIAN

#' Fit TPD using the EMG function and get AIC
#'
#' @param d A thermal performance dataset with "t" and "p" as temperature and performance.
#'
#' @return An nls-format model object
#'
#' @examples
#'
#' @export

fit_emg <- function(d){

  # Set starting values for param estimation depending on the carachteristics of the data
  s_s <- max(d$p, na.rm = T) # Starting value for the scale parameter is Pmax
  s_a <- mean(d$t[d$p == max(d$p,na.rm = TRUE)]) #bStarting value for a paramter is Topt
  s_b <- (max(d$t,na.rm = TRUE) - min(d$t, na.rm = TRUE))/2 # Starting value for b paramter is Tbr
  s_c <- 2 # Set as recommended by rTPD package
  starting_values <- c(s = s_s, a = s_a, b = s_b, c = s_c)

  # Set lower limit values for param estimation.
  low_s <- 1 # Lower limit is 1, there cannot be negative TPCs
  low_a <- min(d$t, na.rm = T) # Lower limit is CTmin
  low_b <- 0.1 # Should be zero but just to keep it positive
  low_c <- 0.01 # Should be zero but jist to keep it positive
  lower_values <- c(low_s, low_a, low_b, low_c)

  # Set upper limit values for param estimation
  up_s <- 100 # Rare to find something higher than this in TPCs
  up_a <- max(d$t, na.rm = T) # Upper limit is CTmax
  up_b <- (max(d$t,na.rm = TRUE) - min(d$t, na.rm = TRUE))*10
  up_c <- 10
  upper_values <- c(up_s, up_a, up_b, up_c)

  # Fitting model using non-linear least squares
  fit <- nlsLM(p ~ pred_emg(x=t, s, a, b, c),
             data = d,
             start = starting_values,
             lower =  lower_values,
             upper = upper_values,
             algorithm = "port")

  return(fit)
}


### * FIT WEIBULL

#' Fit Weibull
#'
#' @description Fits a Non-linear Least Squares(nls) Model using the Weibull function. Extracted from rTPC and Weibull 1995
#'
#' @param d A thermal performance dataset with "t" and "p" as temperature and performance.
#'
#' @return An nls-format model object
#'
#' @examples
#'
#' @export

fit_weibull <- function(d){

  # Set starting values for param estimation depending on the carachteristics of the data
  s_s <- max(d$p, na.rm = T) # Starting value for the scale parameter is Pmax
  s_a <- mean(d$t[d$p == max(d$p,na.rm = TRUE)]) #bStarting value for a paramter is Topt
  s_b <- max(d$t,na.rm = TRUE) - min(d$t, na.rm = TRUE) # Starting value for b paramter is Tbr
  s_c <- 4 # Set as recommended by rTPD package
  starting_values <- c(s = s_s, a = s_a, b = s_b, c = s_c)

  # Set lower limit values for param estimation.
  low_s <- 1 # Lower limit is 1, there cannot be negative TPCs
  low_a <- min(d$t, na.rm = T) # Lower limit is CTmin
  low_b <- 0.1 # Should be zero but just to keep it positive
  low_c <- 0.1 # Should be zero but jist to keep it positive
  lower_values <- c(low_s, low_a, low_b, low_c)

  # Set upper limit values for param estimation
  up_s <- 100 # Rare to find something higher than this in TPCs
  up_a <- max(d$t, na.rm = T) # Upper limit is CTmax
  up_b <- Inf
  up_c <- Inf
  upper_values <- c(up_s, up_a, up_b, up_c)

  # Fitting model using non-linear least squares
  fit <- nlsLM(p ~ pred_weibull(x=t, s, a, b, c),
             data = d,
             start = starting_values,
             lower =  lower_values,
             upper = upper_values,
             algorithm = "port")

  return(fit)
}