Learn R Programming

ROptEst (version 1.0)

roptest: Optimally robust estimation

Description

Function to compute optimally robust estimates for L2-differentiable parametric families via k-step construction.

Usage

roptest(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est, neighbor = ContNeighborhood(), risk = asMSE(), steps = 1L, distance = CvMDist, startPar = NULL, verbose = NULL, OptOrIter = "iterate", useLast = getRobAStBaseOption("kStepUseLast"), withUpdateInKer = getRobAStBaseOption("withUpdateInKer"), IC.UpdateInKer = getRobAStBaseOption("IC.UpdateInKer"), withICList = getRobAStBaseOption("withICList"), withPICList = getRobAStBaseOption("withPICList"), na.rm = TRUE, initial.est.ArgList, ..., withLogScale = TRUE, ..withCheck = FALSE, withTimings = FALSE, withMDE = NULL, withEvalAsVar = NULL) roptest.old(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est, neighbor = ContNeighborhood(), risk = asMSE(), steps = 1L, distance = CvMDist, startPar = NULL, verbose = NULL, OptOrIter = "iterate", useLast = getRobAStBaseOption("kStepUseLast"), withUpdateInKer = getRobAStBaseOption("withUpdateInKer"), IC.UpdateInKer = getRobAStBaseOption("IC.UpdateInKer"), withICList = getRobAStBaseOption("withICList"), withPICList = getRobAStBaseOption("withPICList"), na.rm = TRUE, initial.est.ArgList, ..., withLogScale = TRUE)

Arguments

