These functions are used to detemine which models we can propose moves to, given the current model in the MCMC algorithm, and the principle of marginality.
add_term(curr.index, data, maximal.mod)
drop_term(curr.index, data, maximal.mod)
prop_mod(curr.index,data,maximal.mod,null.move.prob=0.5)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.
An object of class "data.frame" containing the variables in the model.
An object of class "glm" giving the fit of the maximal model.
An optional scalar argument giving the probability of performing a null move, i.e. proposing a move to the current model. The default value is 0.5.
The functions add_term and drop_term will output a character
vector containing the names of terms that can be dropped.
The function prop_mod will return a list with the following components.
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.
A character string which will be one of c("null","drop","add")
depending on the type of move proposed.
If type is not "null", then total.choices will be scalar giving the total number of non-null moves available. If type is equal to "null", then total.choices will be 0.
A scalar giving the probability of a null move.
In the reversible jump algorithm we propose a move to a model given the current model. The function prop_mod
implements a scheme whereby only local proposals are made, i.e. either a term is added or dropped. These
types of move are called birth and death moves, respectively, by Forster et al (2012).
When a term is either added or dropped, we preserve the principle of marginality, e.g. we can only propose to add a three-way interaction if all the possible two-way interactions between the three factors are included in the present model.
The functions add_term and drop_term determine which terms can be added or dropped whilst preserving the
principle of marginality.
The function prop_mod will call add_term and drop_term thus determining which terms can be added or
dropped. With probability null.move.prob it will choose to remain in the current model; otherwise it
will choose one of the possible terms to add or drop.
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.
# 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.
add_term(curr.index=curr.index,data=AOH,maximal.mod=maximal.mod)
## See what terms we can add - will get:
#[1] "alc:hyp" "alc:obe"
drop_term(curr.index=curr.index,data=AOH,maximal.mod=maximal.mod)
## See what terms we can drop - will get:
#[1] "hyp:obe"
set.seed(4)
## Set the seed for reproducibility.
prop_mod(curr.index=curr.index,data=AOH,maximal.mod=maximal.mod)
## Propose a model. Will be a drop move, proposing the independence model by
## dropping the hyp:obe interaction. The total.choices object is 3, i.e. one
## drop move and two add moves. Specifically:
#$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
#
#$type
#[1] "drop"
#
#$total.choices
#[1] 3
#
#$null.move.prob
#[1] 0.5
# }
Run the code above in your browser using DataLab