Drawing path diagrams - options part 1

One of my main projects at Charles University involves path analysis. A problem I came up against in my PhD work was with finding a time-efficient way to draw multiple path diagrams. I didn't put too much effort into finding a solution back then and just used PowerPoint and Inkscape to draw diagrams. However, when I had to make small modifications to the diagrams and produce numerous diagrams this method quickly became overly time consuming and annoying.

My current project is going to demand even more complex diagrams. I want to find a better method.

Some initial reading has led me to the idea of doing-away with the manual drawing methods and I now think that a command line programming method could be better in the long run. Although there is usually a steep learning curve at the start, which is a bit off-putting.  Luckily, I can already use Latex and R so it's not too intimidating and I'll give it a try.

Here is a summary of what I've learned today:

This post from Andrew Wheeler on stack exchange is where I started off and I found it helpful because he describes the exact issues I've been having. The post steered me towards thinking that the Tikz/pgf drawing library in Latex will be the way to go rather than Graphviz. As I saw form another post, Graphviz has trouble drawing curved arrows and also can't handle some of the mathematical notation one might need.

I also rejected the idea of using the in-built visualisation functions from the sem package in R because I'm not planning on using sem to run the analysis in R [I want to use some method for confirmatory path analysis - to be decided if it's going Bayesian or not. For now, I need it to be phylogenetic, so thinking phylogenetic confirmatory path analysis].

Here is a basic diagram. It's only taken me a couple of hours to get to this stage. It's not pretty yet and doesn't have all my variables, but it's an ok start. Mostly, it's adapted from this code by Ivan Griffin.

The document class doesn't actually convert it into a .png. A thread says that this will be possible eventually. But having it like this does at least produce a cropped .pdf document.

Latex code:


\documentclass[convert={density=300,size=1080x800,outext=.png}]{standalone}
\usepackage{tikz} % Load tikz package
\usetikzlibrary{fit,positioning,calc,backgrounds}
\begin{document}
\centering
\begin{tikzpicture} % Encloses drawings in tikz env.

% Styles for states, and state edges
\tikzstyle{state} = [draw, very thick, fill=blue!10, rectangle, minimum height=3em, minimum width=7em, node distance=8em, font={\sffamily\bfseries}]
\tikzstyle{stateEdgePortion} = [black,thick];
\tikzstyle{stateEdge} = [stateEdgePortion,->];
\tikzstyle{edgeLabel} = [pos=0.5, text centered, font={\sffamily\small}];
\tikzstyle{main}=[circle, minimum size = 10mm, thick, draw =red!80, node distance = 16mm]
\tikzstyle{connect}=[-latex, thick]
\tikzstyle{box}=[rectangle, draw=green!100]

% Position the nodes (boxes)
\node[state, name=mrt] {MRT};
\node[state, name=gridcells, below of=mrt, left of=mrt, xshift=-2em] {Grid cells};
\node[state, name=propagule, below of=gridcells] {Propagule length};
\node[state, name=dna, below of=propagule, right of=propagule, xshift=2em] {DNA};
\node[state, name=states, below of=gridcells, right of=gridcells, xshift=20em, node distance=4em] {States};

% Connect the nodes (boxes) via edges (arrows)
\draw ($(propagule.north) + (-.0em,0)$)
edge[stateEdge] node[edgeLabel, xshift=-0em]{\emph{}}
($(gridcells.south) + (-.0em,0)$);
\draw ($(gridcells.north) + (.0em,0)$)
edge[stateEdge, bend left=22.5] node[edgeLabel, xshift=-0em]{\emph{}}
($(mrt.west) + (.0em,0)$);
\draw ($(dna.west) + (-0em,0)$)
edge[stateEdge, bend left=22.5] node[edgeLabel, xshift=-0em, yshift=0em]{}
($(propagule.south) + (0,0em)$);
\draw ($(mrt.east) + (0em,0)$)
edge[stateEdge] node[edgeLabel, xshift=0em, yshift=0em]{\emph{}}
($(states.north) + (0,0em)$);

\end{tikzpicture}
\end{document}
% note - compiled with pdflatex

The importance of meta-data

Creating and maintaining accurate meta-data for a database is important, but often overlooked.

Meta-data is especially important when working on large databases with many collaborators. You may know what you've done, but to others it's not always obvious. Equally, come back to your database after a few months away and you may not remember what you did, what unit the data is in, or it's source.

I am working on a large collaborative project and received a database with relatively little meta-data. I am now spending the day deciphering what each variable is and will probably have to consult my collaborators for more information at some point.

