propagate (version 1.3-0)

predictNLS: Confidence/prediction intervals for (weighted) nonlinear models based on uncertainty propagation

Description

A function for calculating confidence/prediction intervals of (weighted) nonlinear models for the supplied or new predictor values, 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 (correlated) errors.

Usage

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

Value

A list with the following items:

summary: a data frame with one row per predictor value, containing the prop and sim results of propagate: the mean and uncertainty estimates from first-/second-order Taylor expansion with the \(t_{n-p}\)-based interval (columns 1-6), and from Monte Carlo simulation with the probabilistically symmetric Monte Carlo coverage interval (columns 7-12).

prop: the complete output from propagate for each row of newdata.

Arguments

model

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

newdata

a data frame with new predictor values, with one column per predictor variable of model, named as in the model formula (in any order). See predict.nls and 'Examples'. If omitted, the model's predictor values are employed.

newerror

optional uncertainties of the predictor values: a data frame of standard uncertainties with the same column names and number of rows as newdata (independent predictors), or, for correlated predictors, a covariance matrix with row and column names as in newdata (used for all rows), or a list of such matrices, one per row of newdata. See 'Examples'.

interval

a character string: "confidence" for confidence intervals of the fitted curve, "prediction" for prediction intervals of a new observation, or "none". With "none", the results are those of "confidence".

alpha

a numeric value: \(1 - \alpha\) is the confidence level of the intervals (the coverage probability of the Monte Carlo coverage interval).

...

other parameters to be supplied to propagate.

Author

Andrej-Nikolai Spiess

Details

Taylor expansion. Calculation of the propagated variance \(\sigma_y^2 = \nabla \mathbf{\Sigma} \nabla^T\) is called the "Delta Method" and is widely applied in NLS fitting. It is based on a first-order Taylor expansion and thus assumes that \(f\) is linear in \(\beta\) around \(\hat{\beta}\). The second-order approach of propagate partially corrects for this restriction by using a second-order polynomial.

Confidence and prediction intervals are calculated as $$\hat{y} \pm t_{n-p}(1 - \alpha/2)\,\sigma_y \quad (1), \qquad \hat{y} \pm t_{n-p}(1 - \alpha/2)\sqrt{\sigma_y^2 + \textcolor{red}{\sigma_r^2}} \quad (2),$$ with \(n\) the number of observations, \(p\) the number of coefficients and the residual variance $$\textcolor{red}{\sigma_r^2} = \frac{\sum_{i=1}^n w_i (y_i - \hat{y}_i)^2}{n - p} \quad (3)$$ (\(w_i = 1\) for unweighted fits). The coefficients and the residual term share this variance estimate and are therefore one component with \(n - p\) degrees of freedom for the effective degrees of freedom (Willink 2007, Result 1; passed to propagate as ws.group). Without predictor errors, the coverage factor is thus exactly \(t_{n-p}(1 - \alpha/2)\), the classical delta-method interval; uncertainties of the predictors from newerror (Gaussian, \(\nu = \infty\)) increase the effective degrees of freedom according to the generalized Welch-Satterthwaite formula (see WelchSatter). A df.tot passed in ... replaces these effective degrees of freedom, i.e. it only sets the coverage factor of the Taylor expansion interval; it does not change the Monte Carlo simulation.

