Non linear regression with nls()

November 18, 2018

  R regression nls
  stats

Kevin Cazelles  

          


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[1] - (cr[1] - 10)*exp(-cr[2]*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>   [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#R>   [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#R>   [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#R>  [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#R>  
#R>  attached base packages:
#R>  [1] stats     graphics  grDevices utils     datasets  methods   base     
#R>  
#R>  other attached packages:
#R>  [1] magrittr_2.0.3         inSilecoRef_0.0.1.9000
#R>  
#R>  loaded via a namespace (and not attached):
#R>   [1] tidyselect_1.2.0  xfun_0.35         bslib_0.4.1       vctrs_0.5.1      
#R>   [5] generics_0.1.3    miniUI_0.1.1.1    htmltools_0.5.4   yaml_2.3.6       
#R>   [9] utf8_1.2.2        rlang_1.0.6       later_1.3.0       pillar_1.8.1     
#R>  [13] jquerylib_0.1.4   httpcode_0.3.0    glue_1.6.2        DBI_1.1.3        
#R>  [17] lifecycle_1.0.3   plyr_1.8.8        stringr_1.5.0     blogdown_1.15    
#R>  [21] htmlwidgets_1.5.4 evaluate_0.18     knitr_1.41        fastmap_1.1.0    
#R>  [25] httpuv_1.6.6      curl_4.3.3        fansi_1.0.3       highr_0.9        
#R>  [29] Rcpp_1.0.9        xtable_1.8-4      promises_1.2.0.1  backports_1.4.1  
#R>  [33] DT_0.26           cachem_1.0.6      jsonlite_1.8.4    rcrossref_1.2.0  
#R>  [37] mime_0.12         digest_0.6.30     stringi_1.7.8     bookdown_0.30    
#R>  [41] dplyr_1.0.10      shiny_1.7.3       bibtex_0.5.0      cli_3.4.1        
#R>  [45] tools_4.2.2       sass_0.4.4        tibble_3.1.8      RefManageR_1.4.0 
#R>  [49] crul_1.3          pkgconfig_2.0.3   ellipsis_0.3.2    xml2_1.3.3       
#R>  [53] timechange_0.1.1  lubridate_1.9.0   assertthat_0.2.1  rmarkdown_2.18   
#R>  [57] httr_1.4.4        R6_2.5.1          compiler_4.2.2

Comment with Disqus