rout_fitter {OptimModel}

Fit Model with ROUT

Description

The rout_fitter method in R fits nonlinear regression using the ROUT method as described in the reference below. The starting point is to fit a robust nonlinear regression approach assuming the Lorentzian distribution. An adaptive method then proceeds. The False Discovery Rate is used to detect outliers and the method fits in an iterative fashion.

The rout_fitter function provides a wrapper algorithm to the optim function in package stats.

Usage

   rout_fitter(theta0 = NULL, f.model, x, y, lbs = FALSE, ntry = 0, tol = 1e-3, 
               Q = 0.01, check.pd.tol = 1e-8, ...)

Arguments

theta0

a numeric vector of starting values. Alternatively, if given as NULL, theta0 can be computed within [rout.fitter()] if a starting values function is supplied as attr(f.model, "start"), as a function of x and y. If theta0 is user supplied, the last entry of theta0 should be log(sigma), where sigma = residual standard deviation. Otherwise, log(sigma) will be estimated and appended to the results from attr(f.model, "start").

f.model

Name of mean model function. See detials below.

x

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

y

a numeric vector for the response variable.

lbs

Parmeter

ntry

Parmeter

tol

Tolerance for optim algorithm.

Q

The test size for ROUT testing.

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

[rout.fitter()] is a wrapper for [optim()], specifically for Cauchy likelihood linear and non-linear regression. The Default algorithm uses method=“BFGS” or “L-BFGS-B”, if lower= and upper= arguments are specified. These can easily be overridden using the “...”.

The assumed model is:

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

After the Cauchy likelihood model is fitted to data, the residuals are interrogated to determine which observations might be outliers. An FDR correction is made to p-values (for outlier testing) through the p.adjust(method="fdr") function of the stats package.

The package supports several mean model functions for the f.model parameter. This includes beta_model, exp_2o_decay, exp_decay_pl, gompertz_model, hill_model, hill_quad_model, hill_switchpoint_model, hill5_model and linear_model.

Value

An object with class “rout_fit” is returned that gives a list with the following components and attributes:

par = estimated Cauchy model coefficients. The last term is log(sigma)

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)

Converge = logical value to indicate hessian is positive definite

call = [rout.fitter()] function call

residuals = model residuals

rsdr = robust standard deviation from ROUT method

sresiduals = residuals/rsdr

outlier = logical vector. TRUE indicates observation is an outlier via hypothesis testing with unadjust p-values.

outlier.adj = logical vector. TRUE indicates observation is an outlier via hypothesis testing with FDR-adjust p-values.

attr(object, "Q") = test size for outlier detection

Author(s)

Steven Novick

References

Motulsky, H.J. and Brown, R.E. (2006) Detecting Outliers When Fitting Data with Nonlinear Regression: A New Method Based on Robust Nonlinear Regression and the False Discovery Rate. BMC Bioinformatics, 7, 123.

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

Examples

set.seed(123L)
x = rep( c(0, 2^(-4:4)), each=4 )
theta = c(0, 100, log(.5), 2)
y = hill_model(theta, x) + rnorm( length(x), sd=2 )

rout_fitter(NULL, hill_model, x=x, y=y)
rout_fitter(c( theta, log(2) ), hill_model, x=x, y=y)

ii = sample( 1:length(y), 2 )
y[ii] = hill_model(theta, x[ii]) + 5.5*2 + rnorm( length(ii), sd=2 )

rout_fitter(c( theta, log(2) ), hill_model, x=x, y=y, Q=0.01)
rout_fitter(c( theta, log(2) ), hill_model, x=x, y=y, Q=0.05)
rout_fitter(c( theta, log(2) ), hill_model, x=x, y=y, Q=0.0001)

## Use optim method="L-BFGS-B"
rout_fitter(NULL, hill_model, x=x, y=y, Q=0.01, lower=c(-2, 95, NA, 0.5), upper=c(5, 110, NA, 4) )