Huber's Least Favourable Distribution Family Function

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

Usage

huber1(llocation = "identitylink", k = 0.862, imethod = 1)
huber2(llocation = "identitylink", lscale = "loglink",
       k = 0.862, imethod = 1, zero = "scale")

Arguments

llocation, lscale

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

k

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 intercept-only.

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.

Value

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

References

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

Author

T. W. Yee. Help was given by Arash Ardalan.

Note

Warning: actually, huber2() may be erroneous since the first derivative is not continuous when there are two parameters to estimate. huber1() is fine in this respect.

The response should be univariate.

Examples

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)
#> Iteration 1: loglikelihood = -67.24962
#> Iteration 2: loglikelihood = -61.145177
#> Iteration 3: loglikelihood = -59.955971
#> Iteration 4: loglikelihood = -59.836191
#> Iteration 5: loglikelihood = -59.82747
#> Iteration 6: loglikelihood = -59.826981
#> Iteration 7: loglikelihood = -59.826955
#> Iteration 8: loglikelihood = -59.826954
#> Iteration 9: loglikelihood = -59.826954
fit.huber1 <- vglm(y ~ x2, huber1(imethod = 3), hdata, trace = TRUE)
#> Iteration 1: loglikelihood = -71.478426
#> Iteration 2: loglikelihood = -61.066324
#> Iteration 3: loglikelihood = -60.041872
#> Iteration 4: loglikelihood = -59.957341
#> Iteration 5: loglikelihood = -59.955772
#> Iteration 6: loglikelihood = -59.955742
#> Iteration 7: loglikelihood = -59.955742

coef(fit.huber2, matrix = TRUE)
#>             location loglink(scale)
#> (Intercept) 1.465782     0.08756576
#> x2          8.624609     0.00000000
summary(fit.huber2)
#> 
#> Call:
#> vglm(formula = y ~ x2, family = huber2(imethod = 3), data = hdata, 
#>     trace = TRUE)
#> 
#> Coefficients: 
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1  1.46578    0.49572   2.957  0.00311 ** 
#> (Intercept):2  0.08757    0.16225   0.540  0.58941    
#> x2             8.62461    0.90399   9.541  < 2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Names of linear predictors: location, loglink(scale)
#> 
#> Log-likelihood: -59.827 on 57 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 9 
#> 


if (FALSE)  # Plot the results
plot(y ~ x2, data = hdata, col = "blue", las = 1)
lines(fitted(fit.huber2) ~ x2, data = hdata, col = "darkgreen", lwd = 2)
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet

fit.lm <- lm(y ~ x2, hdata)  # Compare to a LM:
lines(fitted(fit.lm) ~ x2, data = hdata, col = "lavender", lwd = 3)
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet

# Compare to truth:
lines(coef1 + coef2 * x2 ~ x2, data = hdata, col = "orange",
      lwd = 2, lty = "dashed")
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet

legend("bottomright", legend = c("truth", "huber", "lm"),
       col = c("orange", "darkgreen", "lavender"),
       lty = c("dashed", "solid", "solid"), lwd = c(2, 2, 3))  # \dontrun{}
#> Error in (function (s, units = "user", cex = NULL, font = NULL, vfont = NULL,     ...) {    if (!is.null(vfont))         vfont <- c(typeface = pmatch(vfont[1L], Hershey$typeface),             fontindex = pmatch(vfont[2L], Hershey$fontindex))    .External.graphics(C_strWidth, as.graphicsAnnot(s), pmatch(units,         c("user", "figure", "inches")), cex, font, vfont, ...)})(dots[[1L]][[1L]], cex = dots[[2L]][[1L]], font = dots[[3L]][[1L]],     units = "user"): plot.new has not been called yet