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 simulate_gnls
> 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.simulate.gnls.test <- Sys.info()[["user"]] == "fernandomiguez" && FALSE
> 
> if(run.simulate.gnls.test){
+   
+   set.seed(101)
+   
+   data(Orange)
+   
+   fitg <- gnls(circumference ~ SSlogis(age, Asym, xmid, scal), data = Orange, weights = varPower())
+   
+   system.time(fitg.bt0 <- boot_nlme(fitg, fitted, psim = 0, cores = 4))
+   
+   lwr0.q <- apply(t(fitg.bt0$t), 1, quantile, probs = 0.05, na.rm = TRUE)
+   upr0.q <- apply(t(fitg.bt0$t), 1, quantile, probs = 0.95, na.rm = TRUE)
+  
+   system.time(fitg.bt1 <- boot_nlme(fitg, fitted, psim = 1, cores = 4))
+   
+   lwr1.q <- apply(t(fitg.bt1$t), 1, quantile, probs = 0.05, na.rm = TRUE)
+   upr1.q <- apply(t(fitg.bt1$t), 1, quantile, probs = 0.95, na.rm = TRUE)
+   
+   ggplot() + 
+     geom_point(data = Orange, aes(x = age, y = circumference)) + 
+     geom_line(data = Orange, aes(x = age, y = fitted(fitg))) + 
+     geom_ribbon(aes(x = Orange$age, ymin = lwr0.q, ymax = upr0.q), fill = "blue", alpha = 0.2) + 
+     geom_ribbon(aes(x = Orange$age, ymin = lwr1.q, ymax = upr1.q), fill = "purple", alpha = 0.2)
+   
+   ## What about the model coefficients?
+   system.time(fitg.cfs.bt0 <- boot_nlme(fitg, psim = 0, cores = 4))
+   system.time(fitg.cfs.bt1 <- boot_nlme(fitg, psim = 1, cores = 4)) 
+ 
+   summary(fitg.cfs.bt0)
+   summary(fitg.cfs.bt1)
+   
+   ## Somewhat wider using bootstrap
+   ## Wider when simulating the parameters
+   confint(fitg.cfs.bt0)
+   confint(fitg.cfs.bt1)
+   intervals(fitg)
+   
+   hist(fitg.cfs.bt0)
+   hist(fitg.cfs.bt1)
+   
+   sims0 <- simulate_gnls(fitg, psim = 0) ## Fitted values
+   sims1 <- simulate_gnls(fitg, psim = 1) ## One simulation from the mean or population-level
+   sims2 <- simulate_gnls(fitg, psim = 2) ## One simulation from the individual-level
+   
+   Orange2 <- Orange
+   Orange2$sims0 <- sims0
+   Orange2$sims1 <- sims1
+   Orange2$sims2 <- sims2
+   
+   ggplot() + 
+     geom_point(data = Orange2, aes(x = age, y = circumference)) + 
+     geom_line(data = Orange2, aes(x = age, y = sims0)) + 
+     geom_point(data = Orange2, aes(x = age, y = sims1), color = "red") +
+     geom_point(data = Orange2, aes(x = age, y = sims2), color = "green3")  
+     
+   ## What about a nlme?
+   
+   fitL <- nlsList(circumference ~ SSlogis(age, Asym, xmid, scal), data = Orange)
+   
+   fit.nm <- nlme(fitL, random = pdDiag(Asym + xmid + scal ~ 1))  
+   
+   plot(augPred(fit.nm, level = 0:1))
+   
+   ## Let's look at intervals first
+   system.time(fit.nm.bt <- boot_nlme(fit.nm, cores = 4))
+ 
+   confint(fit.nm.bt)  
+   intervals(fit.nm, which = "fixed")
+   
+   ## Fitted values
+   prdf <- function(x) predict(x, level = 0, newdata = data.frame(age = 50:1600)) 
+   
+   system.time(fit.nm.bt.f0 <- boot_nlme(fit.nm, prdf, psim = 0, cores = 4))
+  
+   lwr0.q <- apply(t(fit.nm.bt.f0$t), 1, quantile, probs = 0.05, na.rm = TRUE)
+   upr0.q <- apply(t(fit.nm.bt.f0$t), 1, quantile, probs = 0.95, na.rm = TRUE)
+  
+   Orange$fttd <- predict(fit.nm, level = 0)
+  
+   system.time(fit.nm.bt.f1 <- boot_nlme(fit.nm, prdf, psim = 1, cores = 4))
+   
+   lwr1.q <- apply(t(fit.nm.bt.f1$t), 1, quantile, probs = 0.05, na.rm = TRUE)
+   upr1.q <- apply(t(fit.nm.bt.f1$t), 1, quantile, probs = 0.95, na.rm = TRUE)
+   
+   ggplot() + 
+     geom_point(data = Orange, aes(x = age, y = circumference)) + 
+     geom_line(data = Orange, aes(x = age, y = fttd)) + 
+     geom_ribbon(aes(x = 50:1600, ymin = lwr0.q, ymax = upr0.q), fill = "blue", alpha = 0.2) + 
+     geom_ribbon(aes(x = 50:1600, ymin = lwr1.q, ymax = upr1.q), fill = "purple", alpha = 0.2)
+     
+   ## Testing the feature of sampling correlated errors
+   ## The way to really test this would be to look at the correlation
+   ## structure in the new simulations. This is not a good dataset for this
+   ## Need to check in Pinheiro and Bates for a good example
+   data("ChickWeight")
+   
+   fitcw <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal), data = ChickWeight, 
+                 correlation = corCAR1(form = ~ Time | Chick),
+                 weights = varPower())
+   
+   fitcw.vc <- var_cov(fitcw) 
+   image(fitcw.vc[,ncol(fitcw.vc):1], xaxt = "n", yaxt = "n")
+   fitcw.vc2 <- fitcw.vc[1:36,1:36]
+   image(log(fitcw.vc2[,ncol(fitcw.vc2):1]), xaxt = "n", yaxt = "n")
+   ## In this case, this is much slower, because of the huge matrix involved
+   fitcw.sim <- simulate_nlme(fitcw, nsim = 100, psim = 2, value = "data.frame")
+  
+   fitcw.sim$Chick_ID <- with(fitcw.sim, paste0(Chick,"_",ii))
+   
+   ## It looks like this approach does not capture the specific
+   ## Chick-level effect, but it does seem to produce
+   ## correlated lines.
+   ## Perhaps I cannot expect to reproduce the original data for this 'gnls'
+   ## object, but rather, I should check if it creates correlated errors
+   ## similar in structure to the ones in the original dataset
+   ggplot(data = fitcw.sim) + 
+     facet_wrap(~ Chick) + 
+     geom_line(aes(x = Time, y = sim.y, color = Chick, group = Chick_ID)) + 
+     geom_point(aes(x = Time, y = weight)) +
+     theme(legend.position = "none")
+ }
> 
> proc.time()
   user  system elapsed 
   1.31    0.21    1.53