iwls_mh: Iterated Weighted Least Square Metropolis Hastings Algorithm

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

View source: R/iwls_mh.R

Description

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.

Usage

1
iwls_mh(curr.y, curr.X, curr.beta, iprior.var)

Arguments

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.

Details

For details of the original algorithm see Gamerman (1997). For its application to log-linear models see Overstall & King (2014), and the references therein.

Value

The function will output a vector of length p giving the new values of the log-linear parameters.

Note

This function will not typically be called by the user.

Author(s)

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

References

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/

Examples

 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

conting documentation built on May 1, 2019, 8:47 p.m.