Beta hook-effect model, gradient, starting values, and back-calculation functions

Description

Five-parameter hook-effect model for dose-response curve fitting

Usage

 
       beta_model(theta, x)

Arguments

Details

The five-parameter Beta model is given by:

y = e_{\min} + e_{\max} \times \exp( log( \beta(\delta_1, \delta_2) ) + \delta_1 \times\log(x) + \delta_2*\log(\text{sc}-x) - (\delta_1+\delta_2)\times\log(\text{sc})

where

\beta(\delta_1, \delta_2) = (\delta_1+\delta_2)^(\delta_1+\delta_2) /(\delta_1^{\delta_1} \times \delta_2^{\delta_2})

and

\text{sc} = \max(x) + \delta_3.

Note that the Beta model depends on the maximum x value. For a particular data set, this may be set by

attr(theta), "maxX") = max(x).

Value

Let N = length(x). Then

• beta_model(theta, x) returns a numeric vector of length N.

• attr(beta_model, "gradient")(theta, x) returns an N x 5 matrix.

• attr(beta_model, "start")(x, y) returns a numeric vector of length 5 with starting values for

  (e_{\min}, e_{\max}, \log(\delta_1), \log(\delta_2), \log(\delta_3)).

• attr(beta_model, "backsolve")(theta, y) returns a numeric vector of length=length(y) with the first x such that beta_model(theta, x)=y.

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123L)
x = rep( c(0, 2^(-4:4)), each=4 )
theta = c(emin=0, emax=115, ldelta1=-1.5, ldelta2=9, ldelta3=11.5)
y = beta_model(theta, x)  + rnorm( length(x), mean=0, sd=1 )

beta_model(theta, x)
attr(beta_model, "gradient")(theta, x)
attr(beta_model, "start")(x, y)

attr(theta, "maxX") = max(x)
attr(beta_model, "backsolve")(theta, 50)

Five-parameter second-order exponential decay, gradient, starting values, and back-calculation functions

Description

Five-parameter second-order exponential decay, gradient, starting values, and back-calculation functions.

Usage

 
        exp_2o_decay(theta, x)

Arguments

Details

The five-parameter exponential decay model is given by:

y = A + B\times P\times \exp(-K1 \times x) + B\times (1 - P)\times \exp(-K2\times x)

The parameter vector is (A, B, k1, k2, p) where A = \min y (min y value), A+B = \max y (max y value), K1 = \exp(k1) which is the shape parameter for first term, K2 = \exp(k2) which is the shape parameter for second term, and P = 1/(1+\exp(p)) which is the proportion of signal from the first term.

Value

Let N = length(x). Then

• exp_2o_decay(theta, x) returns a numeric vector of length N.

• attr(exp_2o_decay, "gradient")(theta, x) returns an N x 5 matrix.

• attr(exp_2o_decay, "start")(x, y) returns a numeric vector of length 5 with starting values for (A, B, k1, k2, p).

• attr(exp_2o_decay, "backsolve")(theta, y) returns a numeric vector of length = length(y).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123L)
x = 2^(-4:4)
theta = c(25, 75, log(3), log(1.2), 1/(1+exp(.7)))
y = exp_2o_decay(theta, x) + rnorm( length(x), mean=0, sd=1 )
attr(exp_2o_decay, "gradient")(theta, x)
attr(exp_2o_decay, "start")(x, y)
attr(exp_2o_decay, "backsolve")(theta, 38)

Three-parameter exponential decay, gradient, starting values, and back-calculation functions

Description

Three-parameter exponential decay, gradient, starting values, and back-calculation functions.

Usage

 
        exp_decay(theta, x)

Arguments

Details

The three-parameter exponential decay model is given by:

y = A + B \times \exp(-Kx).

The parameter vector is (A, B, k) where A =\min y ( minimum y value), A + B = \max y (maximum y value), and K = \exp(k) whichi is the shape parameter.

Value

Let N = length(x). Then

• exp_decay(theta, x) returns a numeric vector of length N.

• attr(exp_decay, "gradient")(theta, x) returns an N x 3 matrix.

• attr(exp_decay, "start")(x, y) returns a numeric vector of length 3 with starting values for (A, B, k).

