huber2: Huber's least favourable distribution family function

Description

M-estimation of the two parameters of Huber's least favourable distribution. The one parameter case is also implemented.

Usage

huber1(llocation = "identity", k = 0.862, imethod = 1)
huber2(llocation = "identity", lscale = "loge",
       k = 0.862, imethod = 1, zero = 2)

Arguments

llocation, lscale

Link functions applied to the location and scale parameters. See Links for more choices.

Tuning constant. See rhuber for more information.

imethod, zero

See CommonVGAMffArguments for information. The default value of zero means the scale parameter is modelled as an intercept-only.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Details

Huber's least favourable distribution family function is popular for resistant/robust regression. The center of the distribution is normal and its tails are double exponential.

By default, the mean is the first linear/additive predictor (returned as the fitted values; this is the location parameter), and the log of the scale parameter is the second linear/additive predictor. The Fisher information matrix is diagonal; Fisher scoring is implemented.

The VGAM family function huber1() estimates only the location parameter. It assumes a scale parameter of unit value.

References

Huber, P. J. and Ronchetti, E. (2009) Robust Statistics, 2nd ed. New York: Wiley.

Examples

Run this code

set.seed(1231); NN <- 30; coef1 <- 1; coef2 <- 10
hdata <- data.frame(x2 = sort(runif(NN)))
hdata <- transform(hdata, y  = rhuber(NN, mu = coef1 + coef2 * x2))

hdata$x2[1] <- 0.0  # Add an outlier
hdata$y[1] <- 10  

fit.huber2 <- vglm(y ~ x2, huber2(imethod = 3), hdata, trace = TRUE)
fit.huber1 <- vglm(y ~ x2, huber1(imethod = 3), hdata, trace = TRUE)

coef(fit.huber2, matrix = TRUE)
summary(fit.huber2)


# Plot the results
plot(y ~ x2, hdata, col = "blue", las = 1)
lines(fitted(fit.huber2) ~ x2, hdata, col = "darkgreen", lwd = 2)

fit.lm <- lm(y ~ x2, hdata)  # Compare to a LM:
lines(fitted(fit.lm) ~ x2, hdata, col = "lavender", lwd = 3)

# Compare to truth:
lines(coef1 + coef2 * x2 ~ x2, hdata, col = "orange", lwd = 2, lty = "dashed")

legend("bottomright", legend = c("truth", "huber", "lm"),
       col = c("orange", "darkgreen", "lavender"),
       lty = c("dashed", "solid", "solid"), lwd = c(2, 2, 3))

Run the code above in your browser using DataLab

Get 50% off unlimited learning