[R-sig-Geo] openbugs on linux

Thu Jul 15 11:34:33 CEST 2010

  Hi, please check this code.
Is a comparison between INLA and OpenBUGS/BRugs, using a quite simple model.
(INLA -Integrated Nested Laplace Approximation- performs Bayesian 
inference without using MCMC)
The dataset of the example is included in the INLA package, so you 
should get it from http://www.r-inla.org/.

Hope it helps.
ƒacu.-

PD: Sorry about the comments in spanish.

library(BRugs)
library(INLA)

data(Epil) # datos en formato adecuado para inla, provenientes del 
paquete INLA
Epil.inla <- Epil
load(file="Epil.RData") # mismos datos en formato adecuado para bugs, 
sacados del WinBUGS
attach(Epil)

# BUGS

epil.model <- function()
{
for(j in 1:N)
{
for(k in 1:T)
{
y[j,k] ~ dpois(mu[j,k])
log(mu[j,k]) <- a0 + a.Base * (l.Base4[j]-bar.l.Base4)
+ a.Trt * (Trt[j]-bar.Trt)
+ a.BT * (BT[j]-bar.BT)
+ a.Age * (l.Age[j]-bar.l.Age)
+ a.V4 * (V4[k]-bar.V4)
+ b1[j] + b[j,k]
b[j,k] ~ dnorm(0.0, tau.b) # Efecto aleatorio individuo*visita
}
b1[j] ~ dnorm(0.0, tau.b1) # Efecto aleatorio del individuo
BT[j] <- Trt[j] * l.Base4[j] # Interacción
l.Base4[j] <- log(Base[j]/4)
l.Age[j] <- log(Age[j])
}
# Medias de covariables
bar.l.Base4 <- mean(l.Base4[])
bar.Trt <- mean(Trt[])
bar.BT <- mean(BT[])
bar.l.Age <- mean(l.Age[])
bar.V4 <- mean(V4[])

# Previas
a0 ~ dnorm(0.0, 1.0E-4)
a.Base ~ dnorm(0.0, 1.0E-4)
a.Trt ~ dnorm(0.0, 1.0E-4)
a.BT ~ dnorm(0.0, 1.0E-4)
a.Age ~ dnorm(0.0, 1.0E-4)
a.V4 ~ dnorm(0.0, 1.0E-4)
# sigma.b1~ dunif(0, 100)
# sigma.b ~ dunif(0, 100)
# tau.b1 <- 1.0/(sigma.b1*sigma.b1)
# tau.b <- 1.0/(sigma.b*sigma.b)
tau.b1 ~ dgamma(1.0E-3,1.0E-3); sigma.b1 <- 1.0 / sqrt(tau.b1)
tau.b ~ dgamma(1.0E-3,1.0E-3); sigma.b <- 1.0/ sqrt(tau.b)

# recalcular la interceptación en la escala original
alpha0 <- a0 - a.Base*bar.l.Base4 - a.Trt*bar.Trt - a.BT*bar.BT - 
a.Age*bar.l.Age - a.V4*bar.V4
}

epil.data <- Epil # lista con N, T y los valores de todas las variables. 
La variable respuesta es bidimensional NxT.

epil.ini <- list(
list(a0 = 1, a.Base = 0, a.Trt = 0, a.BT = 0, a.Age = 0, a.V4 = 0, 
sigma.b1 = 1, sigma.b = 1, b1=rep(0,N), b=matrix(0,N,T)),
list(a0 = 0, a.Base = 1, a.Trt = 1, a.BT = 1, a.Age = 1, a.V4 = 1, 
sigma.b1 =10, sigma.b =10, b1=rep(1,N), b=matrix(1,N,T)),
list(a0 =-1, a.Base = 1, a.Trt =-1, a.BT = 1, a.Age =-1, a.V4 = 1, 
sigma.b1 =10, sigma.b = 1, b1=rep(-1,N),b=matrix(1,N,T)))

# para previas gamma sobre la precisión
epil.ini <- list(
list(a0 = 1, a.Base = 0, a.Trt = 0, a.BT = 0, a.Age = 0, a.V4 = 0, 
tau.b1 = 1, tau.b = 1, b1=rep(0,N), b=matrix(0,N,T)),
list(a0 = 0, a.Base = 1, a.Trt = 1, a.BT = 1, a.Age = 1, a.V4 = 1, 
tau.b1 =10, tau.b =10, b1=rep(1,N), b=matrix(1,N,T)),
list(a0 =-1, a.Base = 1, a.Trt =-1, a.BT = 1, a.Age =-1, a.V4 = 1, 
tau.b1 =10, tau.b = 1, b1=rep(-1,N),b=matrix(1,N,T)))

writeModel(epil.model, file.path(tempdir(), "epilmodel.txt"))
bugsData(epil.data, file.path(tempdir(), "epildata.txt"))
inifiles <- paste(file.path(tempdir(), "epilini"),1:3,".txt", sep="")
bugsInits(epil.ini, numChains = length(epil.ini), fileName=inifiles)

# Parámetros de interés
pts=c("alpha0", "a.Base", "a.Trt", "a.BT", "a.Age", "a.V4", "sigma.b1", 
"sigma.b")

# Inicializaciones
nburn = 5000 # etapa de quemado para llegar a convergencia en la cadena 
MCMC
nsim = 10000 # tamaño de la cadena MCMC buena

