Description Usage Arguments Details Value Note Author(s) References Examples
This function implements one iteration of the Iterated Weight Least Square Metropolis Hastings Algorithm as proposed by Gamerman (1997) for generalised linear models as applied to log-linear models.
1 | iwls_mh(curr.y, curr.X, curr.beta, iprior.var)
|
curr.y |
A vector of length n giving the cell counts. |
curr.X |
An n by p design matrix for the current model, where p is the number of log-linear parameters. |
curr.beta |
A vector of length p giving the current log-linear parameters. |
iprior.var |
A p by p matrix giving the inverse of the prior variance matrix. |
For details of the original algorithm see Gamerman (1997). For its application to log-linear models see Overstall & King (2014), and the references therein.
The function will output a vector of length p giving the new values of the log-linear parameters.
This function will not typically be called by the user.
Antony M. Overstall A.M.Overstall@soton.ac.uk.
Gamerman, D. (1997) Sampling from the posterior distribution in generalised linear mixed models. Statistics and Computing, 7 (1), 57–68.
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/
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | set.seed(1)
## Set seed for reproducibility
data(AOH)
## Load AOH data
maximal.mod<-glm(y~alc+hyp+obe,family=poisson,x=TRUE,contrasts=list(alc="contr.sum",
hyp="contr.sum",obe="contr.sum"),data=AOH)
## Fit independence model to get a design matrix
IP<-t(maximal.mod$x)%*%maximal.mod$x/length(AOH$y)
IP[,1]<-0
IP[1,]<-0
## Set up inverse prior variance matrix under the UIP
## Let the current parameters be the MLE under the independence model
as.vector(coef(maximal.mod))
#[1] 2.89365105 -0.04594959 -0.07192507 0.08971628 -0.50545335 0.00818037
#[7] -0.01636074
## Update parameters using MH algorithm
iwls_mh(curr.y=AOH$y,curr.X=maximal.mod$x,curr.beta=coef(maximal.mod),iprior.var=IP)
## Will get:
#[1] 2.86468919 -0.04218623 -0.16376055 0.21656167 -0.49528676 -0.05026597
#[7] 0.02726671
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.