Hi Peter,
You will probably need to define 'N' as the number of areas in your
problem. For this, you could use
N<-nrow(n...@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...@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
> >
>
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
