## simulation for Poisson GLMM
set.seed(1234)
n <- 20
beta <- c(2, -1)
sigma <- 0.1
alpha <- rnorm(n, 0, sigma)
x <- runif(n)
X <- model.matrix(~x)
linpred <- crossprod(t(X), beta) + alpha
Y <- rpois(n, exp(linpred))
## JAGS model as a function
jfun1 <- function() {
for (i in 1:n) {
Y[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha[i] + inprod(X[i,], beta)
alpha[i] ~ dnorm(0, 1/sigma^2)
}
for (j in 1:np) {
beta[j] ~ dnorm(0, 0.001)
}
sigma ~ dlnorm(0, 0.001)
}
## data
jdata <- list(n = n, Y = Y, X = X, np = NCOL(X))
## inits with latent variable and parameters
ini <- list(alpha=rep(0,n), beta=rep(0, NCOL(X)))
## function to update inits
ifun <- function(model, n.clones) {
list(alpha=dclone(rep(0,n), n.clones),
beta=coef(model)[-length(coef(model))])
}
## iteartive fitting
jmod <- dc.fit(jdata, c("beta", "sigma"), jfun1, ini,
n.clones = 1:5, multiply = "n", unchanged = "np",
initsfun=ifun)
## summary with DC SE and R hat
summary(jmod)
dct <- dctable(jmod)
plot(dct)
## How to use estimates to make priors more informative?
glmm.model.up <- function() {
for (i in 1:n) {
Y[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha[i] + inprod(X[i,], beta[1,])
alpha[i] ~ dnorm(0, 1/sigma^2)
}
for (j in 1:p) {
beta[1,j] ~ dnorm(priors[j,1], priors[j,2])
}
sigma ~ dgamma(priors[(p+1),2], priors[(p+1),1])
}
## function for updating, x is an MCMC object
upfun <- function(x) {
if (missing(x)) {
p <- ncol(X)
return(cbind(c(rep(0, p), 0.001), rep(0.001, p+1)))
} else {
par <- coef(x)
return(cbind(par, rep(0.01, length(par))))
}
}
updat <- list(n = n, Y = Y, X = X, p = ncol(X), priors = upfun())
dcmod <- dc.fit(updat, c("beta", "sigma"), glmm.model.up,
n.clones = 1:5, multiply = "n", unchanged = "p",
update = "priors", updatefun = upfun)
summary(dcmod)
## time series example
## data and model taken from Ponciano et al. 2009
## Ecology 90, 356-362.
paurelia <- c(17,29,39,63,185,258,267,392,510,
570,650,560,575,650,550,480,520,500)
dat <- list(ncl=1, n=length(paurelia), Y=dcdim(data.matrix(paurelia)))
beverton.holt <- function() {
for (k in 1:ncl) {
for(i in 2:(n+1)){
## observations
Y[(i-1), k] ~ dpois(exp(X[i, k]))
## state
X[i, k] ~ dnorm(mu[i, k], 1 / sigma^2)
mu[i, k] <- X[(i-1), k] + log(lambda) - log(1 + beta * exp(X[(i-1), k]))
}
## state at t0
X[1, k] ~ dnorm(mu0, 1 / sigma^2)
}
# Priors on model parameters
beta ~ dlnorm(-1, 1)
sigma ~ dlnorm(0, 1)
tmp ~ dlnorm(0, 1)
lambda <- tmp + 1
mu0 <- log(2) + log(lambda) - log(1 + beta * 2)
}
mod <- dc.fit(dat, c("lambda","beta","sigma"), beverton.holt,
n.clones=c(1, 2, 5, 10), multiply="ncl", unchanged="n")
## compare with results from the paper:
## beta = 0.00235
## lambda = 2.274
## sigma = 0.1274
summary(mod)
## Using WinBUGS/OpenBUGS
library(R2WinBUGS)
data(schools)
dat <- list(J = nrow(schools), y = schools$estimate,
sigma.y = schools$sd)
bugs.model <- function(){
for (j in 1:J){
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}
inits <- function(){
list(theta=rnorm(nrow(schools), 0, 100), mu.theta=rnorm(1, 0, 100),
sigma.theta=runif(1, 0, 100))
}
param <- c("mu.theta", "sigma.theta")
if (.Platform$OS.type == "windows") {
sim2 <- dc.fit(dat, param, bugs.model, n.clones=1:2,
flavour="bugs", program="WinBUGS", multiply="J",
n.iter=2000, n.thin=1)
summary(sim2)
}
sim3 <- dc.fit(dat, param, bugs.model, n.clones=1:2,
flavour="bugs", program="brugs", multiply="J",
n.iter=2000, n.thin=1)
summary(sim3)
library(R2OpenBUGS)
sim4 <- dc.fit(dat, param, bugs.model, n.clones=1:2,
flavour="bugs", program="openbugs", multiply="J",
n.iter=2000, n.thin=1)
summary(sim4)Run the code above in your browser using DataLab