• attr(exp_decay, "backsolve")(theta, y) returns a numeric vector of length=length(y).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123L)
x = 2^(-4:4)
theta = c(25, 75, log(3))
y = exp_decay(theta, x)  + rnorm( length(x), mean=0, sd=1 )
attr(exp_decay, "gradient")(theta, x)
attr(exp_decay, "start")(x, y)
attr(exp_decay, "backsolve")(theta, 38)

Three-parameter exponential decay with initial plateau, gradient, starting values, and back-calculation functions

Description

Three-parameter exponential decay with initial plateau, gradient, starting values, and back-calculation functions.

Usage

 
        exp_decay_pl(theta, x)

Arguments

Details

The three-parameter exponential decay with initial plateau model is given by y = \text{yMax} whenever x\le 0 otherwise

y = \text{yMin} + (\text{yMax}-\text{yMin})\times\exp(-K(x-X0))\text{ if } x > X0,

where X0 = \exp(x0) is an inflection point between plateau and exponential decay curve, \text{yMin} = \min y (min response), \text{yMax} = \max y (maximum response), and K = \exp(k) is the shape parameter.

Value

Let N = length(x). Then

• exp_decay_pl(theta, x) returns a numeric vector of length N.

• attr(exp_decay_pl, "gradient")(theta, x) returns an N x 4 matrix.

• attr(exp_decay_pl, "start")(x, y) returns a numeric vector of length 4 with starting values for (x0, yMax, yMin, k).

• attr(exp_decay_pl, "backsolve")(theta, y) returns a numeric vector of length=length(y).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(100)
x = 2^(-4:4)
theta = c(0.4, 75, 10, log(3))
y = exp_decay_pl(theta, x)  + rnorm( length(x), mean=0, sd=1 )
attr(exp_decay_pl, "gradient")(theta, x)
attr(exp_decay_pl, "start")(x, y)
attr(exp_decay_pl, "backsolve")(theta, 38)

Compute derivative with respect to parameters

Description

Compute derivative with respect to parameters.

Usage

f2djac(Func, theta, ...)

Arguments

Value

Returns an n x p matrix of derivatives with respect to theta. Computes \frac{\delta \, Func( \theta, ...)}{ \delta \, \theta}, where \theta = theta

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

    f = function(theta, x){ theta[1] + (theta[2]-theta[1])/(1 + (x/theta[3])^theta[4]) }
    theta0 = c(0, 100, .5, 2)
    x = 0:10
    f2djac(f, theta0, x=x)

Extract data object from an optim fit

Description

Extract data object from an optim_fit object.

Usage

  getData(object)

Arguments

Value

Returns a data frame with elements x and y.

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

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 )
fit = optim_fit(c(0, 100, .5, 1), f.model=hill_model, x=x, y=y)
d=getData(fit)

Compute standard error for a function of model parameter estimates

Description

Compute standard error for a function of model parameter estimates via the delta method.

Usage

get_se_func(object, Func, ..., level=0.95)

Arguments

Details

Func is of the form function(theta, ...). For example,

Func = function(theta, x){ exp(theta[1])*log(x)/theta[2] }

Value

Returns a data.frame with a single row for the estimated Func call (Est), its standard error (SE), and a confidence interval (lower, upper).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

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 )
fit = optim_fit(theta, hill_model, x=x, y=y)
  
## Get SE for IC20 and IC40
ic.z = function(theta, z){  attr(hill_model, "backsolve")(theta, z)  }
get_se_func(object=fit, Func=ic.z, z=20)
get_se_func(object=fit, Func=ic.z, z=40)
  

Four-parameter Gompertz model, gradient, starting values, and back-calculation functions

Description

Four-parameter Gompertz model, gradient, starting values, and back-calculation functions.

Usage

 
        gompertz_model(theta, x)

Arguments

Details

The four parameter Gompertz model is given by:

y = A + (B-A)\times\exp( -\exp( m(x-\text{offset}) ) )\text{, where}

A = \min y (minimum y value), A+(B-A)\exp(-\exp( -m*\text{offset} )) is the maximum y value, m is the shape parameter, and offset shifts the curve, relative to the concentration x.

Value

Let N = length(x). Then

