[R-sig-Geo] BYM model in R-INLA - how to specify priors individually

James Rooney ROONEYJ4 at tcd.ie
Mon May 5 11:03:28 CEST 2014


Hi all,

I am attempting to implement the Besag-York-Mollie model in R-INLA. Specifically to recreate the WinBUGS model specified in ASDAR fig 11.10.

However I cannot figure out how to specify the spatial and random priors to be different as they are in the WinBUGS model fig11.10.

i.e. from the WinBUGS model fig 11.10 priors are specified as follows:
precu ~dgamma(0.001,0.001)
precv ~dgamma(0.1,0.1)

sigmau <-1/precu
sigmav <-1/precv

I need to do this in the INLA format. I have been playing with an example from the R-INLA webpage also using the North Carolina SIDS dataset - the code is listed below. This uses the default priors and initial values for the spatial and random parts of the INLA model. I would like to specify those priors differently for each component - equivalent to how they are specified in the above WinBugs model.

I'd be greatful for any assistance!
Kind regards,
James



#Load libraries
library(spdep)
library(maptools)
library(INLA)

#Load data from a shapefile included in the spdep package
nc.sids <- readShapePoly(system.file("etc/shapes/sids.shp", package="spdep")[1])

#Create adjacency matrix
nc.nb <- poly2nb(nc.sids)

#Compute expted number of cases
nc.sids$EXP<-nc.sids$BIR74*sum(nc.sids$SID74)/sum(nc.sids$BIR74)

#Compute proportion of non-white births
nc.sids$NWPROP<-nc.sids$NWBIR74/nc.sids$BIR74

#Convert the adjacency matrix into a file in the INLA format
nb2INLA("nc.adj", nc.nb)

#Create areas IDs to match the values in nc.adj
nc.sids$ID<-1:100

#Besag model with random spatial effect (i.e. BYM model)
m2<-inla(SID74~NWPROP+f(nc.sids$ID, model="besag", graph="nc.adj")+f(ID, model="iid"),
         family="poisson", E=nc.sids$EXP, data=as.data.frame(nc.sids),
         control.predictor=list(compute=TRUE))

#Get realtive risk estimates
nc.sids$RR2<-m2$summary.fitted.values[,1]

#Display relative risk estimates to show that both examples fit the same model
spplot(nc.sids, c( "RR2"))


More information about the R-sig-Geo mailing list