Often it's easiest to add another worksheet tab into the excel file if you're using that, titled 'meta-data', rather than listing meta-data in a separate document that could easily become detached from the database.

Here, then, is a basic outline of what meta-data should always be included. 

1. A few lines describing who created the database and when, who it was received from if it was emailed to you, and what it is about.

2. List of variable names as they appear in the raw database; a brief description of the variable; the coding or units; the type of variable (cat, cont, binary, percentage, etc); if it's a response or explanatory variable; the shortened name used for the variable in any R scripts.

I often find that this sort of meta-data is useful when writing a paper and rarely a waste of time to compile, because you almost have a ready-to-go table that could be added to the paper or supplemental information about your data. When I read papers with lots of variables I find such tables very helpful.

Growing your own creativity

Since starting to seriously use Twitter a couple of weeks ago to follow like-minded researchers and tweet my own interests,  I've begun to feel better informed about ecology and also more creative. These outcomes alone are excellent reasons to be using social media and I'm glad I've taken that step.

Perhaps I shouldn't admit this - but coming up with creative ideas for research is something that I've found challenging. Useful creativity requires a deep knowledge of your subject area as well as the ideas to advance it in a novel way.

Why has this been a problem for me? When I was younger I enjoyed and even excelled at traditionally "creative" activities such as creative writing, music and art. Yet somehow through university I lost that natural creative feeling. I suspect it was because I felt I should focus all my energy on learning, understanding and memorising in order to get good grades at the expense of taking time out to really think about things.

I've heard people say that your PhD is a time to really think about an area in depth and to enjoy it because you'll never get this chance again...implying perhaps that this is a time to be creative. However, I didn't find that this was how I tackled my PhD. With the pressure of having to complete it in three years due to funding constraints (everyone gets this one), coupled with the (totally correct) expectation that I should publish parts of my thesis before submission, I found myself feeling rather like I was on a treadmill.

I felt I had no time for "extraneous" things, which in my mind then included Twitter, blogging and time out to relax and contemplate. I used to stuff a sandwich into my face at my desk instead of taking a break and would make furtive dashes to Coffee Culture to grab a flat white, which I would get to take-away, because of the need to get back to work.

This probably squashed any creativity in me and it's totally my own doing!

I've decided to try and nurture my creativity and original thought. Even just admitting I need to do this has freed me to think in a less goal-orientated way for at least a little bit of the day.

Watch this space.

Logistic regression in BUGS: model fit and performance using Chi-square and deviance

