Bayesian Analysis of Incomplete Contingency Tables
bict.fit(priornum, missing1, missing2, maximal.mod, IP, eta.hat, ini.index, ini.beta, ini.sig, ini.y0, iters, save, name, null.move.prob, a, b, progress)
priornum |
A numeric scalar indicating which prior is to be used: 1 = |
missing1 |
A vector of the same length as the number of missing cell counts giving the row numbers of
the data.frame in |
missing2 |
A vector of the same length as the number of censored cell counts giving the row numbers of
the data.frame in |
maximal.mod |
An object of class |
IP |
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. |
ini.y0 |
A numeric vector giving the starting values of the missing and censored cell entries 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 |
name |
A prefix to the external files saved if the argument |
null.move.prob |
A scalar argument giving the probability of performing a null move, i.e. proposing a move to the current model. |
a |
The shape hyperparameter of the Sabanes-Bove & Held prior, see Overstall & King (2014). |
b |
The scale hyperparameter of the Sabanes-Bove & Held prior, see Overstall & King (2014). |
progress |
Logical argument. If |
The function will return a list with the following components.
BETA |
An |
MODEL |
A vector of length |
SIG |
A vector of length |
Y0 |
An |
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. |
This function will not typically be called by the user.
Antony M. Overstall A.M.Overstall@soton.ac.uk.
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/
data(spina) ## Load spina data. spina$z<-spina$y spina$z[is.na(spina$y)]<-0 ## Define a new variable in spina data.frame which is equal to y except where ## y is NA, in which case z=0. This is just so we can fit maximal model to the ## complete contingency table. maximal.mod<-glm(formula=z~(S1+S2+S3+eth)^2,data=spina,x=TRUE,y=TRUE, contrasts=list(S1="contr.sum",S2="contr.sum",S3="contr.sum", eth="contr.sum")) ## Fit maximal model to complete contingency table. curr.index<-formula2index(big.X=maximal.mod$x,formula=z~S1+S2+S3+eth,data=spina) ## Set up binary vector for independence model. 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[!is.na(spina$y),],prior="UIP", y=maximal.mod$y[!is.na(spina$y)],IP=IP) ## Find the posterior mode under the maximal model fitted to observed cell counts. 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<-bict.fit(priornum=1, missing1=(1:length(maximal.mod$y))[is.na(spina$y)], missing2=NULL,maximal.mod=maximal.mod, IP=IP, eta.hat=eta.hat, ini.index=curr.index, ini.beta=bmod[curr.index==1], ini.sig=1, ini.y0=c(500,200,20),iters=10, save=0, name=NULL, null.move.prob=0.5, a=0.001, b=0.001, progress = FALSE) ## Run for 10 iterations starting at model defined by curr.index. test1$MODEL ## Look at sampled model indicators. Should be: # [1] "7e00" "7e00" "7e00" "7e00" "7e00" "7e00" "7e00" "7e00" "7f00" "7f00" model2index(test1$MODEL,dig=15) ## Convert these to binary indicators of the log-linear parameters. ## Will get: # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7e00 1 1 1 1 1 1 0 0 0 0 0 #7f00 1 1 1 1 1 1 1 0 0 0 0 #7f00 1 1 1 1 1 1 1 0 0 0 0 # [,12] [,13] [,14] [,15] #7e00 0 0 0 0 #7e00 0 0 0 0 #7e00 0 0 0 0 #7e00 0 0 0 0 #7e00 0 0 0 0 #7e00 0 0 0 0 #7e00 0 0 0 0 #7e00 0 0 0 0 #7f00 0 0 0 0 #7f00 0 0 0 0
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.