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 of a bad model vs. uncertainty vs. model averaging
> require(nlraa)
Loading required package: nlraa
> require(car)
Loading required package: car
Loading required package: carData
> require(ggplot2)
Loading required package: ggplot2
> 
> run.predict.nls <- Sys.info()[["user"]] == "fernandomiguez" && FALSE
> 
> if(run.predict.nls){
+  
+   data(barley, package = "nlraa")
+   
+   ggplot(data = barley, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     xlab("NF (g/m2)") + ylab("Yield (g/m2)") + 
+     ggtitle("Barley yield response to N fertilizer")
+   
+   ## This is not a 'good' model but we'll go with it
+   fm.LP <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley)
+   
+   sim.LP <- simulate_nls(fm.LP, nsim = 1e3)
+   
+   ## Does predict work for a single model?
+   prd.LP <- predict_nls(fm.LP)
+   prd.LP.ci <- predict_nls(fm.LP, interval = "confidence")
+   prd.LP.pi <- predict_nls(fm.LP, interval = "prediction")
+   
+   ggplot(data = barley, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     geom_line(aes(y = fitted(fm.LP))) + 
+     geom_vline(xintercept = coef(fm.LP)[3]) + 
+     xlab("NF (g/m2)") + ylab("Yield (g/m2)") + 
+     ggtitle("Linear-plateau fit with break-point")
+   
+   barleyA <- cbind(barley, summary_simulate(sim.LP, probs = c(0.05, 0.95)))
+   
+   fm.LP.bt <- boot_nls(fm.LP) ## Bootstrap
+   fm.LP.bt.ci <- confint(fm.LP.bt) ## Bootstrap CI
+   fm.LP.ci <- confint(fm.LP) ## Profiled CI
+   
+   ggplot(data = barleyA, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     geom_line(aes(y = fitted(fm.LP))) + 
+     geom_ribbon(aes(ymin = Q5, ymax = Q95), alpha = 0.3, fill = "purple") + 
+     geom_vline(xintercept = fm.LP.bt$t0[3]) +
+     geom_errorbarh(aes(y = 100, xmin = fm.LP.ci[3,1], xmax = fm.LP.ci[3,2], 
+                        color = "profiled"), color = "blue") + 
+     geom_errorbarh(aes(y = 50, xmin = fm.LP.bt.ci[3,1], xmax = fm.LP.bt.ci[3,2], 
+                        color = "bootstrap"), color = "purple") + 
+     geom_text(aes(x = 13, y = 100, label = "profiled"), color = "blue") +  
+     geom_text(aes(x = 13, y = 50, label = "bootstrap"), color = "purple") +
+     xlab("NF (g/m2)") + ylab("Yield (g/m2)") + 
+     ggtitle("90% uncertainty bands and intervals for the break-point")
+   
+   ## What if we fit several models?
+   fm.L <- lm(yield ~ NF, data = barley)
+   fm.Q <- lm(yield ~ NF + I(NF^2), data = barley)
+   fm.A <- nls(yield ~ SSasymp(NF, Asym, R0, lrc), data = barley)
+   fm.BL <- nls(yield ~ SSblin(NF, a, b, xs, c), data = barley)
+   
+   print(IC_tab(fm.L, fm.Q, fm.A, fm.LP, fm.BL), digits = 2)
+   
+   ggplot(data = barley, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     geom_line(aes(y = fitted(fm.L), color = "Linear")) +
+     geom_line(aes(y = fitted(fm.Q), color = "Quadratic")) +
+     geom_line(aes(y = fitted(fm.A), color = "Asymptotic")) +  
+     geom_line(aes(y = fitted(fm.LP), color = "Linear-plateau")) + 
+     geom_line(aes(y = fitted(fm.BL), color = "Bi-linear")) + 
+     xlab("NF (g/m2)") + ylab("Yield (g/m2)") + 
+     ggtitle("Different model fits")
+   
+   ## Each model prediction is weighted using the AIC values
+   prd <- predict_nls(fm.L, fm.Q, fm.A, fm.LP, fm.BL)
+   prdc <- predict_nls(fm.L, fm.Q, fm.A, fm.LP, fm.BL, interval = "confidence")
+   prdp <- predict_nls(fm.L, fm.Q, fm.A, fm.LP, fm.BL, interval = "prediction")
+   
+   ggplot(data = barley, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     geom_line(aes(y = fitted(fm.L), color = "Linear")) +
+     geom_line(aes(y = fitted(fm.Q), color = "Quadratic")) +
+     geom_line(aes(y = fitted(fm.A), color = "Asymptotic")) +  
+     geom_line(aes(y = fitted(fm.LP), color = "Linear-plateau")) + 
+     geom_line(aes(y = fitted(fm.BL), color = "Bi-linear")) + 
+     geom_line(aes(y = prd, color = "Avg. Model"), size = 1.2, color = "black") + 
+     xlab("NF (g/m2)") + ylab("Yield (g/m2)") + 
+     ggtitle("Different model fits and average model weighted by AIC")
+   
+   ggplot(data = barley, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     geom_line(aes(y = fitted(fm.L), color = "Linear")) +
+     geom_line(aes(y = fitted(fm.Q), color = "Quadratic")) +
+     geom_line(aes(y = fitted(fm.A), color = "Asymptotic")) +  
+     geom_line(aes(y = fitted(fm.LP), color = "Linear-plateau")) + 
+     geom_line(aes(y = fitted(fm.BL), color = "Bi-linear")) + 
+     geom_line(aes(y = prd, color = "Avg. Model"), size = 1.2, color = "black") + 
+     geom_ribbon(aes(ymin = prdc[,3], ymax = prdc[,4]), 
+                 fill = "purple", alpha = 0.3) + 
+     geom_ribbon(aes(ymin = prdp[,3], ymax = prdp[,4]), 
+                 fill = "purple", alpha = 0.1) + 
+     xlab("NF (g/m2)") + ylab("Yield (g/m2)") + 
+     ggtitle("Model fits, 90% uncertainty bands for confidence and prediction")
+ 
+   ## Do GAMs work?
+   require(mgcv)
+   
+   fm.L <- lm(yield ~ NF, data = barley)
+   fm.Q <- lm(yield ~ NF + I(NF^2), data = barley)
+   fm.C <- lm(yield ~ NF + I(NF^2) + I(NF^3), data = barley)
+   fm.A <- nls(yield ~ SSasymp(NF, Asym, R0, lrc), data = barley)
+   fm.LP <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley)
+   fm.G <- gam(yield ~ NF + s(NF, k = 3), data = barley)
+ 
+   fm.Gs <- simulate_lm(fm.G, nsim = 1e3)
+   fm.Gss <- summary_simulate(fm.Gs, probs = c(0.05, 0.95))
+   barleyAS <- cbind(barley, fm.Gss)
+   
+   ## The default predict method for GAMs does not produce intervals
+   ## But we can generate them
+   fm.Gp <- predict(fm.G, se.fit = TRUE)
+   qnt <- qt(0.05, 72)
+   fm.Gpd <- data.frame(prd = fm.Gp$fit, 
+                        lwr = fm.Gp$fit + qnt * fm.Gp$se.fit,
+                        upr = fm.Gp$fit - qnt * fm.Gp$se.fit)
+   ## These intervals are almost exactly the same as the ones
+   ## obtained through simulation
+   
+   print(IC_tab(fm.L, fm.Q, fm.C, fm.A, fm.LP, fm.G), digits = 2)
+   
+   fm.prd <- predict_nls(fm.L, fm.Q, fm.C, fm.A, fm.LP, fm.G)
+   
+   ggplot(data = barleyAS, aes(x = NF, y = yield)) + 
+     geom_point() + 
+     geom_line(aes(y = fitted(fm.G), color = "gam")) + 
+     geom_line(aes(y = fitted(fm.C), color = "cubic")) + 
+     geom_line(aes(y = Estimate, color = "simulate_lm")) + 
+     geom_line(aes(y = fm.prd, color = "Avg. Model")) + 
+     geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "purple", alpha = 0.3) +
+     ggtitle("90% bands based on simulation")
+ 
+ }
> 
> 
> proc.time()
   user  system elapsed 
   1.29    0.26    1.54