Editing Software Plugins and Codes (section)

=R=
==JAGS Statistics==
Analysis of Bayesian models using Markov Chain Monte Carlo.

===Overview===

JAGS is a replacement for WinBUGS that can be run from within R.  It is cross-platform and will run on any system (Windows, OS-X, Linux).

===Installation===

====JAGS====
You can get get JAGS from here: http://www-fis.iarc.fr/~martyn/software/jags/

====R====

You can get R from here: http://www.r-project.org/ 

====runjags====

You will also need the runjags package from R.  This is the webpage for it: http://cran.stat.sfu.ca/web/packages/runjags/index.html .  There is a very good PDF manual there, with examples.  (The help is also available from within R).

Note that there is an alternate package rjags, which I haven't tried.  I'm not sure which is better.

Installation of runjags is very easy.  From within R, go to the menu bar and select Packages&Data -> Package Installer.  Click Get List, then enter runjags in the search field.  Select it and click Install Selected.  You will also need to install the package coda.

Now go to Packages&Data ->Package Manager, scroll down and find coda and runjags and make sure that they are both selected.  Note that if you click on runjags here will will get it's on-line help, which can be useful.

===Reading a spreadsheet into R===

Refer to http://cran.r-project.org/doc/manuals/R-intro.html#The-read_002etable_0028_0029-function .  There is also a whole import/export manual for R at http://cran.r-project.org/doc/manuals/R-data.html

The first step is to export your spreadsheet as a text file from your spreadsheet program.  You might want to copy only the relevant columns into a new spreadsheet.  

If you are using Excel, I believe that you can export as a text file, although I am not certain of this.

I was using Apple's Numbers.  In Numbers, select the spreadsheet or the columns that you want to Export, then start up TextEdit.  Make sure that the new blank text document is plain text (Format -> Make Plain Test), then Paste in the spreadsheet.  This resulted in a text file in exactly the right format for importing into R, like this

 cage	fnbmd10
 20	0.72
 20.1	0.74
 20.1	1.02
 20.2	0.97
 20.2	0.84
 20.4	0.69
 20.4	0.98
 20.5	
 20.6	0.82

It may happen that there are some missing data in your columns (as in the above example).  This command (from the command line) will copy all lines from file fnbmd10.data to file fnbmd10_trimmed.data that have exactly two fields, which is what we want:

 awk '(NF==2) {print $0}' fnbmd10.data > fnbmd10_trimmed.data

If you have a column for sex, and you want to split it into two files, one for each sex, you can do this:

 awk '($2==0) {print $0}' input.txt > output_w.txt
 awk '($2==1) {print $0}' input.txt > output_m.txt

In this example, the sex is the second field (hence $2).  You will need to manually add the column headings again after doing this.

Once your text file looks correct, you can load it into R with the following command (from within R obviously):

 dataframe <- read.table("fnbmd10_trimmed.data", header=TRUE)

I recommend at this point using some of R's graphical capabilities to plot and examine the data, just to make sure it seems to be reasonable.

 plot (dataframe$cage, dataframe$fnbmd10)

====Tip: Setting the Working Directory in R====

 setwd("/Users/eric/Documents/change_point")

====Tip: Converting line endings====

In case it's required, this can be done from a command prompt in OS-X or Linux:

 cat input.txt | tr '\r' '\n' > output.txt

===Model file===

Here is the JAGS model I used, findcp.model .  I suggest thinking carefully about what distributions, with what values, are most appropriate for b.intercept, b.slope and tau as this will affect how the solution converges.

<source lang=text>
#############################################
# Obtain slopes for females                 #
# 25 to 45 for the BMD value.               #
# Determine the cut point from 25 to 45 .   #
#############################################

var
  n,           # Read from data file; number of subjects
  age[n],      # Read from data file; age of subject
  y[n],        # Read from data file; BMD of subject
  b.cutoff,    # change point age
  b.intercept, # constant value of BMD function
  b.slope,     # slope of BMD function
  mu[n],       # expected mean BMD of subject
  tau,         # tau of the normal distribution of BMD
  sigma;       # std. dev. of the normal distribution of BMD

