conting (version 1.7)

RJ_update: Reversible Jump Algorithm

Description

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.

Usage

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)

Arguments

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.

Value

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.

Details

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.

References

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/

Examples

# NOT RUN {
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

# }