Hi, folks,
This reminds me a question about car.normal and car.proper:
Since Holly puts intercept (alpha) in the model, so the sum of 'EPS' must be zero. I check the manual of WinBugs and car.normal already assumes this while car.proper does not. In the code, we can see that Holly set the EPS.mean=0 from 1 to N, but this should not be the case of car.normal. Of course, the sum of N zeroes are zero, but these N numbers could not all be zeroes necessarily. So, I guess, the estimates for these two methods should not be the same.
Any thoughts?
Wu
________________________________
From: Holly G <[log in to unmask]>
To: [log in to unmask]
Sent: Tuesday, July 31, 2012 1:29 PM
Subject: [BUGS] car.normal runs yet car.proper won’t compile
Hello list,
I can get the car.normal to run on my dataset (N= 8715),
but I run into an error at the compiling stage using car.proper. On one hand,
this doesn’t appear to be a memory issue, because the BUGS process seems to
only take a max of 28 MB before crashing. On the other hand, the car.proper
model runs on a subset of 1200 datapoints but produces the error for a subset
of 1500 datapoints and above. The prompt that arises at initiation of
compilation is:
model is syntactically correct
data loaded
Error in samplesSetLastChain(numChains) :
it is required to
have 1 <= last <= nchains
Calls: BRugsFit -> modelCompile -> samplesSetLastChain
If you do not think that this is a memory issue, have I
written a bug into either the car.proper vectors or the model itself?
Those scripts are below.
Thank you,
Holly
#####car.proper vectors
N<- nrow(FallData)
library(spdep)
junk=dnearneigh(as.matrix(FallData[,c("x","y")]),
d1=0, d2=4)
#anything less than d2=4 will produce empty neighbour
sets
#spatial weights output for car.proper:
NeighbSpatWeights <- nb2listw(junk)
InfoWeighted <- listw2WB(NeighbSpatWeights)
str(InfoWeighted)
#List of 3
# $ adj : int
[1:163488] 2 3 4 5 1 3 4 5 6 1 ...
# $ weights: num [1:163488] 0.25 0.25 0.25 0.25 0.2 ...
# $ num : int
[1:8715] 4 5 6 7 8 8 8 8 8 8 ...
sum(InfoWeighted$num==0)
#0
InfoWeighted$M <- 1/InfoWeighted$num
str(InfoWeighted$M)
# num [1:8715] 0.25 0.2 0.167 0.143 0.125 ...
#num =
InfoWeighted$num,
#C =
InfoWeighted$weights,
#M =
InfoWeighted$M,
#adj =
InfoWeighted$adj
##########ZIP GLMM CAR
model
{
#Priors
for(i in 1:7) { beta[i] ~ dnorm(0, 0.001) }
for(i in 1:NYdate)
{ a[i] ~ dnorm(0, tau.Ydate) } #RE
alpha ~ dnorm(0,
0.001) #intercept? for offset
psi ~ dunif(0, 1)
#ZI
tau.Ydate <- 1
/ (sigma.Ydate * sigma.Ydate) #RE
sigma.Ydate ~
dunif(0,10) #RE
#CAR precision
count prior
SP.tau<-1/(SP.sig*SP.sig)
SP.sig~dunif(0,
5)
#CAR spatial
correlation count prior
rho ~ dunif(0, 1)
#CAR count
EPS[1:N] ~
car.proper(EPS.mean[], C[], adj[], num[], M[], SP.tau, rho)
for (i in 1:N) {
W[i] ~
dbern(psi) #ZI
#Poisson part
(count)
CRATcnt[i] ~
dpois(eff.mu[i]) #note eff
eff.mu[i] <- W[i] * mu[i] #ZI
log(mu[i]) <- alpha
+ LogTransAreaSqM[i] +
beta[1]*STscale[i]
+
beta[2]*SalScale[i]
+
beta[3]*DistScale[i]
+
beta[4]*NbthScale[i]
+
beta[5]*TnDensScale[i]
+
beta[6]*DlDensScale[i]
+
beta[7]*WhlScale[i]
+
a[YdateNew[i]]
+
EPS[i]
EPS.mean[i] <-
0 #CAR
}
}
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask] To mail the BUGS list, mail to [log in to unmask] Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask] If this fails, mail [log in to unmask], NOT the whole list
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