Augmented model for prediction intervals. The new observation is \(y = f(x, \beta) + \varepsilon\) (4), with a residual term \(\varepsilon\) of expectation 0 and variance \(\textcolor{red}{\sigma_r^2}\). It enters as an additional variable rv of the expression, and the covariance matrix \(\mathbf{\Sigma}\) of the \(p\) coefficients is augmented to a \((p+1) \times (p+1)\) matrix with \(\mathbf{\Sigma}_{p+1,p+1} = \textcolor{red}{\sigma_r^2}\). For two coefficients with variances \(\sigma_1^2, \sigma_2^2\) and covariance \(\sigma_{12}\): $$\left[\frac{\partial f}{\partial \beta_1}\; \frac{\partial f}{\partial \beta_2}\; \textcolor{red}{1}\right] \left[ \begin{array}{ccc} \sigma_1^2 & \sigma_{12} & 0 \\ \sigma_{12} & \sigma_2^2 & 0 \\ 0 & 0 & \textcolor{red}{\sigma_r^2} \end{array} \right] \left[ \begin{array}{c} \frac{\partial f}{\partial \beta_1} \\ \frac{\partial f}{\partial \beta_2} \\ \textcolor{red}{1} \end{array} \right] = \left(\frac{\partial f}{\partial \beta_1}\right)^2\sigma_1^2 + 2 \frac{\partial f}{\partial \beta_1} \frac{\partial f}{\partial \beta_2} \sigma_{12} + \left(\frac{\partial f}{\partial \beta_2}\right)^2\sigma_2^2 + \textcolor{red}{\sigma_r^2} \equiv \sigma_y^2 + \textcolor{red}{\sigma_r^2},$$ which then goes into (2). The prediction interval assumes that the new observation has unit weight.

Monte Carlo simulation. The coefficients \(\hat{\beta}\) and the residual term \(\varepsilon\) all depend on the same estimated residual variance \(\textcolor{red}{\sigma_r^2}\). Their joint distribution is therefore a multivariate \(t\)-distribution with \(\nu = n - p\) degrees of freedom and one common \(\chi^2_\nu\) variable, not a set of independent or copula-joined \(t\)-distributions: $$(\beta, \varepsilon) \sim t_{n-p}\!\left((\hat{\beta}, 0),\; \left[\begin{array}{cc} \mathbf{\Sigma} & 0 \\ 0 & \textcolor{red}{\sigma_r^2} \end{array}\right]\right).$$ This is the posterior distribution for a model that is linear in the coefficients, with Gaussian errors and a non-informative prior, and it has the same structure as the multivariate \(t\)-distribution of JCGM 102:2011, 5.3.2. predictNLS samples it with propagate(..., type = "stat", mc.margins = "mvt"), with \(\nu = n - p\) in the 3rd row of the input data. The uncertainties of the predictor values from newerror are Type B uncertainties independent of the fit (degrees of freedom Inf); the predictor values are sampled from a (multivariate) Gaussian distribution with the covariance matrix given in newerror, or, if a data frame of standard uncertainties is given, from independent Gaussian distributions.

For a model that is linear in the coefficients and without predictor errors, the Taylor expansion reproduces the classical intervals of predict.lm exactly, and the Monte Carlo simulation within the Monte Carlo error: \(\hat{y} \pm t_{n-p}(1 - \alpha/2)\,\sigma_y\) for confidence and \(\hat{y} \pm t_{n-p}(1 - \alpha/2)\sqrt{\sigma_y^2 + \sigma_r^2}\) for prediction intervals (see 'Examples'). For nonlinear models, the Monte Carlo intervals additionally account for the nonlinearity of \(f\).

Predictor errors (newerror). newerror can be given in three forms:
1) a data frame (or a matrix without matching row names) of standard uncertainties, with the predictor names as column names and one row per row of newdata: the predictor errors are independent.
2) a covariance matrix (variances on the diagonal, covariances off the diagonal, not standard uncertainties) whose row and column names are the predictor names, in any order: it is used for all rows of newdata.
3) a list of such covariance matrices, one per row of newdata.
In forms 2) and 3), the matrices must be symmetric. The predictor covariance matrix is added as a separate block to the covariance matrix of the coefficients (see mixCov), i.e. the predictor errors are uncorrelated with the coefficients.

References

Nonlinear Regression.
Seber GAF & Wild CJ.
John Wiley & Sons, New York (1989).

Nonlinear Regression Analysis and Its Applications.
Bates DM & Watts DG.
John Wiley & Sons, New York (1988).

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

