if (FALSE) {
## simple regression example from the JAGS manual
jfun <- function() {
for (i in 1:N) {
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + beta * (x[i] - x.bar)
}
x.bar <- mean(x)
alpha ~ dnorm(0.0, 1.0E-4)
beta ~ dnorm(0.0, 1.0E-4)
sigma <- 1.0/sqrt(tau)
tau ~ dgamma(1.0E-3, 1.0E-3)
}
## data generation
set.seed(1234)
N <- 100
alpha <- 1
beta <- -1
sigma <- 0.5
x <- runif(N)
linpred <- crossprod(t(model.matrix(~x)), c(alpha, beta))
Y <- rnorm(N, mean = linpred, sd = sigma)
## list of data for the model
jdata <- list(N = N, Y = Y, x = x)
## what to monitor
jpara <- c("alpha", "beta", "sigma")
## write model onto hard drive
jmodnam <- write.jags.model(jfun)
## fit the model
regmod <- jags.fit(jdata, jpara, jmodnam, n.chains = 3)
## cleanup
clean.jags.model(jmodnam)
## model summary
summary(regmod)
}
## let's customize this model
jfun2 <- structure(
c(" model { ",
" for (i in 1:n) { ",
" Y[i] ~ dpois(lambda[i]) ",
" Y[i] <- alpha[i] + inprod(X[i,], beta[1,]) ",
" log(lambda[i]) <- alpha[i] + inprod(X[i,], beta[1,]) ",
" alpha[i] ~ dnorm(0, 1/sigma^2) ",
" } ",
" for (j in 1:np) { ",
" beta[1,j] ~ dnorm(0, 0.001) ",
" } ",
" sigma ~ dlnorm(0, 0.001) ",
" } "),
class = "custommodel")
custommodel(jfun2)
## GLMM
custommodel(jfun2, 4)
## LM
custommodel(jfun2, c(3,5))
## deparse when print
print(custommodel(jfun2), deparse=TRUE)
Run the code above in your browser using DataLab