• gompertz_model(theta, x) returns a numeric vector of length N.

• gompertz_model(hill_model, "gradient")(theta, x) returns an N x 4 matrix.

• attr(gompertz_model, "start")(x, y) returns a numeric vector of length 4 with starting values for (A, B, m, offset).

• attr(gompertz_model, "backsolve")(theta, y) returns a numeric vector of length=length(y).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(100)
x = rep( c(0, 2^(-4:4)), each=4 )
theta = c(0, 100, log(.5), 2)
y = gompertz_model(theta, x)  + rnorm( length(x), mean=0, sd=1 )
attr(gompertz_model, "gradient")(theta, x)
attr(gompertz_model, "start")(x, y)
attr(gompertz_model, "backsolve")(theta, 50)

Five-parameter Hill model, gradient, starting values, and back-calculation functions

Description

Five-parameter Hill model, gradient, starting values, and back-calculation functions.

Usage

 
        hill5_model(theta, x)

Arguments

Details

The five parameter Hill model is given by:

y = e_{\min} + \frac{e_{\max}-e+{\min}}{ 1 + \exp( m\log(x) - m\text{ log.ic50}) )^{\exp(\text{log.sym})}}

e_{\min} = \min y (minimum y value), e_{\max} = \max y (maximum y value), \text{log.ic50} = \log( \text{ic50} ), m is the shape parameter, and \text{log.sym} = \log( \text{symmetry parameter}).

Note: ic50 is defined such that hill5_model(theta, ic50) = e_{\min}+(e_{\max}-e_{\min})/2^{\exp(\text{log.sym})}

Value

Let N = length(x). Then

• hill5_model(theta, x) returns a numeric vector of length N.

• attr(hill5_model, "gradient")(theta, x) returns an N x 5 matrix.

• attr(hill5_model, "start")(x, y) returns a numeric vector of length 5 with starting values for (e_{\min}, e_{\max}, \text{log.ic50}, m, \text{log.sym}).

• attr(hill5_model, "backsolve")(theta, y) returns a numeric vector of length=length(y).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123L)
x = rep( c(0, 2^(-4:4)), each=4 )
theta = c(0, 100, log(.5), 2, log(10))
y = hill5_model(theta, x)  + rnorm( length(x), mean=0, sd=1 )
attr(hill5_model, "gradient")(theta, x)
attr(hill5_model, "start")(x, y)
attr(hill5_model, "backsolve")(theta, 50)

Four-parameter Hill model, gradient, starting values, and back-calculation functions

Description

Four-parameter Hill model, gradient, starting values, and back-calculation functions.

Usage

 
        hill_model(theta, x)

Arguments

Details

The four parameter Hill model is given by:

y = e_{\min} + \frac{(e_{\max}-e_{\min})}{( 1 + \exp( m\log(x) - m*\text{lec50} ) )}\text{, where }

e_{\min} = \min y (minimum y value), e_{\max} = \max y (maximum y value), \text{lec50} = \log( \text{ec5} ), and m is the shape parameter. Note: ec50 is defined such that hill.model(theta, ec50) = .5*( emin+ emax ).

Value

Let N = length(x). Then

• hill_model(theta, x) returns a numeric vector of length N.

• attr(hill_model, "gradient")(theta, x) returns an N x 4 matrix.

• attr(hill_model, "start")(x, y) returns a numeric vector of length 4 with starting values for (e_{\min}, e_{\max}, \text{lec50}, m).

• attr(hill_model, "backsolve")(theta, y) returns a numeric vector of length=length(y).

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

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), mean=0, sd=1 )
attr(hill_model, "gradient")(theta, x)
attr(hill_model, "start")(x, y)
attr(hill_model, "backsolve")(theta, 50)

Five-parameter Hill model with quadratic component, gradient, starting values, and back-calculation functions

Description

Five-parameter Hill model with quadratic component, gradient, starting values, and back-calculation functions.

Usage

 
        hill_quad_model(theta, x)

Arguments

Details

The five parameter Hill model with quadratic component is given by:

y = A + \frac{B-A}{( 1 + \exp( -(a + bz + cz^2) ) )}\text{, where }z = \log(x)

A =\min y ( minimum y value), B = \max y (maximum y value), (a, b, c) are quadratic parameters for \log(x).

