huber.RdM-estimation of the two parameters of Huber's least favourable distribution. The one parameter case is also implemented.
huber1(llocation = "identitylink", k = 0.862, imethod = 1)
huber2(llocation = "identitylink", lscale = "loglink",
k = 0.862, imethod = 1, zero = "scale")Link functions applied to the location and scale parameters.
See Links for more choices.
Tuning constant.
See rhuber for more information.
See CommonVGAMffArguments for information.
The default value of zero means the scale parameter is
modelled as intercept-only.
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.
An object of class "vglmff"
(see vglmff-class).
The object is used by modelling functions
such as vglm,
and vgam.
Huber, P. J. and Ronchetti, E. (2009). Robust Statistics, 2nd ed. New York: Wiley.
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.
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
#>
#> No Hauck-Donner effect found in any of the estimates
#>
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