[R-sig-Geo] BYM Model Problem

[Peter Larson]

Virgilio,

It appears to be working now. Thank you very much for your help.

Pete



On 2010/05/31 16:21, Virgilio Gomez Rubio wrote:
> 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
>>>
>>>        
>>      
>