Bugs Group,
Thanks for all the responses to my question about using the car.normal
function to have a random area effect that is simultaneously
correlated in space and time. I will summarize the responses below.
I have been working on this for some time and couldn’t quite figure it
out. Many people pointed out that I should look at the book “Disease
Mapping with WinBUGS and MLwiN.” I have looked at this book but not
for a long time. Many people pointed out the model of Bernardinelli
(1995) and of Waller (1997) discussed in this book (starting on page
128). Bernardinelli’s model is a “random slope” type of model where
the time trend is linear but areaspecific with some spatial
correlation determined by the inputs to car.normal. Waller’s model is
a nested type model that has different areaspecific spatially
correlated effects for each separate year in the time series, but no
temporal correlation. This allows more freedom for the time trend as
it does not force a linear trend with time (pointed out to me by
Rodney Sparapani). Anyhow, it turns out that I had been fitting a
model similar to Waller’s model but realized the correlation across
time was important. It also turns out
that I had fit a “random slope” type model like Bernardinelli long
ago but had big stability problems with the random effect variance and
was getting very unrealistic (I thought) maps. Hence, my plea for
help for the Bugs group.
I had been trying to figure out an indexing scheme with the
adjacency matrix input to the car.normal but couldn't quite get my
head around it. Rodney Sparapani pointed out the trick that made
solving this problem quite easy (see below). I ended up writing an R
function that takes the spatial vectors num[] and adj[] and the number
of time points T, then outputs the temp.num[] and temp.adj[] vector
that refer to the number of neighbors and the index which refers to
the adjacent neighbors in spacetime. The way that I have it now, any
given position in spacetime is correlated to all the neighbors in the
3x3x3 spacetime cube with equal weight. This seems like a reasonable
and simple first step for my situation, but I’m not convinced that
equal weights for the whole cube are the best a priori expectation.
I’m dealing with fairly sparse data in a smallscale disease mapping
problem, so using the whole cube as neighbors (as I see it) helps deal
with the sparseness of the da
ta.
I have included the R code for the function that I describe.
I’m new to R programming and there are probably more stylish ways to
do the programming. I would be interested in people’s comments on the
program, on its use, and on any problems that are discovered. I have
tested this program on small (2x2 and 3x3) spatial maps, which I found
to work fine. The function also works well and fast for my space
time: 3600 locations X 4 years.
When I use the spacetime correlated random effect to make my
map in WinBUGS, the iterations go fairly fast and convergence for the
randomeffect variance and fixed effect parameters are good (at least
as good as the Waller model). The map produced is also much more
reasonable given my prior beliefs based on past explorations of the
data and theoretical expectations. I should also mention that I am
not fitting the same Poisson disease mapping models as Bernardinelli
and Waller. Instead of the Poisson distribution as the likelihood, I
am using the Bernoulli, and instead of the loglink function, I am
using the complimentaryloglog link function. (Just for those who
are interested.)
I would appreciate any general comments on the use of this
Spacetime correlation effect for disease mapping as I have described
it. I know I haven’t been very thorough, but the details will have to
wait until the manuscript is finished. Thanks for all the help!
Here is the R function:
######################################################################
spacetime < function(num,adj,T)
{
if(sum(num)!=length(adj))
return('Error: num[] vector does not match adj[]
vector');
sumnum < cumsum(num);
N < length(num);
temp_adj < c();
temp_num < c();
cnter < 1;
for(i in 1:T)
{
for(j in 1:N)
{ ######Make temp_num vector######
if(i==1i==T)
{
temp_num[cnter] < num[j]*2+1;
}
else
{
temp_num[cnter] < num[j]*3+2;
}
cnter < cnter+1;
#####Make temp_adj vector####
if(i==1)
{
if(j==1)from < 1;
if(j>1)from < sumnum[j1]+1;
to < sumnum[j];
part < c(adj[from:to],adj[from:to]
+N,j+N);
}
else if(i==T)
{
if(j==1)from < 1;
if(j>1)from < sumnum[j1]+1;
to < sumnum[j];
part < c(adj[from:to]+N*(T1),
adj[from:to]+N*(T2),j+N*(T2));
}
else
{
if(j==1)from < 1;
if(j>1)from < sumnum[j1]+1;
to < sumnum[j];
part < c(adj[from:to]+N*(i1),
adj[from:to]+N*(i2),
adj[from:to]+N*i,j+N*(i2),j+N*i);
}
temp_adj < c(temp_adj,part);
}
}
result < list(temp_num,temp_adj);
return(result);
#####################################################################
Here is the original request:
Original Message
From: (The BUGS software mailing list) [mailto:[log in to unmask]] On
Behalf
Of ERIK E OSNAS
Sent: Tuesday, July 18, 2006 11:20 AM
To: [log in to unmask]
Subject: spacetime correlation in WinBUGS
Bugs Group,
I am interested in using the car.normal function in WinBUGS for
mapping an emerging infectious disease in a deer population. I have
successfully used this function to implement the Besag, York et al.
model of spatial correlation across adjacent locations. However, I
have several years of data and am interested in modeling the dynamics
across years as well as spatial locations. I would like to be able to
implement both space and time correlations where the randomeffect for
any area is correlated to its adjoining spatial neighbors and to its
temporal neighbors (including itself) both one time step before and
after. Is it possible to use car.normal in WinBUGS to do this? I
have looked at many examples but have only found examples for temporal
or spatial correlation but not simultaneous space and time
correlations. If this is not possible in WinBUGS, can someone please
point me to a resource where this could be modeled using MCMC.
Thanks,
Erik
Erik Osnas
Assistant Scientist
Department of Wildlife Ecology
University of Wisconsin
218 Russell Labs
1630 Linden Dr.
Madison, WI 53706
6082621984
Here is the email correspondence (chronological order):
From "Best, Nicky G" <[log in to unmask]>
Sent Monday, July 24, 2006 7:51 am
To ERIK E OSNAS <[log in to unmask]>
Cc
Bcc
Subject RE: spacetime correlation in WinBUGS
Attachments winmail.dat
5K
Dear Erik
We have recently published an article on spacetime modelling of two
related diseases (all done in WinBUGS) which might help you. See
Richardson S; Abellan JJ; Best N. (2006) Bayesian spatiotemporal
analysis of joint patterns of male and female lung cancer risks in
Yorkshire (UK) STAT METHODS MED RES. 15: 385407
Nicky
From Sujit Ghosh <[log in to unmask]>
Sent Monday, July 24, 2006 10:47 am
To ERIK E OSNAS <[log in to unmask]>
Cc Lee Hyeyoung <[log in to unmask]>
Bcc
Subject Re: spacetime correlation in WinBUGS
Hi Erik,
We have a paper on using spatial DLM to model spatiotemporal data and
we have implemented it in winbugs using a reparamerization approach.
You
may find more details in our tech report #2587 available at the
following site: http://www.stat.ncsu.edu/library/mimeo.html
Hope it helps. If you want to collaborate and need hep with winbugs
programming, let me know.
SG
From [log in to unmask]
Sent Monday, July 24, 2006 10:49 am
To ERIK E OSNAS <[log in to unmask]>
Cc
Bcc
Subject Re: [BUGS] spacetime correlation in WinBUGS
Hi Erik,
Hope all's well in the CWD modeling world. I don't have a good answer
for you, however I can offer a bit of advice on incorporating space
and time simultaneously. I found (Farnsworth et al. 2006, Eco. Apps 16
(3)) that incorporating two random effects in the form of spatial and
nonspatial heterogeneity was, not surprisingly, sensitive to the
choice of priors on the precisions and that the choice of gridcell
size had a large impact on subsequent inference, which I imagine may
be the same for your work, unless of course you're only interested in
a single spatial scale (e.g., management units). I imagine that
incorporating temporal correlation into a spatial model will be quite
data hungry, although for a large enough unit size and given the
amount of data you're working with you might be able to partition out
signal from noise. I'm sure you're well aware of all this. Anyway,
if you can post to the group with useful replies I, and others I'm
sure, will appreciate it. One
person who comes to mind is Chris Wikle's work on house finch
invasion across the U.S. If I remember correctly he used integro
difference models in a hierarchical framework to look at space time
invasion dynamics. Beyond that, you might try looking at Alan
Gelfand's recent work on modeling plant diversity in South Africa,
although I can't recall if he considers both space and time
simultaneously. You're squarely in an area that is at the cutting
edge of research, which of course makes your efforts challenging, at
least by my standards. I think some folks have used RJMCMC
cheers,
matt
Matthew Farnsworth, Ph.D.
Research Biologist
USDA/APHIS/WS
National Wildlife Research Center
4101 LaPorte Ave.
Fort Collins, CO. 80521
(970) 2666129
(970) 2666138 (fax)
From Rudy Banerjee <[log in to unmask]>
Sent Monday, July 24, 2006 5:03 pm
To [log in to unmask]
Cc
Bcc
Subject Re: [BUGS] spacetime correlation in WinBUGS
Hi Eric,
Have you tried the Bernardinelli et al (1995) model. You can find
examples of space time models based on car.normal autocorr (see below).
'Disease Mapping with WinBUGS and MLWin' by Lawson et al. 2003 (pages
128133) shows some examples. I have implemented spacetime models with
my data successfully. I will be happy to share the code...
Best of luck,
Rudy
model>
model
{
for (i in 1:m)
{
for (k in 1:T)
{
# Poisson likelihood for observed counts; y is the observed
value of the dependent variable and e is the expected value
y[i,k]~dpois(mu[i,k])
log(mu[i,k])<log(e[i,k])+alpha+ alpha1*variable1+...(your other
variables)..+
+u[i]+v[i]+beta*t[k]+delta[i]*t[k] #u is the structured
spatial heterogeneity or spatial autocorr and v is the unstrsuctured or
random effects, beta is the time random component and delta is the
space
time interaction effect.
# Relative Risk in each area and period of time
theta[i,k]<exp(alpha
+ alpha1*variable1[i,k] +...
+u[i]+v[i]+beta*t[k]+delta[i]*t[k])
}
TT[i]<exp(beta+delta[i])
}
# CAR prior distribution for spatial correlated heterogeneity
u[1:m]~car.normal(adj[],weights[],num[],tau.u)
delta[1:m]~car.normal(adj[],weights[],num[],tau.delta)
# Prior distributions for the Uncorrelated Heterogeneity
for(i in 1:m)
{
v[i]~dnorm(0,tau.v)
}
# Weights
for(k in 1:sumNumNeigh)
{
weights[k]<1
}
alpha~dflat()
alpha1~dflat()
.....
# Hyperprior distributions on inverse variance parameter of random
effects
beta~dnorm(0,1.0E5)
tau.v~dgamma(0.5,0.0005)
tau.u~dgamma(0.5,0.0005)
tau.delta~dgamma(0.5,0.0005)
sigma.v < 1/sqrt(tau.v)
sigma.u < 1/sqrt(tau.u)
sigma.delta < 1/sqrt(tau.delta)
}
NOTE: theta[i,k] the Relative Risk and TT will give the time period
RR.
DATA>
list(y=structure(.Data=c(),.Dim=c(#of areal units,# of time units)),
m = 304, T=5,t=c(1,2,3,4,5),
sumNumNeigh = 1860,
num = c(),
adj = c())
INITS>
list(alpha=0,
alpha1=0,
...,
beta=0,
tau.v=1,
tau.u=1,
tau.delta=1,
u=c(), # zeores
delta=c(), # zeores
v=c()) # zeores
From Bridget Freisthler <[log in to unmask]>
Sent Monday, July 24, 2006 6:12 pm
To [log in to unmask]
Cc
Bcc
Subject spacetime correlation in WinBUGS
In case no one else has responded, I have seen two different ways of
modeling spacetime variation. Details can be found in:
Bernardinelli, L., Clayton, D., Pascutto, C., Montomoli, C.,
Ghislandi, M., & Songini, M. (1995). Bayesiananalysis of spacetime
variation in disease risk. Statistics in Medicine, 14(2122), 243343.
Waller, L.A., Carlin, B.P., Xia, H., and Gelfand, A.
(1997). "Hierarchical spatiotemporal mapping of disease rates".
Journal of the American Statistical Association 92, 607617.
Best,
Bridget
Bridget Freisthler, Ph.D., Assistant Professor
Department of Social Welfare, UCLA School of Public Affairs
3250 Public Policy Building, Box 951656, Los Angeles, CA 90095
(Phone) 310/2061602 (Fax) 310/2067564
(email) [log in to unmask]
From Rodney Sparapani <[log in to unmask]>
Sent Tuesday, July 25, 2006 4:12 pm
To ERIK E OSNAS <[log in to unmask]>
Cc
Bcc
Subject Re: spacetime correlation
Hi Erik:
This sounds like a very interesting project that I'd love to hear more
about.
I think this can be done, but you'd need to do some bookkeeping.
Suppose
that we have N areas and T times, letting M=TxN. If N=4 and T=2, then
it
should be simple enough to display if adjacency is only nearest
neighbors and
the weights are all 1:
S[8] ~ car.normal(adj[], weight[], num[], tau);
Areas
1 2
3 4
Index Area Time adj[1] adj[2] adj[3] num[]
1 1 1 2 3
0 2
2 2 1 1 4
0 2
3 3 1 1 4
0 2
4 4 1 2 3
0 2
5 1 2 1 6
7 3
6 2 2 2 5
8 3
7 3 2 3 5
8 3
8 4 2 4 6
7 3
Rodney
From Rodney Sparapani <[log in to unmask]>
Sent Wednesday, July 26, 2006 12:15 pm
To ERIK E OSNAS <[log in to unmask]>
Cc
Bcc
Subject Re: spacetime correlation
ERIK E OSNAS wrote:
> Thanks for your response. I have been trying to think of an
indexing
> scheme like what you propose, but I haven't been able to get my head
> around it.
>
> I do think that there will have to be N x T "locations" in the
> adjacency matrix and I have been able to code that in for the num[]
> vector. The adj[] vector is more problematic. Your example helps.
I
> will have to think more about it to see if it solves my problem.
> The "trick" you propose is to refer to an adjacent neighbor in the
new
> NT index rather than to the original index of N areas. That makes
> sense to me. If I can write reasonable code for that problem and it
> works, I'll let you know.
>
> Do you think that explicitly coding the spatial and temporal
neighbors
> into the car.normal function is different than using a model of
> Bernardinelli (1995)? The model is:
>
> x(ij) = a + u(i) + B t(j) + v(i) t(j)
>
> i is an index for spatial location, j is an index for time t, a is
the
> average effect (intercept), u(i) is the car.normal random effect (+
an
> unstructured effect as in the BYM model that I left out for
> simplicity), B is a regression effect for time, and v(i) is a area
> specific car.normal "random slope" for time.
>
> Is this different than the indexing model that I asked my original
> question about? I think so, but I'm not sure exactly how to think
> about this. It seem to me that the Bernardinelli model allows the
> spatially correlated time effect v(i) to be independent of the
> spatially correlated areaspecific random effect u(i); whereas, the
> model that I propose forces simultaneous correlation across time and
> space for each location.
>
> Thanks for any help,
>
> Erik
>
Hi Erik:
Well, I haven't seen your model, but I suspect the two models are
different. The Bernardinelli
model has a random slope which seems to imply that there is a linear
relationship with time.
These are called "growth models" and I don't care for them. They make
sense only
in circumstances where a specific relationship with time makes sense:
although, that sounds like it is
obvious, these models are often misused, especially in psychiatric
clinical trials that I am familiar
with. So, in my opinion, your model makes alot more sense. Of
course,
you can try both models
and see which one is a better fit for your data. However, it is nice
to
have a hypothesis that one
will be better than another before trying it out. I think this would
lead to a nice test via a
Bayes Factor.
Rodney
From Darren Mayne <[log in to unmask]>
Sent Wednesday, July 26, 2006 8:50 pm
To ERIK E OSNAS <[log in to unmask]>
Cc
Bcc
Subject RE: spacetime correlation in WinBUGS
Hi Erik
I am certain you have already received plenty of advise from people
much
more qualified than me; however, you can find examples of space x time
models of the formulation you describe in Disease Mapping with WinBUGS
and
MLwiN (Lawson et al, 2003, pp 128133). The programs for implementing
the
examples in this book are also available from
http://www.sph.sc.edu/alawson/. I have not implemented these models
myself;
however, a college has a paper that may be of interest: Spatiotemporal
analysis of acute admissions for Ischemic Heart Disease in NSW,
Australia.
Burden S, GuhaS, Morgan G, Ryan L, Sparks R, Young L. Environmental and
Ecological Statistics 2005;12(4):427448.
Cheers
Darren
Erik Osnas
Assistant Scientist
Department of Wildlife Ecology
University of Wisconsin
218 Russell Labs
1630 Linden Dr.
Madison, WI 53706
6082621984

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
