Running 'test_CV_response.R' [0s]
  Running 'test_boot_nlme.R' [166s]
  Running 'test_environments.R' [148s]
  Running 'test_nlsLMList.R' [4s]
  Running 'test_predict_gam.R' [2s]
  Running 'test_predict_nlme.R' [158s]
  Running 'test_predict_nls.R' [3s]
  Running 'test_predict_summary_simulate.R' [2s]
  Running 'test_predict_varFunc.R' [2s]
  Running 'test_sim_boot_lfmc.R' [172s]
  Running 'test_simulate_gam.R' [2s]
  Running 'test_simulate_gls.R' [156s]
  Running 'test_simulate_gnls.R' [2s]
  Running 'test_simulate_lm.R' [162s]
  Running 'test_simulate_lme.R' [2s]
  Running 'test_simulate_nlme.R' [160s]
  Running 'test_var_cov.R' [2s]
  Running 'tests_SSall.R' [167s]
  Running 'tests_SSother.R' [162s]
Running the tests in 'tests/test_boot_nlme.R' failed.
Complete output:
  > ## Test for boot_nlme
  > require(nlme)
  Loading required package: nlme
  > require(nlraa)
  Loading required package: nlraa
  > require(car)
  Loading required package: car
  Loading required package: carData
  > 
  > ## I only run this on my computer and occasionally
  > run.boot.test <- Sys.info()[["user"]] == "fernandomiguez" && FALSE
  > 
  > data(barley, package = "nlraa")
  > set.seed(101)
  > ## Does Boot produce 'good' confidence intervals?
  > 
  > if(run.boot.test){
  +   
  +   fit.nls <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley)
  +   
  +   ## Profiled confidence intervals
  +   confint(fit.nls)
  + 
  +   ## Does bootstrap as implemented in car underestimates confidence intervals?  
  +   fit.nls.bt <- Boot(fit.nls)
  +   ## It is not possible to parallelize the previous code, why?
  +   
  +   fit.gnls <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley)
  +   
  +   intervals(fit.gnls)
  +   
  +   ## This takes: 9.3 seconds
  +   system.time(fit.gnls.bt <- boot_nlme(fit.gnls, R = 999, cores = 4))
  +   
  +   ## This takes also <10 seconds 
  +   fit.gnls.bt2 <- boot_nlme(fit.gnls, R = 999, psim = 0, cores = 4)
  +   
  +   ## Compare confidence intervals
  +   confint(fit.nls) ## nls profile
  +   confint(fit.nls.bt) ## Boot
  +   confint(fit.gnls.bt) ## boot_nlme psim = 1
  +   confint(fit.gnls.bt2) ## boot_nlme psim = 0
  +   ## These intervals do not agree exaclty because different assumptions 
  +   ## are made and I'm running bootstrap a relatively small number of times
  + }
  > 
  > if(run.boot.test){
  +   
  +   ## Simple test for barley and an object of class 'gnls'
  +   ## Simplify the dataset to make the set up simpler
  +   barley2 <- subset(barley, year < 1974)
  + 
  +   fit.lp.gnls2 <- gnls(yield ~ SSlinp(NF, a, b, xs), data = barley2)
  + 
  +   intervals(fit.lp.gnls2)
  + 
  +   ## Compare this to the bootstrapping approach
  +   ## This take about ~13 seconds
  +   system.time(fit.lp.gnls2.bt <- boot_nlme(fit.lp.gnls2, R = 2000, cores = 4))
  + 
  +   summary(fit.lp.gnls2.bt) ## Bias is low, which is good
  + 
  +   confint(fit.lp.gnls2.bt, type = "perc")
  + 
  +   hist(fit.lp.gnls2.bt, 1, ci = "perc")
  +   hist(fit.lp.gnls2.bt, 2, ci = "perc")
  +   hist(fit.lp.gnls2.bt, 3, ci = "perc")
  + 
  +   ## Testing the bootstrap function with an object which
  +   ## contains factors
  +   
  +   barley2$year.f <- as.factor(barley2$year)
  + 
  +   cfs <- coef(fit.lp.gnls2)
  + 
  +   fit.lp.gnls3 <- update(fit.lp.gnls2, 
  +                          params = list(a + b + xs ~ year.f),
  +                          start = c(cfs[1], 0, 0, 0, 
  +                                    cfs[2], 0, 0, 0,
  +                                    cfs[3], 0, 0, 0))
  + 
  +   intervals(fit.lp.gnls3)
  + 
  +   ## This takes 13 seconds (Mac), not bad. Windows (21 seconds)
  +   system.time(fit.lp.gnls3.bt <- boot_nlme(fit.lp.gnls3, R = 3000, cores = 4))
  + 
  +   summary(fit.lp.gnls3.bt)
  +   
  +   confint(fit.lp.gnls3.bt, type = "perc")
  +   
  +   hist(fit.lp.gnls3.bt, 1, ci = "perc")
  +   hist(fit.lp.gnls3.bt, 2, ci = "perc")
  +   hist(fit.lp.gnls3.bt, 3, ci = "perc")
  + 
  +   ## Testing the function for 'nlme'. 
  +   barley$year.f <- as.factor(barley$year)
  + 
  +   barleyG <- groupedData(yield ~ NF | year.f, data = barley)
  + 
  +   fitL.bar <- nlsList(yield ~ SSlinp(NF, a, b, xs), data = barleyG)
  + 
  +   fit.bar.nlme <- nlme(fitL.bar, random = pdDiag(a + b + xs ~ 1))
  + 
  +   plot(augPred(fit.bar.nlme, level = 0:1))
  +   ## Confidence intervals of the model fixed parameters
  +   intervals(fit.bar.nlme, which = "fixed")
  + 
  +   ## Bootstrap for the asymptote
  +   fna <- function(x) fixef(x)[1] + fixef(x)[2] * fixef(x)[3]
  + 
  +   ## This takes much longer... 215-237 seconds on my laptop
  +   system.time(fit.bar.nlme.bt <- boot_nlme(fit.bar.nlme, f = fna, R = 2000, cores = 4))
  + 
  +   confint(fit.bar.nlme.bt, type = "perc")
  + 
  +   ## This is good because the estimate on the original data is
  +   ## 348.9912, which is in the middle of the confidence interval
  +   hist(fit.bar.nlme.bt, 1, ci = "perc")
  + 
  + }
  > 
  > ## Trying to fit the model using different approaches
  > ## This time with the SSasymp function
  > 
  > if(run.boot.test){
  +   
  +   fmL1 <- nlsList(yield ~ SSasymp(NF, Asym, R0, lrc), data = barleyG)
  +   
  +   fm1 <- nlme(fmL1)
  +   
  +   fm2 <- update(fm1, random = pdDiag(Asym + R0 + lrc ~ 1))
  +   ## Arguably the second model is better
  +   
  +   ## fm3 <- nlmer(yield ~ SSasymp(NF, Asym, R0, lrc) ~ (Asym | year.f) + (R0 | year.f), data = barley,
  +   ##             start = getInitial(yield ~ SSasymp(NF, Asym, R0, lrc), data = barley))
  +   ## This does not fail now, but it did before. 
  +   
  +   ## This is very slow... 17 seconds for 10 attempts
  +   ## system.time(fm1.bt <- boot_nlme(fm1, R = 10, cores = 4))
  +   ## What about the second model?
  +   ## This one takes about 39 sec (Mac and Windows?)
  +   system.time(fm2.bt <- boot_nlme(fm2, R = 1000, cores = 4))
  +   
  +   summary(fm2.bt)
  +   
  +   confint(fm2.bt)
  +   
  +   hist(fm2.bt, 1)
  +   hist(fm2.bt, 2)
  +   hist(fm2.bt, 3)
  +   
  +   ## I feel validated! These intervals are very, very close to the ones
  +   ## obtained through normal approximation
  +   ## I wrote the previous statement a long time aga, I don't know
  +   ## if it is still true (2020-5-22)
  + 
  +   ## rstanarm failed with this example. Chains did not converge (2020-05-22)
  +   ## fm3 <- stan_nlmer(yield ~ SSasymp(NF, Asym, R0, lrc) ~ Asym + R0 + lrc | year.f, 
  +   ##                  data = barley, cores = 2)  
  +   
  +   ## What is the impact of psim = 2 vs. psim = 3 on bootstrapping?
  +   data(barley, package = "nlraa")
  +   barley$year.f <- as.factor(barley$year)
  +   barleyG <- groupedData(yield ~ NF | year.f, data = barley)
  +   fmL1 <- nlsList(yield ~ SSasymp(NF, Asym, R0, lrc), data = barleyG)
  +   fm1 <- nlme(fmL1)
  +   fm2 <- update(fm1, random = pdDiag(Asym + R0 + lrc ~ 1))
  +   fm3 <- update(fm1, random = pdDiag(Asym + R0 ~ 1))
  + 
  +   anova(fm2, fm3)
  +   ## Bootstrap the random effects
  +   vcov_est <- function(x) diag(var_cov(x, type = "random"))
  +   
  +   system.time(fm3.vc.bt2 <- boot_nlme(fm3, vcov_est, psim = 2, cores = 4))
  +   
  +   fm3.vc.bt2
  +   
  +   hist(fm3.vc.bt2, 1, ci = "perc")
  +   hist(fm3.vc.bt2, 2, ci = "perc")
  +   
  +   system.time(fm3.vc.bt3 <- boot_nlme(fm3, vcov_est, psim = 3, cores = 4))
  + 
  +   fm3.vc.bt3
  +   
  +   hist(fm3.vc.bt3, 1, ci = "perc")
  +   hist(fm3.vc.bt3, 2, ci = "perc")
  +   
  +   ## Does this show that this method results in biased estimates of the
  +   ## variance components? Bootstrapping has substantial bias
  +   ## Is it because I attempt to fix the variance (I use empirical = TRUE in MASS::rmvnorm)
  +   ## I could either set empirical to FALSE or simulate new variance components...
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.15    1.48 
Running the tests in 'tests/test_environments.R' failed.
Complete output:
  > ## Test and learn scoping and environments
  > ## I need to refine this and add more notes related to whay I learned
  > 
  > run.test.environments <- FALSE
  > 
  > if(run.test.environments){
  + 
  +   ## Scoping, frames and environments
  +   g <- function(a, b, c = 0) a + b + c
  +   print(parent.frame())  
  +   f1 <- function(){
  +     x <- 1
  +     print(ls())
  +     print(ls(envir = parent.frame()))
  +     print(parent.frame())
  +     f2 <- function(){
  +       w <- 2
  +       print(ls())
  +       print(ls(envir = parent.frame()))
  +       print(parent.frame())
  +       ans1 <- g(a = w, eval(x))
  +       f3 <- function(){
  +         print(ls(envir = parent.frame()))  
  +         print(parent.frame())
  +         print(eval(x))
  +         z <- eval(ans1) + eval(w) + eval(x)
  +       }
  +       ans1 <- f3()
  +     }
  +     ans <- f2()
  +     ans
  +   }
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.04    0.20 
Running the tests in 'tests/test_predict_nlme.R' failed.
Complete output:
  > ## Example
  > require(ggplot2)
  Loading required package: ggplot2
  > require(nlme)
  Loading required package: nlme
  > require(nlraa)
  Loading required package: nlraa
  > 
  > ## I do not want to run this on CRAN servers
  > run.predict.nlme <- Sys.info()[["user"]] == "fernandomiguez" && FALSE
  > 
  > if(run.predict.nlme){
  +   
  +   data(Orange)
  +   
  +   ## All models should be fitted using Maximum Likelihood
  +   fm.L <- nlme(circumference ~ SSlogis(age, Asym, xmid, scal), 
  +                random = pdDiag(Asym + xmid + scal ~ 1), 
  +                method = "ML", data = Orange)
  +   fm.G <- nlme(circumference ~ SSgompertz(age, Asym, b2, b3), 
  +                random = pdDiag(Asym + b2 + b3 ~ 1), 
  +                method = "ML", data = Orange)
  +   fm.F <- nlme(circumference ~ SSfpl(age, A, B, xmid, scal), 
  +                random = pdDiag(A + B + xmid + scal ~ 1), 
  +                method = "ML", data = Orange)
  +   fm.B <- nlme(circumference ~ SSbg4rp(age, w.max, lt.e, ldtm, ldtb), 
  +                random = pdDiag(w.max + lt.e + ldtm + ldtb ~ 1), 
  +                method = "ML", data = Orange)
  +   
  +   ## Print the table with weights
  +   IC_tab(fm.L, fm.G, fm.F, fm.B)
  +   
  +   ## Each model prediction is weighted according to their AIC values
  +   prd <- predict_nlme(fm.L, fm.G, fm.F, fm.B)
  +   
  +   ggplot(data = Orange, aes(x = age, y = circumference)) + 
  +     geom_point() + 
  +     geom_line(aes(y = predict(fm.L, level = 0), color = "Logistic")) +
  +     geom_line(aes(y = predict(fm.G, level = 0), color = "Gompertz")) +
  +     geom_line(aes(y = predict(fm.F, level = 0), color = "4P-Logistic")) +  
  +     geom_line(aes(y = predict(fm.B, level = 0), color = "Beta")) +
  +     geom_line(aes(y = prd, color = "Avg. Model"), size = 1.2)
  +   
  +   ## Testing confidence intervals
  +   prdc <- predict_nlme(fm.L, fm.G, fm.F, fm.B, interval = "conf")
  +   OrangeA <- cbind(Orange, prdc)
  +   
  +   ggplot(data = OrangeA, aes(x = age, y = circumference)) + 
  +     geom_point() + 
  +     geom_line(aes(y = predict(fm.L, level = 0), color = "Logistic")) +
  +     geom_line(aes(y = predict(fm.G, level = 0), color = "Gompertz")) +
  +     geom_line(aes(y = predict(fm.F, level = 0), color = "4P-Logistic")) +  
  +     geom_line(aes(y = predict(fm.B, level = 0), color = "Beta")) +
  +     geom_line(aes(y = prd, color = "Avg. Model"), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3)
  +   
  +   ## Confidence at the level of individual or Tree in this case
  +   prdc1 <- predict_nlme(fm.L, fm.G, fm.F, fm.B, interval = "conf", plevel = 1, level = 0.9)
  +   OrangeAC <- cbind(Orange, prdc1)
  +   
  +   ggplot(data = OrangeAC, aes(x = age, y = circumference, color = Tree)) + 
  +     facet_wrap(~ Tree, ncol = 5) + 
  +     geom_point() + 
  +     geom_line(aes(y = Estimate), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95, color = NULL, fill = Tree), alpha = 0.3) + 
  +     ggtitle("90% Confidence bands at the level of 'Tree'")
  +   
  +   ## Predictions only make sense at the level of individual or Tree in this case
  +   prdp <- predict_nlme(fm.L, fm.G, fm.F, fm.B, interval = "pred", plevel = 1, level = 0.9)
  +   OrangeAP <- cbind(Orange, prdp)
  +   
  +   ggplot(data = OrangeAP, aes(x = age, y = circumference, color = Tree)) + 
  +     facet_wrap(~ Tree, ncol = 5) + 
  +     geom_point() + 
  +     geom_line(aes(y = Estimate), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95, color = NULL, fill = Tree), alpha = 0.3) + 
  +     ggtitle("90% Prediciton bands at the level of 'Tree'")
  +   
  +   ## Testing new-prediction
  +   prdnp <- predict_nlme(fm.L, interval = "new-prediction", plevel = 1)
  +   OrangeANP <- cbind(Orange, prdnp)
  +   
  +   ## These should really be unlabal
  +   OrangeANP$lbs <- with(OrangeANP, paste("New", Tree, sep = "_"))
  +   ggplot(data = OrangeANP, aes(x = age, y = circumference)) + 
  +     facet_wrap(~ lbs) + 
  +     geom_line(aes(y = Estimate, group = Tree)) + 
  +     geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.5) + 
  +     ggtitle("95% Prediction Band for a NEW Tree")
  +   
  +   ## Example using Soybean
  +   data(Soybean)
  +   
  +   ## All models should be fitted using Maximum Likelihood
  +   fms.L <- nlme(weight ~ SSlogis(Time, Asym, xmid, scal), 
  +                 random = pdDiag(Asym + xmid ~ 1), 
  +                 weights = varPower(),
  +                 method = "ML", data = Soybean)
  +   fms.G <- nlme(weight ~ SSgompertz(Time, Asym, b2, b3), 
  +                 random = pdDiag(Asym + b2 ~ 1), 
  +                 weights = varPower(),
  +                 method = "ML", data = Soybean)
  +   fms.B <- nlme(weight ~ SSbg4rp(Time, w.max, lt.e, ldtm, ldtb), 
  +                 random = pdDiag(w.max + lt.e ~ 1), 
  +                 weights = varPower(),
  +                 method = "ML", data = Soybean) 
  +   fms.W <- nlme(weight ~ SSweibull(Time, Asym, Drop, lrc, pwr),
  +                  random = pdDiag(Asym + lrc ~ 1),
  +                  weights = varPower(),
  +                  method = "ML", data = Soybean)
  +   fms.E <- nlme(weight ~ SSexplin(Time, cm, rm, tb),
  +                 random = pdDiag(cm + rm ~ 1),
  +                 weights = varPower(),
  +                 method = "ML", data = Soybean)
  +   fms.T <- nlme(weight ~ SStrlin(Time, a, b, xs1, c, xs2, d),
  +                 random = pdDiag(b + xs1 + c + d ~ 1),
  +                 weights = varPower(),
  +                 method = "ML", data = Soybean)
  +   
  +   IC_tab(fms.L, fms.G, fms.B, fms.W, fms.E, fms.T)
  +   
  +   ## The beta growth function is the best one by far
  +   ## No need to perform averaging
  +   prdsc <- predict_nlme(fms.B, interval = "conf", level = 0.9)
  +   SoybeanAC <- cbind(Soybean, prdsc)
  +   
  +   ggplot(data = SoybeanAC, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = Estimate), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "purple", alpha = 0.3) + 
  +     ggtitle("90% Confidence bands ")
  + 
  +   prdsp <- predict_nlme(fms.B, interval = "pred", level = 0.9, plevel = 1)
  +   SoybeanAP <- cbind(Soybean, prdsp)
  +   
  +   ggplot(data = SoybeanAP, aes(x = Time, y = weight)) + 
  +     facet_wrap(~ Plot) + 
  +     geom_point() + 
  +     geom_line(aes(y = Estimate), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "purple", alpha = 0.3) + 
  +     ggtitle("90% Prediction bands for each plot ")
  +   
  +   ## Example with ChickWeight
  +   data(ChickWeight)
  +   
  +   fmc.L <- nlme(weight ~ SSlogis(Time, Asym, xmid, scal), 
  +                 random = pdDiag(Asym + xmid ~ 1), 
  +                 method = "ML", data = ChickWeight)
  +   
  +   fmc.G <- nlme(weight ~ SSgompertz(Time, Asym, b2, b3), 
  +                 random = pdDiag(Asym + b2 ~ 1), 
  +                 method = "ML", data = ChickWeight)
  +   
  +   fmc.F <- nlme(weight ~ SSfpl(Time, A, B, xmid, scal), 
  +                 random = pdDiag(A + B + xmid ~ 1), 
  +                 method = "ML", data = ChickWeight)
  + 
  +   fmc.B <- nlme(weight ~ SSbg4rp(Time, w.max, lt.e, ldtm, ldtb), 
  +                 random = pdDiag(w.max + lt.e ~ 1), 
  +                 method = "ML", data = ChickWeight)
  +   
  +   ## FPL FTW  
  +   IC_tab(fmc.L, fmc.G, fmc.F, fmc.B)
  + 
  +   prdsch <- predict_nlme(fmc.F, interval = "conf", level = 0.90)
  +   ChickWeightAC <- cbind(ChickWeight, prdsch)  
  +   
  +   ggplot(data = ChickWeightAC, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = Estimate), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "purple", alpha = 0.3) + 
  +     ggtitle("90% Confidence bands")
  + 
  +   prdschp <- predict_nlme(fmc.F, interval = "pred", level = 0.99, plevel = 1)
  +   ChickWeightAP <- cbind(ChickWeight, prdschp)  
  +   
  +   ggplot(data = ChickWeightAP, aes(x = Time, y = weight)) + 
  +     facet_wrap(~Chick) + 
  +     geom_point() + 
  +     geom_line(aes(y = Estimate), size = 1.2) + 
  +     geom_ribbon(aes(ymin = Q0.5, ymax = Q99.5), fill = "purple", alpha = 0.3) + 
  +     ggtitle("99% Prediction bands")
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     2.03    0.34    2.31 
Running the tests in 'tests/test_sim_boot_lfmc.R' failed.
Complete output:
  > ## For this test I'm planning to use the lfmc dataset
  > ## The goals is to test both simulate and bootstrap
  > require(nlme)
  Loading required package: nlme
  > require(nlraa)
  Loading required package: nlraa
  > require(ggplot2)
  Loading required package: ggplot2
  > 
  > run.sim.boot.lfmc <- Sys.info()[["user"]] == "fernandomiguez" && FALSE
  > 
  > if(run.sim.boot.lfmc){
  +   
  +   data(lfmc, package = "nlraa")
  +   
  +   lfmc <- droplevels(subset(lfmc, leaf.type != "Grass E"))
  + 
  +   lfmcG <- groupedData(lfmc ~ time | group, data = lfmc)
  +   
  +   fitL <- nlsList(lfmc ~ SSdlf(time, asym, a2, xmid, scal), data = lfmcG)  
  + 
  +   plot(fitL)
  +   plot(intervals(fitL))
  + 
  +   fm <- nlme(fitL, random = pdDiag(asym + a2 + xmid + scal ~ 1))
  +   
  +   fm1 <- update(fm, random = pdDiag(asym + xmid ~ 1))
  +   ## I'll ignore site for this example
  +   fxf <- fixef(fm1)
  +   
  +   fm2 <- update(fm1, fixed = list(asym + a2 + xmid + scal ~ leaf.type),
  +                 weights = varPower(),
  +                 start = c(fxf[1], rep(0, 2), fxf[2], rep(0, 2), fxf[3], rep(0, 2), fxf[4], rep(0, 2)))
  +   
  +   anova(fm2)
  +   
  +   ## Simulate from this model 
  +   ## psim = 1, level = 0
  +   ## This is at the level of the mean function for each level of factor 'leaf.type'
  +   sim10 <- simulate_nlme(fm2, nsim = 100, level = 0, value = "data.frame") 
  + 
  +   ggplot(data = sim10, aes(x = time, y = lfmc)) + 
  +     facet_wrap(~ leaf.type) + 
  +     geom_line(aes(x = time, y = sim.y, group = ii, color = leaf.type)) + 
  +     geom_point()
  + 
  +   ## This is at the level of each 'group' 
  +   sim11 <- simulate_nlme(fm2, nsim = 100, level = 1, value = "data.frame") 
  +   
  +   sim11$ID <- with(sim11, paste0(ii,"_",group))
  + 
  +   ## This is interesting because it shows why we need bootstrapping, right?                   
  +   ggplot(data = sim11, aes(x = time, y = lfmc)) + 
  +     facet_wrap(~ leaf.type) + 
  +     geom_line(aes(x = time, y = sim.y, group = ID, color = leaf.type)) + 
  +     geom_point()
  + 
  +   ## This is at the level of observation 
  +   ## We can't go deeper than level 1 (fm2$dims$Q)
  +   sim21 <- simulate_nlme(fm2, nsim = 100, psim = 2, level = 1, value = "data.frame") 
  +   
  +   ggplot(data = sim21, aes(x = time, y = lfmc)) + 
  +     facet_wrap(~ leaf.type) + 
  +     geom_point(aes(x = time, y = sim.y, color = leaf.type)) + 
  +     geom_point()
  +   
  +   ## Bootstrapping section
  +   prd_fun <- function(x) predict(x, level = 0)
  +   ## This takes about 111 seconds
  +   system.time(bprd <- boot_nlme(fm2, prd_fun, cores = 4))
  +   
  +   lfmcG$lwr.q <- apply(t(bprd$t), 1, quantile, probs = 0.05, na.rm = TRUE)
  +   lfmcG$upr.q <- apply(t(bprd$t), 1, quantile, probs = 0.95, na.rm = TRUE)
  + 
  +   lfmcG$prd <- predict(fm2, level = 0)
  +   
  +   ggplot(data = lfmcG, aes(x = time, y = lfmc)) + 
  +     facet_wrap(~ leaf.type) + 
  +     geom_ribbon(aes(x = time, ymin = lwr.q, ymax = upr.q, fill = leaf.type), alpha = 0.5) + 
  +     geom_point() + 
  +     geom_line(aes(y = prd)) + 
  +     ggtitle("Bootstrapped fits at the fixed level of species")
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.93    0.29    2.18 
Running the tests in 'tests/test_simulate_gls.R' failed.
Complete output:
  > ## test simulate_lme and simulate_gls
  > require(nlme)
  Loading required package: nlme
  > require(nlraa)
  Loading required package: nlraa
  > require(ggplot2)
  Loading required package: ggplot2
  > 
  > run.simulate.gls.test <- Sys.info()[["user"]] == "fernandomiguez"
  > 
  > if(run.simulate.gls.test){
  +   
  +   data(Orthodont)
  +   
  +   fm1 <- gls(distance ~ age, data = Orthodont,
  +              correlation = corCAR1(form = ~ age | Subject)) ## Sub by default
  +   
  +   ## Visualize correlation
  +   fm1.vc <- var_cov(fm1)
  +   fm1.vc2 <- fm1.vc[1:28, 1:28]
  +   image(log(fm1.vc2[,ncol(fm1.vc2):1]), 
  +         main = "Covariance (residual or marginal)",
  +         xaxt = "n", yaxt = "n")
  +   
  +   ## auto-correlation
  +   plot(ACF(fm1))
  +   
  +   sim1 <- simulate_lme(fm1, nsim = 2e2, value = "data.frame")
  + 
  +   sim1$Sub_ii <- paste0(sim1$Subject,"_",sim1$ii)
  +   
  +   ggplot(data = sim1, aes(age, distance )) + 
  +     geom_line(aes(y = sim.y, group = Sub_ii), color = "grey",  
  +               alpha = 0.3) + 
  +     geom_point() + ggtitle("Orthodont")
  +   
  +   ## Trying the ChickWeight dataset
  +   data(ChickWeight)
  +   
  +   fm2 <- lme(weight ~ Time + I(Time^2) + I(Time^3), 
  +              data = ChickWeight, weights = varPower(),
  +              random = ~ Time)
  +   
  +   ## Level = 0 is at the population-level           
  +   sim2 <- simulate_lme(fm2, nsim = 200, level = 0, value = "data.frame")
  +   
  +   sim2$Chick_ii <- paste0(sim2$Chick,"_",sim2$ii)
  +   
  +   ggplot(data = sim2, aes(Time, weight)) + 
  +     geom_line(aes(y = sim.y, group = Chick_ii), color = "grey",  
  +               alpha = 0.3) + 
  +     geom_point()
  +   
  +   ## Test with the Orange dataset
  +   data("Orange")
  +   
  +   fm3 <- lme(circumference ~ age, data = Orange,
  +              random = ~ 1 | Tree)
  +   
  +   ## Small data set, no issues
  +   plot(fm3)
  +   
  +   sim3 <- simulate_lme(fm3, nsim = 200, value = "data.frame")
  +   
  +   sim3$Tree_ii <- paste0(sim3$Tree,"_",sim3$ii)
  +   
  +   ggplot(data = sim3, aes(age, circumference)) + 
  +     geom_line(aes(y = sim.y, group = Tree_ii), color = "grey",  
  +               alpha = 0.3) + 
  +     geom_point()
  +   
  +   ## Testing scoping issues
  +   ## Script from bug report
  +   simulation_summary_fixed <- function(duration_per_period = 6L) {
  + 
  +     treatment <- as.factor(c(0,1))
  +     time <- 1:20
  +     cluster <- as.factor(c("a", "b", "c"))
  +     SWdesign <- expand.grid(treatment = treatment, time = time, cluster = cluster)
  +     SWdesign$id <- with(SWdesign, paste0(cluster, "_", treatment))
  +     SWdesign$outcome <- with(SWdesign, as.numeric(treatment) + 
  +                                (1 + as.numeric(treatment)) * time + 
  +                                scale(as.numeric(cluster)) * 0.5 + 
  +                                rnorm(nrow(SWdesign)))
  +     
  +     cor_init <- corAR1(form= ~ time|id )
  +     cor_init <- Initialize(cor_init, data = SWdesign )
  +     
  +     SWdesign$outcome <- rnorm( nrow(SWdesign) )
  +     
  +     local.gls <- gls( outcome ~ treatment*time, correlation=cor_init, data=SWdesign, method="ML") 
  +     local.lme <- lme( outcome ~ treatment*time, random = ~ 1 | cluster, 
  +                       data=SWdesign, method="ML") 
  + 
  +     ans1 <- simulate_lme(local.lme, psim=2, data = SWdesign)
  +     ans2 <- simulate_gls(local.gls, psim=2, data = SWdesign)
  +   }
  +   
  +   ## If this runs without errors it means the above works
  +   simulation_summary_fixed()
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.92    0.31    2.21 
Running the tests in 'tests/test_simulate_lm.R' failed.
Complete output:
  > ## Testing the levels of simulation
  > require(ggplot2)
  Loading required package: ggplot2
  > require(nlraa)
  Loading required package: nlraa
  > 
  > run.simulte.lm.test <- Sys.info()[["user"]] == "fernandomiguez"
  > 
  > if(run.simulte.lm.test){
  +   
  +   data(Orange)
  +   
  +   ## Simple simulation
  +   fit1 <- lm(circumference ~ age, data = Orange)
  +   sims1 <- simulate_lm(fit1, nsim = 100, value = "data.frame")
  +    
  +   ggplot(data = sims1) + 
  +     geom_line(aes(x = age, y = sim.y, group = ii), 
  +                color = "gray", alpha = 0.5) + 
  +     geom_point(aes(x = age, y = circumference)) + 
  +     geom_smooth(aes(x = age, y = circumference), method = "lm", 
  +                 linetype = 2, se = FALSE)
  +   
  +   prd0 <- cbind(Orange, as.data.frame(predict(fit1, interval = "confidence")))
  +   
  +   ggplot() + 
  +     geom_ribbon(data = prd0, aes(x = age, ymin = lwr, ymax = upr),
  +                 fill = "purple", alpha = 0.5) + 
  +     geom_point(data = Orange, aes(age, circumference))
  +     
  +   ## Testing quadratic
  +   fit2 <- lm(circumference ~ age + I(age^2), data = Orange)
  +   sims2 <- simulate_lm(fit2, nsim = 100, value = "data.frame")
  +   
  +   ggplot(data = sims2) + 
  +     geom_line(aes(x = age, y = sim.y, group = ii), 
  +               color = "gray", alpha = 0.5) + 
  +     geom_point(aes(x = age, y = circumference)) 
  +   
  +   ## Testing psim == 2
  +   sims3 <- simulate_lm(fit1, psim = 2, nsim = 100, value = "data.frame")
  +   
  +   ggplot(data = sims3) + 
  +     geom_line(aes(x = age, y = sim.y, group = ii), 
  +               color = "gray", alpha = 0.5) + 
  +     geom_point(aes(x = age, y = circumference)) 
  +   
  +   ## Testing psim == 3
  +   ## This level is very similar to the 'prediction' below
  +   sims4 <- simulate_lm(fit1, psim = 3, nsim = 100, value = "data.frame")
  +   
  +   sims4a <- aggregate(sim.y ~ age, data = sims4, FUN = quantile,
  +                       probs = 0.05)
  +   sims4b <- aggregate(sim.y ~ age, data = sims4, FUN = quantile,
  +                       probs = 0.95)
  +   sims4c <- merge(sims4a, sims4b, by = "age")
  +   
  +   names(sims4c) <- c("age","lwr","upr")
  +   
  +   ggplot() + 
  +     geom_ribbon(data = sims4c,
  +                 aes(x = age, ymin = lwr, ymax = upr), 
  +                 fill = "blue2", alpha = 0.3) + 
  +     geom_point(data = Orange, aes(x = age, y = circumference)) 
  + 
  +   prd <- predict(fit1, interval = "prediction")  
  +   prdd <- cbind(Orange, prd)
  +     
  +   ggplot() + 
  +     geom_ribbon(data = prdd, aes(x = age, ymin = lwr, ymax = upr), 
  +               fill = "purple", alpha = 0.3) + 
  +     geom_point(data = Orange, aes(x = age, y = circumference)) 
  +   
  +   ## I fitted this model with brms
  +   ## library(brms)
  +   ## obr <- brm(circumference ~ age, data = Orange)
  +   ## This is the same as interval 'confidence'
  +   ## plot(conditional_effects(obr), points = T)
  +   ## This is the same as interval 'prediction'
  +   ## plot(conditional_effects(obr, method = "posterior_predict"), points = T)
  +  
  +   ## Trying wild residuals
  +   ## Testing psim == 2
  +   sims3 <- simulate_lm(fit1, psim = 2, nsim = 100, value = "data.frame", resid.type = "wild")
  +   
  +   sims3$Tree_ID <- with(sims3, paste0(Tree,"_",ii))
  +   
  +   ggplot(data = sims3) + 
  +     geom_line(aes(x = age, y = sim.y, group = Tree_ID), 
  +               color = "gray", alpha = 0.5) + 
  +     geom_point(aes(x = age, y = circumference)) 
  +   
  +   ## Testing scoping issues
  +   ## This version does not work
  +   # f1 <- function(){
  +   #   dat <- data.frame(x = rnorm(10), y = rnorm(10))
  +   #   fm00 <- lm(y ~ x, data = dat)
  +   #   ans <- simulate_lm(fm00)
  +   #   ans
  +   # }
  +   
  +   f1 <- function(){
  +     dat <- data.frame(x = rnorm(10), y = rnorm(10))
  +     fm00 <- lm(y ~ x, data = dat)
  +     ans <- simulate_lm(fm00, data = dat)
  +     ans
  +   }
  +   ## Well, this seems to work because 'eval' is not inside a function
  +   ## such as 'getData'???
  +   res1 <- f1()
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.81    0.34    2.14 
Running the tests in 'tests/test_simulate_nlme.R' failed.
Complete output:
  > ## Testing simulate_nlme levels
  > require(nlme)
  Loading required package: nlme
  > require(nlraa)
  Loading required package: nlraa
  > require(ggplot2)
  Loading required package: ggplot2
  > require(car)
  Loading required package: car
  Loading required package: carData
  > 
  > run.test.simulate.nlme <- Sys.info()[["user"]] == "fernandomiguez" && FALSE
  > 
  > if(run.test.simulate.nlme){
  + 
  +   data(Orange)
  + 
  +   fitL <- nlsList(circumference ~ SSlogis(age, Asym, xmid, scal), data = Orange)
  + 
  +   fmm <- nlme(fitL, random = pdDiag(Asym + xmid + scal ~ 1))
  + 
  +   ## This first simulation is at the level of population (fixed, level = 0)
  +   ## It also means that psim = 0, so there is no sampling from mvrnorm
  +   sim00 <- simulate_nlme(fmm, nsim = 100, psim = 0, level = 0, value = "data.frame")
  +   ## Level one means that we simlate at the level of individual Tree (or BLUP in this case)
  +   ## But psim = 0, so this is the same as predict(x, level = 0)
  +   sim01 <- simulate_nlme(fmm, nsim = 100, psim = 0, level = 1, value = "data.frame")
  + 
  +   ## What if I simulate?
  +   sim10 <- simulate_nlme(fmm, nsim = 100, psim = 1, level = 0, value = "data.frame")
  +   sim11 <- simulate_nlme(fmm, nsim = 100, psim = 1, level = 1, value = "data.frame")
  + 
  +   ## Level 2?
  +   ## Not sure if it is correct, but it seems to work as intended
  +   ## The line below does not make sense because we are adding error
  +   ## at the population level
  +   ## sim20 <- simulate_nlme(fmm, nsim = 100, psim = 2, level = 0)
  +   sim21 <- simulate_nlme(fmm, nsim = 100, psim = 2, level = 1, value = "data.frame")
  + 
  +   ggplot(data = sim00, aes(x = age, y = circumference)) + 
  +     geom_point() + 
  +     geom_line(aes(y = sim.y, group = ii)) + 
  +     ggtitle("psim = 0, level = 0")
  + 
  +   ggplot(data = sim01) + 
  +     geom_point(aes(x = age, y = circumference, color = Tree)) + 
  +     geom_line(aes(y = sim.y, x = age, color = Tree)) + 
  +     theme(legend.position = "none") + 
  +     ggtitle("psim = 0, level = 1")
  +     
  +   ggplot(data = sim10) + 
  +     geom_line(aes(y = sim.y, x = age, group = ii), color = "gray", alpha = 0.4) + 
  +     geom_point(aes(x = age, y = circumference)) + 
  +     theme(legend.position = "none") + 
  +     ggtitle("psim = 1, level = 0")
  + 
  +   sim11$Tree_ID <- with(sim11, paste0(Tree,"_",ii))
  +   
  +   ggplot(data = sim11) + 
  +     facet_wrap(~ Tree) + 
  +     geom_line(aes(y = sim.y, x = age, color = Tree, group = Tree_ID), alpha = 0.4) + 
  +     geom_point(aes(x = age, y = circumference)) + 
  +     theme(legend.position = "none") + 
  +     ggtitle("psim = 1, level = 1")
  +   
  +   ## Test with Soybean example
  +   data(Soybean)
  +   
  +   fmsL <- nlsList(SSlogis, Soybean)
  +   fmsm <- nlme(fmsL)  
  + 
  +   fxf <- fixef(fmsm)
  +   fmsm2 <- update(fmsm, fixed = list(Asym + xmid + scal ~ Variety),
  +                   start = c(fxf[1], 0, fxf[2], 0, fxf[3], 0),
  +                   random = pdDiag(Asym + xmid + scal ~ 1))
  +   
  +   fmsm.s <- simulate_nlme(fmsm2, nsim = 100, value = "data.frame")
  + 
  +   fmsm.s$ID <- with(fmsm.s, paste0(ii,"_",Plot))
  +   
  +   ggplot(data = fmsm.s, aes(x = Time, y = weight)) + 
  +     facet_grid(Variety  ~ Year) + 
  +     geom_line(aes(x = Time, y = sim.y, group = ID, color = Variety), alpha = 0.5) + 
  +     geom_point() 
  + 
  +   ## Test with CO2 example
  +   data(CO2)
  +   
  +   fmL <- nlsList(SSasymp, CO2)
  +   fmcm <- nlme(fmL, random = pdDiag(Asym + R0 ~ 1))
  +   
  +   fxf <- fixef(fmcm)
  +   fmcm2 <- update(fmcm, fixed = list(Asym + R0 + lrc ~ Treatment),
  +                   start = c(fxf[1], 0, fxf[2], 0, fxf[3], 0))
  + 
  +   fmcm2.s <- simulate_nlme(fmcm2, nsim = 100, value = "data.frame")
  +  
  +   fmcm2.s$ID <- with(fmcm2.s, paste0(ii,"_",Plant))
  +   
  +   ggplot(data = fmcm2.s, aes(x = conc, y = uptake)) + 
  +     facet_grid(Type ~ Treatment) + 
  +     geom_line(aes(x = conc, y = sim.y, group = ID, color = Treatment), alpha = 0.5) + 
  +     geom_point() 
  +   
  +   ## What about bootstrap?
  +   system.time(fmcm2.bt <- boot_nlme(fmcm2, cores = 4)) ## 35 seconds on my laptop MacBook Pro
  +   
  +   hist(fmcm2.bt, 2) 
  + 
  +   ## predictions
  +   ## Incorporate the fixed effect of type
  +   fxf2 <- fixef(fmcm2)
  +   fmcm3 <- update(fmcm2, fixed = list(Asym + R0 + lrc ~ Treatment + Type),
  +                   start = c(fxf2[1:2], 0, fxf2[3:4], 0, fxf2[5:6], 0),
  +                   weights = varPower())  
  +   
  +   fxf3 <- fixef(fmcm3)
  +   fmcm4 <- update(fmcm3, fixed = list(Asym + R0 + lrc ~ Treatment + Type + Treatment:Type),
  +                   start = c(fxf3[1:3], 0, fxf3[4:6], 0, fxf3[7:9], 0))
  +   
  +   prd_fun <- function(x) predict(x, level = 0)
  +   
  +   fmcm4.bt2 <- boot_nlme(fmcm4, prd_fun, cores = 4)
  +   
  +   CO2$lwr <- apply(t(fmcm4.bt2$t), 1, quantile, probs = 0.05, na.rm = TRUE)
  +   CO2$upr <- apply(t(fmcm4.bt2$t), 1, quantile, probs = 0.95, na.rm = TRUE)
  + 
  +   CO2$ID <- with(CO2, paste0(Treatment,"_",Type))
  +   
  +   ggplot(data = CO2, aes(x = conc, y = uptake)) + 
  +     facet_grid( Type ~ Treatment) + 
  +     geom_ribbon(aes(x = conc, ymin = lwr, ymax = upr, group = ID, fill = Treatment), alpha = 0.5) + 
  +     geom_point() 
  +   
  +   ## Trying the nested approach
  +   fmcm4 <- update(fmcm2, random = list(Type = pdDiag(Asym + R0 ~ 1),
  +                                        Plot = list(Asym + R0 ~ 1)),
  +                   groups = ~Type/Plant)
  +   
  +   prd0 <- predict(fmcm4, newdata = CO2, level = 1)
  + 
  +   fmcm4.s <- simulate_nlme(fmcm4, nsim = 100, value = "data.frame")
  +   ## This shows that while predict.nlme fails, simulate_nlme works 
  +   ## (or seems to work)
  +   
  +   fmcm4.s$ID <- with(fmcm4.s, paste0(ii,"_",Plant))
  +   
  +   ggplot(data = fmcm4.s, aes(x = conc, y = uptake)) + 
  +     facet_grid( Type ~ Treatment) + 
  +     geom_line(aes(x = conc, y = sim.y, group = ID, color = Treatment), alpha = 0.5) + 
  +     geom_point() 
  +   
  +   ## GAMS don't quite get it right
  +   library(mgcv)
  +   CO2$id <- as.factor(with(CO2, paste0(Type, "_", Treatment)))
  +   fit.gam <- gam(uptake ~ Type*Treatment + s(conc, k = 3, by = id), data = CO2)
  +   # 
  +   ggplot(data = CO2, aes(x = conc, y = uptake)) + 
  +      facet_grid( Type ~ Treatment) + 
  +      geom_line(aes(y = fitted(fit.gam))) + 
  +      geom_point()
  +   
  +   ## Compared to gnls model
  +   fit.lis <- nlsLMList(uptake ~ SSnrh(conc, asym, phi, theta, rd), data = CO2)
  +   
  +   plot(fit.lis)
  +   
  +   fmco2 <- nlme(fit.lis, random = pdDiag(asym + phi + rd ~ 1))
  +   
  +   fmco2.2 <- update(fmco2, fixed = asym + phi + theta + rd ~ Type,
  +                     start = c(fixef(fmco2)[1], 0, 
  +                               fixef(fmco2)[2], 0, 
  +                               fixef(fmco2)[3], 0,
  +                               fixef(fmco2)[4], 0))
  + 
  +   fmco2.3 <- update(fmco2, fixed = asym + phi + theta + rd ~ Type + Treatment,
  +                     start = c(fixef(fmco2.2)[1:2], 0, 
  +                               fixef(fmco2.2)[3:4], 0, 
  +                               fixef(fmco2.2)[5:6], 0,
  +                               fixef(fmco2.2)[7:8], 0))
  +   
  +   ## Turns out this model has too many parameters and they are not stable
  +   ## for this reason the predictions are all over the place
  +   fmco2.4 <- update(fmco2, fixed = asym + phi + theta + rd ~ Type + Treatment + Type:Treatment,
  +                     start = c(fixef(fmco2.3)[1:3], 0, 
  +                               fixef(fmco2.3)[4:6], 0, 
  +                               fixef(fmco2.3)[7:9], 0,
  +                               fixef(fmco2.3)[10:12], 0))
  +   
  +   prd <- predict_nlme(fmco2.4, interval = "conf", level = 0.9)
  +   CO2A <- cbind(CO2, prd)
  +   
  +   CO2AS <- subset(CO2A, Type == "Quebec")
  +   ggplot(data = CO2AS, aes(x = conc, y = uptake)) + 
  +     facet_grid( Type ~ Treatment) + 
  +     geom_line(aes(y = Estimate)) +
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "purple", alpha = 0.3) + 
  +     geom_point()
  +   
  +   ## A simpler model is better
  +   fmco2.5 <- update(fmco2.4, fixed = list(asym + phi ~ Type + Treatment + Type:Treatment,
  +                                         theta + rd ~ 1),
  +                     start = c(fixef(fmco2.3)[1:3], 0, 
  +                               fixef(fmco2.3)[4:6], 0, 
  +                               fixef(fmco2.3)[7],
  +                               fixef(fmco2.3)[10]))
  +   
  +   prd <- predict_nlme(fmco2.5, interval = "conf", level = 0.9)
  +   CO2A <- cbind(CO2, prd)
  +   
  +   ggplot(data = CO2A, aes(x = conc, y = uptake)) + 
  +     facet_grid( Type ~ Treatment) + 
  +     geom_line(aes(y = Estimate)) +
  +     geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "purple", alpha = 0.3) + 
  +     geom_point()
  +   
  +   fit.gnls <- gnls(uptake ~ SSnrh(conc, asym, phi, theta, rd), data = CO2,
  +                    params = list(asym + phi ~ Type*Treatment, theta + rd ~ 1),
  +                    start = c(35.35, 0, 0, 0, 0.1577, 0, 0, 0, 0.955, 2.35))
  + 
  +   ## This shows that a 'good' nls model is better than a 'GAM' but 'gams'
  +   ## are much, much more flexible
  +   IC_tab(fit.gnls, fit.gam)
  +   
  +   ## Testing the implementation of psim = 3 for nlme objects
  +   ## Picking up the previous Soybean model
  +   data(Soybean)
  +   
  +   fmsL <- nlsList(SSlogis, Soybean)
  +   fmsm <- nlme(fmsL)  
  +   
  +   fxf <- fixef(fmsm)
  +   fmsm2 <- update(fmsm, fixed = list(Asym + xmid + scal ~ Variety),
  +                   start = c(fxf[1], 0, fxf[2], 0, fxf[3], 0),
  +                   random = pdDiag(Asym + xmid + scal ~ 1))
  +   
  +   ## Ok, this seems to work
  +   sim2.psim2 <- simulate_nlme(fmsm2, nsim = 1, 
  +                               psim = 2, level = 1,
  +                               value = "data.frame")
  +   
  +   sim2.psim2$id <- with(sim2.psim2, paste(Plot, Variety, sep = "_"))
  +   
  +   ggplot(data = sim2.psim2, aes(x = Time, y = weight)) + 
  +     geom_line(aes(x = Time, y = sim.y, group = id), color = "gray", alpha = 0.4) + 
  +     geom_point() 
  +   
  +   ## psim = 3 is not specific for a given subject
  +   sim2.psim3 <- simulate_nlme(fmsm2, nsim = 10, 
  +                               psim = 3, level = 1,
  +                               value = "data.frame")
  + 
  +   sim2.psim3$id <- with(sim2.psim3, paste(Plot, Variety, ii, sep = "_"))
  +   
  +   ggplot(data = sim2.psim3, aes(x = Time, y = weight)) + 
  +     geom_line(aes(x = Time, y = sim.y, group = id), color = "gray", alpha = 0.3) + 
  +     geom_point() 
  +   
  +   ## Since psim = 3 is not specific to an individual subject
  +   ## we should not expect for the lines to vary around the observed data
  +   ## in each individual panel. If that is what you want then
  +   ## you need to use psim = 2
  +   ggplot(data = sim2.psim3, aes(x = Time, y = weight)) + 
  +     facet_wrap(~ Variety * Plot) + 
  +     geom_line(aes(x = Time, y = sim.y, group = ii), color = "gray", alpha = 0.3) + 
  +     geom_point() 
  +   
  +   fmsm3 <- update(fmsm2, random = list(Year = pdDiag(Asym + xmid + scal ~ 1),
  +                                        Plot = pdDiag(Asym + xmid + scal ~ 1)),
  +                   groups = ~Year/Plot)
  +   
  +   ## I won't check if this gives you the right answer at the moment, but it runs
  +   sim3.psim3 <- simulate_nlme(fmsm3, nsim = 10, psim = 3, value = "data.frame")
  +   
  +   sim3.psim3$id <- with(sim3.psim3, paste(Plot, Variety, Year, ii, sep = "_"))
  +   
  +   ggplot(data = sim3.psim3, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(x = Time, y = sim.y, group = id), alpha = 0.1)
  +   
  +   ggplot(data = sim3.psim3, aes(x = Time, y = weight)) + 
  +     facet_wrap(~ Year * Variety) + 
  +     geom_point() + 
  +     geom_line(aes(x = Time, y = sim.y, group = id), alpha = 0.1)
  + 
  +   varT <- var(Soybean$weight)  
  +   
  +   sim3.psim2.2 <- simulate_nlme(fmsm3, nsim = 2e3, psim = 2)
  +   sim3.psim3.2 <- simulate_nlme(fmsm3, nsim = 2e3, psim = 3)
  +   
  +   varT.psim2 <- apply(sim3.psim2.2, 2, var, na.rm = TRUE)  
  +   varT.psim3 <- apply(sim3.psim3.2, 2, var, na.rm = TRUE)  
  +   
  +   ggplot() + 
  +     geom_density(aes(x = varT.psim2, color = "psim = 2")) + 
  +     geom_density(aes(x = varT.psim3, color = "psim = 3")) + 
  +     geom_vline(xintercept = varT)
  +   ## Unexpectedly, this shows that the variance for psim = 3
  +   ## is lower than the variance for psim2, perhaps I need to do more simulations?
  +   ## Try a different dataset
  +   
  +   data(Orange)
  +   fitL <- nlsList(circumference ~ SSlogis(age, Asym, xmid, scal), data = Orange)
  +   fmm <- nlme(fitL, random = pdDiag(Asym + xmid + scal ~ 1))
  +   
  +   varT <- var(Orange$circumference)
  +   fmm.sim.psim2 <- simulate_nlme(fmm, nsim = 1e3, psim = 2)  
  +   fmm.sim.psim3 <- simulate_nlme(fmm, nsim = 1e3, psim = 3)  
  +  
  +   varT.psim2 <- apply(fmm.sim.psim2, 2, var, na.rm = TRUE)  
  +   varT.psim3 <- apply(fmm.sim.psim3, 2, var, na.rm = TRUE)  
  +   
  +   ggplot() + 
  +     geom_density(aes(x = varT.psim2, color = "psim = 2")) + 
  +     geom_density(aes(x = varT.psim3, color = "psim = 3")) + 
  +     geom_vline(xintercept = varT) 
  +   
  +   ## Prediction Interval for the mean of a NEW Tree
  +   nprd <- simulate_nlme(fmm, nsim = 1e3, psim = 3, value = "data.frame")  
  +   
  +   ## First step: remove the observation-level variation
  +   nprdA1 <- aggregate(sim.y ~ ii + age + Tree, data = nprd, FUN = mean)
  +   
  +   nprdA <- aggregate(sim.y ~ age, data = nprdA1, FUN = median)
  +   nprdA$lwr <- aggregate(sim.y ~ age, data = nprdA1, FUN = quantile, probs = 0.05)$sim.y
  +   nprdA$upr <- aggregate(sim.y ~ age, data = nprdA1, FUN = quantile, probs = 0.95)$sim.y
  + 
  +   ggplot() + 
  +     geom_point(data = Orange, aes(x = age, y = circumference, color = Tree)) + 
  +     geom_line(data = nprdA, aes(x = age, y = sim.y)) + 
  +     geom_ribbon(data = nprdA, aes(x = age, ymin = lwr, ymax = upr), alpha = 0.3) + 
  +     ggtitle("90% Prediction Interval for a NEW Tree regression")
  + }
  > 
  > proc.time()
     user  system elapsed 
     2.07    0.40    2.42 
Running the tests in 'tests/tests_SSall.R' failed.
Complete output:
  > ## 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 25.102151 25.43329
  ldt   4.571701  4.95934
  > 
  > ## 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 
     3.53    0.31    3.84 
Running the tests in 'tests/tests_SSother.R' failed.
Complete output:
  > ## Running other tests to see if the code actually works
  > 
  > test.other.examples <- Sys.info()[["user"]] == "fernandomiguez"
  > 
  > if(test.other.examples){
  + 
  +   require(nlraa)
  +   require(ggplot2)
  +   require(segmented)
  +   require(minpack.lm)
  +   require(car)
  +   require(nlstools)
  +   
  +   data(plant)
  +   
  +   ## For this first case I need to use nlsLM instead
  +   fit.rwc <- nlsLM(y ~ SSblin(time, a, b, xs, c), data = plant, subset = group == "RWC")
  +   fit.rkw <- nls(y ~ SSblin(time, a, b, xs, c), data = plant, subset = group == "RKW")
  +   fit.rkv <- nls(y ~ SSblin(time, a, b, xs, c), data = plant, subset = group == "RKV")
  +   
  +   fit.lm <- lm(y ~ time, data = subset(plant, group == "RWC"))
  +   fit.seg <- segmented(fit.lm)
  +   ## When comparing the outputs they differ in the interpretation of the
  +   ## second coefficient for SSblin it is for (x - xs) and for 'segmented'
  +   ## it must be with 'x' as a predictor, but I'm not sure
  +   
  +   ## plot for RWC
  +   ggplot(data = subset(plant, group == "RWC"), aes(x = time, y = y)) + 
  +     geom_point() + geom_line(aes(y = fitted(fit.rwc)))
  +   ## plot for RKW
  +   ggplot(data = subset(plant, group == "RKW"), aes(x = time, y = y)) + 
  +     geom_point() + geom_line(aes(y = fitted(fit.rkw)))
  +   ## plot for RKV
  +   ggplot(data = subset(plant, group == "RKV"), aes(x = time, y = y)) + 
  +     geom_point() + geom_line(aes(y = fitted(fit.rkv)))
  +   
  +   ## Data 'down'
  +   ## data(down)
  +   ## fit <- nls(births ~ SSricker(age, a, b), data = down)
  +   
  +   ## ggplot(data = down, aes(age, births)) + 
  +     ## geom_point() + 
  +     ## geom_line(aes(y = fitted(fit)))
  +   
  +   ## Testing with the swpg data
  +   data(swpg)
  +   ## blinear
  +   fit1 <- nls(lfgr ~ SSblin(ftsw, a, b, xs, c), data = swpg)
  +   ggplot(data = swpg, aes(x = ftsw, y = lfgr)) + 
  +          geom_point() + 
  +          geom_line(aes(y = fitted(fit1)))
  +   ## linear plateau
  +   fit2 <- nls(lfgr ~ SSlinp(ftsw, a, b, xs), data = swpg)
  +   ggplot(data = swpg, aes(x = ftsw, y = lfgr)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit2)))
  +   
  +   ## Some tests with datasets from package nlme
  +   require(nlme)
  +   data(Relaxin)
  +   fit.nlis <- nlsList(cAMP ~ SSquadp(conc, a, b, c, xs), data = Relaxin)
  +   fit.nlme <- nlme(fit.nlis, random = pdDiag(b + xs ~ 1))
  +   
  +   data(Fatigue)
  +   fit.nlis <- nlsList(relLength ~ SSexpf(cycles, a, c), data = Fatigue)
  +   fit.nlme <- nlme(fit.nlis, random = pdDiag(a + c ~ 1))
  +   
  +   data(Cefamandole)
  +   fit.nlis <- nlsList(conc ~ SSricker(Time, a, b), data = Cefamandole)
  +   fit.nlme <- nlme(fit.nlis, random = pdDiag(a + b ~ 1))
  +   
  +   data(Soybean)
  +   ## Default model
  +   fit.nlis <- nlsList(weight ~ SSbgf(Time, w.max, t.e, t.m), data = Soybean)
  +   ## The previous results in 7 warnings
  +   fit.soy.1 <- nlme(fit.nlis, random = pdDiag(w.max ~ 1))
  +   fit.nlis <- nlsList(weight ~ SSbgrp(Time, w.max, lt.e, ldt), data = Soybean)
  +   ## This results in just one
  +   fit.soy.2 <- nlme(fit.nlis, random = pdDiag(w.max + lt.e ~ 1))
  +   ## Four parameter
  +   fit.nlis <- nlsList(weight ~ SSbgf4(Time, w.max, t.e, t.m, t.b), data = Soybean)
  +   fit.soy.3 <- nlme(fit.nlis, random = pdDiag(w.max  ~ 1))
  +   ## Four parameter - reparameterized
  +   fit.nlis <- nlsList(weight ~ SSbg4rp(Time, w.max, lt.e, ldtm, ldtb), data = Soybean)
  +   fit.soy.4 <- nlme(fit.nlis, random = pdDiag(w.max + lt.e ~ 1))
  +   anova(fit.soy.1, fit.soy.2, fit.soy.3, fit.soy.4)
  +   ## The second model is better, but also the reparameterized versions
  +   ## allow for including the lt.e random effect which improves the fit
  +   
  +   ## Testing the behavior of the ratio function
  +   
  +   ## Changing a
  +   xx <- 1:100
  +   y1 <- ratio(xx, 1, 1, 1, 1)
  +   y2 <- ratio(xx, 0.5, 1, 1, 1)
  +   dat <- data.frame(xx = xx, y1 = y1, y2 = y2)
  +   ggplot(data = dat) + 
  +     geom_point(aes(x = xx, y = y1), color = "red") +
  +     geom_point(aes(x = xx, y = y2), color = "blue")
  +   ## If I decrease a, y decreases
  +   ## It makes sense because a is in the numerator
  +   
  +   ## Changing b
  +   y1 <- ratio(xx, 1, 1, 1, 1)
  +   y2 <- ratio(xx, 1, 0.5, 1, 1)
  +   dat <- data.frame(xx = xx, y1 = y1, y2 = y2)
  +   ggplot(data = dat) + 
  +     geom_point(aes(x = xx, y = y1), color = "red") +
  +     geom_point(aes(x = xx, y = y2), color = "blue")
  +   ## If I decrease b, y increases
  +   ## It makes sense because b is in the denominator
  +   ## When c and d = 1, then a/b is the asymptote
  +   
  +   ## Changing c
  +   y1 <- ratio(xx, 1, 1, 1, 1)
  +   y2 <- ratio(xx, 1, 1, 0.5, 1)
  +   dat <- data.frame(xx = xx, y1 = y1, y2 = y2)
  +   ggplot(data = dat) + 
  +     geom_point(aes(x = xx, y = y1), color = "red") +
  +     geom_point(aes(x = xx, y = y2), color = "blue")
  +   ## This turns the equation into an exponential decay
  +   
  +   ## Changing d, 1
  +   y1 <- ratio(xx, 1, 1, 1.1, 1.3)
  +   y2 <- ratio(xx, 1, 1, 1.1, 1.31)
  +   dat <- data.frame(xx = xx, y1 = y1, y2 = y2)
  +   ggplot(data = dat) + 
  +     geom_line(aes(x = xx, y = y1), color = "red") +
  +     geom_line(aes(x = xx, y = y2), color = "blue")
  +   ## Changing d, 2
  +   y1 <- ratio(xx, 1, 0.5, 1.1, 1.3)
  +   y2 <- ratio(xx, 1, 0.5, 1.1, 1.31)
  +   dat <- data.frame(xx = xx, y1 = y1, y2 = y2)
  +   ggplot(data = dat) + 
  +     geom_line(aes(x = xx, y = y1), color = "red") +
  +     geom_line(aes(x = xx, y = y2), color = "blue")
  +   
  +   ## Testing SSratio
  +   data(Theoph)
  +   Theoph.10 <- subset(Theoph, Subject == 10)
  +   fit.10 <- nlsLM(conc ~ SSratio(Time, a, b, c, d), data = Theoph.10)
  +   dat <- data.frame(x = Theoph.10$Time, y = Theoph.10$conc)
  +   ggplot(data = dat, aes(x = x, y = y)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit.10)))
  +   
  +   ## Dataset ChickWeight is an interesting one to analyze
  +   data(ChickWeight)
  +   chick.11 <- subset(ChickWeight, Chick == 11)
  +   fit.11.b <- nls(weight ~ bgf2(Time, w.max, w.b = 43, t.e, t.m, t.b = 0), data = chick.11, start = list(w.max = 170, t.e = 17, t.m = 10))
  +   fit.11.l <- nls(weight ~ SSfpl(Time, A, B, xmid, scal), data = chick.11)
  +   anova(fit.11.b, fit.11.l) ## RSS is lower for bgf2
  +   AIC(fit.11.l, fit.11.b) ## AIC is lower for bgf2
  +   
  +   ggplot(data = chick.11, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit.11.b)))
  +   
  +   ggplot(data = chick.11, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit.11.l)))
  +   
  +   chick.43 <- subset(ChickWeight, Chick == 43)
  +   fit.43.b <- nls(weight ~ SSbgf(Time, w.max, t.e, t.m), data = chick.43)
  +   fit.43.brp <- nls(weight ~ SSbgrp(Time, w.max, t.e, t.m), data = chick.43, subset = Time > 0)
  +   fit.43.b4rp <- nls(weight ~ SSbg4rp(Time, w.max, lt.e, ldtm, ldtb), data = chick.43)
  +   
  +   ggplot(data = chick.43, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit.43.b)))
  +   
  +   ggplot(data = subset(chick.43, Time > 0), aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit.43.brp)))
  +   
  +   ggplot(data = chick.43, aes(x = Time, y = weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit.43.b4rp)))
  +   
  +   ## Looking at dataset 'Indometh'
  +   data(Indometh)
  +   
  +   fit <- nlsLM(conc ~ SSratio(time, a, b, c, d), data = Indometh)
  +   
  +   ggplot(data = Indometh, aes(x = time, y = conc)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   fit.nlis <- nlsLMList(conc ~ SSratio(time, a, b, c, d), data = Indometh)  
  +   ## Can't fit a reasonable model using nlme...
  +   
  +   ## Looking at 'Orange' dataset
  +   data(Orange)
  +   
  +   fit <- nlsList(circumference ~ SSbg4rp(age, w.max, lt.e, ldtm , ldtb), data = Orange)
  +   
  +   fit.nlme2 <- nlme(fit, random = list(w.max  ~ 1))
  +   
  +   plot(augPred(fit.nlme2, level = 0:1))
  +   
  +   ## Looking at dataset 'pressure'
  +   data(pressure)
  +   
  +   fit <- nls(pressure ~ SSexpf(temperature, a, c), data = pressure)
  +   
  +   ggplot(data = pressure, aes(x = temperature, y = pressure)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   ## Looking at dataset 'Loblloly'
  +   data(Loblolly)
  +   
  +   fit <- nlsLM(height ~ SSratio(age, a, b, c, d), data = Loblolly)
  +   
  +   ggplot(data = Loblolly, aes(x = age, y = height)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   fit <- nlsLM(height ~ SSblin(age, a, b, xs, c), data = Loblolly)
  +   
  +   ggplot(data = Loblolly, aes(x = age, y = height)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   ## Should create some examples from pacakge 'nlstools'
  +   require(nlstools)
  +   
  +   data(O2K)
  +   fit <- nls(VO2 ~ SSpquad(t, a, xs, b, c), data = O2K)
  +   
  +   plotfit(fit)
  +   
  +   ggplot(data = O2K, aes(x = t, y = VO2)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   data(vmkm)
  +   
  +   fit <- nlsLM(v ~ SSblin(S, a, b, xs, c), data = vmkm)
  +   
  +   ggplot(data = vmkm, aes(x = S, y = v)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   ## Let's review the datasets in NISTnls
  +   library(NISTnls)
  +   
  +   ## Bennett5 - not interested
  +   data(Chwirut1)
  +   fit <- nls(y ~ SSexpfp(x, a, c, d), data = Chwirut1)
  +   ggplot(data = Chwirut1, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   data(Chwirut2)
  +   fit <- nls(y ~ SSexpfp(x, a, c, d), data = Chwirut2)
  +   ggplot(data = Chwirut2, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   ## DanielWood - not interested
  +   ## ENSO - not interested
  +   data(Eckerle4)
  +   fit <- nls(y ~ SSbell(x, ymax, a, b, xc), data = Eckerle4)
  +   ggplot(data = Eckerle4, aes(x, y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   ## Gauss1, Gauss2, Gauss3 - nlraa cannot handle this!
  +   ## should use splines or gams instead
  +   
  +   data(Hahn1)
  +   fit <- nlsLM(y ~ SSratio(x, a, b, c, d), data = Hahn1, control = list(maxiter = 1e3))
  +   ggplot(data = Hahn1, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   ## Almost perfect fit, but there are complains
  +   fit2 <- nlsLM(y ~ SSratio(x, a, b, c, d), data = Hahn1, start = coef(fit))
  +   ggplot(data = Hahn1, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit2)))
  +   
  +   ## Kirby2 - maybe, but there is not much of an error
  +   ## Lanczos1, 2, and 3 - not too interesting
  +   data(MGH09)
  +   fit <- nls(y ~ SSblin(x, a, b, xs, c), data = MGH09)
  +   ggplot(data = MGH09, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   ## MGH10 - not interested
  +   ## MGH17 - challenging, but do not have a good candidate
  +   ## Misra1a, b, c and d - not interested
  +   ## Nelson - not sure what this one is about
  +   ## Ratkowsky2 and 3 - look like good candidates for logistic or Gompertz
  +   data(Ratkowsky2)
  +   fit1 <- nls(y ~ SSlogis(x, Asym, xmid, scal), data = Ratkowsky2)
  +   ggplot(data = Ratkowsky2, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit1)))
  +   
  +   ## What about a 5 parameter logistic?
  +   fit2 <- nlsLM(y ~ SSlogis5(x, asym1, asym2, xmid, iscal, theta), data = Ratkowsky2, control = list(maxiter = 1e3))
  +   ggplot(data = Ratkowsky2, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit2)))
  +   ## Comparing the models
  +   anova(fit1, fit2)
  +   ## fit2.bt <- Boot(fit2) It takes too long to run
  +   ## hist(fit2.bt)
  +   ## These shows that the parameters are not well constrained under this model
  +   
  +   data("Ratkowsky3")
  +   fit1 <- nls(y ~ SSlogis(x, asym, xmid, scal), data = Ratkowsky3)
  +   ggplot(data = Ratkowsky3, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit1)))
  +   
  +   fit2 <- nlsLM(y ~ SSlogis5(x, asym1, asym2, xmid, iscal, theta), data = Ratkowsky3, control = list(maxiter = 1e3))
  +   fit3 <- nls(y ~ SSgompertz(x, a, b, c), data = Ratkowsky3)
  +   ggplot(data = Ratkowsky3, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit2)))
  +   ggplot(data = Ratkowsky3, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit3)))
  +   
  +   ## Five-parameter logistic fits a little bit better
  +   anova(fit1, fit2, fit3)
  +   ## fit2.bt <- Boot(fit2) takes too long to run
  +   ## hist(fit2.bt)
  +   
  +   data(Roszman1)
  +   fit <- nls(y ~ SSexpf(x, a, c), data = Roszman1)
  +   ggplot(data = Roszman1, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   ## Blinear is not actually a terrible fit
  +   fit <- nls(y ~ SSblin(x, a, b, xs, c), data = Roszman1)
  +   ggplot(data = Roszman1, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   data(Thurber)
  +   ## Several possible models
  +   fit <- nls(y ~ SSfpl(x, a, b, c, d), data = Thurber)
  +   ggplot(data = Thurber, aes(x = x, y = y)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   ## Another idea is to show and compare the car::Boot and nlstools:nlsBoot functions
  +   ## These could be very useful, especially when the confidence interval function 'confint'
  +   ## fails
  +   
  +   ## Datasets from 'nlstools'
  +   data(L.minor)
  +   ## Blinear
  +   fit1 <- nls(rate ~ SSblin(conc, a, b, xs1, c), data = L.minor)
  +   ## MicMen
  +   fit2 <- nls(rate ~ SSmicmen(conc, Vm, K), data = L.minor)
  +   ## rational with minpack.lm
  +   fit3 <- nlsLM(rate ~ SSratio(conc, a, b, c, d), data = L.minor)
  +   
  +   ggplot(data = L.minor, aes(conc, rate)) + geom_point() + 
  +     geom_line(aes(y = fitted(fit1), color = "blinear")) + 
  +     geom_line(aes(y = fitted(fit2), color = "MicMen")) +
  +     geom_line(aes(y = fitted(fit3), color = "rational"))
  +   
  +   ## Oxygen uptake
  +   data("O2K")
  +   fit <- nls(VO2 ~ SSpquad(t, a, xs, b, c), data = O2K)
  +   ggplot(data = O2K, aes(x = t, y = VO2)) + geom_point() + geom_line(aes(y = fitted(fit)))
  +   
  +   ## nlstool::vmkm
  +   data(vmkm)
  +   fit1 <- nls(v ~ SSmicmen(S, Vm, K), data = vmkm)
  +   ggplot(data = vmkm, aes(x = S, y = v)) + geom_point() + geom_line(aes(y = fitted(fit1)))
  +   
  +   fit2 <- nls(v ~ SSasymp(S, Asym, R0, lrc), data = vmkm)
  +   ggplot(data = vmkm, aes(x = S, y = v)) + geom_point() + geom_line(aes(y = fitted(fit2)))
  +   
  +   ## Testing function SShill3
  +   ## With Soybean
  +   fit <- nls(weight ~ SShill3(Time, Ka, n, a), data = Soybean)
  +   
  +   ggplot(data = Soybean, aes(Time, weight)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   ## With Indometh
  +   fit <- nls(conc ~ SShill3(time, Ka, n, a), data = Indometh)
  +   
  +   ggplot(data = Indometh, aes(x = time, y = conc)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   ## With Loblolly
  +   fit <- nls(height ~ SShill3(age, Ka, n, a), data = Loblolly)
  +   
  +   ggplot(data = Loblolly, aes(x = age, y = height)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit)))
  +   
  +   ## With Thurber
  +   Thurber$x2 <- Thurber$x + 6
  +   fit <- nlsLM(y ~ SShill3(x2, Ka, n, a), data = Thurber)
  + 
  +   ggplot(data = Thurber, aes(x = x2, y = y)) + 
  +     geom_point() + 
  +     geom_line(aes(y = fitted(fit))) 
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.51    0.03    0.53