R Under development (unstable) (2023-12-18 r85702 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## Testing all SS functions in nlraa package
> require(nlraa)
Loading required package: nlraa
> 
> ## 1. SSbgf
> ## This function won't be tested here as it is used extensively in the 
> ## vignette nlraa::nlraa-AgronJ-paper
> ## However, in the future I will create a version reparameterized in 
> ## terms of unconstrained parameters, because the condition t.m < t.e 
> ## is not guranteed
> 
> ## 2. SSbgf4 
> ## It would also be beneficial to reparameterize this function
> ## for routine work in terms of unconstrained parameters
> data(sm)
> #' ## Let's just pick one crop
> sm2 <- subset(sm, Crop == "M")
> ## For this particular problem it is easier to 'fix' t.b and w.b
> fit <- nls(Yield ~ bgf2(DOY, w.max, w.b = 0, t.e, t.m, t.b = 141), 
+             data = sm2, start = list(w.max = 16, t.e= 240, t.m = 200))
> 
> ## 3. SSbgrp
> x <- 1:30
> y <- bgrp(x, 20, log(25), log(5)) + rnorm(30, 0, 1)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSbgrp(x, w.max, lt.e, ldt), data = dat)
> exp(confint(fit)[2:3,])
Waiting for profiling to be done...
          2.5%     97.5%
lt.e 24.894853 25.289427
ldt   4.803164  5.283775
> 
> ## 4. SSbg4rp
> set.seed(1234)
> x <- 1:100
> y <- bg4rp(x, 20, log(70), log(30), log(20)) + rnorm(100, 0, 1)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSbg4rp(x, w.max, lt.e, ldtm, ldtb), data = dat)
> exp(coef(fit))[-1]
    lt.e     ldtm     ldtb 
70.32360 30.02220 17.75343 
> 
> ## 5. SSdlf 
> ## Extended example in vignette 'nlraa-Oddi-LFMC'
> 
> ## 6. SSricker
> set.seed(123)
> x <- 1:30
> y <- 30 * x * exp(-0.3 * x) + rnorm(30, 0, 0.25)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSricker(x, a, b), data = dat)
> confint(fit)
Waiting for profiling to be done...
        2.5%      97.5%
a 29.7588789 30.4225878
b  0.2982059  0.3020252
> 
> ## 7. SSprofd
> ## I'm not a huge fan of this function as I think the dlf is more stable
> set.seed(1234)
> x <- 1:10
> y <- profd(x, 0.3, 0.05, 0.5, 4) + rnorm(10, 0, 0.01)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSprofd(x, a, b, c, d), data = dat)
> confint(fit, level = 0.9)
Waiting for profiling to be done...
          5%         95%
a 0.26478673  0.32478558
b 0.01251864  0.06535237
c 0.30950669  0.81949619
d 1.92550894 15.01686244
> 
> ## 8. SSnrh
> set.seed(1234)
> x <- seq(0, 2000, 100)
> y <- nrh(x, 35, 0.04, 0.83, 2) + rnorm(length(x), 0, 0.5)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSnrh(x, asym, phi, theta, rd), data = dat)
> confint(fit)
Waiting for profiling to be done...
             2.5%       97.5%
asym  34.29435363 40.14274335
phi    0.03732059  0.04944883
theta  0.57992430  0.86057599
rd     1.50629012  3.27600106
> 
> ## 9. SSlinp
> set.seed(123)
> x <- 1:30
> y <- linp(x, 0, 1, 20) + rnorm(30, 0, 0.5)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSlinp(x, a, b, xs), data = dat)
> confint(fit)
Waiting for profiling to be done...
         2.5%      97.5%
a  -0.3479829  0.6014031
b   0.9547236  1.0376608
xs 19.1980347 20.3174589
> 
> ## 10. SSplin
> set.seed(123)
> x <- 1:30
> y <- plin(x, 10, 20, 1) + rnorm(30, 0, 0.5)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSplin(x, a, xs, b), data = dat)
> confint(fit)
Waiting for profiling to be done...
         2.5%     97.5%
a   9.8543980 10.287226
xs 19.9897839 21.260686
b   0.9660968  1.187966
> 
> ## 11. SSquadp
> set.seed(123)
> x <- 1:30
> y <- quadp(x, 5, 1.7, -0.04, 20) + rnorm(30, 0, 0.6)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSquadp(x, a, b, c, xs), data = dat, algorithm = "port")
> ## Using port because default does not work
> summary(fit)

