hdeffsev.RdComputes the severity of the Hauck-Donner effect for each regression coefficient of a fitted VGLM.
A fitted vglm object,
although not all VGAM family functions
will work, e.g.,
GAITD regression.
Alternatively, use wsdm.
Fed into wsdm.
Fed into wsdm.
Character vector of descriptors, plus
the last value for initialization.
Usually users should not assign anything to
this argument.
Used in conjunction with lookup.
Numeric, thresholds used for assigning
severity.table values
based on WSDM statistics.
This look-up table should be sorted
and the first element equal to 0.
Usually users should not assign anything to
this argument,
else be careful when using it.
Logical, transform WSDM before comparisons?
Applies to certain links only
(and if object was inputted), currently
loglink and
cauchitlink.
This is because it is possible to slightly
improve
the calibration between WSDM statistics
and the look-up table.
In particular,
a square root transformation shrinks
certain values towards unity.
The WSDM statistics can be inputted directly into the function here.
This function is intended to replace all
previous code for measuring HDE severity.
In particular,
hdeffsev0 and
hdeffsev2 are old and are
not recommended.
Details behind this function spring from
wsdm.
By default this function
(hdeffsev)
returns a labelled vector with
names names(coef(object)) and
elements selected from
severity.table.
Yee, T. W. (2022). On the Hauck-Donner effect in Wald tests: Detection, tipping points and parameter space characterization, Journal of the American Statistical Association, 117, 1763–1774. doi:10.1080/01621459.2021.1886936 .
For VGAM version 1.1-13,
hdeffsev() was renamed to hdeffsev0(),
hdeffsev2() to hdeffsev2() [no change],
and hdeffsev() is new and based on wsdm(vglmfit).
This function has not been tested
extensively and the thresholds may change
slightly in the future.
Improvements are intended.
The function was written specifically for
binomialff, but they should work
for almost all other family functions.
example(genpoisson0)
#>
#> gnpss0> gdata <- data.frame(x2 = runif(nn <- 500))
#>
#> gnpss0> gdata <- transform(gdata, y1 = rgenpois0(nn, theta = exp(2 + x2),
#> gnpss0+ logitlink(1, inverse = TRUE)))
#>
#> gnpss0> gfit0 <- vglm(y1 ~ x2, genpoisson0, data = gdata, trace = TRUE)
#> Iteration 1: loglikelihood = -2294.282
#> Iteration 2: loglikelihood = -2258.2725
#> Iteration 3: loglikelihood = -2255.347
#> Iteration 4: loglikelihood = -2255.3314
#> Iteration 5: loglikelihood = -2255.3314
#>
#> gnpss0> coef(gfit0, matrix = TRUE)
#> loglink(theta) logitlink(lambda)
#> (Intercept) 1.988143 1.021921
#> x2 1.041146 0.000000
#>
#> gnpss0> summary(gfit0)
#>
#> Call:
#> vglm(formula = y1 ~ x2, family = genpoisson0, data = gdata, trace = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 1.98814 0.05230 38.02 <2e-16 ***
#> (Intercept):2 1.02192 0.05200 19.65 <2e-16 ***
#> x2 1.04115 0.07552 13.79 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: loglink(theta), logitlink(lambda)
#>
#> Log-likelihood: -2255.331 on 997 degrees of freedom
#>
#> Number of Fisher scoring iterations: 5
#>
#> No Hauck-Donner effect found in any of the estimates
#>
summary(gfit0, wsdm = TRUE)
#>
#> Call:
#> vglm(formula = y1 ~ x2, family = genpoisson0, data = gdata, trace = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|) WSDM
#> (Intercept):1 1.98814 0.05230 38.02 <2e-16 0.57 ***
#> (Intercept):2 1.02192 0.05200 19.65 <2e-16 0.48 ***
#> x2 1.04115 0.07552 13.79 <2e-16 0.46 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: loglink(theta), logitlink(lambda)
#>
#> Log-likelihood: -2255.331 on 997 degrees of freedom
#>
#> Number of Fisher scoring iterations: 5
#>
#> No Hauck-Donner effect found in any of the estimates
#>
#> Max-WSDM (excluding intercepts): 0.463
#>
hdeffsev(gfit0)
#> (Intercept):1 (Intercept):2 x2
#> "Weak" "None" "Faint"