Having narrowly escaped being an April-fools baby, this week saw me celebrating my birthday (no, I won't say which one!). The great thing about being an Easter birthday is that I can legitimately scoff a lot of chocolate.

Project-wise, I've managed to code a full two path models using my PhD data and test the independence claims in each model, as well as the overall model, following some methods in Clough (2012). This has been quite an achievement for me. The models run in R, though I can't say whether they're totally correct or not yet. This exercise has been useful because I've had to code logistic (Bernoulli), Poisson and regular Normal linear regression models in BUGS. I've also had to tackle some issue around data standardisation prior to modelling.

Another issue I wanted to resolve was how to calculate something akin to Bayesian p-values, or a measure of the fit, for a BUGS logistic regression model. Here is some code provided by Richard Duncan that does just this using a Chi-square goodness of fit tests and deviance measure. However, I haven't managed to get this working on my own data yet.

# load packages

  library(rjags)

  library(R2jags)

  library(mcmcplots)

  load.module("glm")          # loads the glm module for JAGS, which

                              # may help with convergence

  

# set working directory

setwd("c:\\pgrad\\kirsty mcgregor")

# generate some data

  unlogit <- function(x) exp(x)/(1+exp(x))

  

# an explanatory variable

  x <- seq(1, 10, 0.1)

# probability as a function of x on the logit scale  

  p <- unlogit(-1 + 0.3*x)

  plot(p ~ x)

  

# now generate some bernoulli data using this probability info

  y <- rbinom(length(p), 1, p)

  y

  

# logistic regression model in R

  summary(m1 <- glm(y ~ x, family=binomial))

  

  N <- length(y)

  

# in JAGS

  mod <- "model

  {

  for(i in 1:N) {

    y[i] ~ dbern(p[i])

    logit(p[i]) <- b0 + b1*x[i]

# calculate goodness of fit statistic for logistic regression model using data

    predicted[i] <- p[i]

    res.y[i] <- ((y[i] - predicted[i]) / sqrt(predicted[i]*(1-predicted[i])))^2

# calculate goodness of fit statistic for logistic regression model using new predicted data

    y.rep[i] ~ dbern(p[i])

    res.y.rep[i] <- ((y.rep[i] - predicted[i]) / sqrt(predicted[i]*(1-predicted[i])))^2

  }

  fit <- sum(res.y[])           # test statistic for data   

  fit.rep <- sum(res.y.rep[])   # test statistic for new predicted data   

  test <- step(fit.rep - fit)   # Test whether new data set more extreme

  bpvalue <- mean(test)   # Bayesian p-value 

  #priors

  b0 ~ dnorm(0, 0.0001)

  b1 ~ dnorm(0, 0.0001)

}"

  # write model

  write(mod, "modelRD.txt")

  

  set.seed(rnbinom(1, mu=200, size=1))

  mod <- jags(model = "modelRD.txt",

              data = list(N=N, y=y, x=x),

              inits = function() list(b0=rnorm(1), b1=rnorm(1)), 

              param = c("b0", "b1", "bpvalue"),

              n.chains = 3,

              n.iter =20000,

              n.burnin = 10000)

  

# put output into mcmc form and plot chains  

  out <- as.mcmc(mod)

  mcmcplot(out, parms=c("b0", "b1"))

# save and view output  

  all.sum <- mod$BUGSoutput$summary

  all.sum

Reading this week:

Clough, Y. 2012. A generalized approach to modeling and estimating indirect effects in ecology. Ecology 93:1809–1815. http://dx.doi.org/10.1890/11-1899.1

Imai, K., Keele, L., and Yamamoto, T. 2010. Identification, inference and sensitivity analysis for causal mediation effects. Statistical Science 25: 51–71. http://dx.doi.org/10.1214/10-STS321

Morin, L., Paini, D.R., Randall, R.P. 2013. Can global weed assemblages be used to predict future weeds? PLoS ONE 8(2): e55547.  http://dx.doi.org/10.1371/journal.pone.0055547

Bechara, F.C., Reis, A., Bourscheid, K., Vieira, N.K., and Trentin, B.E. 2013. Reproductive biology and early establishment of Pinus elliottii var. elliottii in Brazilian sandy coastal plain vegetation: implications for biological invasion. Scientia Agricola 70: 88-92. Available at: http://www.scielo.br/pdf/sa/v70n2/05.pdf

First steps into Bayesian method for confirmatory path analysis: changing WinBUGS implementation to Jags

After a reasonably productive week where my understanding of some issues has improved a lot, I find I've got that Friday feeling...I'm thinking slowly and can't figure out simple problems (such as how to do a simple Poisson regression model in BUGS with my pine data). 

Here, then, is a small part of my work this week: converting model and code intended for WinBUGS into something that'll run in Jags. And hopefully, getting the same answer. I'm at the very start of really understanding how this method for d-sep confirmatory path analysis works, so, baby steps!

Required reading: 

Shipley, B. 2009. Confirmatory path analysis in a generalized multilevel context. Ecology 90:363–368. http://dx.doi.org/10.1890/08-1034.1

Clough, Y. 2012. A generalized approach to modeling and estimating indirect effects in ecology. Ecology 93:1809–1815. http://dx.doi.org/10.1890/11-1899.1

Some of the following code and all the data to run the example is available online in the supplemental information in Clough (2012). I have modified the WinBUGS section to include R scripts that specify and write the model file (although the appropriate model file is also supplied in the supplemental information of the article).

WinBUGS version:


# Load library
  library(R2WinBUGS)

#Set working directory
  setwd("C:\\Users\\kirsty.mcgregor\\wd")

# Data preparation
data.plot=read.table("data.plot.txt",h=T,dec=".")
data.subplot=read.table("data.subplot.txt",h=T,dec=".")
data.tree=read.table("data.tree.txt",h=T,dec=".")

# Defining tree-level variable vectors
Pc=data.tree$Pc
Hs=data.tree$Hs
Cc=data.tree$Cc
Npodh=data.tree$npodh
subplot=data.tree$subplot

# Defining subplot-level variable vectors
Nfert=data.subplot$Nfert
plot=data.subplot$plot

# Defining plot-level variable vectors
temp=data.plot$temp
age=data.plot$age

# Number of observations at each level
n.tree=430
n.subplot=86
n.plot=43

# Specify the model in BUGS language
mod1 <- " model{

for (i in 1:n.plot){

temp[i] ~ dnorm(temp.hat[i],tau)
temp.hat[i] <- a.temp+b.age_temp*age[i]

}

tau <- 1/sigma2
sigma2 <- pow(sigma,2)
sigma ~ dunif(0, 100)
a.temp ~ dnorm(0,1.0E-6)
b.age_temp ~ dnorm(0,1.0E-6)
}"

