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.

> ## 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 
   1.34    0.20    1.53