Formula: y ~ SSquadp(x, a, b, c, xs)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
a   5.217705   0.504927  10.334 1.07e-10 ***
b   1.654013   0.136705  12.099 3.49e-12 ***
c  -0.036514   0.007817  -4.671 8.01e-05 ***
xs 16.759684   1.166384  14.369 7.02e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5908 on 26 degrees of freedom

Algorithm "port", convergence message: both X-convergence and relative convergence (5)

> ## It's strange but confint will return NAs unless level is 0.5
> confint(fit, level = 0.5)
Waiting for profiling to be done...
           25%         75%
a   4.86984193  5.56309409
b   1.56049856  1.74703911
c  -0.04170704 -0.03116638
xs 16.10215516 18.10527059
> 
> ## 12. SSpquad
> set.seed(12345)
> x <- 1:40
> y <- pquad(x, 5, 20, 1.7, -0.04) + rnorm(40, 0, 0.6)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSpquad(x, a, xs, b, c), data = dat)
> confint(fit)
Waiting for profiling to be done...
          2.5%       97.5%
a   4.75353507  5.33189468
xs 19.59114287 20.65057362
b   1.62375057  1.95410208
c  -0.05397683 -0.03487333
> 
> ## 13. SSblin
> set.seed(1234)
> x <- 1:30
> y <- blin(x, 0, 0.75, 15, 1.75) + rnorm(30, 0, 0.5)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSblin(x, a, b, xs, c), data = dat)
> confint(fit)
Waiting for profiling to be done...
         2.5%      97.5%
a  -0.6619779  0.4050879
b   0.6757369  0.8041629
xs 13.7666329 15.4484429
c   1.6706631  1.7726210
> 
> ## 14. SSexpf
> set.seed(1234)
> x <- 1:15
> y <- expf(x, 10, -0.3) + rnorm(15, 0, 0.2)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSexpf(x, a, c), data = dat)
> confint(fit)
Waiting for profiling to be done...
       2.5%      97.5%
a  9.330309 10.5981097
c -0.329983 -0.2832748
> 
> ## 15. SSexpfp
> set.seed(12345)
> x <- 1:30
> y <- expfp(x, 10, 0.1, 15) + rnorm(30, 0, 1.5)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSexpfp(x, a, c, xs), data = dat)
> confint(fit)
Waiting for profiling to be done...
          2.5%      97.5%
a   8.74904806 10.8792954
c   0.09302385  0.1130315
xs 14.30133599 15.4435476
> 
> ## 16. SSpexpf
> set.seed(1234)
> x <- 1:30
> y <- pexpf(x, 20, 15, -0.2) + rnorm(30, 0, 1)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSpexpf(x, a, xs, c), data = dat)
> 
> ## 17. SSbell
> set.seed(1234)
> x <- 1:20
> y <- bell(x, 8, -0.0314, 0.000317, 13) + rnorm(length(x), 0, 0.5)
> dat <- data.frame(x = x, y = y)
> fit <- nls(y ~ SSbell(x, asym, a, b, xc), data = dat)
> confint(fit)
Waiting for profiling to be done...
              2.5%        97.5%
asym  7.2302693788  8.301058072
a    -0.0369890550 -0.026000855
b    -0.0007388103  0.001929091
xc   12.5419350314 13.606979346
> 
> ## 18. SSratio
> require(minpack.lm)
Loading required package: minpack.lm
> set.seed(1234)
> x <- 1:100
> y <- ratio(x, 1, 0.5, 1, 1.5) + rnorm(length(x), 0, 0.025)
> dat <- data.frame(x = x, y = y)
> fit <- nlsLM(y ~ SSratio(x, a, b, c, d), data = dat)
> 
> ## Testing nlsLMList
> require(nlme)
Loading required package: nlme
> data(Orange)
> fit.nlis.o <- nlsLMList(circumference ~ SSlogis(age, asym, xmid, scal), data = Orange)
> data(Soybean)
> fit.nlis1.s <- nlsLMList(weight ~ SSbgf(Time, w.max, t.e, t.m), data = Soybean)
Warning message:
In nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower,  :
  lmdif: info = -1. Number of iterations has reached `maxiter' == 50.

> fit.nlis2.s <- nlsLMList(weight ~ SSbgrp(Time, w.max, lt.e, ldt), data = Soybean)
Warning message:
In nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower,  :
  lmdif: info = -1. Number of iterations has reached `maxiter' == 50.

> fit.nlme <- nlme(fit.nlis2.s, random = pdDiag(w.max + lt.e ~ 1))
> ##plot(augPred(fit.nlme, level=0:1))
> 
> 
> 
> proc.time()
   user  system elapsed 
   2.53    0.20    2.67