Notes:

1. If c = 0, this model is equivalent to the four-parameter Hill model (hill.model).

2. The ic50 is defined such that a + bz + cz^2 = 0. If the roots of the quadratic equation are real, then the ic50 is given by \tfrac{-b \pm\sqrt{b^2 - 4ac }}{2a}.

Value

Let N = length(x). Then

• hill_quad_model(theta, x) returns a numeric vector of length N.

• attr(hill_quad_model, "gradient")(theta, x) returns an N x 5 matrix.

• attr(hill_quad_model, "start")(x, y) returns a numeric vector of length 5 with starting values for (A, B, a, b, c).

  If the quadratic roots are real-valued, attr(hill_quad_model, "backsolve")(theta, y) returns a numeric vector of length=2.

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123L)
x = rep( c(0, 2^(-4:4)), each=3 )      ## Dose
theta = c(0, 100, 2, 1, -0.5)          ## Model parameters
y = hill_quad_model(theta, x) + rnorm( length(x), mean=0, sd=5 )

## Generate data
hill_quad_model(theta, x)
attr(hill_quad_model, "gradient")(theta, x)
attr(hill_quad_model, "start")(x, y)
attr(hill_quad_model, "backsolve")(theta, 50)

Five-parameter Hill model with switch point component, gradient, starting values, and back-calculation functions

Description

Five-parameter Hill model with switch point component, gradient, starting values, and back-calculation functions.

Usage

 
        hill_switchpoint_model(theta, x)

Arguments

Details

The five parameter Hill model with switch point component is given by:

y = e_{\min} + \frac{e_{\max}-e_{\min}}{1 + \exp( m\times f(\exp(\text{lsp}), x)\times(\log(x) - \text{log.ic50}) )}\text{, where}

e_{\min} = \min y (minimum y value), e_{\max} = \max y (maximum y value), \text{log.ic50} = \log( \text{ic50} ), m is the shape parameter, and f(s, x) is the switch point function.

The function f(s, x) = (s-x)/(s+x) = 2/(1+x/s) - 1. This function is constrained to be between -1 and +1 with s > 0.

Notes:

1. At x = 0, f(s, x) = 1, which reduces to hill_model(theta[1:4], 0).

2. The hill_switchpoint_model() is more flexible compared to hill_quad_model().

3. When the data does not contain a switchpoint, then lsp should be a large value so that f(\exp(\text{lsp}), x) will be near 1 for all x.

Value

Let N = length(x). Then

• hill_switchpoint_model(theta, x) returns a numeric vector of length N.

• attr(hill_switchpoint_model, "gradient")(theta, x) returns an N x 5 matrix.

• attr(hill_switchpoint_model, "start")(x, y) returns a numeric vector of length 5 with starting values for (e_{\min}, e_{\max}, \text{lec50}, m, \text{lsp}).

  Because hill_switchpoint_model() can be fitted to biphasic data with a hook-effect, attr(hill_switchpoint_model, "backsolve")(theta, y0) returns the first x that satisfies y0=hill_switchpoint_model(theta, x)

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123L)
x = rep( c(0, 2^(-4:4)), each=3 )      ## Dose
## Create model with no switchpoint term
theta = c(0, 100, log(.5), 2)
y = hill_model(theta, x) + rnorm( length(x), mean=0, sd=5 )

## fit0 and fit return roughly the same r-squared and sigma values.
## Note that BIC(fit0) < BIC(fit), as it should be.
fit0 = optim_fit(NULL, hill_model, x=x, y=y)
fit = optim_fit(c(coef(fit0), lsp=0), hill_switchpoint_model, x=x, y=y)
fit = optim_fit(NULL, hill_switchpoint_model, x=x, y=y) 

## Generate data from hill.quad.model() with a biphasic (hook-effect) profile
set.seed(123L)
theta = c(0, 100, 2, 1, -0.5)          ## Model parameters
y = hill_quad_model(theta, x) + rnorm( length(x), mean=0, sd=5 )

## fit.qm and fit.sp return nearly identical fits
fit.qm = optim_fit(theta, hill_quad_model, x=x, y=y)  
fit.sp = optim_fit(NULL, hill_switchpoint_model, x=x, y=y, ntry=50)  

