[R-sig-Geo] BYM Model Problem

[Virgilio Gomez Rubio]

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
> >    
>