aboutsummaryrefslogtreecommitdiff
path: root/README.md
blob: ff6f638355c1f841711074063c44c96868330e41 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
# jagsUI: Run JAGS from R

<!-- badges: start -->
[![CRAN status](https://www.r-pkg.org/badges/version/jagsUI)](https://cran.r-project.org/web/packages/jagsUI/index.html)
[![R build status](https://github.com/kenkellner/jagsUI/workflows/R-CMD-check/badge.svg)](https://github.com/kenkellner/jagsUI/actions)
<!-- badges: end -->

This package runs `JAGS` (Just Another Gibbs Sampler) analyses from within `R`. It acts as a wrapper and alternative interface for the functions in the `rjags` package and adds some custom output and graphical options. It also makes running chains in parallel quick and easy.

## Installation

You can install the package from [CRAN](https://cran.r-project.org/web/packages/jagsUI/index.html), or get the development version from Github:

```r
devtools::install_github('kenkellner/jagsUI')
```

You will also need to separately install JAGS, which you can download [here](https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/).

## Example


```r
library(jagsUI)
```

Format data:


```r
jags_data <- list(
  gnp = longley$GNP,
  employed = longley$Employed,
  n = nrow(longley)
)
```

Write BUGS model file:


```r
modfile <- tempfile()
writeLines("
model{

  # Likelihood
  for (i in 1:n){ 
    # Model data
    employed[i] ~ dnorm(mu[i], tau)
    # Calculate linear predictor
    mu[i] <- alpha + beta*gnp[i]
  }
    
  # Priors
  alpha ~ dnorm(0, 0.00001)
  beta ~ dnorm(0, 0.00001)
  sigma ~ dunif(0,1000)
  tau <- pow(sigma,-2)

}
", con=modfile)
```

Set initial values and parameters to save:


```r
inits <- function(){  
  list(alpha=rnorm(1,0,1),
       beta=rnorm(1,0,1),
       sigma=runif(1,0,3)
  )  
}

params <- c('alpha','beta','sigma')
```

Run JAGS:


```r
out <- jags(data = jags_data,
            inits = inits,
            parameters.to.save = params,
            model.file = modfile,
            n.chains = 3,
            n.adapt = 100,
            n.iter = 1000,
            n.burnin = 500,
            n.thin = 2)
```

```
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 3
##    Total graph size: 74
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 500 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
```

View output:


```r
out
```

```
## JAGS output for model '/tmp/Rtmp9aOSoz/file15fda23cf96c2', generated by jagsUI.
## Estimates based on 3 chains of 1000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 500 iterations and thin rate = 2,
## yielding 750 total samples from the joint posterior. 
## MCMC ran for 0.001 minutes at time 2024-01-21 17:29:30.558591.
## 
##            mean    sd   2.5%    50%  97.5% overlap0 f  Rhat n.eff
## alpha    51.911 0.739 50.539 51.911 53.402    FALSE 1 0.999   750
## beta      0.035 0.002  0.031  0.035  0.038    FALSE 1 1.000   750
## sigma     0.715 0.148  0.487  0.699  1.095    FALSE 1 1.006   535
## deviance 33.247 2.798 30.050 32.466 40.263    FALSE 1 1.002   750
## 
## Successful convergence based on Rhat values (all < 1.1). 
## Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## For each parameter, n.eff is a crude measure of effective sample size. 
## 
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 3.9 and DIC = 37.161 
## DIC is an estimate of expected predictive error (lower is better).
```

## Acknowledgments

* Martyn Plummer, developer of the excellent JAGS software package and the `rjags` R package.
* Andrew Gelman, Sibylle Sturtz, Uwe Ligges, Yu-Sung Su, and Masanao Yajima, developers of the `R2WinBUGS` and `R2jags` packages on which the package was originally based.
* Robert Swihart, Marc Kery, Jerome Guelat, Michael Schaub, and Mike Meredith who tested and provided helpful suggestions and improvements for the package.