nlpredci

Nonlinear regression prediction confidence intervals

Syntax

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Covar',CovB)

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Covar',CovB,Name,Value)

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Jacobian',J)

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Jacobian',J,Name,Value)

Description

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Covar',CovB) returns predictions, Ypred, and 95% confidence interval half-widths, delta, for the nonlinear regression model modelfun at input values X. Before calling nlpredci, use nlinfit to fit modelfun and get the estimated coefficients, beta, residuals, R, and variance-covariance matrix, CovB.

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Covar',CovB,Name,Value) uses additional options specified by one or more name-value pair arguments.

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Jacobian',J) returns predictions, Ypred, and 95% confidence interval half-widths, delta, for the nonlinear regression model modelfun at input values X. Before calling nlpredci, use nlinfit to fit modelfun and get the estimated coefficients, beta, residuals, R, and Jacobian, J.

If you use a robust option with nlinfit, then you should use the Covar syntax rather than the Jacobian syntax. The variance-covariance matrix, CovB, is required to properly take the robust fitting into account.

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Jacobian',J,Name,Value) uses additional options specified by one or more name-value pair arguments.

Examples

Confidence Interval for Nonlinear Regression Curve

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Fit the Hougen-Watson model to the rate data using the initial values in beta0.

[beta,R,J] = nlinfit(X,y,@hougen,beta0);

Obtain the predicted response and 95% confidence interval half-width for the value of the curve at average reactant levels.

[ypred,delta] = nlpredci(@hougen,mean(X),beta,R,'Jacobian',J)

ypred = 
5.4622

delta = 
0.1921

Compute the 95% confidence interval for the value of the curve.

[ypred-delta,ypred+delta]

ans = 1×2

    5.2702    5.6543

Prediction Interval for New Observation

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Fit the Hougen-Watson model to the rate data using the initial values in beta0.

[beta,R,J] = nlinfit(X,y,@hougen,beta0);

Obtain the predicted response and 95% prediction interval half-width for a new observation with reactant levels [100,100,100].

[ypred,delta] = nlpredci(@hougen,[100,100,100],beta,R,'Jacobian',J,...
                         'PredOpt','observation')

ypred = 
1.8346

delta = 
0.5101

Compute the 95% prediction interval for the new observation.

[ypred-delta,ypred+delta]

ans = 1×2

    1.3245    2.3447

Simultaneous Confidence Intervals for Robust Fit Curve

Generate sample data from the nonlinear regression model $$y=b_1+b_2\cdot exp\{b_{3}x\}+ϵ$$, where $$b_1$$, $$b_2$$, and $$b_3$$ are coefficients, and the error term is normally distributed with mean 0 and standard deviation 0.5.

modelfun = @(b,x)(b(1)+b(2)*exp(-b(3)*x));

rng('default') % for reproducibility
b = [1;3;2];
x = exprnd(2,100,1);
y = modelfun(b,x) + normrnd(0,0.5,100,1);

Fit the nonlinear model using robust fitting options.

opts = statset('nlinfit');
opts.RobustWgtFun = 'bisquare';
beta0 = [2;2;2];
[beta,R,J,CovB,MSE] = nlinfit(x,y,modelfun,beta0,opts);

Plot the fitted regression model and simultaneous 95% confidence bounds.

xrange = min(x):.01:max(x);
[ypred,delta] = nlpredci(modelfun,xrange,beta,R,'Covar',CovB,...
                         'MSE',MSE,'SimOpt','on');
lower = ypred - delta;
upper = ypred + delta;

figure()
plot(x,y,'ko') % observed data
hold on
plot(xrange,ypred,'k','LineWidth',2)
plot(xrange,[lower;upper],'r--','LineWidth',1.5)

Confidence Interval Using Observation Weights

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Specify a function handle for observation weights, then fit the Hougen-Watson model to the rate data using the specified observation weights function.

a = 1; b = 1;
weights = @(yhat) 1./((a + b*abs(yhat)).^2);
[beta,R,J,CovB] = nlinfit(X,y,@hougen,beta0,'Weights',weights);

Compute the 95% prediction interval for a new observation with reactant levels [100,100,100] using the observation weight function.

[ypred,delta] = nlpredci(@hougen,[100,100,100],beta,R,'Jacobian',J,...
                         'PredOpt','observation','Weights',weights);
[ypred-delta,ypred+delta]

ans = 1×2

    1.5264    2.1033

Confidence Interval Using Nonconstant Error Model

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Fit the Hougen-Watson model to the rate data using the combined error variance model.

