# Non linear regression with nls()

#### November 18, 2018

R regression nls stats   ## Context

Last week, I was discussing about how to use nls() for a specific model with one of my colleague and I ended creating a piece of code to show what I was talking about! Even though there are many posts exploring nls() in more depth that I did (for instance this post on datascienceplus by Lionel Herzog), I thought I could share these lines of command here!

Basically, we were talking about a model where the temperature ($T$) follows a saturation curve starting from 10°C at t=0 (so T(0) = 10) and plateauing at $$T_{\inf}$$.

$$T(t) = T_{\inf} - (T_{\inf} - T_0)\exp(-kt)$$

## Goal and data

The goal here is to use nls() (Nonlinear Least Squares) to find $$k$$ and $$T_{inf}$$. For the sack of clarity, I simulate the data, i.e. I use the saturation curve with known parameter values, then I add some noise (here a white noise):

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  library(magrittr) # Parameters ## known T0 = 10 ## the ones we are looking for k = 0.2 Tinf = 20 # Simulate data ## time seqt <- seq(0, 50, .25) ## create a data frame simdata <- cbind( seqt = seqt, sim = Tinf - (Tinf- T0) * exp(-k * seqt) + .5 * rnorm(length(seqt)) ) %>% as.data.frame head(simdata) #R> seqt sim #R> 1 0.00 9.644506 #R> 2 0.25 11.062452 #R> 3 0.50 10.653373 #R> 4 0.75 10.816932 #R> 5 1.00 11.442564 #R> 6 1.25 12.968414 

## Use nls()

Now I call nls() to fit the data:

 1  res <- nls(sim ~ Tinf - (Tinf - 10)*exp(-k*seqt), simdata, list(Tinf = 1, k = .1)) 

All the information needed are stored in res and display via the print method:

  1 2 3 4 5 6 7 8 9 10  res #R> Nonlinear regression model #R> model: sim ~ Tinf - (Tinf - 10) * exp(-k * seqt) #R> data: simdata #R> Tinf k #R> 20.0955 0.2001 #R> residual sum-of-squares: 54.95 #R> #R> Number of iterations to convergence: 9 #R> Achieved convergence tolerance: 9.882e-06 

Let’s draw a quick plot:

 1 2 3 4 5  ## get the coefficients values cr <- coef(res) fitC <- function(x) cr - (cr - 10)*exp(-cr*x) plot(simdata[,1], simdata[,2], xlab = "Time (min)", ylab = "Temperature (°C)") curve(fitC, 0, 50, add = TRUE, col = "#2f85a4", lwd = 3) #### That’s all folks!

Session info
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35  sessionInfo() #R> R version 4.3.2 (2023-10-31) #R> Platform: x86_64-pc-linux-gnu (64-bit) #R> Running under: Ubuntu 22.04.3 LTS #R> #R> Matrix products: default #R> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 #R> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0 #R> #R> locale: #R>  LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 #R>  LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C #R>  LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C #R> #R> time zone: UTC #R> tzcode source: system (glibc) #R> #R> attached base packages: #R>  stats graphics grDevices utils datasets methods base #R> #R> other attached packages: #R>  magrittr_2.0.3 inSilecoRef_0.1.1 #R> #R> loaded via a namespace (and not attached): #R>  sass_0.4.7 utf8_1.2.4 generics_0.1.3 xml2_1.3.5 blogdown_1.18 #R>  stringi_1.8.1 httpcode_0.3.0 digest_0.6.33 evaluate_0.23 bookdown_0.36 #R>  fastmap_1.1.1 plyr_1.8.9 jsonlite_1.8.7 backports_1.4.1 crul_1.4.0 #R>  promises_1.2.1 fansi_1.0.5 jquerylib_0.1.4 bibtex_0.5.1 cli_3.6.1 #R>  shiny_1.8.0 rlang_1.1.2 ellipsis_0.3.2 cachem_1.0.8 yaml_2.3.7 #R>  tools_4.3.2 dplyr_1.1.4 httpuv_1.6.12 DT_0.30 rcrossref_1.2.0 #R>  curl_5.1.0 vctrs_0.6.4 R6_2.5.1 mime_0.12 lifecycle_1.0.4 #R>  stringr_1.5.1 fs_1.6.3 htmlwidgets_1.6.2 miniUI_0.1.1.1 pkgconfig_2.0.3 #R>  bslib_0.5.1 pillar_1.9.0 later_1.3.1 glue_1.6.2 Rcpp_1.0.11 #R>  highr_0.10 xfun_0.41 tibble_3.2.1 tidyselect_1.2.0 knitr_1.45 #R>  xtable_1.8-4 htmltools_0.5.7 rmarkdown_2.25 compiler_4.3.2