Computes Monte Carlo exact P-values for general log-linear models.

Share:

Description

This function computes Monte Carlo estimates of conditional P-values for goodness of fit tests for general log-linear models.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
mcexact(formula,
        data,
        stat = gof,
        dens = hyper,
        nosim = 10 ^ 3,
        method = "bab",
        savechain = FALSE,
        tdf = 3,
        maxiter = nosim,
        p = NULL,
        batchsize = NULL)

build.mcx.obj(formula,
        data,
        stat = gof,
        dens = hyper,
        nosim = 10 ^ 3,
        method = "bab",
        savechain = FALSE,
        tdf = 3,
        maxiter = nosim,
        p = NULL,
        batchsize = NULL)

Arguments

formula

Null model formula specified as in glm

data

Data frame

stat

The test statistic, a function of the form function(y, mu.hat) where y is the observed and mu.hat are the fitted values. Current default gof is a bivariate function of the deviance and the Pearson chi-squared.

dens

The target density on the log scale up to a constant of proportionallity. A function of the form function(y). Current default is (proportional to) the log of the generalized hypergeometric density.

nosim

Desired number of simulations.

method

Possibly two values, the importance sampling method of Booth and Butler, method = "bab" or the MCMC approach of Caffo and Booth method = "cab".

savechain

If TRUE saves the values of the chain.

tdf

A tuning parameter

maxiter

For method = "bab" number of iterations is different from the number of simulations. maxiter is a bound on the total number of iterations.

p

A tuning parameter for method = "cab".

batchsize

Required batchsizes for method = "cab".

Value

Returns a list of class either "bab" or "cab" depending on method. The list contains all of the inputs plus all required information to resume the simulation. Generic functions print and summary format the output while update can be used to resume simulations. mcexact is the front end while build.mcx.obj simply builds the basic object that mcexact applies to. simulate.conditional generates a matrix of simulated tables.

Author(s)

Brian Caffo

References

Booth and Butler (1999), "An importance sampling algorithm for exact conditional tests in log-linear models", Biometrika 86: 321-332.

Caffo and Booth (2001). "A Markov Chain Monte Carlo Algorithm for Approximating Exact Conditional Probabilities", The Journal of Computational and Graphical Statistics 10: 730-45.

http://www.biostat.jhsph.edu/~bcaffo/downloads.htm

See Also

fisher.test

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
#library(mcexact)
set.seed(1)

#importance sampling
data(residence.dat)
mcx <- mcexact(y ~ res.1985 + res.1980 + factor(sym.pair), data = residence.dat) 
summary(mcx)

#mcmc
data(pathologist.dat)
mcx <- mcexact(y ~ factor(A) + factor(B) + I(A * B),
               data = pathologist.dat,
               method = "cab",
               p = .5,
               nosim = 10 ^ 4,
               batchsize = 100)
summary(mcx)