nls.multstart_rmse_rsquared.R

# nls.multstart example calculating rsquared and RMSE

# load packages
library(nls.multstart)
library(soilphysics)
library(modelr)

# load in data

data("Chlorella_TRC")
Chlorella_TRC_test <- Chlorella_TRC[Chlorella_TRC$curve_id == 1,]

# run nls_multstart()

# define the Sharpe-Schoolfield equation
schoolfield_high <- function(lnc, E, Eh, Th, temp, Tc) {
  Tc <- 273.15 + Tc
  k <- 8.62e-5
  boltzmann.term <- lnc + log(exp(E/k*(1/Tc - 1/temp)))
  inactivation.term <- log(1/(1 + exp(Eh/k*(1/Th - 1/temp))))
  return(boltzmann.term + inactivation.term)
}

fits <- nls_multstart(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20),
                      data = Chlorella_TRC_test,
                      iter = 500,
                      start_lower = c(lnc=-10, E=0.1, Eh=0.5, Th=285),
                      start_upper = c(lnc=10, E=2, Eh=5, Th=330),
                      lower = c(lnc=-10, E=0, Eh=0, Th=0),
                      supp_errors = 'Y')

#-------------------#
# calculate RMSE ####
#-------------------# 

# https://stackoverflow.com/questions/43123462/how-to-obtain-rmse-out-of-lm-result

# method 1:
sqrt(mean(residuals(fits)^2))

# method 2:
# residual sum of squares
RSS <- c(crossprod(residuals(fits)))

# Mean squared error
MSE <- RSS / length(residuals(fits))

# Root MSE:
RMSE <- sqrt(MSE)     
RMSE

# method 3: use function from modelr
modelr::rmse(fits, Chlorella_TRC_test)

#------------------------#
# calculate pseudo r2 ####
#------------------------#

# i dont like this because its been shown that it is not a good indicator of fit

# discussion of this here: https://stackoverflow.com/questions/14530770/calculating-r2-for-a-nonlinear-least-squares-fit

# can use functions available in various packages already
soilphysics::Rsq(fits)
modelr::rsquare(fits, Chlorella_TRC_test)


# https://gist.github.com/e6a4c7e7fff6c37a761261131df8b6ca