optim_fit {OptimModel}

Fit Model with optim

Description

Fit nonlinear model using the optim function in the stats library. This defaults to Ordinary Least Squares (OLS) The other options are Iterative Reweighted Least Squares (IRWLS), and Maximum Likelihood Estimator (MLE).

Usage

optim_fit( theta0, f.model, gr.model=NULL, x, y, wts,
            fit.method=c("ols", "irwls", "mle"),
            var.method=c("hessian", "normal", "robust"),
            phi0=NULL, phi.fixed=TRUE, conf.level = 0.95, tol = 1e-3, 
            n.start=1000, ntry=1, start.method=c("fixed", "random"), 
            until.converge=TRUE, check.pd.tol = 1e-8, ...)


robust_fit(theta0, f.model, gr.model=NULL, x, y, wts=c("huber", "tukey"), 
    var.method=c("hessian", "normal", "robust"), conf.level=.95, tol=1e-3, ...)             

Arguments

theta0

starting values. Alternatively, if given as NULL, theta0 can be computed within optim.fit() if a starting values function is supplied as attr(f.model, "start"), as a function of x and y.

f.model

Name of mean model function.

gr.model

If specified, name of the partial derivative of f.model with respect to its parameter argument. If not specified, f2djac will approximate the derivative. Alternatively, the gradient of f.model can be specified as attr(f.model, "gradient")

x

Explanatory variable(s). Can be vector, matrix, or data.frame

y

Response variable.

fit.method

"ols" for ordinary least squares, "irwls" for iterative re-weighted least squares, "mle" for maximum likelihood.

wts

For optim.fit(), can be a numeric vector or a function. Functions supplied in the library are weights_varIdent, weights_tukey_bw, weights_huber, weights_varExp, weights_varPower, and weights_varConstPower. For robust_fit(), choices are a character string of "huber" for weights_huber and "tukey" for weights_tukey_bw

var.method

Method to compute variance-covariance matrix of estimated model parameters. Choices are "hessian" to use the hessian inverse, "normal" to use the so-called 'folk-lore' theorem estimator, and "robust" to use a sandwich-variance estimator. When fit.method = "mle", var.method = "hessian" and cannot be overridden.

phi0

Not meaningful for fit.method = "ols". Starting value(s) for variance parameters (for weights).

phi.fixed

Not meaningful for fit.method = "ols". If set to TRUE, the variance parameter(s) remain fixed at the given starting value, phi0. Otherwise, the variance parameter(s) are estimated.

conf.level

Confidence level of estimated parameters.

tol

Tolerance for optim algorithm.

n.start

Number of starting values to generate (see details).

ntry

Maximum number of calls to optim.fit().

start.method

Parameter

until.converge

Logical (TRUE/FALSE) indicating when algorithm should stop.

check.pd.tol

absolute tolarence for determing whether a matrix is positive definite. Refer to test_fit.

...

Other arguments to passed to optim. See ?optim. For example, lower=, upper=, method=

Details

optim_fit() is a wrapper for stats::optim(), specifically for non-linear regression. The Default algorithm is ordinary least squares (ols) using method="BFGS" or "L-BFGS-B", if lower= and upper= are specified. These can easily be overridden.

The assumed model is:

y = \text{f.model}(\theta, x) + g(\theta, \phi, x)\sigma\epsilon\text{, where }\epsilon\sim N(0, 1).

Usually the function g(\cdot) = 1.

With the exception of weights.tukey.bw and weights.huber, the weights functions are equivalent to g(\theta, \phi, x).

Note that robust_fit() is a convenience function to simplify the model call with fit.method = "irwls", phi0 = NULL, and phi.fixed = TRUE.

Algorithms:

1. OLS. Minimize sum( (y - f.model(theta, x))^2 )

2. IRWLS. Minimize sum( g(theta, phi, x)*(y - f.model(theta, x))^2 ), where g(theta, phi, x) act as weights. See section on Variance functions below for more details on g(\cdot).

3. MLE. Minimize the -log(Likelihood). See section on Variance functions below for more details on g(\cdot).

Variance functions:

Weights are given by some variance function. Some common variance functions are supplied.

See weights_varIdent, weights_varExp, weights_varPower, weights_varConstPower, weights_tukey_bw, weights_huber.

User-specified variance functions can be provided for weighted regression.

Value

The returned object is a list with the following components and attributes:

coefficients = estimated model coefficients

value, counts, convergence = returns from optim()

message = character, indicating problem if any. otherwise=NULL

hessian = hessian matrix of the objective function (e.g., sum of squares)

fit.method = chosen fit.method (e.g., "ols")

var.method = chosen var.method (e.g., "hessian")

call = optim.fit() function call

fitted, residuals = model mean and model residuals

r.squared, bic = model statistics

df = error degrees of freedom = N - p, where N = # of observations and p = # of parameters

dims = list containing the values of N and p

sigma = estimated standard deviation parameter

varBeta = estimated variance-covariance matrix for the coefficients

beta = data.frame summary of the fit

attr(object, "weights") = weights

attr(object, "w.func") = weights model for the variance

attr(object, "var.param") = variance parameter values

attr(object, "converge.pls") = logical indicating if IRWLS algorithm converged.

Author(s)

Steven Novick

See Also

optim_fit, rout_outlier_test, beta_model, exp_2o_decay, exp_decay_pl, gompertz_model, hill_model, hill5_model, hill_quad_model, hill_switchpoint_model, linear_model, weights_huber, weights_tukey_bw, weights_varConstPower, weights_varExp, weights_varIdent, weights_varPower

Examples

set.seed(123L)
x = rep( c(0, 2^(-4:4)), each=4 )
theta = c(0, 100, log(.5), 2)
y1 = hill_model(theta, x) + rnorm( length(x), sd=2 )
y2 = hill_model(theta, x) + rnorm( length(x), sd=.1*hill_model(theta, x) )
wts = runif( length(y1) )
fit1=optim_fit(theta, hill_model, x=x, y=y1)
fit2=optim_fit(theta, hill_model, x=x, y=y1, wts=weights_varIdent)
fit3=optim_fit(theta, hill_model, x=x, y=y1, wts=wts)
fit4=optim_fit(theta, hill_model, attr(hill_model, "gradient"), x=x, y=y1, wts=wts)

fit5=optim_fit(NULL, hill_model, x=x, y=y2, wts=weights_varPower, fit.method="irwls")
fit6=optim_fit(theta, hill_model, x=x, y=y2, wts=weights_varPower, fit.method="mle")

fit7=optim_fit(theta, hill_model, x=x, y=y2, wts=weights_varPower, fit.method="irwls", 
               phi0=0.5, phi.fixed=FALSE)
fit8=optim_fit(theta, hill_model, x=x, y=y2, wts=weights_varPower, fit.method="mle", 
              phi0=0.5, phi.fixed=FALSE)

fit9a=optim_fit(theta, hill_model, x=x, y=y1, wts=weights_tukey_bw, fit.method="irwls", 
             phi0=4.685, phi.fixed=TRUE)

fit9b=robust_fit(theta, hill_model, x=x, y=y1, wts="tukey")