x
sample
L2Fam
object of class "L2ParamFamily"
eps
positive real (0 < eps
eps.lower
positive real (0 <= eps.lower <= eps.upper): lower bound for the amount of gross errors. See details below.
eps.upper
positive real (eps.lower <= eps.upper
fsCor
positive real: factor used to correct the neighborhood radius; see details.
initial.est
initial estimate for unknown parameter. If missing minimum distance estimator is computed.
neighbor
object of class "UncondNeighborhood"
risk
object of class "RiskType"
steps
positive integer: number of steps used for k-steps construction
distance
distance function
startPar
initial information used by optimize resp. optim; i.e; if (total) parameter is of length 1, startPar is a search interval, else it is an initial parameter value; if NULL slot startPar of ParamFamily is used to produce it; in the multivariate case, startPar may also be of class Estimate, in which case slot untransformed.estimate is used.
verbose
logical: if TRUE, some messages are printed
useLast
which parameter estimate (initial estimate or k-step estimate) shall be used to fill the slots pIC, asvar and asbias of the return value.
OptOrIter
character; which method to be used for determining Lagrange multipliers A and a: if (partially) matched to "optimize", getLagrangeMultByOptim is used; otherwise: by default, or if matched to "iterate" or to "doubleiterate", getLagrangeMultByIter is used. More specifically, when using getLagrangeMultByIter, and if argument risk is of class "asGRisk", by default and if matched to "iterate" we use only one (inner) iteration, if matched to "doubleiterate" we use up to Maxiter (inner) iterations.
withUpdateInKer
if there is a non-trivial trafo in the model with matrix $D$, shall the parameter be updated on $ker(D)$?
IC.UpdateInKer
if there is a non-trivial trafo in the model with matrix $D$, the IC to be used for this; if NULL the result of getboundedIC(L2Fam,D) is taken; this IC will then be projected onto $ker(D)$.
withPICList
logical: shall slot pICList of return value be filled?
withICList
logical: shall slot ICList of return value be filled?
na.rm
logical: if TRUE, the estimator is evaluated at complete.cases(x).
initial.est.ArgList
a list of arguments to be given to argument start if the latter is a function; this list by default already starts with two unnamed items, the sample x, and the model L2Fam.
...
further arguments
withLogScale
logical; shall a scale component (if existing and found with name scalename) be computed on log-scale and backtransformed afterwards? This avoids crossing 0.
..withCheck
logical: if TRUE, debugging info is issued.
withTimings
logical: if TRUE, separate (and aggregate) timings for the three steps evaluating the starting value, finding the starting influence curve, and evaluating the k-step estimator is issued.
withMDE
logical or NULL: Shall a minimum distance estimator be used as starting estimator---in addition to the function given in slot startPar of the L2 family? If NULL (default), the content of slot .withMDE in the L2 family is used instead to take this decision.
withEvalAsVar
logical or NULL: if TRUE (default), tells R to evaluate the asymptotic variance or if FALSE just to produces a call to do so. If withEvalAsVar is NULL (default), the content of slot .withEvalAsVar in the L2 family is used instead to take this decision.

Value

"kStepEstimate". In addition, it has an attribute "timings" where computation time is stored.

Details

Computes the optimally robust estimator for a given L2 differentiable parametric family. The computation uses a k-step construction with an appropriate initial estimate; cf. also kStepEstimator. Valid candidates are e.g. Kolmogorov(-Smirnov) or von Mises minimum distance estimators (default); cf. Rieder (1994) and Kohl (2005).

Before package version 0.9, this computation was done with the code of function roptest.old (with the same formals). From package version 0.9 on, this function uses the modularized function robest internally.

If the amount of gross errors (contamination) is known, it can be specified by eps. The radius of the corresponding infinitesimal contamination neighborhood is obtained by multiplying eps by the square root of the sample size.

If the amount of gross errors (contamination) is unknown, try to find a rough estimate for the amount of gross errors, such that it lies between eps.lower and eps.upper.

In case eps.lower is specified and eps.upper is missing, eps.upper is set to 0.5. In case eps.upper is specified and eps.lower is missing, eps.lower is set to 0.

If neither eps nor eps.lower and/or eps.upper is specified, eps.lower and eps.upper are set to 0 and 0.5, respectively.

If eps is missing, the radius-minimax estimator in sense of Rieder et al. (2001, 2008), respectively Section 2.2 of Kohl (2005) is returned.

Finite-sample and higher order results suggest that the asymptotically optimal procedure is to liberal. Using fsCor the radius can be modified - as a rule enlarged - to obtain a more conservative estimate. In case of normal location and scale there is function finiteSampleCorrection which returns a finite-sample corrected (enlarged) radius based on the results of large Monte-Carlo studies.

The default value of argument useLast is set by the global option kStepUseLast which by default is set to FALSE. In case of general models useLast remains unchanged during the computations. However, if slot CallL2Fam of IC generates an object of class "L2GroupParamFamily" the value of useLast is changed to TRUE. Explicitly setting useLast to TRUE should be done with care as in this situation the influence curve is re-computed using the value of the one-step estimate which may take quite a long time depending on the model.

If useLast is set to TRUE the computation of asvar, asbias and IC is based on the k-step estimate.

References

Kohl, M. (2005) Numerical Contributions to the Asymptotic Theory of Robustness. Bayreuth: Dissertation.

Kohl, M. and Ruckdeschel, P. (2010): R package distrMod: Object-Oriented Implementation of Probability Models. J. Statist. Softw. 35(10), 1--27 Kohl, M. and Ruckdeschel, P., and Rieder, H. (2010): Infinitesimally Robust Estimation in General Smoothly Parametrized Models. Stat. Methods Appl., 19, 333--354.

Rieder, H. (1994) Robust Asymptotic Statistics. New York: Springer.

Rieder, H., Kohl, M. and Ruckdeschel, P. (2008) The Costs of not Knowing the Radius. Statistical Methods and Applications 17(1) 13-40.

Rieder, H., Kohl, M. and Ruckdeschel, P. (2001) The Costs of not Knowing the Radius. Appeared as discussion paper Nr. 81. SFB 373 (Quantification and Simulation of Economic Processes), Humboldt University, Berlin; also available under www.uni-bayreuth.de/departments/math/org/mathe7/RIEDER/pubs/RR.pdf

See Also

roblox, L2ParamFamily-class UncondNeighborhood-class, RiskType-class

Examples

Run this code
## Don't run to reduce check time on CRAN
## Not run: 
# #############################
# ## 1. Binomial data
# #############################
# ## generate a sample of contaminated data
# ind <- rbinom(100, size=1, prob=0.05) 
# x <- rbinom(100, size=25, prob=(1-ind)*0.25 + ind*0.9)
# 
# ## ML-estimate
# MLest <- MLEstimator(x, BinomFamily(size = 25))
# estimate(MLest)
# confint(MLest)
# 
# ## compute optimally robust estimator (known contamination)
# robest1 <- roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3)
# robest1.0 <- roptest.old(x, BinomFamily(size = 25), eps = 0.05, steps = 3)
# identical(robest1,robest1.0)
# estimate(robest1)
# confint(robest1, method = symmetricBias())
# ## neglecting bias
# confint(robest1)
# plot(pIC(robest1))
# tmp <- qqplot(x, robest1, cex.pch=1.5, exp.cex2.pch = -.25,
#               exp.fadcol.pch = .55, jit.fac=.9)
# 
# ## compute optimally robust estimator (unknown contamination)
# robest2 <- roptest(x, BinomFamily(size = 25), eps.lower = 0, eps.upper = 0.2, steps = 3)
# estimate(robest2)
# confint(robest2, method = symmetricBias())
# plot(pIC(robest2))
# 
# ## total variation neighborhoods (known deviation)
# robest3 <- roptest(x, BinomFamily(size = 25), eps = 0.025, 
#                    neighbor = TotalVarNeighborhood(), steps = 3)
# estimate(robest3)
# confint(robest3, method = symmetricBias())
# plot(pIC(robest3))
# 
# ## total variation neighborhoods (unknown deviation)
# robest4 <- roptest(x, BinomFamily(size = 25), eps.lower = 0, eps.upper = 0.1, 
#                    neighbor = TotalVarNeighborhood(), steps = 3)
# estimate(robest4)
# confint(robest4, method = symmetricBias())
# plot(pIC(robest4))
# 
# #############################
# ## 2. Poisson data
# #############################
# ## Example: Rutherford-Geiger (1910); cf. Feller~(1968), Section VI.7 (a)
# x <- c(rep(0, 57), rep(1, 203), rep(2, 383), rep(3, 525), rep(4, 532), 
#        rep(5, 408), rep(6, 273), rep(7, 139), rep(8, 45), rep(9, 27), 
#        rep(10, 10), rep(11, 4), rep(12, 0), rep(13, 1), rep(14, 1))
# 
# ## ML-estimate
# MLest <- MLEstimator(x, PoisFamily())
# estimate(MLest)
# confint(MLest)
# 
# ## compute optimally robust estimator (unknown contamination)
# robest <- roptest(x, PoisFamily(), eps.upper = 0.1, steps = 3)
# estimate(robest)
# confint(robest, symmetricBias())
# 
# plot(pIC(robest))
# tmp <- qqplot(x, robest, cex.pch=1.5, exp.cex2.pch = -.25,
#               exp.fadcol.pch = .55, jit.fac=.9)
#  
# ## total variation neighborhoods (unknown deviation)
# robest1 <- roptest(x, PoisFamily(), eps.upper = 0.05, 
#                   neighbor = TotalVarNeighborhood(), steps = 3)
# estimate(robest1)
# confint(robest1, symmetricBias())
# plot(pIC(robest1))
# ## End(Not run)

#############################
## 3. Normal (Gaussian) location and scale
#############################
## 24 determinations of copper in wholemeal flour
library(MASS)
data(chem)
plot(chem, main = "copper in wholemeal flour", pch = 20)

## ML-estimate
MLest <- MLEstimator(chem, NormLocationScaleFamily())
estimate(MLest)
confint(MLest)

## Don't run to reduce check time on CRAN
## Not run: 
# ## compute optimally robust estimator (known contamination)
# ## takes some time -> you can use package RobLox for normal 
# ## location and scale which is optimized for speed
# robest <- roptest(chem, NormLocationScaleFamily(), eps = 0.05, steps = 3)
# estimate(robest)
# confint(robest, symmetricBias())
# plot(pIC(robest))
# ## plot of relative and absolute information; cf. Kohl (2005)
# infoPlot(pIC(robest))
# 
# tmp <- qqplot(chem, robest, cex.pch=1.5, exp.cex2.pch = -.25,
#               exp.fadcol.pch = .55, withLab = TRUE, which.Order=1:4,
#               exp.cex2.lbl = .12,exp.fadcol.lbl = .45,
#               nosym.pCI = TRUE, adj.lbl=c(1.7,.2),
#               exact.pCI = FALSE, log ="xy")
# 
# ## finite-sample correction
# if(require(RobLox)){
#     n <- length(chem)
#     r <- 0.05*sqrt(n)
#     r.fi <- finiteSampleCorrection(n = n, r = r)
#     fsCor <- r.fi/r
#     robest <- roptest(chem, NormLocationScaleFamily(), eps = 0.05, 
#                       fsCor = fsCor, steps = 3)
#     estimate(robest)
# }
# 
# ## compute optimally robust estimator (unknown contamination)
# ## takes some time -> use package RobLox!
# robest1 <- roptest(chem, NormLocationScaleFamily(), eps.lower = 0.05, 
#                    eps.upper = 0.1, steps = 3)
# estimate(robest1)
# confint(robest1, symmetricBias())
# plot(pIC(robest1))
# ## plot of relative and absolute information; cf. Kohl (2005)
# infoPlot(pIC(robest1))
# ## End(Not run)

Run the code above in your browser using DataLab