plot(log(x+0.01), y)
lines(log(x+0.01), fitted(fit.qm))
lines(log(x+0.01), fitted(fit.sp), col="red")

## Generate data from hill.switchback.model()
set.seed(123)
theta = c(0, 100, log(0.25), -3, -2)
y = hill_switchpoint_model(theta, x) + rnorm( length(x), mean=0, sd=5 )
plot( log(x+0.01), y )

## Note that this model cannot be fitted by hill.quad.model()
fit = optim_fit(NULL, hill_switchpoint_model, x=x, y=y, ntry=50, 
       start.method="fixed", until.converge=FALSE)
pred = predict(fit, x=exp(seq(log(0.01), log(16), length=50)), interval='confidence')

plot(log(x+0.01), y, main="Fitted curve with 95% confidence bands")
lines(log(pred[,'x']+0.01), pred[,'y.hat'], col='black')
lines(log(pred[,'x']+0.01), pred[,'lower'], col='red', lty=2)
lines(log(pred[,'x']+0.01), pred[,'upper'], col='red', lty=2)


## Other functions
hill_switchpoint_model(theta, x)
attr(hill_switchpoint_model, "gradient")(theta, x)
attr(hill_switchpoint_model, "start")(x, y)
attr(hill_switchpoint_model, "backsolve")(theta, 50)

Linear model, gradient, and starting values

Description

Linear model, gradient, and starting values.

Usage

 
        linear_model(theta, x)

Arguments

Details

The linear model is given by:

y = x * \theta \text{, where}

x is a matrix, possibly generated from model.matrix() \theta is a vector of linear parameters

Value

Let N = nrow(x). Then

• linear_model(theta, x) returns a numeric vector of length N.

• attr(linear_model, "gradient")(theta, x) returns x.

• attr(linear_model, "start")(x, y) returns solve(t(x) * x) * t(x) * y

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Examples

set.seed(123)
d = data.frame( Group=factor(rep(LETTERS[1:3], each=5)), age=rnorm(15, mean=20, sd=3) )
d$y = c(80, 100, 120)[unclass(d$Group)] - 3*d$age + rnorm(nrow(d), mean=0, sd=5)

X = model.matrix(~Group+age, data=d)  ## This is the "x" in linear.model()
theta = c(80, 20, 40, -3)   ## Intercept, effect for B, effect for C, slope for age
linear_model(theta, x=X)
attr(linear_model, "gradient")(theta, x=X)
attr(linear_model, "start")(x=X, y=d$y)

Negative log-likelihood function for Cauchy regression

Description

The negative log-likelihood function for Cauchy regression, for use with rout_fitter. Usually not called by the end user.

Usage

 
       nlogLik_cauchy(theta, x, y, f.model, lbs)

Arguments

Details

The function provides the negative log-likelihood for Cauchy regression

Let mu = f.model(theta[1:(p-1)], x) and sigma = exp(theta[p]), where p = number of parameters in theta.

The Cauchy likelihood is

L = \prod \frac{1}{\pi \sigma} ( 1 + ( \frac{y_i - \mu_i}{\sigma} )^2 )^{-1}

.

The function returns \log(L).

If lbs == TRUE, then \mu is replaced with \log(mu).

Value

Returns a single numerical value.

Author(s)

Steven Novick

See Also

rout_fitter

Examples

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

theta1 = c(theta, lsigma=log(2))
nlogLik_cauchy(theta1, x=x, y=y, f.model=hill_model, lbs=FALSE)

  ## Cauchy regression via maximum likelihood
optim( theta1, nlogLik_cauchy, x=x, y=y, f.model=hill_model, lbs=FALSE )

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

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

Weight functions for optim_fit

Description

Weight functions for optim_fit. May be used when fit.method=="irwls" or fit.method=="mle". Generally, not called by the user.

Usage

 
       weights_varIdent(phi, mu)
       weights_varExp(phi, mu)
       weights_varPower(phi, mu)
       weights_varConstPower(phi, mu)
       weights_tukey_bw(phi = 4.685, resid)
       weights_huber(phi=1.345, resid)

Arguments

Details