Least-squares analysis of data with uncertainty in x and y: A Monte Carlo methods comparison.
Tellinghuisen J.
Chemometrics and Intelligent Laboratory Systems (2010), 103: 160-169.

A generalization of the Welch-Satterthwaite formula for use with correlated uncertainty components.
Willink R.
Metrologia (2007), 44: 340-349.

Evaluation of measurement data - Supplement 2 to the "Guide to the expression of uncertainty in measurement" - Extension to any number of output quantities.
JCGM 102:2011. Joint Committee for Guides in Metrology (2011).
tools:::Rd_expr_doi("10.59161/JCGM102-2011")

Examples

## In these examples, 'nsim' is 1E4 or 1E5 to save
## Rcmd check time (CRAN). It is advocated
## to use at least 'nsim = 1E6' though...

## 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), nsim = 10000)
PRED1 <- predict(fm3DNase1, newdata = data.frame(conc = 2))
PROP1$summary
PRED1
## => mean.1 equal to PRED1

## Using a single predictor value with error.
PROP2 <- predictNLS(fm3DNase1, newdata = data.frame(conc = 2), 
                    newerror = data.frame(conc = 0.5), nsim = 10000)
PROP2$summary

## Two predictors with correlated errors (r = 0.9)
set.seed(1); x <- runif(30, 1, 10); z <- runif(30, 0, 2)
y <- 5 * exp(-0.3 * x) + 2 * z^1.5 + rnorm(30, 0, 0.05)
fm2 <- nls(y ~ a * exp(-b * x) + g * z^d, start = list(a = 4, b = 0.2, g = 1, d = 1))
S <- matrix(c(0.2^2, 0.9 * 0.2 * 0.05, 0.9 * 0.2 * 0.05, 0.05^2), 2,
            dimnames = list(c("x", "z"), c("x", "z")))
## 'S' is a covariance matrix (variances, not standard uncertainties)
## whose dimnames are the predictor names.
predictNLS(fm2, newdata = data.frame(x = 5, z = 1), newerror = S, nsim = 1E4)$summary
## A list of covariance matrices, one per row of 'newdata'.
predictNLS(fm2, newdata = data.frame(x = c(5, 6), z = c(1, 1)),
           newerror = list(S, 2 * S), nsim = 1E4)$summary

# \donttest{
## Using a sequence of predictor values without error.
CONC <- seq(1, 12, by = 1)
PROP3 <- predictNLS(fm3DNase1, newdata = data.frame(conc = CONC), nsim = 10000)
PRED3 <- predict(fm3DNase1, newdata = data.frame(conc = CONC))
PROP3$summary
PRED3
## => mean.1 equal to PRED3

## Plot mean 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, PROP3$summary[, 1], col = 2, pch = 16)
lines(CONC, PROP3$summary[, 5], col = 2)
lines(CONC, PROP3$summary[, 6], col = 2)
lines(CONC, PROP3$summary[, 11], col = 4)
lines(CONC, PROP3$summary[, 12], col = 4)
# } 

## Gauge: for a model that is linear in the coefficients,
## the Taylor and Monte Carlo intervals equal those of predict.lm.
set.seed(70)
xx <- 1:15
yy <- 2 + 0.5 * xx + rnorm(15, 0, 0.4)
fitN <- nls(yy ~ a + b * xx, start = list(a = 1, b = 1))
fitL <- lm(yy ~ xx)
PN <- predictNLS(fitN, newdata = data.frame(xx = 7.3), interval = "prediction", nsim = 1E5)
PN$summary[, c(5, 6, 11, 12)]
predict(fitL, newdata = data.frame(xx = 7.3), interval = "prediction")
## => Taylor interval (columns 5, 6) equal to lwr/upr (4.806, 6.373),
##    Monte Carlo interval (columns 11, 12) within the Monte Carlo error