[beta,R,J,CovB,MSE,S] = nlinfit(X,y,@hougen,beta0,'ErrorModel','combined');

Compute the 95% prediction interval for a new observation with reactant levels [100,100,100] using the fitted error variance model.

[ypred,delta] = nlpredci(@hougen,[100,100,100],beta,R,'Jacobian',J,...
                         'PredOpt','observation','ErrorModelInfo',S);
[ypred-delta,ypred+delta]

ans = 1×2

    1.3245    2.3447

Input Arguments

modelfun — Nonlinear regression model function

function handle

Nonlinear regression model function, specified as a function handle. modelfun must accept two input arguments, a coefficient vector and an array X—in that order—and return a vector of fitted response values.

For example, to specify the hougen nonlinear regression function, use the function handle @hougen.

Data Types: function_handle

X — Input values for predictions

matrix

Input values for predictions, specified as a matrix. nlpredci makes a prediction for the covariates in each row of X. There should be a column in X for each coefficient in the model.

Data Types: single | double

beta — Estimated regression coefficients

vector returned by nlinfit

Estimated regression coefficients, specified as the vector of fitted coefficients returned by a previous call to nlinfit.

Data Types: single | double

R — Residuals

vector returned by nlinfit

Residuals for the fitted modelfun, specified as the vector of residuals returned by a previous call to nlinfit.

CovB — Estimated variance-covariance matrix

matrix returned by nlinfit

Estimated variance-covariance matrix for the fitted coefficients, beta, specified as the variance-covariance matrix returned by a previous call to nlinfit.

J — Estimated Jacobian

matrix returned by nlinfit

Estimated Jacobian of the nonlinear regression model, modelfun, specified as the Jacobian matrix returned by a previous call to nlinfit.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'Alpha',0.1,'PredOpt','observation' specifies 90% prediction intervals for new observations.

Alpha — Significance level

0.05 (default) | scalar value in the range (0,1)

Significance level for the confidence interval, specified as the comma-separated pair consisting of 'Alpha' and a scalar value in the range (0,1). If Alpha has value α, then nlpredci returns intervals with 100×(1–α)% confidence level.

The default confidence level is 95% (α = 0.05).

Example: 'Alpha',0.1

Data Types: single | double

ErrorModelInfo — Information about error model fit

structure returned by nlinfit

Information about the error model fit, specified as the comma-separated pair consisting of 'ErrorModelInfo' and a structure returned by a previous call to nlinfit.

ErrorModelInfo only has an effect on the returned prediction interval when PredOpt has the value 'observation'. If you do not use ErrorModelInfo, then nlpredci assumes the error variance model is 'constant'.

The error model structure returned by nlinfit has the following fields:

ErrorModel	Chosen error model
ErrorParameters	Estimated error parameters
ErrorVariance	Function handle that accepts an N-by-p matrix, X, and returns an N-by-1 vector of error variances using the estimated error model
MSE	Mean squared error
ScheffeSimPred	Scheffé parameter for simultaneous prediction intervals when using the estimated error model
WeightFunction	Logical with value true if you used a custom weight function previously in nlinfit
FixedWeights	Logical with value true if you used fixed weights previously in nlinfit
RobustWeightFunction	Logical with value true if you used robust fitting previously in nlinfit

MSE — Mean squared error

MSE returned by nlinfit

Mean squared error (MSE) for the fitted nonlinear regression model, specified as the comma-separated pair consisting of 'MSE' and the MSE value returned by a previous call to nlinfit.

If you use a robust option with nlinfit, then you must specify the MSE when predicting new observations to properly take the robust fitting into account. If you do not specify the MSE, then nlpredci computes the MSE from the residuals, R, and does not take the robust fitting into account.

For example, if mse is the MSE value returned by nlinfit, then you can specify 'MSE',mse.

Data Types: single | double

PredOpt — Prediction interval to compute

'curve' (default) | 'observation'

Prediction interval to compute, specified as the comma-separated pair consisting of 'PredOpt' and either 'curve' or 'observation'.

• If you specify the value 'curve', then nlpredci returns confidence intervals for the estimated curve (function value) at the observations X.

• If you specify the value 'observation', then nlpredci returns prediction intervals for new observations at X.

If you specify 'observation' after using a robust option with nlinfit, then you must also specify a value for MSE to provide the robust estimate of the mean squared error.

Example: 'PredOpt','observation'

SimOpt — Indicator for specifying simultaneous bounds

'off' (default) | 'on'

