BYM Model Problem

Mon, May 31, 2010 1:21 PM

Hi Peter,

You will probably need to define 'N' as the number of areas in your
problem. For this, you could use 

N<-nrow(nc at data)

This is the most likely cause for your problem, as I mentioned in my
previous e-mail


Also, please use these two lines to get the neighbour list in WB format:

nc.neig<-poly2nb(nc)
nc.nb <- nb2WB(nc.neig)


I believe that your code will work with these changes.

Best,

Virgilio

El lun, 31-05-2010 a las 16:19 -0400, Peter Larson escribi?:

Hello,

I can get the example in the book to work fine, but it will not work 
with my data. I have loaded a shp file in using the code below but don't 
know if I can attach it to the whole list or not.

What could be going wrong?

Thank you very much for your help,

Pete


nc <- readShapePoly("LynchingHotSpots.shp",  proj4string=llCRS)

statesWanted <- c("Alabama","Mississippi","Louisiana","Georgia","Kentucky",
"North Carolina","South Carolina","Florida","Arkansas","Tennessee")

nc <- nc[nc at data$STATE_NAME %in% statesWanted,]

#names(nc)#Variables in the dataset

nc$Observed<-nc$Lynchings
nc$Population<-nc$POP1990   #Population at risk
r<-sum(nc$Observed)/sum(nc$Population)
nc$Expected<-nc$Population*r

#Computed Standardised Mortality Ratio
nc$SMR<-nc$Observed/nc$Expected

#######BESAG

nc.nb <- nb2WB(nc)
nc$nwprop <- nc$Lynchings/nc$POP1990
d<-list(N=N, observed=nc$Observed, expected=nc$Expected,
   nonwhite=nc$nwprop,#log(nwprop/(1-nwprop)),
   adj=nc.nb$adj,  weights=nc.nb$weights, num=nc.nb$num)

dwoutcov<-list(N=N, observed=nc$Observed, expected=nc$Expected,
   adj=nc.nb$adj,  weights=nc.nb$weights, num=nc.nb$num)

inits<-list(u=rep(0,N), v=rep(0,N), alpha=0, beta=0, precu=.001, 
precv=.001)
#inits$v[d$num==0]<-NA


#### Winbugs ####

  bymmodelfile<-paste(getwd(), "/besag.txt", sep="")
  wdir<-paste(getwd(), "/BYM", sep="")
  if(!file.exists(wdir)){dir.create(wdir)}
  BugsDir <- "C:/Users/Pete/Downloads/WinBUGS14"
  MCMCres<- bugs(data=d, inits=list(inits),
  working.directory=wdir,
  parameters.to.save=c("theta", "alpha", "beta", "u", "v", "sigmau", 
"sigmav"),
  n.chains=1, n.iter=30000, n.burnin=20000, n.thin=10,
  model.file=bymmodelfile,
  bugs.directory=BugsDir,
  WINEPATH="/usr/bin/winepath")
      save(file="BYM.RData", list=c("d", "inits", "MCMCres") )
#Load the data obtained by running WinBUGS in Windows
nc$BYMmean<-MCMCres$mean$theta
#nc$BYMmedian<-MCMCres$median$theta
nc$BYMumean<-MCMCres$mean$u
#nc$BYMumedian<-MCMCres$median$u
nc$BYMvmean<-MCMCres$mean$v
#nc$BYMvmedian<-MCMCres$median$v
On 2010/05/31 14:20, Virgilio Gomez Rubio wrote:

On 2010/05/31 14:20, Virgilio Gomez Rubio wrote:

Dear Peter,

I am working through the examples in Applied Spatial Data Analysis with
R using my own data.

When I attempt to run the Besage-York-Mollie model in Chapter 11, I get
an index out of range error in WinBugs.

How should I go about finding where the problem is?

PLease, could you provide more information? In particular, are you
trying to reproduce the example in the book or adapting it to your own
data set? In that case, could you provide the code and data that you are
using (in a private e-mail, if you prefer)?

The error usually happens when there is some mismatch between the index
in the loop exceeds some vector size(i.e., looping over 20 areas when
there are just 15).

Best,

Virgilio

BYM Model Problem

Thread (8 messages)