propagate (version 1.0-4)

predictNLS: Confidence intervals for nonlinear models based on uncertainty propagation

Description

A function for calculating confidence intervals for the fitted values of nonlinear models by using first-/second-order Taylor expansion and Monte Carlo simulation. This approach can be used to construct more realistic error estimates and confidence/prediction intervals for nonlinear models than what is possible with only a simple linearization (first-order Taylor expansion) approach. Another application is when there is an "error in x" setup with uncertainties in the predictor variable (See 'Examples'). This function will also work in the presence of multiple predictors with/without errors.

Usage

predictNLS(model, newdata, interval = c("confidence", "prediction", "none"),
           alpha = 0.05, ...)

Arguments

model

a model obtained from nls or nlsLM (package 'minpack.lm').

newdata

A named list or data frame in which to look for variables with which to predict. The first \(n\) columns must contain the predictor values, the following \(n\) columns can contain errors. See predict.nls and 'Examples'.

interval

A character string indicating if confidence/prediction intervals on the mean responses are to be calculated or not.

alpha

the 1 - confidence level.

...

other parameters to be supplied to propagate.

Value

A list with the following items: summary: The mean/error estimates obtained from first-/second-order Taylor expansion and Monte Carlo simulation, together with calculated confidence/prediction intervals based on asymptotic normality. prop: the complete output from propagate for each value in newdata.

Details

Calculation of the propagated uncertainty \(\sigma_y^2\) using \(\nabla_x C_x \nabla_x^T\) is called the "Delta Method" and is widely applied in NLS fitting. However, this method is based on first-order Taylor expansion and thus assummes linearity around \(f(x)\). The second-order approach as implemented in the propagate function can partially correct for this restriction by using a second-order polynomial around \(f(x)\). Confidence/prediction intervals are calculated in a usual way using \(t(1 - \frac{\alpha}{2}, \nu) \cdot \sigma_{prop}\) or \(t(1 - \frac{\alpha}{2}, \nu) \cdot \sqrt{\sigma_{prop}^2 + rv}\), respectively, where \(rv\) = the residual variance \(\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{n - \nu}\). If errors are supplied to the predictor values in newdata, for \(n\) predictor values they have to be in the next \(n\) columns, i.e. predictors in [, 1:3], and errors in [, 4:6].

References

Nonlinear Regression. Seber GAF & Wild CJ. John Wiley & Sons; 1ed, 2003.

Nonlinear Regression Analysis and its Applications. Bates DM & Watts DG. Wiley-Interscience; 1ed, 2007.

Statistical Error Propagation. Tellinghuisen J. J. Phys. Chem. A (2001), 47: 3917-3921.

Least-squares analysis of data with uncertainty in x and y: A Monte Carlo methods comparison. Tellinghuisen J. Chemometr Intell Lab (2010), 47: 160-169.

From the author's blog: http://rmazing.wordpress.com/2013/08/14/predictnls-part-1-monte-carlo-simulation-confidence-intervals-for-nls-models/ http://rmazing.wordpress.com/2013/08/26/predictnls-part-2-taylor-approximation-confidence-intervals-for-nls-models/

Examples

# NOT RUN {
## Example from ?nls.
DNase1 <- subset(DNase, Run == 1)
fm3DNase1 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),
                 data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))

## Using a single predictor value without error.
PROP1 <- predictNLS(fm3DNase1, newdata = data.frame(conc = 2))
PRED1 <- predict(fm3DNase1, newdata = data.frame(conc = 2))
PROP1$summary
PRED1
## => Prop.Mean.1 equal to PRED1

# }
# NOT RUN {
## Using a sequence of predictor values without error.
CONC <- seq(1, 12, by = 1)
PROP2 <- predictNLS(fm3DNase1, newdata = data.frame(conc = CONC))
PRED2 <- predict(fm3DNase1, newdata = data.frame(conc = CONC))
PROP2$summary
PRED2
## => Prop.Mean.1 equal to PRED2

## Using a sequence of predictor values with error.
DAT <- data.frame(conc = CONC, error = rnorm(12, 0, 0.1))
PROP3 <- predictNLS(fm3DNase1, newdata = DAT)
PRED3 <- predict(fm3DNase1, newdata = DAT)
PROP3$summary
PRED3
## => Prop.Mean.1 equal to PRED3

## Plot predicted and confidence values from 
## first-/second-order Taylor expansion 
## and Monte Carlo simulation.
plot(DNase1$conc, DNase1$density)
lines(DNase1$conc, fitted(fm3DNase1), lwd = 2, col = 1)
points(CONC, PROP2$summary[, 1], col = 2, pch = 16)
lines(CONC, PROP2$summary[, 5], col = 2)
lines(CONC, PROP2$summary[, 6], col = 2)
# }
# NOT RUN {
## Using multiple predictor values
## 1: Setup of response values
## with gaussian error of 10%.
x <- seq(1, 10, by = 0.01)
y <- seq(10, 1, by = -0.01)
a <- 2
b <- 5
c <- 10
z <- a * exp(b * x)^sin(y/c)
z <- z + sapply(z, function(x) rnorm(1, x, 0.10 * x))
## 2: Fit 'nls' model.
MOD <- nls(z ~ a * exp(b * x)^sin(y/c), 
             start = list(a = 2, b = 5, c = 10))
## 3: newdata without errors and prediction.
DAT1 <- data.frame(x = 4, y = 3)
PROP4 <- predictNLS(MOD, newdata = DAT1)
PROP4$summary
## 4: newdata with errors and prediction.
DAT2 <- data.frame(x = 4, y = 3, error.x = 0.2, error.y = 0.1)
PROP5 <- predictNLS(MOD, newdata = DAT2)
PROP5$summary
# }