• weights_varIdent returns a vector of ones.

• weights_varExp returns \exp(\phi*\mu)

• weights_varPower returns |\mu|^{\phi}

• weights_varConstPower returns \phi_1 + |\mu|^{\phi_2} where \phi_i =\phi[i]

• weights_tukey_bw is a Tukey bi-weight function. Let

  r=\tfrac{|\text{resid}|}{\text{mad}(\text{resid},\text{center}=\text{TRUE})}.

  Then this function returns

  \left(1-\left(\tfrac{r}{\phi}\right)^2\right)^2\text{ whenever } r <= \phi\text{ and }0{ o.w.}

  For this the user should use phi.fixed=TRUE in the optim_fit function.

• weights_huber is a Huber weighting function that returns \min(1, \phi/r), where r = |\text{resid}|/\text{sig} and \text{sig} = \text{mad}(\text{resid}, \text{center} = \text{TRUE}). For this the user should use phi.fixed = TRUE in the optim_fit function.

Value

A vector of numeric weights.

Author(s)

Steven Novick

See Also

optim_fit, rout_fitter

Predicted values for optim.fit objects

Description

Provides predicted values, standard errors, confidence intervals and prediction intervals for optim_fit objects.

Usage

 
    ## S3 method for class 'optim_fit'
predict(object, x, se.fit=FALSE, 
    interval=c("none", "confidence", "prediction"), K = 1, level = 0.95,...)

Arguments

Value

Returns a vector (if interval='none'). Otherwise returns a data.frame with possible columns 'x', 'y.hat', 'se.fit', 'lower', and 'upper'.

Author(s)

Steven Novick

See Also

optim_fit

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 )


fit1=optim_fit(theta, hill_model, x=x, y=y1)
fitted(fit1)
predict(fit1)
predict(fit1, x=x)
predict(fit1, x=seq(0, 1, by=.1), se.fit=TRUE)
predict(fit1, x=seq(0, 1, by=.1), interval="conf")
predict(fit1, x=seq(0, 1, by=.1), interval="pred")

Prints optim_fit objects

Description

Provides a nice printed summary of optim_fit objects.

Usage

 
  ## S3 method for class 'optim_fit'
print(x, digits=4,...)

Arguments

Value

No Return Value.

Author(s)

Steven Novick

See Also

optim_fit

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 )

fit1=optim_fit(theta, hill_model, x=x, y=y1)
print(fit1)
fit1

Residuals for optim.fit objects

Description

Provides raw and studentized residuals for optim_fit objects.

Usage

 
    ## S3 method for class 'optim_fit'
residuals(object, type=c("raw", "studentized"), ...)

Arguments

Value

Returns a numeric vector.

Author(s)

Steven Novick

See Also

optim_fit

Examples

set.seed(123)

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 )

fit1=optim_fit(theta, hill_model, x=x, y=y1)
residuals(fit1)
residuals(fit1, type="s")

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

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

ROUT Outlier Testing

Description

Perform ROUT outlier testing on rout.fitter object.

Usage

rout_outlier_test(fit, Q = 0.01)

Arguments

Details

rout_outlier_test() is typically called from rout_fitter(), but may also be called directly by the user.

Value

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

See Also

rout_fitter

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 )

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

fit = rout_fitter(c( theta, log(2) ), hill_model, x=x, y=y, Q=0.01)
rout_outlier_test(fit, Q=0.001)

Test Fit Parameters

Description

Test if estimated parameters optimize the regression system (i.e., minimize sums of squares, maximize likelihood).

Usage

 
       test_fit(obj, check.pd.tol = 1e-8)

Arguments

Details

The function checks if optim convergence has been reached and also checks if the cholesky decompoisition of the Hessian matrix is positive definite. The latter is an indication that optimization has been reached. Sometimes the chol decomposition check doesn't work and to enforce that constriant we use the check.pd.tol to make sure all the eigenvalues are larger than this minimum threshhold.

Value

Returns a TRUE or FALSE as to whether or not hessian component of the object is Positive Definite.

Author(s)

Steven Novick

See Also

optim_fit

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 )
wts = runif( length(y1) )
fit1=optim_fit(theta, hill_model, x=x, y=y1)
test_fit(fit1)