Example of Logistic Regression
library ("R2jags")
############################# define the model for jags #######################
regmodel <- '
model{
beta0 ~ dnorm (0, 0.000001)
beta1 ~ dnorm (0, 0.01)
beta2 ~ dnorm (0, 0.01)
for (i in 1:n)
{
y[i] ~ dbern (ilogit(beta0 + beta1 * x1[i] + beta2 * x2[i]))
}
}
'
write (regmodel, file = "regmodel.bug")
##################### use jags to simulate MCMC ###############################
n <-50
ilogit <- function (x) 1/(1 + exp (-x))
x1 <- rnorm (n, 0, 2)
x2 <- rnorm (n, 0, 2)
truebeta <- c(0,2,2)
y <- rbinom (n, 1, ilogit( truebeta[1] + truebeta[2]*x1 + truebeta[3]*x2))
data <- list (x1 = x1, x2 = x2, y = y, n = n)
# look at the data
plot (x1,x2,col=y+1, pch = y+1)
abline(a=-truebeta[1]/truebeta[3],-truebeta[2]/truebeta[3])
# define a function to generate initial values for parameters
inits <- function ()
{
list (beta0 = 0, beta1 = 0, beta2 = 0)
}
# call jags to simulate MCMC
fitj <- jags (model.file = "regmodel.bug",
data = data, inits = inits,
parameters = c("beta0", "beta1","beta2"),
n.chains = 4, n.thin = 1, n.burnin = 1000,
n.iter = 10000)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 357
##
## Initializing model
traceplot (fitj,ask=FALSE)
## Inference for Bugs model at "regmodel.bug", fit using jags,
## 4 chains, each with 10000 iterations (first 1000 discarded)
## n.sims = 36000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta0 -0.335 0.584 -1.526 -0.715 -0.324 0.058 0.780 1.001 11000
## beta1 3.702 1.193 1.756 2.816 3.569 4.458 6.343 1.006 470
## beta2 2.531 0.758 1.263 1.976 2.465 3.011 4.185 1.005 570
## deviance 27.201 2.553 24.205 25.302 26.565 28.460 33.655 1.004 910
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 3.2 and DIC = 30.4
## DIC is an estimate of expected predictive error (lower deviance is better).
fitj.mcmc <- as.mcmc (fitj)
fitj.matrix <- as.matrix (fitj.mcmc)
colnames (fitj.matrix)
## [1] "beta0" "beta1" "beta2" "deviance"
# sample from posterior of (beta0, beta1)
plot (fitj.matrix[seq (1000,10000, by = 10),c("beta0", "beta1")], type = "p")
# sample from posterior of (beta1, beta2)
plot (fitj.matrix[seq (1000,10000, by = 10),c("beta1", "beta2")], type = "p")
# density of beta0
plot(density(fitj.matrix[,"beta0"]))
# density of beta1
plot(density(fitj.matrix[,"beta1"]))
# density of beta2
plot(density(fitj.matrix[,"beta2"]))
#Find Psuedo MLE
imle <- which.min(fitj.matrix[,"deviance"])
mlebeta <- fitj.matrix[imle,]; mlebeta
## beta0 beta1 beta2 deviance
## -0.2057432 2.8282113 1.9416948 23.9518673
# look classification boundaries
plot (x1,x2,col=y+1, pch = y+1)
# MCMC Estimates
for (i in seq (1000,5000, by = 50)) {
abline (a = -fitj.matrix[i,"beta0"]/fitj.matrix[i,"beta2"],
b = -fitj.matrix[i,"beta1"]/fitj.matrix[i,"beta2"],
col = "grey")
}
# True Values
abline(a=-truebeta[1]/truebeta[3],b=-truebeta[2]/truebeta[3],
col="black", lwd=3)
abline(a=-mlebeta[1]/mlebeta[3],b=-mlebeta[2]/mlebeta[3],
col="red", lwd=3)
points (x1,x2,col=y+1, pch = y+1)
legend("bottomleft", legend=c("TRUE", "MLE", "MCMC"),
col = c("black", "red", "grey"), lty = c(1,1,1))
# look at classification boundary coefficient
plot(density(fitj.matrix[,"beta1"]/fitj.matrix[,"beta2"]))
abline(v=truebeta[2]/truebeta[3], col="black")
abline(v=mlebeta[2]/mlebeta[3], col="red")
# look at classification sharpness coefficient
# plot(density(fitj.matrix[,"beta1"]^2+fitj.matrix[,"beta2"]^2))
# abline(v=truebeta[2]^2+truebeta[3]^2, col="black")
# abline(v=mlebeta[2]^2+mlebeta[3]^2, col="red")
Example of Mixture Model
library ("R2jags")
library ("MASS")
# user should provide ifold
mixmodel <- '
var
N, # number of observations
K, # number of components
y[N], # observations
Z[N], # true groups (labelled 1,2)
mu[K], # means of two groups
tau.mu, # prior precision for mu
tau[K], # sampling precision
sigma[K], # sampling standard deviation
P[K], # proportion in first group
alpha[K], # prior parameters for proportions
ind[N,K], # indicator function for group membership
tol[K], # total number of observations in each group
Itot[K]; # Censored indicator for tot (1 if tot[i] > 0; 0 otherwise)
model {
for (k in 1:K){
mu[k] ~ dnorm(20, 1.0E-4);
tau[k] ~ dgamma(0.01, 0.01 * 20);
sigma[k] <- 1/sqrt(tau[k]);
}
P[] ~ ddirch (alpha[]); # prior for mixing proportion
for (i in 1:N){
y[i] ~ dnorm(mu[Z[i]],tau[Z[i]]);
Z[i] ~ dcat(P[]);
for (k in 1:K) {
ind[i,k] <- (Z[i] == k);
}
}
for (k in 1:K) {
tot[k] <- sum(ind[,k]);
Itot[k] ~ dinterval(tot[k], 0);
}
for (i in 1:N){
logp.tr [i] <- log (dnorm (y[i], mu[Z[i]], tau[Z[i]]));
}
}
'
write (mixmodel, file = "mixmodel.jags")
K <- 2 ## number of components
y <- sample(c(rnorm (30, mean = 20), rnorm (70, mean=23)))
hist (y, nclass=15)
N <- length (y)
Itot <- rep (1, K)
alpha <- rep (1, K)
eyesdata <- list ( y = y, K = K, N = N, Itot = Itot, alpha = alpha)
inits <- function ()
{
list (#mu = rnorm (K, mean (y),20),
#tau = rep(1/var (y),K),
Z = sample (1:K, N, replace = T) )
}
mixmodel.jags <- jags (data = eyesdata, inits = inits, model.file = "mixmodel.jags",
parameters = c("P", "mu","sigma"),
n.thin = 5, n.iter = 5000, n.chain = 1)
## Compiling model graph
## Declaring variables
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 102
## Unobserved stochastic nodes: 105
## Total graph size: 827
##
## Initializing model
mixmodel.mcmc <- as.mcmc (mixmodel.jags)
summary (mixmodel.mcmc)
##
## Iterations = 2501:4996
## Thinning interval = 5
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## deviance 329.7696 30.4981 1.363918 4.03457
## mu[1] 21.9249 1.3915 0.062232 0.46314
## mu[2] 21.6052 1.4140 0.063236 0.36169
## P[1] 0.4584 0.2146 0.009596 0.04756
## P[2] 0.5416 0.2146 0.009596 0.04756
## sigma[1] 1.1667 0.3981 0.017804 0.07738
## sigma[2] 1.4571 0.7701 0.034438 0.06901
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## deviance 282.15039 303.3144 327.3941 352.1799 391.6427
## mu[1] 19.44722 20.4506 22.7037 22.9821 23.4917
## mu[2] 19.42948 20.4956 21.7819 22.8431 23.2793
## P[1] 0.09423 0.2741 0.4459 0.6432 0.8323
## P[2] 0.16768 0.3568 0.5541 0.7259 0.9058
## sigma[1] 0.48728 0.9220 1.0835 1.3649 1.9653
## sigma[2] 0.86026 1.0646 1.3215 1.7419 2.3469
traceplot (mixmodel.mcmc)
