[R-sig-Geo] error message using BRugs package - multivariate distribution must have more than one component

[Philip Lane]

Hi there,

I am trying to compare two hierarchical models  - one with a spatial term (using a car.normal function) and one without the term. I am trying to do this using the BRugs package in R.

I can get the non-spatial model to run in R using the BRugs package, but when
I try to add the CAR term to model the possible spatial association, I get an error message
multivariate distribution must have more than one component. 

I can get both models to work in WinBugs itself, I just want to call them from R, without
going through the CODA package to bring in the parameters I am interested in.

I am hoping to get a lead when the asdar book comes out next month, but I thought I might be fortunate to get some help from the list as well. I have attached the two models, data files, and R commands to the end of this email. 

One more thing - the error message has changed with the release of R 2.7.1 - when running the R code for the model with the CAR term, the model will not compile, with the error message
"logical expression contains too many operators" - I can work around this by using an older version of R, but I'd rather not.

Any help would be appreciated, as well as any tips or suggestions as to other approaches I could consider.

Cheers,


Philip.





#R code for non spatial lga analysis works
# lgas nonspatial using bugs command
# setwd("C:/Rwork/brugsfit/lga_sp")
model.file <- ("lga_ns3.bug")
#file.show(model.file) 
lgans <- read.table("lgas_ns_data2.txt", header=TRUE)
N <- nrow(lgans)
dlga <- lgans$dlga
drugs <- lgans$drugs
maliciouspd <- lgans$maliciouspd
personalviolence <- lgans$personalviolence
assault <- lgans$assault
theft <- lgans$theft

data <- list("N", "dlga", "theft", "assault", "drugs", "personalviolence", "maliciouspd")
inits <- function(){ 
list(mulga = rnorm(N, 0, 100), taulga = rnorm(1, 0, 100), 
mulgaerr = rnorm(N, 0, 100)) } 
parameters <- c("mulga", "theft", "beta0", "maliciouspd", "drugs", "assault", "personalviolence")
lgans1.sim <- bugs(data, inits, parameters, model.file,	
	n.chains=3, n.iter=9000, 
	program="OpenBUGS",working.directory="C:/Rwork/brugsfit/lga/")


# where the file 'lga_ns3.bug' is the following
model
	{
		
                for (i in 1 : N) {
		dlga[i] ~ dnorm(mulga[i], taulga)
		mulga[i] <- beta0 + theft[i] + drugs[i] + assault[i] + personalviolence[i] + maliciouspd[i] + mulgaerror[i]
		mulgaerr[i] ~ dunif(0,10)		
		mulgaerror[i] <-  pow(mulgaerr[i], -2) 
		}
	beta0 ~  dnorm(0, 0.00001)
	taulga ~ dunif(0,10)
	sigma <- sqrt(1 /  taulga) 
}

# and the file 'lgas_ns_data2.txt' is the following
dlga	assault	personalviolence	theft	drugs	maliciouspd
-1.308690894	-1.902542937	3.108024792	-1.908211249	-1.236057587	-2.029887038
0.192953031	-0.904543369	-0.168098513	-0.882832347	-0.771429408	-1.272027964
-0.742608421	-0.865276115	-0.977619891	-0.836548675	-0.485808832	-0.700348878
-0.830741022	0.989810005	0.406419865	1.189752862	2.548549645	1.020948188
0.948859566	0.288408929	-0.161014806	0.653297195	-0.694454275	0.31069764
0.511586279	0.871016672	-0.033402987	1.562855456	0.078507236	0.420099579
-1.257845163	-1.098227294	-0.759990924	-0.86449525	-0.562803883	-0.754887637
-1.356146909	-0.011072993	-0.256779458	-0.66885098	-0.217201319	-0.396837197
0.226850185	1.660304858	0.146284286	0.389135584	1.405311729	1.374473092
1.433588869	0.268211733	-0.193027654	0.522926518	-0.396927871	-0.170071478
0.331931362	0.841021781	-0.423071952	0.86942456	-0.052876539	0.608529601
0.426843394	-0.415227213	-0.49603507	-0.344025594	0.397721409	0.497197973
1.423419723	0.278115943	-0.191687687	0.317571921	-0.012530303	1.092114118


