Function to Minimize Jaeckel's Dispersion Function

Description

Uses the built-in function optim to minimize Jaeckel's dispersion function. Alternates between CG and BFGS steps and initially takes a number of IRLS steps.

Usage

jaeckel(x, y, beta0 = rq(y ~ x)$coef[2:(ncol(x) + 1)], scores = Rfit::wscores,...)

Arguments

Details

This function is meant to mimic the minization algorithm implemented in RGLM (Kapenga, et. al. 1988). The main loop of the function alternates between CG and BFGS estimation methods. To speed up convergence, it first takes a number of iterated reweighted least squares (IRLS) steps which generally gets close to the solution in a small number of steps. Using IRLS for rank regression was first considered by Cheng and Hettmansperger (1983). See also Sievers and Abebe (2004).

Value

Results of optim are returned.

Author(s)

John Kloke kloke@biostat.wisc.edu

References

Cheng, K. S. and Hettmansperger, T. P. (1983), Weighted Least-Squares Rank Regression, Communications in Statistics, Part A - Theory and Methods, 12, 1069-1086.

Hettmansperger, T.P. and McKean J.W. (2011), Robust Nonparametric Statistical Methods, 2nd ed., New York: Chapman-Hall.

Jaeckel, L. A. (1972), Estimating regression coefficients by minimizing the dispersion of residuals. Annals of Mathematical Statistics, 43, 1449 - 1458.

Kapenga, J. A., McKean, J. W., and Vidmar, T. J. (1988), RGLM: Users Manual, Statist. Assoc. Short Course on Robust Statistical Procedures for the Analysis of Linear and Nonlinear Models, New Orleans.

Sievers, J. and Abebe, A. (2004), Rank Estimation of Regression Coefficients Using Iterated Reweighted Least Squares, Journal of Statistical Computation and Simulation, 74, 821-831.

See Also

optim, rfit

Examples

##  This is a internal function.  See rfit for user-level examples.

## The function is currently defined as
function (x, y, beta0 = rq(y ~ x - 1)$coef, scores = wscores, 
    maxiter = 100, irls0 = 10, BFGS0 = 20, stepCG = 5, stepBFGS = 2) 
{
    x <- x - outer(rep(1, nrow(x)), apply(x, 2, mean))
    beta0 <- irls(x, y, beta0, max.iter = irls0)
    if (BFGS0 < 1) 
        BFGS0 <- 1
    fit <- optim(beta0, disp, method = "BFGS", x = x, y = y, 
        scores = scores, control = list(maxit = BFGS0))
    iter <- 0
    while (fit$convergence && iter < maxiter) {
        iter <- iter + 1
        fit <- optim(fit$par, disp, method = "CG", x = x, y = y, 
            scores = scores, control = list(maxit = stepCG))
        fit <- optim(fit$par, disp, method = "BFGS", x = x, y = y, 
            scores = scores, control = list(maxit = stepBFGS))
    }
    optim(fit$par, disp, method = "BFGS", x = x, y = y, scores = scores)
  }

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Rfit)
Loading required package: quantreg
Loading required package: SparseM

Attaching package: 'SparseM'

The following object is masked from 'package:base':

    backsolve

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Rfit/jaeckel.Rd_%03d_medium.png", width=480, height=480)
> ### Name: jaeckel
> ### Title: Function to Minimize Jaeckel's Dispersion Function
> ### Aliases: jaeckel
> 
> ### ** Examples
> 
> ##  This is a internal function.  See rfit for user-level examples.
> 
> ## The function is currently defined as
> function (x, y, beta0 = rq(y ~ x - 1)$coef, scores = wscores, 
+     maxiter = 100, irls0 = 10, BFGS0 = 20, stepCG = 5, stepBFGS = 2) 
+ {
+     x <- x - outer(rep(1, nrow(x)), apply(x, 2, mean))
+     beta0 <- irls(x, y, beta0, max.iter = irls0)
+     if (BFGS0 < 1) 
+         BFGS0 <- 1
+     fit <- optim(beta0, disp, method = "BFGS", x = x, y = y, 
+         scores = scores, control = list(maxit = BFGS0))
+     iter <- 0
+     while (fit$convergence && iter < maxiter) {
+         iter <- iter + 1
+         fit <- optim(fit$par, disp, method = "CG", x = x, y = y, 
+             scores = scores, control = list(maxit = stepCG))
+         fit <- optim(fit$par, disp, method = "BFGS", x = x, y = y, 
+             scores = scores, control = list(maxit = stepBFGS))
+     }
+     optim(fit$par, disp, method = "BFGS", x = x, y = y, scores = scores)
+   }
function (x, y, beta0 = rq(y ~ x - 1)$coef, scores = wscores, 
    maxiter = 100, irls0 = 10, BFGS0 = 20, stepCG = 5, stepBFGS = 2) 
{
    x <- x - outer(rep(1, nrow(x)), apply(x, 2, mean))
    beta0 <- irls(x, y, beta0, max.iter = irls0)
    if (BFGS0 < 1) 
        BFGS0 <- 1
    fit <- optim(beta0, disp, method = "BFGS", x = x, y = y, 
        scores = scores, control = list(maxit = BFGS0))
    iter <- 0
    while (fit$convergence && iter < maxiter) {
        iter <- iter + 1
        fit <- optim(fit$par, disp, method = "CG", x = x, y = y, 
            scores = scores, control = list(maxit = stepCG))
        fit <- optim(fit$par, disp, method = "BFGS", x = x, y = y, 
            scores = scores, control = list(maxit = stepBFGS))
    }
    optim(fit$par, disp, method = "BFGS", x = x, y = y, scores = scores)
}
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>