Indicator for specifying simultaneous bounds, specified as the comma-separated pair consisting of 'SimOpt' and either 'off' or 'on'. Use the value 'off' to compute nonsimultaneous bounds, and 'on' for simultaneous bounds.

Weights — Observation weights

vector | function handle

Observation weights, specified as the comma-separated pair consisting of 'Weights' and a vector of positive scalar values or a function handle. The default is no weights.

• If you specify a vector of weights, then it must have the same number of elements as the number of observations (rows) in X.

• If you specify a function handle for the weights, then it must accept a vector of predicted response values as input, and return a vector of real positive weights as output.

Given weights, W, nlpredci estimates the error variance at observation i by mse*(1/W(i)), where mse is the mean squared error value specified using MSE.

Example: 'Weights',@WFun

Data Types: double | single | function_handle

Output Arguments

Ypred — Predicted responses

vector

Predicted responses, returned as a vector with the same number of rows as X.

delta — Confidence interval half-widths

vector

Confidence interval half-widths, returned as a vector with the same number of rows as X. By default, delta contains the half-widths for nonsimultaneous 95% confidence intervals for modelfun at the observations in X. You can compute the lower and upper bounds of the confidence intervals as Ypred-delta and Ypred+delta, respectively.

If 'PredOpt' has value 'observation', then delta contains the half-widths for prediction intervals of new observations at the values in X.

More About

Confidence Intervals for Estimable Predictions

When the estimated model Jacobian is not of full rank, then it might not be possible to construct sensible confidence intervals at all prediction points. In this case, nlpredci still tries to construct confidence intervals for any estimable prediction points.

For example, suppose you fit the linear function $$f(x_i,\beta )={\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+{\beta }_3{x}_{i3}$$ at the points in the design matrix

$$X=\left(\begin{array}{ccc}1& 1& 0\\ 1& 1& 0\\ 1& 1& 0\\ 1& 0& 1\\ 1& 0& 1\\ 1& 0& 1\end{array}\right).$$

The estimated Jacobian at the values in $$X$$ is the design matrix itself, $$J=X.$$ Thus, the Jacobian is not of full rank:

rng('default') % For reproducibility
y = randn(6,1);

linfun = @(b,x) x*b;
beta0 = [1;1;1];
X = [repmat([1 1 0],3,1); repmat([1 0 1],3,1)];   

[beta,R,J]  = nlinfit(X,y,linfun,beta0);

Warning: The Jacobian at the solution is ill-conditioned, and
some model parameters may not be estimated well (they are not
identifiable).  Use caution in making predictions. 
> In nlinfit at 283 

In this example, nlpredci can only compute prediction intervals at points that satisfy the linear relationship

$$x_{i1}=x_{i2}+x_{i3}.$$

If you try to compute confidence intervals for predictions at nonidentifiable points, nlpredci returns NaN for the corresponding interval half-widths:

xpred = [1 1 1;0 1 -1;2 1 1];
[ypred,delta] = nlpredci(linfun,xpred,beta,R,'Jacobian',J)

ypred =

   -0.0035
    0.0798
   -0.0047


delta =

       NaN
    3.8102
    3.8102

Here, the first element of delta is NaN because the first row in xpred does not satisfy the required linear dependence, and is therefore not an estimable contrast.

Tips

• To compute confidence intervals for complex parameters or data, you need to split the problem into its real and imaginary parts. When calling nlinfit:

  1. Define your parameter vector beta as the concatenation of the real and imaginary parts of the original parameter vector.

  2. Concatenate the real and imaginary parts of the response vector Y as a single vector.

  3. Modify your model function modelfun to accept X and the purely real parameter vector, and return a concatenation of the real and imaginary parts of the fitted values.

  With the problem formulated this way, nlinfit computes real estimates, and confidence intervals are feasible.

Algorithms

• nlpredci treats NaN values in the residuals, R, or the Jacobian, J, as missing values, and ignores the corresponding observations.

• If the Jacobian, J, does not have full column rank, then some of the model parameters might be nonidentifiable. In this case, nlpredci tries to construct confidence intervals for estimable predictions, and returns NaN for those that are not.

References

[1] Lane, T. P. and W. H. DuMouchel. “Simultaneous Confidence Intervals in Multiple Regression.” The American Statistician. Vol. 48, No. 4, 1994, pp. 315–321.

[2] Seber, G. A. F., and C. J. Wild. Nonlinear Regression. Hoboken, NJ: Wiley-Interscience, 2003.

Version History

Introduced before R2006a