# Simulación
fit<-BRugsFit("epilmodel.txt",data="epildata.txt",inits=inifiles, 
,numChains = length(epil.ini), parametersToSave=pts,nBurnin = nburn, 
nIter = nsim, working.directory=tempdir())

# Timing
# 5000 updates took 79 s (1'19'')
# 10000 updates took 173 s (2'53'')

save(fit, file="brugs.fit.unif.RData")
save(fit, file="brugs.fit.gamma.RData")
load(file="brugs.fit.unif.RData")
load(file="brugs.fit.gamma.RData")

# Review
fit
s.h <- samplesHistory(c('*'))
s.d <- samplesDensity("*")

# INLA
detach(Epil)
attach(Epil.inla)

formula <- y ~ log(Base/4) + Trt + I(Trt*log(Base/4)) + log(Age) + V4 + 
f(Ind, model = "iid") + f(rand, model = "iid")
model=inla(formula,family="poisson")
# 
model=inla(formula,family="poisson",keep.data.files=TRUE,control.fixed=list(param=c(0.8,0.02)))
# hyper<-hyperpar.inla(model)
# con todas las opciones:
# model=inla(formula,family="poisson", control.compute=list(dic=T, 
mlik=T), control.predictor=list(return.marginals=T))

# Timing
# Prepare inla file.....nr=2
# ..done
# Run inla...
# Wall-clock time used on [Model.ini] using max_threads=[2]
# Preparations : 0.125 seconds
# Approx inference: 1.688 seconds
# Output : 0.406 seconds
# ---------------------------------
# Total : 2.219 seconds

# Comparación inferencia
par(mfrow=c(4,2))
plotDensity("alpha0", ask=F)
lines(model$marginals.fixed$intercept)
plotDensity("a.Base", ask=F)
lines(model$marginals.fixed$"log(Base/4)")
plotDensity("a.Trt", ask=F)
lines(model$marginals.fixed$Trt)
plotDensity("a.BT", ask=F)
lines(model$marginals.fixed$"I(Trt * log(Base/4))")
plotDensity("a.Age", ask=F)
lines(model$marginals.fixed$"log(Age)")
plotDensity("a.V4", ask=F)
lines(model$marginals.fixed$V4)
# como prec=1/sigmasq, -> fsigma(1/sqrt(y)) = 2*y^(1.5)*fprec(y)
plotDensity("sigma.b1", ask=F)
lines(cbind(1/sqrt(model$marginals.hyperpar$"Precision for ind"[,1]),
2*model$marginals.hyperpar$"Precision for 
ind"[,1]^(3/2)*model$marginals.hyperpar$"Precision for ind"[,2]))
plotDensity("sigma.b", ask=F)
lines(cbind(1/sqrt(model$marginals.hyperpar$"Precision for rand"[,1]),
2*model$marginals.hyperpar$"Precision for 
rand"[,1]^(3/2)*model$marginals.hyperpar$"Precision for rand"[,2]))

## more acurate hyperparameters?
#samplesDensity("sigma.b1",mfrow=c(1,1))
#lines(cbind(1/sqrt(hyper$marginals$"Precision for ind"[,1]),
# 2*hyper$marginals$"Precision for 
ind"[,1]^(3/2)*hyper$marginals$"Precision for ind"[,2]))
#
#samplesDensity("sigma.b",mfrow=c(1,1))
#lines(cbind(1/sqrt(hyper$marginals$"Precision for rand"[,1]),
# 2*hyper$marginals$"Precision for 
rand"[,1]^(3/2)*hyper$marginals$"Precision for rand"[,2]))

# Comparación predicciones inla

esp <- model$summary.fitted.values[,"mean"] # valor esperado
esp <- cbind(esp, qpois(.05, esp), qpois(.95, esp)) # IC90% centrado
out <- which(y>esp[,3] | y<esp[,2]) # 40 62 63 221
n = Epil$N*Epil$T

plot(sqrt(y), pch=19, cex=.5, col='navyblue') # sqrt(observaciones) y 
valores esperados
points(sqrt(esp[,1]), lwd=1);
segments(1:n, sqrt(esp[,2]), 1:n, sqrt(esp[,3]), col='gray')
points(out,sqrt(y[out]), pch=19, cex=.5, col='red') # outliers

#res.inla <- (y-esp[,1])/sqrt(esp[,1]) # residuos de Pearson
#plot(res.inla,pch=19, cex=.5, col='red')
#qqnorm(res.inla); qqline(res.inla);
#shapiro.test(res.inla)

# Bondad de ajuste de modelos (mismo modelo, distintos métodos)

fit$DIC # BRugs
# Dbar Dhat DIC pD
#y 1037 916.7 1157 120.2
model$dic # inla
#mean of the deviance: 1039.8321
#deviance of the mean: 921.4959
#effective number of parameters: 118.3363
#dic: 1158.1684

On 15/07/10 07:45, Sathyanarayan Anand wrote:
> Hi,  Can someone please point me towards running OpenBUGS on linux? I have copied the files as instructed on the OpenBUGS website but I can find any examples on how to write a script file to execute my model. I have the model, data and inits as text files. If possible, I'd prefer to run this through R on linux using the BRugs package, but I'm not sure if the package is supported on Linux. Thanks,Sathya.		 		 	   		
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>