Computes the WSDM statistic for each regression coefficient of a fitted VGLM.
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, ...)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.
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
A fitted vglm object.
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.
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 5).
Thus it takes twice as long 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.
Numeric, positive and of unit length;
multiplier for \(h\), i.e.,
relative to hdiff.
It is sometimes easier specifying a multiplier,
instead of the actual value, of hdiff.
So mux.hdiff = 2 will double
hdiff.
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.
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)).
Logical, use hdeff?
Some of the computation can take advantage
of this function, so this is optional.
Actually unimplemented currently.
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.
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.
Logical; return the derivatives?
If TRUE then a list is returned.
Logical; treat \(h\) as fixed?
If FALSE then hdiff is
multiplied by coef(object) so that
\(h\) is more relative than absolute.
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.
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.
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.
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 the (effective) max-WSDM.
Unused just now.
Thomas W. Yee.
Use with caution. This function 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.
This function, which is currently experimental
and very rough-and-ready,
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).
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. tools:::Rd_expr_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.
summaryvglm,
hdeffsev,
hdeff,
cops,
niters.
# 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)
summary(ppom.wine2, wsdm = TRUE)
max(wsdm(ppom.wine2, with1 = FALSE))
Run the code above in your browser using DataLab