 # 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.666048 #R> 2 0.25 12.027247 #R> 3 0.50 10.286901 #R> 4 0.75 11.167841 #R> 5 1.00 11.852154 #R> 6 1.25 12.653063 

## 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> 19.9187 0.2094 #R> residual sum-of-squares: 52.4 #R> #R> Number of iterations to convergence: 9 #R> Achieved convergence tolerance: 1.953e-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 36 37  sessionInfo() #R> R version 4.2.2 (2022-10-31) #R> Platform: x86_64-pc-linux-gnu (64-bit) #R> Running under: Ubuntu 22.04.1 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 #R> #R> locale: #R>  LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 #R>  LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 #R>  LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C #R>  LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C #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.0.1.9000 #R> #R> loaded via a namespace (and not attached): #R>  tidyselect_1.2.0 xfun_0.35 bslib_0.4.1 vctrs_0.5.1 #R>  generics_0.1.3 miniUI_0.1.1.1 htmltools_0.5.4 yaml_2.3.6 #R>  utf8_1.2.2 rlang_1.0.6 later_1.3.0 pillar_1.8.1 #R>  jquerylib_0.1.4 httpcode_0.3.0 glue_1.6.2 DBI_1.1.3 #R>  lifecycle_1.0.3 plyr_1.8.8 stringr_1.5.0 blogdown_1.15 #R>  htmlwidgets_1.5.4 evaluate_0.18 knitr_1.41 fastmap_1.1.0 #R>  httpuv_1.6.6 curl_4.3.3 fansi_1.0.3 highr_0.9 #R>  Rcpp_1.0.9 xtable_1.8-4 promises_1.2.0.1 backports_1.4.1 #R>  DT_0.26 cachem_1.0.6 jsonlite_1.8.4 rcrossref_1.2.0 #R>  mime_0.12 digest_0.6.30 stringi_1.7.8 bookdown_0.30 #R>  dplyr_1.0.10 shiny_1.7.3 bibtex_0.5.0 cli_3.4.1 #R>  tools_4.2.2 sass_0.4.4 tibble_3.1.8 RefManageR_1.4.0 #R>  crul_1.3 pkgconfig_2.0.3 ellipsis_0.3.2 xml2_1.3.3 #R>  timechange_0.1.1 lubridate_1.9.0 assertthat_0.2.1 rmarkdown_2.18 #R>  httr_1.4.4 R6_2.5.1 compiler_4.2.2 

Comment with Disqus