nlreg {nlreg}

Fit a Nonlinear Heteroscedastic Model via Maximum Likelihood

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

nlreg(formula, weights = NULL, data = sys.frame(sys.parent()), start, 
      offset = NULL, subset = NULL, 
      control = list(x.tol = 1e-06, rel.tol = 1e-06, 
                     step.min = 1/2048, maxit = 100), trace = FALSE, 
      hoa = FALSE)

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 \sigma^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 \sigma^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 \sigma^2 in Var(error) = s^2 V(predictor) is included by default. In order to work with a real value, \sigma^2 is estimated on the logarithmic scale, that is, the model is reparametrized into \log(\sigma^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 \sigma^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

nlreg.object, nls

Examples

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) )