bcct.fit: Bayesian Analysis of Complete Contingency Tables

Description

These functions are the workhorses behind bcct, bcctu, bcctsubset, and bcctsubsetu.

Usage

bcct.fit(priornum, maximal.mod, IP, eta.hat, ini.index, ini.beta, ini.sig, 
iters, save, name, null.move.prob, a, b, progress)
bcctsubset.fit(priornum, subset.index, maximal.mod, IP, eta.hat, ini.index, 
ini.beta, ini.sig, iters, save, name, null.move.prob, a, b, progress)

Arguments

priornum

A numeric scalar indicating which prior is to be used: 1 = "UIP", 2 = "SBH".

subset.index

A matrix where each row gives the index for each model in the subset of models under consideration.

maximal.mod

An object of class "glm" giving the fit of the maximal model.

A p by p matrix giving the inverse of the prior scale matrix for the maximal model.

eta.hat

A vector of length n (number of cells) giving the posterior mode of the linear predictor under the maximal model.

ini.index

A binary vector, of the same length as the number of log-linear parameters in the maximal model, indicating which parameters are present in the initial model.

ini.beta

A numeric vector giving the starting values of the log-linear parameters for the MCMC algorithm.

ini.sig

A numeric scalar giving the starting value of sigma^2 for the MCMC algorithm.

iters

The number of iterations of the MCMC algorithm to peform.

save

If positive, the function will save the MCMC output to external text files every save iterations. If zero , the function will not save the MCMC output to external files.

name

A prefix to the external files saved if the argument save is positive. If NULL, then the external files will have no prefix.

null.move.prob

A scalar argument giving the probability of performing a null move, i.e. proposing a move to the current model.

The shape hyperparameter of the Sabanes-Bove & Held prior, see Overstall & King (2014).

The scale hyperparameter of the Sabanes-Bove & Held prior, see Overstall & King (2014).

progress

Logical argument. If TRUE, then a progress bar will be displayed.

Value

The function will return a list with the following components:

BETA

An iters by p matrix containing the sampled values of the log-linear parameters, where p is the number of log-linear parameters in the maximal model. For elements of this matrix which correspond to a log-linear parameter which is not present for the current model a zero is returned.

MODEL

A vector of length iters giving the sampled model indicators in hexadecimal form.

SIG

A vector of length iters giving the sampled values for sigma^2 under the Sabanes-Bove & Held prior. If the unit information prior is used then the components of this vector will be one.

rj_acc

A binary vector of the same length as the number of reversible jump moves attempted. A 0 indicates that the proposal was rejected, and a 1 that the proposal was accepted.

mh_acc

A binary vector of the same length as the number of Metropolis-Hastings moves attempted. A 0 indicates that the proposal was rejected, and a 1 that the proposal was accepted.

References

Overstall, A.M. & King, R. (2014) conting: An R package for Bayesian analysis of complete and incomplete contingency tables. Journal of Statistical Software, 58 (7), 1--27. http://www.jstatsoft.org/v58/i07/

Examples

Run this code

# NOT RUN {
data(AOH)
## Load the AOH data.

maximal.mod<-glm(formula=y~(alc+hyp+obe)^3,data=AOH,x=TRUE,y=TRUE,
contrasts=list(alc="contr.sum",hyp="contr.sum",obe="contr.sum"))
## Set up the maximal model which in this case is the saturated 
## model.

curr.index<-formula2index(big.X=maximal.mod$x,formula=y~alc+hyp+obe+hyp:obe,data=AOH)
## Set up the binary vector for the model containing all main effects and the 
## hyp:obe interaction.

IP<-t(maximal.mod$x)%*%maximal.mod$x/length(maximal.mod$y)
IP[,1]<-0
IP[1,]<-0
## Set up the inverse scale matrix for the prior distribution under
## the maximal model. 

bmod<-beta_mode(X=maximal.mod$x,prior="UIP",y=maximal.mod$y,IP=IP)
## Find the posterior mode under the maximal model

eta.hat<-as.vector(maximal.mod$x%*%bmod)
## Find the posterior mode of the linear predictor
## under the maximal model.

set.seed(1)
## Set seed for reproducibility

test1<-bcct.fit(priornum=1, maximal.mod=maximal.mod, IP=IP, eta.hat=eta.hat, 
ini.index=curr.index, ini.beta=bmod[curr.index==1], ini.sig=1, iters=5, save=0, 
name=NULL,null.move.prob=0.5, a=0.001, b=0.001, progress=TRUE) 
## Run for 5 iterations starting at model defined by curr.index.

test1$MODEL
## Look at sampled model indicators. Should be:
## [1] "fe00c0" "fe0000" "fe0000" "fe0000" "fe0000"

model2index(test1$MODEL,dig=24)
## Convert these to binary indicators of the log-linear parameters.
## Will get: 

#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
#fe00c0    1    1    1    1    1    1    1    0    0     0     0     0     0
#fe0000    1    1    1    1    1    1    1    0    0     0     0     0     0
#fe0000    1    1    1    1    1    1    1    0    0     0     0     0     0
#fe0000    1    1    1    1    1    1    1    0    0     0     0     0     0
#fe0000    1    1    1    1    1    1    1    0    0     0     0     0     0
#       [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
#fe00c0     0     0     0     1     1     0     0     0     0     0     0
#fe0000     0     0     0     0     0     0     0     0     0     0     0
#fe0000     0     0     0     0     0     0     0     0     0     0     0
#fe0000     0     0     0     0     0     0     0     0     0     0     0
#fe0000     0     0     0     0     0     0     0     0     0     0     0

## See how the hyp:obe interactions in columns 17 and 18 gets dropped after 
## the 1st iteration. 

# }