# but R code for spatial analysis falls over
setwd("C:/Rwork/brugsfit/lga_sp")
model.file <- ("lgatry2.bug")
#file.show(model.file) 
lgasp <- read.table("lgasp_data2.txt", header=TRUE)
N <- nrow(lgasp)
dlga <- lgasp$dlga
drugs <- lgasp$drugs
maliciouspd <- lgasp$maliciouspd
personalviolence <- lgasp$personalviolence
assault <- lgasp$assault
theft <- lgasp$theft
num=c(4,6,3,3,2,5,3,2,3,4,4,6,5)
adj=c(5,6,10,13,
	3,4,6,11,12,13,
	4,12,2,
	3,11,2,
	1,10,
	1,2,10,11,13,
	8,9,12,
	7,12,
	7,12,13,
	1,5,6,11,
	2,4,6,10,
	2,3,7,8,9,13,
	1,2,6,9,12)


data <- list("N", "dlga", "num", "adj", "theft", "assault", "drugs", "personalviolence", "maliciouspd") 
inits <- function(){ 
list(mulga = rnorm(N, 0, 100), taulga = rnorm(1, 0, 100), b=runif(N,0,10),
mulgaerror = rnorm(N, 0, 100), taulgacar=rnorm(1,0,100)) } 
parameters <- c("mulga", "theft", "maliciouspd", "drugs", "b", "assault", "personalviolence")
lgasp1.sim <- bugs(data, inits, parameters, model.file,debug=TRUE,	
	n.chains=3, n.iter=5000, program="OpenBUGS",
working.directory="C:/Rwork/brugsfit/lga_sp/")



# this is the file 'lgatry2.bug'
model
	{
		
                for (i in 1 : N) {
		dlga[i] ~ dnorm(mulga[i], taulga)
		mulga[i] <- beta0 + theft[i] + b[i] + drugs[i] + assault[i] + 
		personalviolence[i] + maliciouspd[i] + mulgaerror[i]
		mulgaerr[i] ~ dunif(0,10)		
		mulgaerror[i] <-  pow(mulgaerr[i], -2) 
		b[i] ~ car.normal(adj[], weights[], num[], taulgacar)
		}

for(k in 1:50)   {
		weights[k] <- 1}

	beta0 ~  dnorm(0, 0.00001)
	taulga ~ dunif(0,10)
	sigma <- sqrt(1 /  taulga) 
	sigmalgacar ~ dunif(0,10)
	taulgacar<- 1/(sigmalgacar*sigmalgacar)
}


#the file for the data 'lgasp_data2.txt'is the same as for the nonspatial ie
dlga	assault	personalviolence	theft	drugs	maliciouspd
-1.308690894	-1.902542937	3.108024792	-1.908211249	-1.236057587	-2.029887038
0.192953031	-0.904543369	-0.168098513	-0.882832347	-0.771429408	-1.272027964
-0.742608421	-0.865276115	-0.977619891	-0.836548675	-0.485808832	-0.700348878
-0.830741022	0.989810005	0.406419865	1.189752862	2.548549645	1.020948188
0.948859566	0.288408929	-0.161014806	0.653297195	-0.694454275	0.31069764
0.511586279	0.871016672	-0.033402987	1.562855456	0.078507236	0.420099579
-1.257845163	-1.098227294	-0.759990924	-0.86449525	-0.562803883	-0.754887637
-1.356146909	-0.011072993	-0.256779458	-0.66885098	-0.217201319	-0.396837197
0.226850185	1.660304858	0.146284286	0.389135584	1.405311729	1.374473092
1.433588869	0.268211733	-0.193027654	0.522926518	-0.396927871	-0.170071478
0.331931362	0.841021781	-0.423071952	0.86942456	-0.052876539	0.608529601
0.426843394	-0.415227213	-0.49603507	-0.344025594	0.397721409	0.497197973
1.423419723	0.278115943	-0.191687687	0.317571921	-0.012530303	1.092114118