dbd (version 0.0-22)

gof: Goodness of fit test for db and beta binomial distributions.

Description

Either a chi-squared or a Monte Carlo test of goodness of fit of a db distribution.

Usage

gof(object, obsd, ...)
# S3 method for mleDb
gof(object,obsd,...,test=TRUE,MC=FALSE,seed=NULL,
                    nsim=99,maxit=1000,verb=FALSE)
# S3 method for mleBb
gof(object,obsd,...,test=TRUE,MC=FALSE,seed=NULL,
                    nsim=99,maxit=1000,verb=FALSE)

Arguments

object

An object of class "mleDb" or "mleBb" as returned by the function mleDb() or by mleBb().

obsd

The data to which object was fitted.

...

Not used.

test

Logical scalar. Should a hypothesis test be carried out? If test is FALSE then only the test statistic is returned. This argument is present so as to facilitate the calculations used in effecting a Monte Carlo test, by allowing gof() to recursively call itself.

MC

Logical scalar. Should a Monte Carlo test be used rather than a chi squared test?

seed

Integer scalar. The seed for the random number generator used when MC is TRUE. If not supplied, seed is created by sampling one integer from 1:1e5. This argument is ignored if MC is FALSE.

nsim

The number of simulated replicates on which the Monte Carlo test is to be based. Ignored if MC is FALSE.

maxit

Integer scalar. The maximum number of iterations to be undertaken by optim() when fitting models to the simulated data. Ignored if MC is FALSE.

verb

Logical scalar. Should rudimentary “progress reports” be issued during the course of the simulations invoked by the Monte Carlo test procedure? Ignored if MC is FALSE.

Value

A list with components

stat

The test statistic.

pval

The p-value of the test.

degFree

The degrees of freedom of the chi squared test.

The last component is present only if a chi squared test (rather than a Monte Carlo test) is used.

If a chi squared test is used and turns out to be invalid, then the returned value has an attribute "expVals", consisting of the (problematic) expected values.

If a Monte Carlo test is used the returned value has an attribute "seed" which is equal to the seed argument or to the random value selected to replace it if the seed argument was not supplied.

Notes

The Monte Carlo \(p\)-value is calculated as (m+1)/(nsim+1) where m is the number of simulated statistics which greater than or equal to the observed statistic (computed from the “real” data.

The smallest that the Monte Carlo \(p\)-value can be is 1/(nsim + 1), e.g. 0.01 when nsim is 99. For “finer distinctions” you must use larger values of nsim, such as 999 or 9999.

The \(p\)-value is random; if you repeat the test (with the same data) you may well get a different \(p\)-value. Resist the temptation to repeat the test until you get a \(p\)-value that you like!!! This invalidates your inference!

Remark on the Examples

In the Examples, db and beta binomial distributions are fitted to the Parsonnet scores from the cardiacsurgery data set which comes from the spcadjust package. It is not completely clear what the value of ntop (db distribution) or size (beta binomial distribution) should be. The data are not actually counts, and in particular they are not counts of successes out of a given number (“size”) of trials. In the event I chose to use the value 71, the maximium value of the Parsonnet scores, for the value of both ntop and size. This was the value chosen for use as size by Wittenberg (2021) when he fitted a beta binomial distribution to these data.

Details

The function gof() is a generic function with two methods, gof.mleDb() and gof.mleBb().

The test statistic is calculated as $$\sum((O-E)^2/E)$$ where \(O\) means “observed” and \(E\) means “expected”. If the mean of \(E\) is less than 5 or if any of the entries of \(E\) is less than 1, then the chi squared test is invalid and a warning to this effect is issued. In this case the expected values are returned as an attribute of the value returned by gof(). The foregoing applies of course only if a chi squared test (as opposed to a Monte Carlo test) is being used.

The degrees of freedom for the chi squared test are length(E) - 3. The value 3 is equal to 2 (for the number of parameters estimated) plus 1 (for the costraint that the probabilities of the values sum to 1).

If it were actually true that, under the null hypothesis, the observed test statistic and those calculated from simulated data are exchangeable, the Monte Carlo test would be exact. However the real data are distributed as \(f(x,\theta)\) whereas the simulated data are distributed as \(f(x,\hat{\theta})\) where \(\hat{\theta}\) is the estimate of \(\theta\) based on the observed data. Consequently the observed test statistic and simulated test statistics are “not quite” exchangeable. Nevertheless it appears that in practice the Monte Carlo test is very close to being exact.

The meaning of “exact” here is that if the null hypothesis is true then, over the set of instances of collecting the data and simulating the required replicates, the \(p\)-value is uniformly distributed on the set \(\{1/N, 2/N, \ldots, (N-1)/N, 1\}\) where \(N\) is equal to nsim.

References

Philipp Wittenberg (2021). Modeling the patient mix for risk-adjusted CUSUM charts. To appear in Statistical Methods in Medical Research.

Axel Gandy and Jan Terje Kvaloy (2013). Guaranteed conditional performance of control charts via bootstrap methods. Scandinavian Journal of Statistics 40, pp. 647--668. (Reference for spcadjust package.)

See Also

mleDb()

Examples

# NOT RUN {
X <- hmm.discnp::Downloads
f <- mleDb(X,15,TRUE)
tst1 <- gof(f,X) # Gives warning that the chi squared test is invalid.
tst2 <- gof(f,X,MC=TRUE,seed=42)
# The p-value is 0.03 so we reject the adequacy of the fit at the 0.05
# significance level.  Note that the p-value that we get, when the
# random number generator seed is set equal to 42, is very similar in
# value to the p-value (0.0347) from the "invalid" chi squared test.
#
# }
# NOT RUN {
 # Takes too long.
if(requireNamespace("spcadjust")) {
    data("cardiacsurgery", package = "spcadjust")
    xxx  <- cardiacsurgery$Parsonnet
    fit1 <- mleDb(xxx,ntop=71,zeta=TRUE)
    g1   <- gof(fit1,obsd=xxx,MC=TRUE,verb=TRUE,seed=42)
    fit2 <- mleBb(xxx,size=71)
    g2   <- gof(fit2,obsd=xxx,MC=TRUE,verb=TRUE,seed=17)
}
# }