model {
    for (i in 1:n) {
        mu[i] <- b.intercept + b.slope * step(age[i]-b.cutoff) * (age[i]-b.cutoff)
        y[i] ~ dnorm(mu[i],tau)
    }
    b.cutoff ~ dunif(25,55)
    b.intercept ~ dnorm(.9,6)
    b.slope ~ dnorm(0,9)
    tau ~ dgamma(1.0E-3,1.0E-3)
    sigma <- 1.0/sqrt(tau)
}
</source>

====Tip: Using a BUGS style model file instead of a JAGS style model file====

JAGS and BUGS model files are very similar, however, the BUGS model file allows more type of things in the file, such as initialization values, that are not allowed in a JAGS model file (in JAGS these go in a separate file).

runjags can read BUGS-style files directly, without modification, with the R command autorun.jagsfile.  This is instead of autorun.jags.  See Below.

===Running a model in R===

This is pretty much just following the example from the runjags documentation for autorun.jags .

We determine the number of data points and stick our data in an R list format:

 n <- length(dataframe$cage)
 jagsdata <- list (age=dataframe$cage, y=dataframe$fnbmd10, n=n)

We read in the model as a test string.  This worked for me, but I suspect that there is an easier way.

 modelraw <- readBin("findcp.model", raw(), 9000)
 model <- rawToChar(modelraw)

The last argument to readBin is the number of bytes to read.  Just make sure that this is bigger than the file size and you will get the whole file.

Now we run the model:

 results <- autorun.jags(model=model, monitor=c("b.cutoff", "b.slope", "b.intercept"), data=jagsdata)

You can analyze the traceplots of the results to assess convergence:

 results$trace

To get a summary of the monitored variables:

 results$summary

This will look something like this:

<source lang=text>
Iterations = 1:21741
Thinning interval = 2 
Number of chains = 2 
Sample size per chain = 10871 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean        SD  Naive SE Time-series SE
b.cutoff    43.152034 5.9029519 4.003e-02      1.799e-01
b.slope     -0.003935 0.0006817 4.624e-06      1.866e-05
b.intercept  0.832309 0.0102962 6.983e-05      2.305e-04

2. Quantiles for each variable:

                 2.5%       25%       50%       75%     97.5%
b.cutoff    31.641838 39.016150 42.733250 47.586850 53.975075
b.slope     -0.005512 -0.004329 -0.003827 -0.003447 -0.002851
b.intercept  0.813198  0.825117  0.832013  0.839125  0.853183
</source>

===Summary of places to find Documentation===

* The OpenBUGS manual: http://www.openbugs.info/Manuals/Manual.html .
* The JAGS manual : http://www-fis.iarc.fr/~martyn/software/jags/jags_user_manual.pdf .
*The manual for runjags : http://cran.stat.sfu.ca/web/packages/runjags/runjags.pdf .  Note that runjags help is also available on-line in R.  For example, entering help(autorun.jags) will bring up the help on that command.

===OBSOLETE: JAGS run script===

This is not required when running JAGS from within R, which seems to be the best way to run JAGS.

Here is the script I used to run JAGS, fnbmd10.jags.  Of course, you could enter all these commands interactively into JAGS as well.  You probably want to customize this.  For example, I just guessed at the number of iterations.  There are also several output options.  Probably you want to use the coda command to save the monitor values in files that can be read by the coda package of R, and then do some sort of statistical analysis on the convergence.

<source lang=text>
# Read in the model and the data
model in findcp.model
data in fnbmd10.R

# Required stuff
compile
initialize

# Do some iterations
# (The adaptive part of the iterating needs to be done before monitors can be set)
update 500

# Now set some monitors to track parameters over iterations 
monitor b.cutoff, thin(50)
monitor b.intercept, thin(50)
monitor b.slope, thin(50)

# Do lots of iterations
update 10000

# Save parameters and monitors.
# There are other kinds of outputs and files you can write; refer to the
# JAGS documentation
parameters to parameters.out
monitors to monitors.out

# Bye
exit
</source>


[[Image:Redirect.jpg|15px]] [[#top|''Return to Contents'']] or [[Main_Page|''Return to Main Page'']]