5.4 Estimación Bayesiana usando WinBUGS u OPENBUGS dentro de R Como hemos visto en la sección anterior con los dos archivos que genera el BRMUW podemos implementar la Estimación Bayesiana con el programa WinBUGS u OpenBugs. Alternativamente también se puede implementar la estimación bayesiana utilizando interfaces del R con WinBUGS u OpenBUGS. Para esto necesitaremos el archivo de texto original con los datos en columna:
sierr permedyc partco concoca pobrez 1 2 2 1 2 0 0 6 1 2 . . . . . . . . . . . . . . . 1 2 9 0 3 este archivo lo guardaremos con el nombre datos.txt y del archivo de sintaxis del modelo generado en BRMUW deberemos copiar solo la sintaxis del modelo. En el ejemplo abajo implementamos el modelo de regresión binaria Skew logit.
esto es debemos copiar toda la sintaxis antes de “Inits” y guardarla en un archivo, en este ejemplo modelo.txt. Entonces el archivo modelo.txt quedaría como
Luego para implementar la estimación bayesiana en el R seguiremos los siguientes pasos para usar la librería R2WinBUGS. 1. Dentro de R, cargar la librería R2WinBUGS con el siguiente comando que ha sido instalada previamente library(R2WinBUGS) 2. Leer los datos (el archivo datos.txt en este ejemplo esta en la carpeta C:\BRMUW\) datos <- read.table("C:/BRMUW/datos.txt", header=TRUE, sep="", na.strings="NA", dec=".",strip.white=TRUE) 3. Crear una lista que contenga los datos y la información que esta debajo de Data en el archivo generado por el BRMUW
Esto la hacemos con el siguiente comando data<-c(as.list(datos),list(n=1947,k=5)) 4. Crear un programa que genere los valores iniciales.
inits<-function(){list(beta=c(0,0,0,0,0),delta=0.5) 5. Finalmente con el comando bugs se implementa la estimación bayesiana. Aquí explicaremos brevemente la sintaxis del comando Bugs parameters.to.save = es un vector con los nombres de los parámetros del modelo cuyas simulaciones deseamos guardar. model.file = es la dirección donde se encuentra el archivo del modelo. n.chains = es el número de cadenas a ser generadas. n.iter = es el número iteraciones totales de cada cadena. n.burnin = es el número de iteraciones a ser descartadas como burn-in. program = es el programa que se utilizará para implementar la estimación bayesiana. n.burnin = es el número de iteraciones a ser descartadas como burn-in. Luego con el siguiente comando se implementa la estimación bayesiana y las simulaciones son guardadas en el objeto salida. salida<-bugs(data,inits,parameters.to.save=c("beta","lambda"), model.file="C:/BRMUW/modelo.txt", n.chains=1, n.iter=2000, n.burnin=1000,program="OpenBUGS") Si ponemos salida en la línea de comandos del R obtenemos un resumen de la simulación > salida Inference for Bugs model at "C:/Users/vaio/Documents/Jorge/modelo.txt", fit using OpenBUGS, 1 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10 n.sims = 1000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% beta[1] -2.8 0.9 -4.2 -3.5 -2.9 -1.8 -1.5 beta[2] 0.7 0.1 0.5 0.6 0.7 0.7 0.8 beta[3] 0.1 0.0 0.1 0.1 0.1 0.1 0.1 beta[4] -0.2 0.1 -0.3 -0.2 -0.2 -0.2 -0.1 beta[5] 0.5 0.1 0.4 0.4 0.5 0.7 0.8 lambda 0.8 0.9 -0.8 0.0 0.8 1.6 2.3 deviance 4984.5 155.5 4664.7 4862.6 5043.7 5123.0 5151.8
DIC info (using the rule, pD = Dbar-Dhat) pD = -90.0 and DIC = 4894.0 DIC is an estimate of expected predictive error (lower deviance is better). 6. Finalmente, para mayores detalles sobre el comando bugs se puede consultar la ayuda escribiendo en la línea de comandos ?bugs |