nlpredci
Nonlinear regression prediction confidence intervals
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,'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.
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=b1+b2⋅exp{b3x}+ϵ, where b1, b2, and b3 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 Pair Arguments
Specify optional
comma-separated pairs of Name,Value arguments. Name is
the argument name and Value is the corresponding value.
Name must appear inside quotes. You can specify several name and value
pair arguments in any order as
Name1,Value1,...,NameN,ValueN.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:
ErrorModelChosen error model
ErrorParametersEstimated error parameters
ErrorVarianceFunction handle that accepts an N-by-p matrix, X,
and returns an N-by-1 vector of error variances
using the estimated error model
MSEMean squared error
ScheffeSimPredScheffé parameter for simultaneous prediction intervals
when using the estimated error model
WeightFunctionLogical with value true if you used a custom
weight function previously in nlinfit
FixedWeightsLogical with value true if you used fixed
weights previously in nlinfit
RobustWeightFunctionLogical 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(xi,β)=β1xi1+β2xi2+β3xi3 at the points in the design
matrix
X=(110110110101101101).
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
xi1=xi2+xi3.
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.8102Here, 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:
Define your parameter vector beta as
the concatenation of the real and imaginary parts of the original
parameter vector.
Concatenate the real and imaginary parts of the response
vector Y as a single vector.
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.
Introduced before R2006a