The WSDM Function — wsdm • VGAM

Computes the WSDM statistic for each regression coefficient of a fitted VGLM.

Usage

wsdm(object, hdiff = 0.005, retry = TRUE, mux.hdiff = 1,
     maxderiv = 5, theta0 = 0, use.hdeff = FALSE,
     doffset = NULL, subset = NULL,
     derivs.out = FALSE, fixed.hdiff = TRUE,
     eps.wsdm = 0.15, Mux.div = 3, warn.retry = TRUE,
     with1 = TRUE, ...)

Arguments

object

A fitted vglm object.

hdiff

Numeric; the difference \(h\) used for the (original) finite difference approximations to the derivatives of the signed Wald statistics. Required to be positive and of unit length. For example, \(f'(x) = [f(x+h)-f(x)]/h + O(h)\) is used. If NAs are returned then increasing hdiff is often better than decreasing it. And hdiff can be changed via mux.hdiff. See also retry, eps.wsdm and Mux.div.

retry

Logical, compute with two other hdiff values to check that the finite-difference approximations were reasonably accurate? (For example, hdiff is multiplied and divided by 3). Thus it takes thrice times longer to return the answer if this is TRUE. And the original hdiff is used for the vector returned. If the absolute change is more than eps.wsdm then a warning is given.

mux.hdiff

Numeric, positive and of unit length; multiplier for \(h\), i.e., relative to hdiff. It is sometimes easier specifying this multiplier, instead of the actual value of hdiff. So mux.hdiff = 2 will double hdiff, etc.

maxderiv

Numeric, positive, integer-valued and of unit length. The highest order derivative to be computed. This means the highest value that the function will return is maxderiv - 1, i.e., it will be right-censored at that value.

theta0

Numeric; the hypothesized value. The default is appropriate for most symmetric binomialff links, and also for poissonff regression with the natural parameter. Recycled to length length(coef(object)).

use.hdeff

Logical, use hdeff? Some of the computation can take advantage of this function, so this is optional. Actually unimplemented currently.

doffset

Numeric; denominator offset. If inputted and not of sufficient length, then the remaining values are 1s. A value NULL is replaced by 1 unless the appropriate values are stored on object (but set 1 or some vector to override this). In particular, logistic regression has been internally set up so that doffset = c(2.399, 1.667, 2.178, 1.680, 2.2405) hence the WSDM function has continuous first derivatives everywhere except at the origin.

subset

Specify a subset of the regression coefficients? May be numeric, logical or character. Should work if coef(object)[subset] works. The default means subset <- TRUE.

derivs.out

Logical; return the derivatives? If TRUE then a list is returned.

fixed.hdiff

Logical; treat \(h\) as fixed? If FALSE then hdiff is multiplied by coef(object) so that \(h\) is more relative than absolute.

eps.wsdm

Numeric (positive and scalar). Tolerance for computing the WSDM value. Unless the three values are within this quantity of each other, a warning will be issued. It is usually not necessary to compute WSDM statistics very accurately because they are used as a diagnostic, hence this argument should not be too small.

Mux.div

Numeric (\(>1\) and scalar), for perturbing \(h\). If retry then hdiff is both multiplied and divided by Mux.div to give two separate step-sizes to compute the finite-difference approximations. Then a comparison involving eps.wsdm is performed to see if the answers are sufficiently similar.

warn.retry

logical; if retry, give a warning if all three estimates of the WSDM statistics are insufficiently similar? If FALSE then no call to warning will be given. However, see below on the attribute "seems.okay" attached to the answer.

with1

Logical. Include the intercepts? (This is a 1 in the formula language). Since WSDM statistics for the intercepts are less important, it is sometimes a good idea to set with1 = FALSE when computing a model's (effective) max-WSDM.

...

Unused just now.

Details

This function may be used as a diagnostic to see if any regression coefficients are alarmingly close to the parameter space boundary for the regularity conditions to be valid. A zero value denotes the centre of the parameter space (cops; COPS), which can be considered the heart of the interior of the parameter space where the Hauck–Donner effect (HDE) does not occur. A unit value occurs at \(w_1[0]\), the locations where the HDE starts taking place. As the WSDM statistic increases, the estimate is approaching the parameter space boundary, hence standard inference is breaking down and becoming more fraught with various dangers and inaccuracies.

The WSDM (pronounced wisdom) and the WSDM function are invariant to the sample size \(n\) for intercept-only models under random sampling. They are intended to be useful as a regression diagnostic tool for most VGLMs.

In summaryvglm, if the max-WSDM statistic, which is the maximum WSDM over all the regression coefficients bar the intercepts, is greater than 1.3, say, then the model should definitely not be used as it stands. One reason for the HDE is because a covariate is heavily skewed. If so, a suitable transformation can remedy the problem. The HDE may also be caused by complete separation in the covariate space.

Incidentally, another thing to check is the number of Fisher scoring iterations needed for convergence, e.g., any value greater than 10, say, should raise concern. Set trace = TRUE or look at niters(object) or summary(object).

Value

By default this function returns the WSDM statistics as a labelled numeric vector with nonnegative values, i.e., with names names(coef(object)). The attribute (see attr) "seems.okay" will always be attached to the answer and will be FALSE, TRUE, or NA or NULL if uncertain. If FALSE, retry by changing hdiff or mux.hdiff.

The following table is suggested for all link functions except for cauchit:

Range	Comments
[0.0, 0.5)	No HDE. Fine.
[0.5, 0.7)	Faint HDE. A borderline case, approaching something to be concerned with.
[0.7, 1.0)	Weak HDE. Should be of concern. Action is recommended but could possibly be optional.
[1.0, 1.3)	Moderate HDE. Action needed here and beyond.
[1.3, 2.0)	Strong HDE. Action definitely needed for this case.
[2.0, 3.0)	Extreme I HDE. Model should not be used or remedial action urgently needed.
[3.0, 4.0)	Extreme II HDE. Ditto.
[4.0, 5.0)	Extreme III HDE. Ditto.
...	...

This table supersedes the one given in Yee (2022), as this one is totally independent of \(n\) and has several advantages. Consequently, hdeffsev has been rewritten. No more than two or three decimal places should be used because the WSDM statistics are approximated by finite differences and are mainly used as a diagnostic. Probably, for most applications, large WSDM statistics for the intercepts should not be a problem, hence the max-WSDM excludes these in the summaryvglm.

To repeat, being mainly used as a diagnostic, WSDM values need not be computed or stated very accurately. It is suggested that 2 (or maybe 3) decimals places is fine.

If derivs.out = TRUE then higher-order derivatives are returned also within a list().

References

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 .

Yee, T. W. (2025). Mapping the parameter space with the WSDM function: A diagnostic for logistic regression and beyond. In preparation.

Author

Thomas W. Yee.

Warning

This function is currently a little bit experimental and rough-and-ready, hence use with some caution. It has been tested the most for logistic regression, and less so for other VGAM family functions. It will not work for all VGAM family functions.

Note

The results can be sensitive to hdiff so it is recommended that several \(h\) be tried, especially for regression coefficients that are near the parameter space boundary. Hence retry = TRUE is definitely recommended. This function could change in the short future because it is under active development and requires further fine-tuning.

Examples

if (FALSE) # Kosmidis (2014, Table 2), JRSSB 76: 169--196
ppom.wine2 <-
  vglm(cbind(bitter1, bitter2, bitter3, bitter4, bitter5) ~
       temp + contact,
       cumulative(reverse = TRUE,
                  parallel = TRUE ~ contact - 1),
       wine, trace = TRUE)
coef(ppom.wine2, matrix = TRUE)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'coef': object 'ppom.wine2' not found
summary(ppom.wine2, wsdm = TRUE)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'ppom.wine2' not found
max(wsdm(ppom.wine2, with1 = FALSE))
#> Error: object 'ppom.wine2' not found
 # \dontrun{}