nlreg: Fit a Nonlinear Heteroscedastic Model via Maximum Likelihood
In nlreg: Higher Order Inference for Nonlinear Heteroscedastic Models
Description
Usage
Arguments
Details
Value
Side Effects
Note
References
See Also
Examples
Description
Returns an object of class nlreg which represents a nonlinear
heteroscedastic model fit of the data obtained by maximizing the
corresponding likelihood function.
Usage
1
2
3
4
5
Arguments
formula
a formula expression as for other nonlinear regression models, of
the form response ~ f(predictor) where f is a
nonlinear function of the predictor involving a number of
regression coefficients. Only one predictor variable can be
included in the model formula. Missing values are not allowed.
weights
a formula expression of the form ~ V(predictor) where
V is a nonlinear variance function involving the predictor
or some transformation of it, variance parameters and/or
regression coefficients. The error variance nlreg works
with is
Var(error) = s^2 V(predictor)
where the constant term \code{s^2} is included by
default and must
not be specified in the weights argument. The nlreg
routine treats it on the logarithmic scale and assigns to it the
parameter name logs. By default, the error variance is
assumed to be constant.
data
an optional data frame in which to interpret the variables
occurring in the model formula. Missing values are not allowed.
start
a numerical vector containing the starting values that initialize
the iterative estimating procedure. Each component of the vector
must be named and represents one of the parameters included in
the mean and, if defined, variance function. Starting values have
to be supplied for every model parameter, except for the constant
term in the variance function which is included by default in the
model. See the weights argument above.
offset
a numerical vector with a single named element. The name
indicates the parameter of interest which will be fixed to the
value specified. logs is used to identify the constant
term \code{s^2} which is included by default in the
variance function.
subset
expression saying which subset of the rows of the data should be
used in the fit. This can be a logical vector or a numeric vector
indicating which observation numbers are to be included. All
observations are included by default.
control
a list of iteration and algorithmic constants. See the
Details section below for their definition.
trace
logical flag. If TRUE, details of the iterations are
printed.
hoa
logical flag. If TRUE, the first and second derivatives of
the mean and, if defined, variance functions are stored in the
fitted model object. The default is FALSE.
Details
A nonlinear heteroscedastic model representing the relationship
between two scalar quantities is fitted. The response is specified
on the left-hand side of the formula argument. The predictor
appears in the right-hand side of the formula and, if
specified, weights arguments. Only one predictor variable
can be included. Missing values in the data are not allowed.
The fitting criterion is maximum likelihood. The core algorithm
implemented in nlreg alternates minimization of minus the
log likelihood with respect to the regression coefficients and the
variance parameters. The quasi-Newton optimizer
optim is used in both steps. The constant
term \code{s^2} in
Var(error) = s^2 V(predictor) is included by default. In
order to work with a real value,
\code{s^2} is estimated on the logarithmic scale,
that is, the model is reparametrized into
\code{log(s^2)} which gives rise to the
parameter name logs. If the errors are homoscedastic, the
second step is omitted and the algorithm switches automatically to
nls. If the weights argument is
omitted, homoscedasticity of the errors is assumed.
Starting values for all parameters have to be passed through the
start argument except for \code{s^2} for which
the maximum likelihood estimate is available in closed form.
Starting values should be chosen carefully in order to avoid
convergence to a local maximum.
The algorithm iterates until convergence or until the maximum number
of iterations defined by maxit is reached. The stopping rule
considers the relative difference between successive estimates and
the relative increment of the log likelihood. These are governed by
the parameters x.tol and rel.tol/step.min,
respectively.
If convergence has been reached, the results are saved in an object
of class nlreg. The output can be examined by
print and summary.
Components can be extracted using coef,
param, fitted
and residuals. The model fit can be updated using
update. Profile plots and profile pair
sketches are provided by profile, and
contour. Diagnostic plots are obtained from
nlreg.diag.plots or simply
plot.
The details are given in Brazzale (2000, Section 6.3.1).
Value
An object of class nlreg is returned which inherits from
nls. See nlreg.object for
details.
Side Effects
If trace = TRUE, the iteration number and the corresponding
log likelihood are printed.
Note
The arguments hoa and control play an important role.
The first forces the algorithm to save the derivatives of the mean
and variance functions in the fitted model object. This is
imperative if one wants to save execution time, especially for
complex models. Fine-tuning of the control argument which
controls the convergence criteria helps surrounding a well-known
problem in nonlinear regression, that is, convergence failure in
cases where the likelihood is very flat.
If the errors are homoscedastic, the nlreg routine switches
automatically to nls which, although rarely,
dumps because of convergence problems. To avoid this, either
reparametrize the model (see Bates and Watts, 1988) or
choose starting values that are more realistic. This advice also
holds in case of convergence problems for models with non constant
variance function. Use the
trace = TRUE argument to gain insight into what goes on at
the different iteration steps.
The weights argument has a different meaning than in other
model fitting routines such as lm and
glm. It defines the variance function of the
nonlinear model and not a vector of observation weights that are
multiplied into the squared residuals.
References
Bates, D. M. and Watts, D. G. (1988) Nonlinear Regression
Analysis and Its Applications. New York: Wiley.
Brazzale, A. R. (2000) Practical Small-Sample Parametric
Inference. Ph.D. Thesis N. 2230, Department of Mathematics, Swiss
Federal Institute of Technology Lausanne.
Seber, G. A. F. and Wild, C. J. (1989) Nonlinear Regression.
New York: Wiley.
See Also
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
library(boot)
data(calcium)
##
## Homoscedastic model fit
calcium.nl <- nlreg( cal ~ b0*(1-exp(-b1*time)), start = c(b0 = 4, b1 = 0.1),
data = calcium )
##
## Heteroscedastic model fit
calcium.nl <- nlreg( cal ~ b0*(1-exp(-b1*time)), weights = ~ ( 1+time^g )^2,
start = c(b0 = 4, b1 = 0.1, g = 1), data = calcium,
hoa = TRUE)
## or
calcium.nl <- update(calcium.nl, weights = ~ (1+time^g)^2,
start = c(b0 = 4, b1 = 0.1, g = 1), hoa = TRUE )
##
##
## Power-of-X (POX) variance function
##
data(metsulfuron)
metsulfuron.nl <-
nlreg( log(area) ~ log( b1+(b2-b1) / (1+(dose/b4)^b3) ),
weights = ~ ( 1+dose^exp(g) )^2, data = metsulfuron,
start = c(b1 = 138, b2 = 2470, b3 = 2, b4 = 0.07, g = log(0.3)),
hoa = TRUE )
##
##
## Power-of-mean (POM) variance function
##
data(ria)
ria.nl <- nlreg( count ~ b1+(b2-b1) / (1+(conc/b4)^b3),
weights = ~ ( b1+(b2-b1) / (1+(conc/b4)^b3) )^g, data = ria,
start = c(b1 = 1.6, b2 = 20, b3 = 2, b4 = 300, g = 2),
hoa = TRUE, trace = TRUE )
##
##
## Error-in-variables (EIV) variance function
##
data(chlorsulfuron)
options( object.size = 10000000 )
chlorsulfuron.nl <-
nlreg( log(area) ~ log( b1+(b2-b1) / (1+(dose/b4)^b3) ),
weights = ~ ( 1+k*dose^g*(b2-b1)^2/(1+(dose/b4)^b3)^4*b3^2*dose^(2*b3-2)/
b4^(2*b3)/(b1+(b2-b1)/(1+(dose/b4)^b3))^2 ),
start = c(b1 = 2.2, b2 = 1700, b3 = 2.8, b4 = 0.28, g = 2.7, k = 1),
data = chlorsulfuron, hoa = TRUE, trace = TRUE,
control = list(x.tol = 10^-12, rel.tol = 10^-12, step.min = 10^-12) )
Example output
Loading required package: statmod
Loading required package: survival
Package "nlreg" 1.2-2.2 (2019-01-30)
Copyright (C) 2000-2019 R. Bellio & A. R. Brazzale
This is free software, and you are welcome to redistribute
it and/or modify it under the terms of the GNU General
Public License published by the Free Software Foundation.
Package "nlreg" comes with ABSOLUTELY NO WARRANTY.
type `help(package="nlreg")' for summary information
Attaching package: ‘boot’
The following object is masked from ‘package:survival’:
aml
differentiating mean and variance function -- may take a while
differentiating mean and variance function -- may take a while
differentiating mean and variance function -- may take a while
iteration 1 : log likelihood = -3.264223
iteration 2 : log likelihood = -3.263662
iteration 3 : log likelihood = -3.263662
differentiating mean and variance function -- may take a while
iteration 1 : log likelihood = -35.12366
iteration 2 : log likelihood = -35.10621
iteration 3 : log likelihood = -35.10528
iteration 4 : log likelihood = -35.10515
iteration 5 : log likelihood = -35.10513
iteration 6 : log likelihood = -35.10513
iteration 7 : log likelihood = -35.10513
iteration 8 : log likelihood = -35.10513
iteration 9 : log likelihood = -35.10513
iteration 10 : log likelihood = -35.10513
differentiating mean and variance function -- may take a while
nlreg documentation built on May 1, 2019, 9:21 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Side Effects Note References See Also Examples
Description
Returns an object of class nlreg which represents a nonlinear
heteroscedastic model fit of the data obtained by maximizing the
corresponding likelihood function.
Usage
1 2 3 4 5 |
Arguments
formula |
a formula expression as for other nonlinear regression models, of
the form |
weights |
a formula expression of the form
where the constant term \code{s^2} is included by
default and must
not be specified in the |
data |
an optional data frame in which to interpret the variables occurring in the model formula. Missing values are not allowed. |
start |
a numerical vector containing the starting values that initialize
the iterative estimating procedure. Each component of the vector
must be named and represents one of the parameters included in
the mean and, if defined, variance function. Starting values have
to be supplied for every model parameter, except for the constant
term in the variance function which is included by default in the
model. See the |
offset |
a numerical vector with a single named element. The name
indicates the parameter of interest which will be fixed to the
value specified. |
subset |
expression saying which subset of the rows of the data should be used in the fit. This can be a logical vector or a numeric vector indicating which observation numbers are to be included. All observations are included by default. |
control |
a list of iteration and algorithmic constants. See the Details section below for their definition. |
trace |
logical flag. If |
hoa |
logical flag. If |
Details
A nonlinear heteroscedastic model representing the relationship
between two scalar quantities is fitted. The response is specified
on the left-hand side of the formula argument. The predictor
appears in the right-hand side of the formula and, if
specified, weights arguments. Only one predictor variable
can be included. Missing values in the data are not allowed.
The fitting criterion is maximum likelihood. The core algorithm
implemented in nlreg alternates minimization of minus the
log likelihood with respect to the regression coefficients and the
variance parameters. The quasi-Newton optimizer
optim is used in both steps. The constant
term \code{s^2} in
Var(error) = s^2 V(predictor) is included by default. In
order to work with a real value,
\code{s^2} is estimated on the logarithmic scale,
that is, the model is reparametrized into
\code{log(s^2)} which gives rise to the
parameter name logs. If the errors are homoscedastic, the
second step is omitted and the algorithm switches automatically to
nls. If the weights argument is
omitted, homoscedasticity of the errors is assumed.
Starting values for all parameters have to be passed through the
start argument except for \code{s^2} for which
the maximum likelihood estimate is available in closed form.
Starting values should be chosen carefully in order to avoid
convergence to a local maximum.
The algorithm iterates until convergence or until the maximum number
of iterations defined by maxit is reached. The stopping rule
considers the relative difference between successive estimates and
the relative increment of the log likelihood. These are governed by
the parameters x.tol and rel.tol/step.min,
respectively.
If convergence has been reached, the results are saved in an object
of class nlreg. The output can be examined by
print and summary.
Components can be extracted using coef,
param, fitted
and residuals. The model fit can be updated using
update. Profile plots and profile pair
sketches are provided by profile, and
contour. Diagnostic plots are obtained from
nlreg.diag.plots or simply
plot.
The details are given in Brazzale (2000, Section 6.3.1).
Value
An object of class nlreg is returned which inherits from
nls. See nlreg.object for
details.
Side Effects
If trace = TRUE, the iteration number and the corresponding
log likelihood are printed.
Note
The arguments hoa and control play an important role.
The first forces the algorithm to save the derivatives of the mean
and variance functions in the fitted model object. This is
imperative if one wants to save execution time, especially for
complex models. Fine-tuning of the control argument which
controls the convergence criteria helps surrounding a well-known
problem in nonlinear regression, that is, convergence failure in
cases where the likelihood is very flat.
If the errors are homoscedastic, the nlreg routine switches
automatically to nls which, although rarely,
dumps because of convergence problems. To avoid this, either
reparametrize the model (see Bates and Watts, 1988) or
choose starting values that are more realistic. This advice also
holds in case of convergence problems for models with non constant
variance function. Use the
trace = TRUE argument to gain insight into what goes on at
the different iteration steps.
The weights argument has a different meaning than in other
model fitting routines such as lm and
glm. It defines the variance function of the
nonlinear model and not a vector of observation weights that are
multiplied into the squared residuals.
References
Bates, D. M. and Watts, D. G. (1988) Nonlinear Regression Analysis and Its Applications. New York: Wiley.
Brazzale, A. R. (2000) Practical Small-Sample Parametric Inference. Ph.D. Thesis N. 2230, Department of Mathematics, Swiss Federal Institute of Technology Lausanne.
Seber, G. A. F. and Wild, C. J. (1989) Nonlinear Regression. New York: Wiley.
See Also
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | library(boot)
data(calcium)
##
## Homoscedastic model fit
calcium.nl <- nlreg( cal ~ b0*(1-exp(-b1*time)), start = c(b0 = 4, b1 = 0.1),
data = calcium )
##
## Heteroscedastic model fit
calcium.nl <- nlreg( cal ~ b0*(1-exp(-b1*time)), weights = ~ ( 1+time^g )^2,
start = c(b0 = 4, b1 = 0.1, g = 1), data = calcium,
hoa = TRUE)
## or
calcium.nl <- update(calcium.nl, weights = ~ (1+time^g)^2,
start = c(b0 = 4, b1 = 0.1, g = 1), hoa = TRUE )
##
##
## Power-of-X (POX) variance function
##
data(metsulfuron)
metsulfuron.nl <-
nlreg( log(area) ~ log( b1+(b2-b1) / (1+(dose/b4)^b3) ),
weights = ~ ( 1+dose^exp(g) )^2, data = metsulfuron,
start = c(b1 = 138, b2 = 2470, b3 = 2, b4 = 0.07, g = log(0.3)),
hoa = TRUE )
##
##
## Power-of-mean (POM) variance function
##
data(ria)
ria.nl <- nlreg( count ~ b1+(b2-b1) / (1+(conc/b4)^b3),
weights = ~ ( b1+(b2-b1) / (1+(conc/b4)^b3) )^g, data = ria,
start = c(b1 = 1.6, b2 = 20, b3 = 2, b4 = 300, g = 2),
hoa = TRUE, trace = TRUE )
##
##
## Error-in-variables (EIV) variance function
##
data(chlorsulfuron)
options( object.size = 10000000 )
chlorsulfuron.nl <-
nlreg( log(area) ~ log( b1+(b2-b1) / (1+(dose/b4)^b3) ),
weights = ~ ( 1+k*dose^g*(b2-b1)^2/(1+(dose/b4)^b3)^4*b3^2*dose^(2*b3-2)/
b4^(2*b3)/(b1+(b2-b1)/(1+(dose/b4)^b3))^2 ),
start = c(b1 = 2.2, b2 = 1700, b3 = 2.8, b4 = 0.28, g = 2.7, k = 1),
data = chlorsulfuron, hoa = TRUE, trace = TRUE,
control = list(x.tol = 10^-12, rel.tol = 10^-12, step.min = 10^-12) )
|
Example output
Loading required package: statmod
Loading required package: survival
Package "nlreg" 1.2-2.2 (2019-01-30)
Copyright (C) 2000-2019 R. Bellio & A. R. Brazzale
This is free software, and you are welcome to redistribute
it and/or modify it under the terms of the GNU General
Public License published by the Free Software Foundation.
Package "nlreg" comes with ABSOLUTELY NO WARRANTY.
type `help(package="nlreg")' for summary information
Attaching package: ‘boot’
The following object is masked from ‘package:survival’:
aml
differentiating mean and variance function -- may take a while
differentiating mean and variance function -- may take a while
differentiating mean and variance function -- may take a while
iteration 1 : log likelihood = -3.264223
iteration 2 : log likelihood = -3.263662
iteration 3 : log likelihood = -3.263662
differentiating mean and variance function -- may take a while
iteration 1 : log likelihood = -35.12366
iteration 2 : log likelihood = -35.10621
iteration 3 : log likelihood = -35.10528
iteration 4 : log likelihood = -35.10515
iteration 5 : log likelihood = -35.10513
iteration 6 : log likelihood = -35.10513
iteration 7 : log likelihood = -35.10513
iteration 8 : log likelihood = -35.10513
iteration 9 : log likelihood = -35.10513
iteration 10 : log likelihood = -35.10513
differentiating mean and variance function -- may take a while
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.