# Loading package and dataset library("nlstools") set.seed(405) # Define the nonlinear model formulaExp <- as.formula(VO2 ~ (t <= 5.883) * VO2rest + (t > 5.883) * (VO2rest + (VO2peak - VO2rest) * (1 - exp(-(t - 5.883) / mu)))) # Graphical assessment of starting values preview(formulaExp, data = O2K, start = list(VO2rest = 400, VO2peak = 1600, mu = 1)) # figure 2 # Model fit, overview of the results, # plot of the fitted model superimposed to the data O2K.nls1 <- nls(formulaExp, start = list(VO2rest = 400, VO2peak = 1600, mu = 1), data = O2K) overview(O2K.nls1) plotfit(O2K.nls1, smooth = TRUE) # Figure 3 # Residuals O2K.res1 <- nlsResiduals(O2K.nls1) plot(O2K.res1) # Figure 4 test.nlsResiduals(O2K.res1) # Confidence regions O2K.cont1 <- nlsContourRSS(O2K.nls1) plot(O2K.cont1, col = FALSE, nlev = 5) # Figure 5, left panel O2K.conf1 <- nlsConfRegions(O2K.nls1, exp = 2, length = 2000) plot(O2K.conf1, bounds = TRUE) # Figure 5, right panel # Jackknife O2K.jack1 <- nlsJack(O2K.nls1) summary(O2K.jack1) plot(O2K.jack1) # Figure 6, left panel # Bootstrap O2K.boot1 <- nlsBoot(O2K.nls1) summary(O2K.boot1) plot(O2K.boot1, type = "boxplot") # Figure 6, right panel # Comparison jackknife / bootstrap / asymptotic confidence intervals (Figure 7) esti <- summary(O2K.nls1)$parameters[, "Estimate"] ster <- summary(O2K.nls1)$parameters[, "Std. Error"] t95 <- qt(0.975, df = (nrow(O2K) - 3)) binf <- esti - t95 * ster bsup <- esti + t95 * ster x11(height = 5, width = 10) par(mfrow = c(1, 3), mar = c(10,5,5,1)) plot(1:3, c(coef(O2K.nls1)[1], O2K.jack1$estijack[1, 1], O2K.boot1$bootCI[1, 1]), ylim = c(O2K.boot1$bootCI[1, 2]*.96, O2K.boot1$bootCI[1, 3]*1.04), xlim = c(0,4), xaxt = "n", xlab = "", ylab = expression(paste(VO[2], "rest (L)")), main = expression(paste(VO[2], "rest")), pch = 19, cex.axis = 1.5, cex.lab = 1.5, cex.main = 1.5) segments(x0 = 1:3, y0 = c(binf[1], O2K.jack1$jackCI[1, 2], O2K.boot1$bootCI[1, 2]), x1 = 1:3, y1 = c(bsup[1], O2K.jack1$jackCI[1, 3], O2K.boot1$bootCI[1, 3]), lty = 1:3) axis(1, at = 1:3, labels = c("t-based", "Jackknife", "Bootstrap"), las= 3, cex.axis = 1.6) plot(1:3, c(coef(O2K.nls1)[2], O2K.jack1$estijack[2, 1], O2K.boot1$bootCI[2, 1]), ylim = c(O2K.boot1$bootCI[2, 2]*.96, O2K.boot1$bootCI[2, 3]*1.04), xlim = c(0,4), xaxt = "n", xlab = "", ylab = expression(paste(VO[2], "peak (L)")), main = expression(paste(VO[2], "peak")), pch = 19, cex.axis = 1.5, cex.lab = 1.5, cex.main = 1.5); segments(x0 = 1:3, y0 = c(binf[2], O2K.jack1$jackCI[2, 2], O2K.boot1$bootCI[2, 2]), x1 = 1:3, y1 = c(bsup[2], O2K.jack1$jackCI[2, 3], O2K.boot1$bootCI[2, 3]), lty = 1:3) axis(1, at = 1:3, labels = c("t-based", "Jackknife", "Bootstrap"), las = 3, cex.axis = 1.6) plot(1:3, c(coef(O2K.nls1)[3], O2K.jack1$estijack[3, 1], O2K.boot1$bootCI[3,1]), ylim = c(O2K.boot1$bootCI[3, 2]*.92, O2K.boot1$bootCI[3,3]*1.06), xlim = c(0,4), xaxt = "n", xlab = "", ylab = expression(paste(mu, " (L/min)")), main = expression(paste(mu, "")), pch = 19, cex.axis = 1.5, cex.lab = 1.5, cex.main = 1.5) segments(x0 = 1:3, y0 = c(binf[3], O2K.jack1$jackCI[3, 2], O2K.boot1$bootCI[3, 2]), x1 = 1:3, y1 = c(bsup[3], O2K.jack1$jackCI[3, 3], O2K.boot1$bootCI[3, 3]), lty = 1:3) axis(1, at = 1:3, labels = c("t-based", "Jackknife", "Bootstrap"), las = 3, cex.axis = 1.6)