Solved – Geostatistical analysis using spatial.exp in WinBugs

bayesianbugslogisticspatial

I have a logistic regression model in GeoBugs to estimate predictors of prevalence of a disease. Can anyone tell me if there is a simple way to determine the lower and upper bounds of phi for spatial.exp.

This is part of the model:

Infected[i] ~ dbin(p[i], tested[i])
logit(p[i]) <- alpha + beta1*covar1[i] + beta2*covar2[i]  + u[i] 
mu[i] <- 0
u[1:N] ~ spatial.exp(mu[], x[], y[], tau, phi, 1)

This how I have specified the priors:

phi ~ dunif(0.25, 20)

At the moment my selection has been arbitrary. The GeoBugs manual only goes as far as saying that the bounds of the distribution are determined by the data.

What I would like to know is there a dummies version of the best way to select the values.

Best Answer

I have worked this out myself. The lower bound for phi can be estiamted from

-ln(0.5)/(max separating distance between points)

To find the max separating distance I used the following code in R. My data are in a flat file with x and y coords renamed to long and lat respectively:

data <- read.csv(file="file.csv", header=T, sep=",")
coords <- data.frame(data$long,data$lat) 
library(sp)
pointDist <- apply(coords, 1, function(eachPoint) spDistsN1(as.matrix(coords), eachPoint, longlat=TRUE))
distances <- as.vector(pointDist)
max(distances)

Related Solutions

Solved – Sampling variables and calculating likelihood in WinBUGS/OpenBUGS

1) The probabilistic dnorm, dunif etc. functions are just describing the probability distribution which the variable on the left hand side (lhs) is assumed to have. If the variable is a parameter, then it's a prior distribution. If the variable is data, then it's p(data | parameters).

2) The distribution is not what is used for sampling. Don't think that pop.mean ~ dunif(0,5000) means that pop.mean is actually drawn from a Uniform(0,5000) distribution when the code is executed! The sampling is done using various algorithms inside WinBUGS / JAGS. The code describes a statistical model, not an algorithm to be executed.

3) As Tomas points out, the order of statements doesn't matter.

prec <- 1/pop.variance
pop.variance <- pop.sd * pop.sd
pop.sd ~ dunif(0,100)

merely tells the underlying sampler that 1) the prior on pop.sd is a U(0,100) distribution, 2) pop.variance = pop.sd*pop.sd, and 3) prec = 1/pop.variance. You can rearrange those statements however you like without changing their meaning.

Pretty much everything that happens goes on under the hood of WinBUGS / JAGS. The software takes the statistical model, as described by the code, and processes it using whatever algorithms it determines are "best". I'll repeat: the code describes the model, it does not specify an algorithm for processing the model + data. This is the hardest thing to get about WinBUGS / JAGS, or functional programming languages in general (at least, it was for me,) and once you get it, programming becomes much easier.

Solved – Bayesian Analysis in R/WinBUGS and SAS (Proc MCMC)

See if something like this works:

proc mcmc
 data = /* your data here */
 nmc = 20000, nbi = 1000; /* iterations and burn-in */

 /* Specify parameters -- separate statements means they're proposed separately */
 parms a;
 parms b;
 parms theta;

 /* Specify priors */
 hyper a ~ gamma(1.5, scale = .25);
 hyper b ~ gamma(1.5, scale = .25); 
 prior theta ~ gamma(a, scale = b);

 /* Sampling model */
 lam = theta*time;
 model fail ~ poisson(lam);
run;

Best Answer

Related Solutions

Solved – Sampling variables and calculating likelihood in WinBUGS/OpenBUGS

Solved – Bayesian Analysis in R/WinBUGS and SAS (Proc MCMC)

Related Question