RJ_update: Reversible Jump Algorithm
In conting: Bayesian Analysis of Contingency Tables

Description Usage Arguments Details Value Note Author(s) References Examples

These functions implement one iteration of the orthogonal projection reversible jump algorithm for generalised linear models proposed by Forster et al (2012) applied to log-linear models with and without swap moves.

RJ_update_swap(prop.index, curr.index, curr.beta, eta.hat, curr.y, big.X, 
proposal.probs, i.prop.prior.var, i.curr.prior.var)

RJ_update(prop.index, curr.index, curr.beta, eta.hat, curr.y, big.X, 
proposal.probs, i.prop.prior.var, i.curr.prior.var)

`prop.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 proposed model.
`curr.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 current model.
`curr.beta`	A vector of length `sum(curr.index)` giving the log-linear parameters under the current model.
`eta.hat`	A vector of length n (number of cells) giving the posterior mode of the linear predictor under the maximal model.
`curr.y`	A vector of length n giving the cell counts.
`big.X`	The design matrix under the maximal model.
`proposal.probs`	A numeric vector of length 2. The first element gives the probability of proposing a move from the proposed model to the current model. The second element gives the probability of proposing a move from the current model to the proposed model.
`i.prop.prior.var`	A matrix giving the inverse of the prior variance matrix for the log-linear parameters under the proposed model.
`i.curr.prior.var`	A matrix giving the inverse of the prior variance matrix for the log-linear parameters under the current model.

For the original algorithm see Forster et al (2012). For details on its application to log-linear models see Overstall & King (2014), and the references therein.

RJ_update_swap performs birth/death and swap moves whereas RJ_update just performs birth/death moves.

The function will return a list with the following components:

`new.beta`	A vector giving the new log-linear parameters.
`new.index`	A binary vector indicating which log-linear parameters are present in the new model.

This function will not typically be called by the user.

Antony M. Overstall A.M.Overstall@soton.ac.uk.

Forster, J.J., Gill, R.C. & Overstall, A.M. (2012) Reversible jump methods for generalised linear models and generalised linear mixed models. Statistics and Computing, 22 (1), 107–120.

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/

set.seed(4)
## Set seed for reproducibility
data(AOH)
## Load data
maximal.mod<-glm(y~(alc+hyp+obe)^3,family=poisson,x=TRUE,contrasts=list(alc="contr.sum",
hyp="contr.sum",obe="contr.sum"),data=AOH)
## Fit maximal model to get a design matrix

IP<-t(maximal.mod$x)%*%maximal.mod$x/length(AOH$y)
IP[,1]<-0
IP[1,]<-0
## Calculate inverse prior scale matrix under maximal model. Under the UIP this
## is the inverse prior variance matrix. Under the SBH prior, we need to divide
## this matrix by the current value of SIG.

bmod<-beta_mode(X=maximal.mod$x,y=AOH$y,IP=IP)
## Find posterior mode under maximal model with UIP
eta.hat<-as.vector(maximal.mod$x%*%bmod)
## Find posterior mode of linear predictor.

curr.index<-formula2index(big.X=maximal.mod$x,formula=y~alc+hyp+obe+alc:hyp,data=AOH)
## Calculate index for current model including alc:hyp interaction
curr.index
## Print out current index
#[1] 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

pm<-prop_mod(curr.index=curr.index,data=AOH,maximal.mod=maximal.mod)
## Propose a model
p2<-(1-pm$null.move.prob)/pm$total.choices
p2
## Calculate probability of proposing proposed model from current model
#[1] 0.1666667

prop.index<-pm$new.index
prop.index
## Assign and print out proposal index
# [1] 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

dm<-prop_mod(curr.index=prop.index,data=AOH,maximal.mod=maximal.mod,null.move.prob=0)
p1<-(1-pm$null.move.prob)/dm$total.choices
p1
## Calculate probability of proposing current model from proposed model
#[1] 0.1666667

RJ_update(prop.index=prop.index,curr.index=curr.index,
curr.beta=coef(maximal.mod)[curr.index==1],eta.hat=eta.hat,curr.y=AOH$y,big.X=maximal.mod$x,
proposal.probs=c(p1,p2),
i.prop.prior.var=IP[prop.index==1,prop.index==1],
i.curr.prior.var=IP[curr.index==1,curr.index==1])

## Do one iteration of reversible jump algorithm. Will get:

#$new.beta
#(Intercept)        alc1        alc2        alc3        hyp1        obe1 
# 2.87128918 -0.07098006 -0.07221330  0.08748803 -0.51899802 -0.07855115 
#       obe2 
#-0.02474727 
#
#$new.index
# [1] 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0