# Write model to working directory
  write(mod1, "claim1_bugs.txt")



# Model to test independence claim 1
m1.data <-list("n.plot","temp","age")

# Defining the initial values for the MCMC chains
m1.inits <- function (){list(
a.temp=runif(1,-1,1),
b.age_temp=runif(1,-1,1),
sigma=runif(1,1,2))
}

# Defining the parameters for which I want to save the posterior samples
m1.parameters <- c(names(m1.inits()))


m1 <- bugs (m1.data, m1.inits, m1.parameters,  n.chains=3, 
 "claim1_bugs.txt", bugs.directory="c:/Programme/WinBUGS14",
 working.directory=getwd(), clearWD=FALSE, n.iter=50000, n.sims=1000,codaPkg=F,debug=TRUE)

print(m1)

## Define the function to get the largest CI that does not include zero

fun.CI=function(x){
w=sum(x<0)
w2=(1002-(abs(501-w)*2))/1002
return(w2)
}

# Apply the function to the posterior distribution of b.age_temp
  fun.CI(m1$sims.list$b.age_temp)

# This shows the independence claim 1 to be substantiated.


And here is the version I wrote in Jags. Note that the specification of sigma is now sigma2 <-sigma^2:


# Load libraries
  library(rjags)
  library(R2jags)
  load.module("glm")          # loads the glm module for JAGS 
                              # which may help with convergence
  library(mcmcplots)

# Specify the model in BUGS language
mod1 <- " model{

for (i in 1:n.plot){

temp[i]~dnorm(temp.hat[i],tau)
temp.hat[i]<-a.temp+b.age_temp*age[i]
}

tau <-1/sigma2
sigma2 <-sigma^2
sigma ~ dunif(0, 100)
a.temp~dnorm(0,1.0E-6)
b.age_temp~dnorm(0,1.0E-6)

}"

# Write model to working directory
  write(mod1, "claim1_bugs.txt")

# Model to test independence claim 1
m1.inits <- function (){list(
a.temp=runif(1,-1,1),
b.age_temp=runif(1,-1,1),
sigma=runif(1,1,2))
}

# Run model in jags
  set.seed(rnbinom(1, mu=200, size=1))
  mod1 <- jags(model = "claim1_bugs.txt",
              data = list("n.plot","temp","age"),
              inits = m1.inits, 
              param = c(names(m1.inits())),
              n.chains = 3,
              n.iter =50000,
              n.burnin = 25000,
 n.thin=75)
  mod1

# Put output into mcmc form and plot chains and view
  out <- as.mcmc(mod1)
  mcmcplot(out, parms=c(names(m1.inits())))

# Define the function to get the largest CI that does not include
# zero
  fun.CI=function(x){
  w=sum(x<0)
  w2=(1002-(abs(501-w)*2))/1002
  return(w2)
  }

# Summary 
  all.sum <- mod1$BUGSoutput$summary
  all.sum

# Apply the function to the posterior distribution of b.age_temp
  fun.CI(mod1$BUGSoutput$sims.list$b.age_temp)

It seems to work for me and produce a similar result...but is it right?!

R2jags and rjags for Bayesian modelling in R

Directly after handing in my PhD I did what most people probably do: have a little rest from data analysis. The 'rest' continued into the first few months of my postdoc, because my main focus was on data collection to expand the database I'll be working on.

Consequently, the world of Bayesian modelling using R and BUGS has moved on without me. Or, perhaps I was always behind the times! Either way, I now find myself free from data collection and getting back into modelling. 

Previously, I was using BRugs to run models in OpenBUGS from R. My PhD supervisor, Richard Duncan, encouraged me to switch over to using R2jags as the R package of choice. The other package I've looked at is rjags. Switching my code over was relatively straightforward and R2jags seems to run models more quickly than through BRugs. If anyone is using R2WinBUGS or BRugs as the R interface for BUGS, I would encourage you to try switching to a jags package.

The other good package I've discovered is mcmcplots, which has a function that plots the model output in an HTML format using your web browser to display it. This is more convenient than dealing with plots in R itself. 

PhD thesis deposited!