Christopher Gandrud (간드루드 크리스토파)

A Link Between topicmodels LDA and LDAvis

2015-05-08T07:41:00.000+01:00

Carson Sievert and Kenny Shirley have put together the really nice LDAvis R package. It provides a Shiny-based interactive interface for exploring the output from Latent Dirichlet Allocation topic models. If you've never used it, I highly recommend checking out their XKCD example (this paper also has some nice background).

LDAvis doesn't fit topic models, it just visualises the output. As such it is agnostic about what package you use to fit your LDA topic model. They have a useful example of how to use output from the lda package.

I wanted to use LDAvis with output from the topicmodels package. It works really nicely with texts preprocessed using the tm package. The trick is extracting the information LDAvis requires from the model and placing it into a specifically structured JSON formatted object.

To make the conversion from topicmodels output to LDAvis JSON input easier, I created a linking function called topicmodels_json_ldavis. The full function is below. To use it follow these steps:

Create a VCorpus object using the tm package's Corpus function.
Convert this to a document term matrix using DocumentTermMatrix, also from tm.
Run your model using topicmodel's LDA function.
Convert the output into JSON format using topicmodels_json_ldavis. The function requires the output from steps 1-3.
Visualise with LDAvis' serVis.

Set up R/Stan on Amazon EC2

2014-12-10T15:12:00.000+00:00

A few months ago I posted the script that I use to set up my R/JAGS working environment on an Amazon EC2 instance.

Since then I've largely transitioned to using R/Stan to estimate my models. So, I've updated my setup script (see below).

There are a few other changes:

I don't install/use RStudio on Amazon EC2. Instead, I just use R from the terminal. Don't get me wrong, I love RStudio. But since what I'm doing on EC2 is just running simulations (I handle the results on my local machine), RStudio is overkill.
I don't install git anymore. Instead I use source_url (from devtools) and source_data (from repmis) to source scripts from GitHub. Again all of the manipulation I'm doing to these scripts is on my local machine.

Do Political Scientists Care About Effect Sizes: Replication and Type M Errors

2014-10-13T07:53:00.000+01:00

Reproducibility has come a long way in political science. Many major journals now require replication materials be made available either on their websites or some service such as the Dataverse Network. Most of the top journals in political science have formally committed to reproducible research best practices by signing up to the The (DA-RT) Data Access and Research Transparency Joint Statement.

This is certainly progress. But what are political scientists actually supposed to do with this new information? Data and code availability does help avoid effort duplication--researchers don't need to gather data or program statistical procedures that have already been gathered or programmed. It promotes better research habits. It definitely provides ''procedural oversight''. We would be highly suspect of results from authors that were unable or unwilling to produce their code/data.

However, there are lots of problems that data/code availability requirements do not address. Apart from a few journals like Political Science Research and Methods, most journals have no standing policy to check the replication materials' veracity. Reviewers rarely have access to manuscripts' code/data. Even if they did have access to it, few reviewers would be willing or able to undertake the time consuming task of reviewing this material.

Do political science journals care about coding and data errors?

What do we do if someone replicating published research finds clear data or coding errors that have biased the published estimates?

Note that I'm limiting the discussion here to honest mistakes, not active attempts to deceive. We all make these mistakes. To keep it simple, I'm also only talking about clear, knowable, and non-causal coding and data errors.

Probably the most responsible action a journal could take when clear cut coding/data biased results have been found would be to directly adjoin to the original article a note detailing the bias. This way readers will always be aware of the correction and will have the best information possible. This is a more efficient way of getting out corrected information than relying on some probabilistic process where readers may or may not stumble upon the information posted elsewhere.

As far as I know, however, no political science journal has a written procedure (please correct me if I'm wrong) for dealing with this new information. My sense is that there are a series of ad hoc responses that closely correspond to how the bias affects the results:

Statistical significance

The situation where a journal is most likely to do anything is when correcting the bias makes the results no longer statistically significant. This might get a journal to append a note to the original article. But maybe not, they could just ignore it.

Sign

It might be that once the coding/data bias is corrected, the sign of an estimated effect flips--the result of what Andrew Gelman calls Type S errors. I really have no idea what a journal would do in this situation. They might append a note or maybe not.

Magnitude

Perhaps the most likely outcome of correcting honest coding/data bias is that the effect size changes. These errors would be the result of Gelman's Type M errors. My sense (and experience) is that in a context where novelty is greatly privileged over facts journal editors will almost certainly ignore this new information. It will be buried.

Do political scientists care about effect size?

Due to the complexity of what political scientists study, we rarely (perhaps with the exception of specific topics like election forecasting) think that we are very close to estimating a given effect's real magnitude. Most researchers are aiming for statistical significance and a sign that matches their theory.

Does this mean that we don't care about trying to estimate magnitudes as closely as possible?

Looking at political science practice pre-publication, there is a lot of evidence that we do care about Type M errors. Considerable effort is given to finding new estimation methods that produce less biased results. Questions of omitted variable bias are very common at research seminars and in journal reviews. Most researchers do carefully build their data sets and code to minimise coding/data bias. Many of these efforts are focused on the headline stuff--whether or not a given effect is significant and what the direction of the effect is. But, these efforts are also part of a desire to make the most accurate estimate of an effect as possible.

However, the review process and journals' responses to finding Type M errors caused by honest coding/data errors in published findings suggest that perhaps we don't care about effect size. Reviewers almost never look at code and data. Journals (as far as I know, please correct me if I'm wrong) never append information on replications that find Type M errors to original papers.

Prescription

I have a simple prescription for demonstrating that we actually care about estimating accurate effect sizes:

Develop a standard practice of including a short authored write up of the data/code bias with corrected results in the original article's supplementary materials. Append a notice to the article pointing to this.

Doing this would not only give readers more accurate effect size estimates, but also make replication materials more useful.

Standardising the practice of publishing authored notes will incentivise people to use replication materials, find errors, and publicly correct them.

Note to self: brew cleanup r

2014-07-22T10:51:00.000+01:00

Note to self: after updating R with Homebrew remember to cleanup old versions:

brew cleanup r

Otherwise I'm liable to get a segfault. (see also)

Simple script from setting up R, Git, and Jags on Amazon EC2 Ubuntu Instance

2014-06-28T15:38:00.002+01:00

Just wanted to put up the script I've been using to create an Amazon EC2 Ubuntu instance for running RStudio, Git, and Jags. There isn't anything really new in here, but it it has been serving me well.

The script begins after the basic instance has been set up in the Amazon EC2 console (yhat has a nice post on how to do this, though some of their screenshots are a little old). Just SSH into the instance and get started.

Updates to repmis: caching downloaded data and Excel data downloading

2014-05-11T16:51:00.000+01:00

Over the past few months I’ve added a few improvements to the repmis–miscellaneous functions for reproducible research–R package. I just want to briefly highlight two of them:

Caching downloaded data sets.
source_XlsxData for downloading data in Excel formatted files.

Both of these capabilities are in repmis version 0.2.9 and greater.

Caching

When working with data sourced directly from the internet, it can be time consuming (and make the data hoster angry) to repeatedly download the data. So, repmis’s source functions (source_data, source_DropboxData, and source_XlsxData) can now cache a downloaded data set by setting the argument cache = TRUE. For example:

DisData <- source_data("http://bit.ly/156oQ7a", cache = TRUE)

When the function is run again, the data set at http://bit.ly/156oQ7a will be loaded locally, rather than downloaded.

To delete the cached data set, simply run the function again with the argument clearCache = TRUE.

`source_XlsxData`

I recently added the source_XlsxData function to download Excel data sets directly into R. This function works very similarly to the other source functions. There are two differences:

You need to specify the sheet argument. This is either the name of one specific sheet in the downloaded Excel workbook or its number (e.g. the first sheet in the workbook would be sheet = 1).
You can pass other arguments to the read.xlsx function from the xlsx package.

Here’s a simple example:

RRurl <- 'http://www.carmenreinhart.com/user_uploads/data/22_data.xls'

RRData <- source_XlsxData(url = RRurl, sheet = 2, startRow = 5)

startRow = 5 basically drops the first 4 rows of the sheet.

d3Network Plays Nice with Shiny Web Apps

2014-05-09T15:43:00.000+01:00

After some delay (and because of helpful prompting by Giles Heywood and code contributions by John Harrison) d3Network now plays nicely with Shiny web apps. This means you can fully integrate R/D3.js network graphs into your web apps.

Here is what one simple example looks like:

An explanation of the code is here and you can download the app and play with it using:

shiny::runGitHub('d3ShinyExample', 'christophergandrud')

Programmatically download political science data with the psData package

2014-02-23T16:41:00.000+00:00

A lot of progress has been made on improving political scientists’ ability to access data ‘programmatically’, e.g. data can be downloaded with source code R. Packages such as WDI for World Bank Development Indicator and dvn for many data sets stored on the Dataverse Network make it much easier for political scientists to use this data as part of a highly integrated and reproducible workflow.

There are nonetheless still many commonly used political science data sets that aren’t easily accessible to researchers. Recently, I’ve been using the Database of Political Institutions (DPI), Polity IV democracy indicators, and Reinhart and Rogoff’s (2010) financial crisis occurrence data. All three of these data sets are freely available for download online. However, getting them, cleaning them up, and merging them together is kind of a pain. This is especially true for the Reinhart and Rogoff data, which is in 4 Excel files with over 70 individual sheets, one for each country’s data.

Also, I’ve been using variables that are combinations and/or transformations of indicators in regularly updated data sets, but which themselves aren’t regularly updated. In particular, Bueno de Mesquita et al. (2003) devised two variables that they called the ‘winset’ and the ‘selectorate’. These are basically specific combinations of data in DPI and Polity IV. However, the winset and selectorate variables haven’t been updated alongside the yearly updates of DPI and Polity IV.

There are two big problems here:

A lot of time is wasted by political scientists (and their RAs) downloading, cleaning, and transforming these data sets for their own research.
There are many opportunities while doing this work to introduce errors. Imagine the errors that might be introduced and go unnoticed if a copy-and-paste approach is used to merge the 70 Reinhart and Rogoff Excel sheets.

As a solution, I’ve been working on a new R package called psData. This package includes functions that automate the gathering, cleaning, and creation of common political science data and variables. So far (February 2014) it gathers DPI, Polity IV, and Reinhart and Rogoff data, as well as creates winset and selectorate variables. Hopefully the package will save political scientists a lot of time and reduce the number of data management errors.

There certainly could be errors in the way psData gathers data. However, once spotted the errors could be easily reported on the package’s Issues Page. Once fixed, the correction will be spread to all users via a package update.

Types of functions

There are two basic types of functions in psData: Getters and Variable Builders. Getter functions automate the gathering and cleaning of particular data sets so that they can easily be merged with other data. They do not transform the underlying data. Variable Builders use Getters to gather data and then transform it into new variables suggested by the political science literature.

Examples

To download only the polity2 variable from Polity IV:

# Load package
library(psData)

# Download polity2 variable
PolityData <- PolityGet(vars = "polity2")

# Show data
head(PolityData)


##   iso2c     country year polity2
## 1    AF Afghanistan 1800      -6
## 2    AF Afghanistan 1801      -6
## 3    AF Afghanistan 1802      -6
## 4    AF Afghanistan 1803      -6
## 5    AF Afghanistan 1804      -6
## 6    AF Afghanistan 1805      -6

Note that the iso2c variable refers to the ISO two letter country code country ID. This standardised country identifier could be used to easily merge the Polity IV data with another data set. Another country ID type can be selected with the OutCountryID argument. See the package documentation for details.

To create winset (W) and selectorate (ModS) data use the following code:

WinData <- WinsetCreator()

head(WinData)


##    iso2c     country year    W ModS
## 1     AF Afghanistan 1975 0.25    0
## 2     AF Afghanistan 1976 0.25    0
## 3     AF Afghanistan 1977 0.25    0
## 15    AF Afghanistan 1989 0.50    0
## 16    AF Afghanistan 1990 0.50    0
## 17    AF Afghanistan 1991 0.50    0

Install

psData should be on CRAN soon, but while it is in the development stage you can install it with the devtools package:

devtools::install_github('psData', 'christophergandrud')

Suggestions

Please feel free to suggest other data set downloading and variable creating functions. To do this just leave a note on the package’s Issues page or make a pull request with a new function added.

Three Quick and Simple Data Cleaning Helper Functions (December 2013)

2013-12-06T17:43:00.000+00:00

As I go about cleaning and merging data sets with R I often end up creating and using simple functions over and over. When this happens, I stick them in the DataCombine package. This makes it easier for me to remember how to do an operation and others can possibly benefit from simplified and (hopefully) more intuitive code.

I've talked about some of the commands in DataCombine in previous posts. In this post I'll give examples for a few more that I've added over the past couple of months. Note: these examples are based on DataCombine version 0.1.11.

Here is a brief run down of the functions covered in this post:

FindReplace: a function to replace multiple patterns found in a character string column of a data frame.
MoveFront: moves variables to the front of a data frame. This can be useful if you have a data frame with many variables and want to move a variable or variables to the front.
rmExcept: removes all objects from a work space except those specified by the user.

`FindReplace`

Recently I needed to replace many patterns in a column of strings. Here is a short example. Imagine we have a data frame like this:

ABData <- data.frame(a = c("London, UK", "Oxford, UK", "Berlin, DE", "Hamburg, DE", "Oslo, NO"), b = c(8, 0.1, 3, 2, 1))

Ok, now I want to replace the UK and DE parts of the strings with England and Germany. So I create a data frame with two columns. The first records the pattern and the second records what I want to replace the pattern with:

Replaces <- data.frame(from = c("UK", "DE"), to = c("England", "Germany"))

Now I can just use FindReplace to make the replacements all at once:

library(DataCombine)

ABNewDF <- FindReplace(data = ABData, Var = "a", replaceData = Replaces, from = "from", to = "to", exact = FALSE)

# Show changes
ABNewDF

##                  a   b
## 1  London, England 8.0
## 2  Oxford, England 0.1
## 3  Berlin, Germany 3.0
## 4 Hamburg, Germany 2.0
## 5         Oslo, NO 1.0

If you set exact = TRUE then FindReplace will only replace exact pattern matches. Also, you can set vector = TRUE to return only a vector of the column you replaced (the Var column), rather than the whole data frame.

`MoveFront`

On occasion I've wanted to move a few variables to the front of a data frame. The MoveFront function makes this pretty simple. It only has two arguments: data and Var. Data is the data frame and Var is a character vector with the columns I want to move to the front of the data frame in the order that I want them. Here is an example:

# Create dummy data
A <- B <- C <- 1:50
OldOrder <- data.frame(A, B, C)

names(OldOrder)

## [1] "A" "B" "C"

# Move B and A to the front
NewOrder2 <- MoveFront(OldOrder, c("B", "A"))
names(NewOrder2)

## [1] "B" "A" "C"

`rmExcept`

Finally, sometimes I want to clean up my work space and only keep specific objects. I want to remove everything else. This is straightforward with rmExcept. For example:

# Create objects
A <- 1
B <- 2
C <- 3

# Remove all objects except for A
rmExcept("A")

## Removed the following objects:
## ABData, ABNewDF, B, C, NewOrder2, OldOrder, Replaces

# Show workspace
ls()

## [1] "A"

You can set the environment you want to clean up with the envir argument. By default is is your global environment.

Showing results from Cox Proportional Hazard Models in R with simPH

2013-09-02T14:27:00.000+01:00

Update 2 February 2014: A new version of simPH (Version 1.0) will soon be available for download from CRAN. It allows you to plot using points, ribbons, and (new) lines. See the updated package description paper for examples. Note that the ribbons argument will no longer work as in the examples below. Please use type = 'ribbons' (or 'points' or 'lines').

Effectively showing estimates and uncertainty from Cox Proportional Hazard (PH) models, especially for interactive and non-linear effects, can be challenging with currently available software. So, researchers often just simply display a results table. These are pretty useless for Cox PH models. It is difficult to decipher a simple linear variable’s estimated effect and basically impossible to understand time interactions, interactions between variables, and nonlinear effects without the reader further calculating quantities of interest for a variety of fitted values.

So, I’ve been putting together the simPH R package to hopefully make it easier to show results from Cox PH models. The package uses plots of post-estimation simulations (the same idea behind the plotting facilities in the Zelig package) to show Cox PH estimates and the uncertainty surrounding them.

Here I want to give an overview of how to use simPH. First the general process, then a few examples. Have a look at this paper for details about the math behind the graphs.

General Steps

There are three steps to use simPH:

Estimate a Cox PH model in the usual way with the coxph command in the survival package.
Simulate quantities of interest--hazard ratios, first differences, marginal effect, relative hazards, or hazard rates--with the appropriate simPH simulation command.
Plot the simulations with the simGG method.

A Few Examples

Here are some basic examples that illustrate the process and key syntax. Before getting started you can install simPH in the normal R way from CRAN.

Linear effects

Let’s look at a simple example with a linear non-interactive effect. The data set I’m using is from Carpenter (2002). It is included with simPH. See the package documentation for details.

First, let’s estimate a Cox PH model where the event of interest is a drug receiving FDA approval.

# Load packages
library(survival)
library(simPH)

# Load Carpenter (2002) data
data("CarpenterFdaData")

# Estimate survival model
 M1 <- coxph(Surv(acttime, censor) ~ prevgenx + lethal +
            deathrt1 + acutediz + hosp01  + hhosleng +
            mandiz01 + femdiz01 + peddiz01 + orphdum +
            vandavg3 + wpnoavg3 + condavg3 + orderent +
            stafcder, data = CarpenterFdaData)

Now say we want to examine the effect that the number of FDA staff reviewing a drug application has on it being accepted. This variable is called stafcder in the model we just estimated. To do this let’s use simPH to simulate hazard ratios. We will simulate hazard ratios for units \(j\) and \(l\), i.e. \(\frac{h_{j}(t)}{h_{l}(t)}\) using simPH’s coxsimLinear command, because we estimate the effect of the number of staff as linear. In the following code we use the Xj argument to set the \(j\) values. We could use Xl also, but as we don’t coxsimLinear assumes all \(x_{l}\) are 0.

Notice that the argument spin = TRUE. This finds the shortest 95% probability interval of the simulations using the SPIn method. SPIn can be especially useful for showing simulated quantities of interest generated from Cox PH models, because then can often be crowded close to a lower boundary (0 in the case of hazard rates). We should be more interested in the area with the highest probability–most simulated values–rather than an arbitrary central interval. That being said, if spin = FALSE, then we will simply find the central 95% interval of the simulations.

# Simulate and plot Hazard Ratios for stafcder variable
Sim1 <- coxsimLinear(M1, b = "stafcder", 
                    qi = "Hazard Ratio", ci = 0.95,
                    Xj = seq(1237, 1600, by = 2), spin = TRUE)
# Plot
simGG(Sim1)

Notice in the plot that each simulation is plotted as an individual point. These are all of the simulations in the shortest 95% probability interval. Each point has a bit of transparency (they are 90% transparent by default). So the plot is visually weighted; the darker areas of the graph have a higher concentration of the simulations. This gives us a very clear picture of the simulation distribution, i.e. the estimated effect and our uncertainty about it.

If you don’t want to plot every point, you can simply use ribbons showing the constricted 95% and the middle 50% of this interval. To do this simply use the ribbons = TRUE argument.

simGG(Sim1, ribbons = TRUE, alpha = 0.5)

Notice the alpha = 0.5 argument. This increased the transparency of the widest ribbon to 50%.

The syntax we’ve used here is very similar to what we use when we are working with nonlinear effects estimated with polynomials and splines. Post-estimation simulations can be run with the coxsimPoly and coxsimSpline commands. See the simPH documentation for more examples.

Interactions

The syntax for two-way interactions simulated with the coxsimInteract command is a little different. Using the same data, let’s look at how to show results for interactions. In the following model we are interacting two variables lethal and prevgenx. We can think of these variables as X1 and X2, respectively. For interactions it can be useful to examine the marginal effect. To find the marginal effect of a one unit increase in the lethal variable given various values of prevgenx let’s use the following code:

# Estimate the model
M2 <- coxph(Surv(acttime, censor) ~ lethal*prevgenx, data = CarpenterFdaData)
 
# Simulate Marginal Effect of lethal for multiple values of prevgenx
Sim2 <- coxsimInteract(M2, b1 = "lethal", b2 = "prevgenx",
                       qi = "Marginal Effect",
         X2 = seq(2, 115, by = 2), nsim = 1000)

# Plot the results
simGG(Sim2, ribbons = TRUE, alpha = 0.5, xlab = "\nprevgenx",
      ylab = "Marginal effect of lethal\n")

The order of the X1 and X2 variables in the interactions matters. The marginal effect is always calculated for the X1 variable over a range of X2 values.

Notice also that we set the x and y axis labels with the xlab and ylab arguments.

Time-varying Effects

Finally, let’s look at how to use with the coxsimtvc command to show results from effects that we estimate to vary over time. Here we are going to use another data set that is also included with simPH. The event of interest in the following model is when deliberation is stopped on a European Union directive (the model is from Licht (2011)). We will create hazard ratios for the effect that the number of backlogged items (backlog) has on deliberation time. We will estimate the effect as a log-time interaction.

# Load Golub & Steunenberg (2007) data. The data is included with simPH.
data("GolubEUPData")

# Create natural log-time interactions
Golubtvc <- function(x){
  tvc(data = GolubEUPData, b = x, tvar = "end", tfun = "log")
}
GolubEUPData$Lcoop <- Golubtvc("coop")
GolubEUPData$Lqmv <- Golubtvc("qmv")
GolubEUPData$Lbacklog <- Golubtvc("backlog")
GolubEUPData$Lcodec <- Golubtvc("codec")
GolubEUPData$Lqmvpostsea <- Golubtvc("qmvpostsea")
GolubEUPData$Lthatcher <- Golubtvc("thatcher")

# Estimate model
M1 <- coxph(Surv(begin, end, event) ~ qmv + qmvpostsea + qmvpostteu +
              coop + codec + eu9 + eu10 + eu12 + eu15 + thatcher + 
              agenda + backlog + Lqmv + Lqmvpostsea + Lcoop + Lcodec +
              Lthatcher + Lbacklog,
            data = GolubEUPData, ties = "efron")

Much of the first half of the code is dedicated to creating the log-time interactions with the tvc command.

Now we simply create the simulations for a range of values of backlog and plot them. Note in the following code that we tell coxsimtvc the name of both the base backlog variable and its log-time interaction term Lbacklog using the btvc argument. We also need to tell coxsimtvc that it is a log-time interaction with tfun = "log".

# Create simtvc object for relative hazard
Sim2 <- coxsimtvc(obj = M1, b = "backlog", btvc = "Lbacklog",
                  qi = "Relative Hazard", Xj = seq(40, 200, 40),
                  tfun = "log", from = 1200, to = 2000, by = 10,
                  nsim = 500)

# Create relative hazard plot
simGG(Sim2, xlab = "\nTime in Days", ribbons = TRUE,
        leg.name = "Backlogged \n Items", alpha = 0.2)

simGG Plots are ggplot2

Finally, because almost every plot created by simGG is a ggplot2 plot, you can use almost the full range of customisation options that that package allows. See the ggplot2 documentation for many examples.

GitHub renders CSV in the browser, becomes even better for social data set creation

2013-08-23T04:42:00.000+01:00

I've written in a number of places about how GitHub can be a great place to store data. Unlike basically all other web data storage sites (many of which I really like such as Dataverse and FigShare) GitHub enables deep social data set development and fits nicely into a reproducible research workflow with R.

One negative though, especially compared to FigShare, was that there was no easy way to view CSV or TSV data files in the browser. Unless you downloaded the data and opened it in Excel or an R viewer or whatever, you had to look at the raw data file in the browser. It's basically impossible to make sense of a data set of any size like this.

However, from at least today, GitHub now renders the data set in the browser as you would expect. Take a look at their blog post for the details.

Getting Started with Reproducible Research: A chapter from my new book

2013-07-16T00:35:00.000+01:00

This is an abridged excerpt from Chapter 2 of my new book Reproducible Research with R and RStudio. It's published by Chapman & Hall/CRC Press.

You can purchase it on Amazon. "Search inside this book" includes a complete table of contents.

Researchers often start thinking about making their work reproducible near the end of the research process when they write up their results or maybe even later when a journal requires their data and code be made available for publication. Or maybe even later when another researcher asks if they can use the data from a published article to reproduce the findings. By then there may be numerous versions of the data set and records of the analyses stored across multiple folders on the researcher’s computers. It can be difficult and time consuming to sift through these files to create an accurate account of how the results were reached. Waiting until near the end of the research process to start thinking about reproducibility can lead to incomplete documentation that does not give an accurate account of how findings were made. Focusing on reproducibility from the beginning of the process and continuing to follow a few simple guidelines throughout your research can help you avoid these problems. Remember "reproducibility is not an afterthought–it is something that must be built into the project from the beginning" (Donoho, 2010, 386).

This chapter first gives you a brief overview of the reproducible research process: a workflow for reproducible research. Then it covers some of the key guidelines that can help make your research more reproducible.

The Big Picture: A Workflow for Reproducible Research

The three basic stages of a typical computational empirical research project are:

data gathering,
data analysis,
results presentation.

Instead of starting to use the individual tools of reproducible research as soon as you learn them I recommend briefly stepping back and considering how the stages of reproducible research tie together overall. This will make your workflow more coherent from the beginning and save you a lot of backtracking later on. The figure below illustrates the workflow. Notice that most of the arrows connecting the workflow’s parts point in both directions, indicating that you should always be thinking how to make it easier to go backwards through your research, i.e. reproduce it, as well as forwards.

Example Workflow & A Selection of Commands to Tie it Together

Around the edges of the figure are some of the commands that make it easier to go forwards and backwards through the process. These commands tie your research together. For example, you can use API-based R packages to gather data from the internet. You can use R’s merge command to combine data gathered from different sources into one data set. The getURL from R’s RCurl package and the read.table commands can be used to bring this data set into your statistical analyses. The knitr package then ties your analyses into your presentation documents. This includes the code you used, the figures you created, and, with the help of tools such as the xtable package, tables of results. You can even tie multiple presentation documents together. For example, you can access the same figure for use in a LaTeX article and a Markdown created website with the includegraphics and ![]() commands, respectively. This helps you maintain a consistent presentation of results across multiple document types.

Practical Tips for Reproducible Research

Before learning the details of the reproducible research workflow with R and RStudio, it’s useful to cover a few broad tips that will help you organize your research process and put these skills in perspective. The tips are:

Document everything!,
Everything is a (text) file,
All files should be human readable,
Explicitly tie your files together,
Have a plan to organize, store, and make your files available.

Using these tips will help make your computational research really reproducible.

Document everything!

In order to reproduce your research others must be able to know what you did. You have to tell them what you did by documenting as much of your research process as possible. Ideally, you should tell your readers how you gathered your data, analyzed it, and presented the results. Documenting everything is the key to reproducible research and lies behind all of the other tips here.

Everything is a (text) file

Your documentation is stored in files that include data, analysis code, the write up of results, and explanations of these files (e.g. data set codebooks, session info files, and so on). Ideally, you should use the simplest file format possible to store this information. Usually the simplest file format is the humble, but versatile, text file.

Text files are extremely nimble. They can hold your data in, for example, comma-separated values (.csv) format. They can contain your analysis code in .R files. And they can be the basis for your presentations as markup documents like .tex or .md, for LaTeX and Markdown files, respectively. All of these files can be opened by any program that can read text files.

One reason reproducible research is best stored in text files is that this helps future proof your research. Other file formats, like those used by Microsoft Word (.docx) or Excel (.xlsx), change regularly and may not be compatible with future versions of these programs. Text files, on the other hand, can be opened by a very wide range of currently existing programs and, more likely than not, future ones as well. Even if future researchers do not have R or a LaTeX distribution, they will still be able to open your text files and, aided by frequent comments (see below), be able to understand how we conducted your research (Bowers, 2011, 3).

Text files are also very easy to search and manipulate with a wide range of programs–such as R and RStudio–that can find and replace text characters as well as merge and separate files. Finally, text files are easy to version and changes can be tracked using programs such as Git.

All files should be human readable

Treat all of your research files as if someone who has not worked on the project will, in the future, try to understand them. Computer code is a way of communicating with the computer. It is 'machine readable' in that the computer is able to use it to understand what you want to do. However, there is a very good chance that other people (or you six months in the future) will not understand what you were telling the computer. So, you need to make all of your files 'human readable'. To make them human readable, you should comment on your code with the goal of communicating its design and purpose (Wilson et al., 2012). With this in mind it is a good idea to comment frequently (Bowers, 2011, 3) and format your code using a style guide (Nagler, 1995). For especially important pieces of code you should use literate programming–where the source code and the presentation text describing its design and purpose appear in the same document (Knuth, 1992). Doing this will make it very clear to others how you accomplished a piece of research.

Explicitly tie your files together

If everything is just a text file then research projects can be thought of as individual text files that have a relationship with one another. They are tied together. A data file is used as input for an analysis file. The results of an analysis are shown and discussed in a markup file that is used to create a PDF document. Researchers often do not explicitly document the relationships between files that they used in their research. For example, the results of an analysis–a table or figure–may be copied and pasted into a presentation document. It can be very difficult for future researchers to trace the table or figure back to a particular statistical model and a particular data set without clear documentation. Therefore, it is important to make the links between your files explicit.

Tie commands are the most dynamic way to explicitly link your files together (see the figure above for examples). These commands instruct the computer program you are using to use information from another file.

Have a plan to organize, store, & make your files available

Finally, in order for independent researchers to reproduce your work they need to be able access the files that instruct them how to do this. Files also need to be organized so that independent researchers can figure out how they fit together. So, from the beginning of your research process you should have a plan for organizing your files and a way to make them accessible.

One rule of thumb for organizing your research in files is to limit the amount of content any one file has. Files that contain many different operations can be very difficult to navigate, even if they have detailed comments. For example, it would be very difficult to find any particular operation in a file that contained the code used to gather the data, run all of the statistical models, and create the results figures and tables. If you have a hard time finding things in a file you created, think of the difficulties independent researchers will have!

Because we have so many ways to link files together there is really no need to lump many different operations into one file. So, we can make our files modular. One source code file should be used to complete one or just a few tasks. Breaking your operations into discrete parts will also make it easier for you and others to find errors (Nagler, 1995, 490).

Quick and Simple D3 Network Graphs from R

2013-06-09T03:00:00.000+01:00

Sometimes I just want to quickly make a simple D3 JavaScript directed network graph with data in R. Because D3 network graphs can be manipulated in the browser–i.e. nodes can be moved around and highlighted–they're really nice for data exploration. They're also really nice in HTML presentations. So I put together a bare-bones simple function–called d3SimpleNetwork for turning an R data frame into a D3 network graph.

Arguments

By bare-bones I mean other than the arguments indicating the Data data frame, as well as the Source and Target variables it only has three arguments: height, width, and file.

The data frame you use should have two columns that contain the source and target variables. Here's an example using fake data:

Source <- c("A", "A", "A", "A", "B", "B", "C", "C", "D")
Target <- c("B", "C", "D", "J", "E", "F", "G", "H", "I")
NetworkData <- data.frame(Source, Target)

The height and width arguments obviously set the graph's frame height and width. You can tell file the file name to output the graph to. This will create a standalone webpage. If you leave file as NULL, then the graph will be printed to the console. This can be useful if you are creating a document using knitr Markdown ~~or (similarly) slidify~~. Just set the code chunk results='asis and the graph will be rendered in the document.

Example

Here's a simple example. First load d3SimpleNetwork:

# Load packages to download d3SimpleNetwork
library(digest)
library(devtools)

# Download d3SimpleNetwork
source_gist("5734624")

Now just run the function with the example NetworkData from before:

d3SimpleNetwork(NetworkData, height = 300, width = 700)

Click here for the fully manipulable version. If you click on individual nodes they will change colour and become easier to see. In the future I might add more customisability, but I kind of like the function's current simplicity.

Update 12 June 2013: The original d3SimpleNetwork command discussed here doesn't work easily with slidify. I have created a new d3Network R package that does work well with slidify (and other knitr-created HTML slideshows). Use its d3Network command and set the argument iframe = TRUE.

Slide: one function for lag/lead variables in data frames, including time-series cross-sectional data

2013-05-21T11:18:00.000+01:00

I often want to quickly create a lag or lead variable in an R data frame. Sometimes I also want to create the lag or lead variable for different groups in a data frame, for example, if I want to lag GDP for each country in a data frame.

I've found the various R methods for doing this hard to remember and usually need to look at old blog posts. Any time we find ourselves using the same series of codes over and over, it's probably time to put them into a function.

So, I added a new command–slide–to the DataCombine R package (v0.1.5).

Building on the shift function TszKin Julian posted on his blog, slide allows you to slide a variable up by any time unit to create a lead or down to create a lag. It returns the lag/lead variable to a new column in your data frame. It works with both data that has one observed unit and with time-series cross-sectional data.

Note: your data needs to be in ascending time order with equally spaced time increments. For example 1995, 1996, 1997.

Examples

Not Cross-sectional data

Let's create an example data set with three variables:

# Create time variable
Year <- 1980:1999

# Dummy covariates
A <- B <- 1:20

Data1 <- data.frame(Year, A, B)

head(Data1)

##   Year A B
## 1 1980 1 1
## 2 1981 2 2
## 3 1982 3 3
## 4 1983 4 4
## 5 1984 5 5
## 6 1985 6 6

Now let's lag the A variable by one time unit.

library(DataCombine)

DataSlid1 <- slide(Data1, Var = "A", slideBy = -1)

head(DataSlid1)

##   Year A B A-1
## 1 1980 1 1  NA
## 2 1981 2 2   1
## 3 1982 3 3   2
## 4 1983 4 4   3
## 5 1984 5 5   4
## 6 1985 6 6   5

The lag variable is automatically given the name A-1.

To lag a variable (i.e. the lag value at a given time is the value of the non-lagged variable at a time in the past) set the slideBy argument as a negative number. Lead variables, are created by using positive numbers in slideBy. Lead variables at a given time have the value of the non-lead variable from some time in the future.

Time-series Cross-sectional data

Now let's use slide to create a lead variable with time-series cross-sectional data. First create the example data:

# Create time and unit ID variables
Year <- rep(1980:1983, 5)
ID <- sort(rep(seq(1:5), 4))

# Dummy covariates
A <- B <- 1:20

Data2 <- data.frame(Year, ID, A, B)

head(Data2)

##   Year ID A B
## 1 1980  1 1 1
## 2 1981  1 2 2
## 3 1982  1 3 3
## 4 1983  1 4 4
## 5 1980  2 5 5
## 6 1981  2 6 6

Now let's create a two time unit lead variable based on B for each unit identified by ID:

DataSlid2 <- slide(Data2, Var = "B", GroupVar = "ID",
                    slideBy = 2)

head(DataSlid2)

##   Year ID A B B2
## 1 1980  1 1 1  3
## 2 1981  1 2 2  4
## 3 1982  1 3 3 NA
## 4 1983  1 4 4 NA
## 5 1980  2 5 5  7
## 6 1981  2 6 6  8

Hopefully you'll find slide useful in your own data analysis. Any suggestions for improvement are always welcome.

Reinhart & Rogoff: Everyone makes coding mistakes, we need to make it easy to find them + Graphing uncertainty

2013-04-17T07:36:00.000+01:00

You may have already seen a lot written on the replication of Reinhart & Rogoff’s (R & R) much cited 2010 paper done by Herndon, Ash, and Pollin. If you haven’t, here is a round up of some of some of what has been written: Konczal, Yglesias, Krugman, Cowen, Peng, FT Alphaville.

This is an interesting issue for me because it involves three topics I really like: political economy, reproducibility, and communicating uncertainty. Others have already commented on these topics in detail. I just wanted to add to this discussion by (a) talking about how this event highlights a real need for researchers to use systems that make finding and correcting mistakes easy, (b) incentivising mistake finding/correction rather than penalising it, and (c) showing uncertainty.

Systems for Finding and Correcting Mistakes

One of the problems Herndon, Ash, and Pollin found in R&R’s analysis was and Excel coding error. I love to hate on Excel as much as the next R devotee, but I think that is missing the point. The real lesson is not “don’t use Excel” the real lesson is: we all make mistakes.

(Important point: I refer throughout this post to errors caused by coding mistakes rather than purposeful fabrications and falsifications.)

Coding mistakes are an ever present part of our work. The problem is not that we make coding mistakes. Despite our best efforts we always will. The problem is that we often use tools and practices that make it difficult to find and correct our mistakes.

This is where I can get in some Excel hating: tools and practices that make it difficult to find mistakes include binary files (like Excel’s) that can’t be version controlled in a way that fully reveals the research process, not commenting code, not making your data readily available in formats that make replication easy, not having a system for quickly fixing mistakes when they are found. Sorry R users, but the last three are definitely not exclusive to Excel.

It took Herndon, Ash, and Pollin a considerable amount of time to replicate R & R’s findings and therefore find the Excel error. This seems partially because R & R did not make their analysis files readily available (Herndon, Ash, and Pollin had to ask for them). I’m not sure how this error is going to be corrected and documented. But I imagine it will be like most research corrections: kind of on the fly, mostly emailing and reposting.

How big of a detail is this? There is some debate over how big of a problem this mistake is. Roger Peng ends his really nice post:

The vibe on the Internets seems to be that if only this problem had been identified sooner, the world would be a better place. But my cynical mind says, uh, no. You can toss this incident in the very large bucket of papers with some technical errors that are easily fixed. Thankfully, someone found these errors and fixed them, and that’s a good thing. Science moves on.

I agree with most of this paragraph. But, given how important R & R’s finding was to major policy debates it would have been much better if the mistake was caught sooner rather than later. The tools and practices R & R used made it harder to find and correct the mistake, so policymakers were operating with less accurate information for longer.

Solutions: I’ve written in some detail in the most recent issue of The Political Methodologist about how cloud-based version control systems like GitHub can be used to make finding and correcting mistakes easier. Pull requests, for example, are a really nice way to directly suggest corrections.

Incentivising Error Finding and Correction

Going forward I think it will be interesting to see how this incident shapes researchers’ perceived incentives to make their work easily replicable. Replication is an important part of finding the mistakes that everyone makes. If being found to make a coding mistake (not a fabrication) has a negative impact on your academic career then there are incentives to make finding mistakes difficult, by for example making replication difficult. Most papers do not receive nearly as much attention as R & R’s. So, for most researchers making replication difficult will make it pretty unlikely that anyone will replicate your research and you’ll be home free.

I can't send you my data b/c I think you might find out I made an error. #overlyhonestmethods
— Carlisle Rainey (@carlislerainey) January 9, 2013

This is a perverse incentive indeed.

What can we do? Many journals now require replicable code to accompany published articles. This is a good incentive. Maybe we should go further, and somehow directly incentivise the finding and correction of errors in data sets and analysis code. Ideas could include giving more weight to replication studies at hiring and promotion committees. Maybe even allowing these committees to include information on researchers’ GitHub pull requests that meaningfully improve other’s work by correcting mistakes.

This of course might create future perversion incentives to add errors so that they can then be found. I think this is a bit fanciful. There are surely enough negative social incentives (i.e. embarrassment) surrounding making mistakes to prevent this.

Showing Uncertainty

Roger Peng’s post highlighted the issue of graphing uncertainty, but I just wanted to build it out a little further. The interpretation of the correlation R & R’s found between GDP Growth and Government Debt could have been tempered significantly before any mistakes were found by more directly communicating their original uncertainty. In their original paper, they presented the relationship using bar graphs of average and median GDP growth per grouped debt/GDP level:

Beyond showing the mean and median there is basically no indication of the distribution of the data they are from.

Herndon, Ash, and Pollin put together some nice graphs of these distributions (and avoid that thing economists do of using two vertical axis with two different meanings).

Here is one that gets rid of the groups altogether:

If R & R had shown a simple scatter plot like this (though they did exclude some of the higher GDP Growth country-years at the high debt end, so their's would have looked different), it would have been much more difficult to overly interpret the substantive–policy–value of a correlation between GDP/growth and debt/GDP.

Maybe this wouldn’t have actually changed the policy debate that much, As Mark Blyth argues in his recent book on austerity “facts never disconfirm a good ideology” (p. 18). But at least Paul Krugman might not have had to debate debt/GDP cutoff points on CNBC (for example time point 12:40):

P.S. To R & R’s credit, they do often make their data available. Their data has been useful for at least one of my papers. However, it is often available in a format that is hard to use for cross-country statistical analysis, including, I would imagine, their own. Though I have never found any errors in the data, reporting and implementing corrections to this data would be piecemeal at best.

Dropbox & R Data

2013-04-11T23:05:00.000+01:00

I'm always looking for ways to download data from the internet into R. Though I prefer to host and access plain-text data sets (CSV is my personal favourite) from GitHub (see my short paper on the topic) sometimes it's convenient to get data stored on Dropbox.

There has been a change in the way Dropbox URLs work and I just added some functionality to the repmis R package. So I though that I'ld write a quick post on how to directly download data from Dropbox into R.

The download method is different depending on whether or not your plain-text data is in a Dropbox Public folder or not.

Dropbox Public Folder

Dropbox is trying to do away with its public folders. New users need to actively create a Public folder. Regardless, sometimes you may want to download data from one. It used to be that files in Public folders were accessible through non-secure (http) URLs. It's easy to download these into R, just use the read.table command, where the URL is the file name. Dropbox recently changed Public links to be secure (https) URLs. These cannot be accessed with read.table.

Instead you need can use the source_data command from repmis:

FinURL <-"https://dl.dropbox.com/u/12581470/code/Replicability_code/Fin_Trans_Replication_Journal/Data/public.fin.msm.model.csv"

# Download data
FinRegulatorData <- repmis::source_data(FinURL,
                             sep = ",",
                             header = TRUE)

Non-Public Dropbox Folders

Getting data from a non-Public folder into R was a trickier. When you click on a Dropbox-based file's Share Link button you are taken to a secure URL, but not for the file itself. The Dropbox webpage you're taken to is filled with lots of other Dropbox information. I used to think that accessing a plain-text data file embedded in one of these webpages would require some tricky web scrapping. Luckily, today I ran across this blog post by Kay Cichini.

With some modifications I was able to easily create a function that could download data from non-Public Dropbox folders. The source_DropboxData command is in the most recent version of repmis (v0.2.4) is the result. All you need to know is the name of the file you want to download and its Dropbox key. You can find both of these things in the URL for the webpage that appears when you click on Share Link. Here is an example:

https://www.dropbox.com/s/exh4iobbm2p5p1v/fin_research_note.csv

The file name is at the very end (fin_research_note.csv) and the key is the string of letters and numbers in the middle (exh4iobbm2p5p1v). Now we have all of the information we need for source_DropboxData:

FinDataFull <- repmis::source_DropboxData("fin_research_note.csv",
                                  "exh4iobbm2p5p1v",
                                  sep = ",",
                                  header = TRUE)

FillIn: a function for filling in missing data in one data frame with info from another

2013-02-15T11:16:00.000+00:00

Update (10 March 2013): FillIn is now part of the budding DataCombine package.

Sometimes I want to use R to fill in values that are missing in one data frame with values from another. For example, I have data from the World Bank on government deficits. However, there are some country-years with missing data. I gathered data from Eurostat on deficits and want to use this data to fill in some of the values that are missing from my World Bank data.

Doing this is kind of a pain so I created a function that would do it for me. It's called FillIn.

An Example

Here is an example using some fake data. (This example and part of the function was inspired by a Stack Exchange conversation between JD Long and Josh O'Brien.)

First let's make two data frames: one with missing values in a variable called fNA. And a data frame with a more complete variable called fFull.

# Create data set with missing values
naDF <- data.frame(a = sample(c(1,2), 100, rep=TRUE), 
                   b = sample(c(3,4), 100, rep=TRUE), 
                   fNA = sample(c(100, 200, 300, 400, NA), 100, rep=TRUE))
                   
# Created full data set
fillDF <- data.frame(a = c(1,2,1,2), 
                     b = c(3,3,4,4),
                     fFull = c(100, 200, 300, 400))

Now we just enter some information into FillIn about what the data set names are, what variables we want to fill in, and what variables to join the data sets on.

# Fill in missing f's from naDF with values from fillDF
FilledInData <- FillIn(D1 = naDF, D2 = fillDF, 
                       Var1 = "fNA", Var2 = "fFull", KeyVar = c("a", "b"))

## [1] "16 NAs were replaced."
## [1] "The correlation between fNA and fFull is 0.313"

D1 and Var1 are for the data frame and variables you want to fill in. D2 and Var2 are what you want to use to fill them in with. KeyVar specifies what variables you want to use to joint the two data frames.

FillIn lets you know how many missing values it is filling in and what the correlation coefficient is between the two variables you are using. Depending on your missing data issues, this could be an indicator of whether or not Var2 is an appropriate substitute for Var1.

Installation

FillIn is currently available as a GitHub Gist and can be installed with this code:

devtools::source_gist("4959237")

You will need the devtools package to install it. For it to work properly you will also need the data.table package.

The Full Code

InstallOldPackages: a repmis command for installing old R package versions

2013-02-03T09:41:00.000+00:00

A big problem in reproducible research is that software changes. The code you used to do a piece of research may depend on a specific version of software that has since been changed. This is an annoying problem in R because install.packages only installs the most recent version of a package. It can be tedious to collect the old versions.

On Toby Dylan Hocking's suggestion, I added tools to the repmis package so that you can install, load, and cite specific R package versions. It should work for any package version that is stored on the CRAN archive (http://cran.r-project.org).

To only install old package versions use the new repmis command InstallOldPackages. For example:

# Install old versions of the e1071 and gtools packages.

# Create vectors of the package names and versions to install
# Note the names and version numbers must be in the same order
Names <- c("e1071", "gtools")
Vers <- c("1.6", "2.6.1")

# Install old package versions into the default library
InstallOldPackages(pkgs = Names, versions = Vers)

You can also now have LoadandCite install specific package versions:

# Install, load, and cite specific package versions

# Create vectors of the package names and versions to install
# Note the names and version numbers must be in the same order
Names <- c("e1071", "gtools") 
Vers <- c("1.6", "2.6.1")

# Run LoadandCite
LoadandCite(pkgs = Names, versions = Vers, install = TRUE, file = "PackageCites.bib")

See this post for more details on LoadandCite.

Future

I intend to continue improving these capabilities. So please post any suggestions for improvement (or report any bugs) at on the GitHub issues page.

repmis: misc. tools for reproducible research in R

2013-01-31T04:32:00.000+00:00

I've started to put together an R package called repmis. It has miscellaneous tools for reproducible research with R. The idea behind the package is to collate commands that simplify some of the common R code used within knitr-type reproducible research papers.

It's still very much in the early stages of development and has two commands:

LoadandCite: a command to load all of the R packages used in a paper and create a BibTeX file containing citation information for them. It can also install the packages if they are on CRAN.
source_GitHubData: a command for downloading plain-text formatted data stored on GitHub or at any other secure (https) URL.

I've written about why you might want to use source_GitHubData before (see here and here).

You can use LoadandCite in a code chunk near the beginning of a knitr reproducible research document to load all of the R packages you will use in the document and automatically generate a BibTeX file you can draw on to cite them. Here's an example:

# Create vector of package names
PackagesUsed <- c("knitr", "xtable")

# Load and Cite
repmis::LoadandCite(PackagesUsed, file = "PackageCitations.bib")

LoadandCite draws on knitr's write_bib command to create the bibliographies, so each citation is given a BibTeX key like this: R-package_name. For example the key for the xtable package is R-xtable. Be careful to save the citations in a new .bib file, because LoadandCite overwrites existing files.

Citation of R packages is very inconsistent in academic publications. Hopefully by making it easier to cite packages more people will do so.

Install/Constribute

Instructions for how to install repmis are available here.

Please feel free to fork the package and suggest additional commands that could be included.

source_GitHubData: a simple function for downloading data from GitHub into R

2013-01-06T09:42:00.000+00:00

Update 31 January: I've folded source_GitHubData into the repmis packaged. See this post.

Update 7 January 2012: I updated the internal workings of source_GitHubData so that it now relies on httr rather than RCurl. Also it is more directly descended from devtool's source_url command.

This has two advantages.

Shortened URL's can be used instead of the data sets' full GitHub URL,
The ssl.verifypeer issue is resolved. (Though please let me know if you have problems).

The post has been rewritten to reflect these changes.

In previous posts I've discussed how to download data stored in plain-text data files (e.g. CSV, TSV) on GitHub directly into R.

Not sure why it took me so long to get around to this, but I've finally created a little function that simplifies the process of downloading plain-text data from GitHub. It's called source_GitHubData. (The name mimicks the devtools syntax for functions like source_gist and source_url. The function's syntax is actually just a modified version of source_url.)

The function is stored in a GitHub Gist HERE (it's also at the end of this post). You can load it directly into R with devtools' source_gist command.

Here is an example of how to use the function to download the electoral disproportionality data I discussed in an earlier post.


# Load source_GitHubData
library(devtools)

# The functions' gist ID is 4466237
source_gist("4466237")

# Create Disproportionality data UrlAddress object
# Make sure the URL is for the "raw" version of the file
# The URL was shortened using bitly
UrlAddress <- "http://bit.ly/Ss6zDO"

# Download data
Data <- source_GitHubData(url = UrlAddress)

# Show Data variable names
names(Data)

## [1] "country"            "year"               "disproportionality"

There you go.

Note that the the function is set by default to load comma-separated data (CSV). This can easily be changed with the sep argument.

Update to Graphing Non-Proportional Hazards in R

2012-12-30T07:40:00.000+00:00

Update 31 July 2013: I've moved all of the functionality described in this post into an R package called simPH. Have a look. It is much easier to use.

This is a quick update for a previous post on Graphing Non-Proportional Hazards in R.

In the previous post I showed how to simulate and graph 1,000 non-proportional hazard ratios at roughly every point in time across an observation period. In the previous example I kept in simulation outliers. Some people have suggested dropping the top and bottom 2.5 percent of simulated values (i.e. keeping the middle 95 percent).

Luckily this can be accomplished with Hadley Wickham's plyr package and three lines of code. The trick is to use plyr's ddply command to subset the data frame at each point in Time where we simulated values. In the previous example the simulated values were in a variable called HRqmv. In each subset we use the quantile command from base R to create logical variables indicating if a simulation of HRqmv is greater than the 0.975 or less than the 0.025 quantile. Then we simply subset the data frame.

Here's how you do it: after all of the values have been simulated and the data set ordered and right before graphing the results add this code:


# Indicate bottom 2.5% of observations
TVSimPerc <- ddply(TVSim, .(Time), transform, Lower = HRqmv < quantile(HRqmv, c(0.025)))

# Indicate top 2.5% of observations
TVSimPerc <- ddply(TVSimPerc, .(Time), transform, Upper = HRqmv > quantile(HRqmv, c(0.975)))

# Drop simulations outside of the middle 95%
TVSimPerc <- subset(TVSimPerc, Lower == FALSE & Upper == FALSE)

Adding this code to the previous example creates this graph:

Here is the full code to replicate the example:

Timeline Maps with googleVis & Twitter Bootstrap Carousel (& updated Slidify)

2012-11-14T06:54:00.000+00:00

I've wanted to create timeline maps with interactive googleVis Geomaps for a while. These would be a nice way to quickly show the spatial distribution of some data over time.

It turns out that it's pretty easy to do with a plugin for Twitter Bootstrap called Carousel. Carousel is probably intended for regular picture slide shows. But because it can hold iframes, it can pretty much include anything, even interactive maps.

Here is a short slide show with examples and code for how to combine googleVis and Twitter Bootstrap Carousel to create interactive timeline maps.

Note: I used the newest version (0.3.1) of Ramnath Vaidyanathan's Slidify to create the presentation. He is really putting a lot of good work into that package. I especially like the choice to set the default slide framework to Google's I/O 2012 style. It has many features you don't find in other HTML slide frameworks. Particularly useful here, it begins to load iframes when you are on the previous slide rather than waiting until you get to the slide with the iframe.

Code:

Code for the example website with a googleVis Carousel timeline map shown in the presentation.
Code for the presentation on GitHub. (Note: the code I used to 'Slidify' the presentation is virtually identical to Ramnath Vaidyanathan's example.

Update (16 Nov. 2012): Ramnath Vaidyanathan has created a great demonstration of how to use the whisker package to include this timeline map directly into Slidify slides.

Graphing Non-Proportional Hazards in R

2012-10-04T06:03:00.000+01:00

Update 30 July 2013: I've moved all of the functionality described in this post into an R package called simPH. Have a look. It is much easier to use.

Update 30 December 2012: I updated the code HERE so that it keeps only the middle 95 percent of the simulated values.

I really like this article by Amanda Licht in Political Analysis. It gives a lot of information on how to interpret nonproportional hazards and includes some nice graphs.

Her source code is really helpful for learning the nuts and bolts of how to simulate quantities of interests over time. However, it's in Stata code, which doesn't really fit into my R-based workflow at the moment. So I decided to port the code over. This post gives an example of what I did.

What is a non-proportional hazard & why use them?

Here is my motivation for being interested in non-proportional hazards: In a few papers I used Cox Proportional Hazard (PH) models to examine countries' policy adoption decisions. The problem with Cox PH models is that they assume the risk of adopting a policy at each point in time changes proportionally. However, in my papers I had good reasons to suspect that this risk actually changed non-proportionally over time. To overcome these model-assumption-busting problems I included time interactions, a pretty common way of dealing with the problem.

When you include time interactions, the variables' coefficients (\( \beta_{cc} \)) are now a combination of the estimates for the base coeffiecient (\( \beta_{1} \)) plus the time interaction (\( \beta_{2}(t) \)):

\[ \beta_{cc} = \beta_{1} + \beta_{2}(t) \]

You'll probably notice that the combined coefficient is different at each point in time \( t \).

See Amanda's original article for a more detailed discussion of these issues.

The Problem

The problem I had before was conveying the substantive effect of the combined coefficient and our uncertainty about it over time. Licht's article solves these problems with Stata. How do we do it in R?

Non-proportional Relative Hazards in R

To demonstrate how to achieve these goals I'll recreate part of a figure from Licht's paper using R. This figure uses data/methods from an earlier paper by Golub & Steunenberg (2007). You can download the data here (it's called GolubEUPdata.tab). My goal is a figure plotting simulated values of the relative non-proportional hazard for the variable qmv between times 80 and 2000.

So, this is the figure I created:

The full R script for doing this is HERE.

Step-by-step

These are the main steps I used to create the figure:

First, run the Cox PH model with time interactions using the coxph command from the survival package. Assign the output of this model to a new object (in my example I call the output object M1).
Extract the coefficeient and variance-covariance matrices from the output object like this:

# Create coefficient matrix
Coef <- matrix(M1$coefficients)

# Create variance-covariance matrix
VC <- vcov(M1)

Using the rmultnorm command from the MSBVAR package, simulate coefficients.

Drawn <- rmultnorm(n = 1000, mu = Coef, vmat = VC)

Clean up the data. For this example we only need two of the columns of simulated coefficients (numbers 1 & 13). We also need to create a simulation ID variable which will be important later.

# Keep qmv and Lqmv
Drawn <- data.frame(Drawn[, c(1, 13)])

# Create Merge Variable (Simulation Number)
Drawn$ID <- 1:1000

Now that we have the simulated coefficients, we can calculate the simulated combined coefficient for qmv. This involves first creating a vector with the times over which we want to calculate the coefficients. We then combine this with the simulated data as in the equation above.

# Create Combined time interactionsCoefficient
TVSim <- outer(Drawn[,2], Range)
TVSim <- data.frame(melt(TVSim))
TVSim <- rename(TVSim, c(X1 = "ID", X2 = "Time", value = "TVR"))

# Merge in the non-time interacted coefficient and combine
TVSim <- merge(Drawn, TVSim, by = "ID")
TVSim$CCqmv <- TVSim$qmv + TVSim$TVR

Now we can exponentiate the coefficient to get the relative hazard

# Create Combined Relative Hazard
TVSim$HRqmv <- exp(TVSim$CCqmv)

Finally, after ordering the data by time, we can plot our simulated relative hazards over time.

# Order Variables
TVSim <- TVSim[order(TVSim$Time),]

# Graph Simulated Combined Hazard Ratios
ggplot(TVSim, aes(Time, HRqmv)) +
  geom_point(shape = 21, alpha = I(0.01), colour = "#FA9FB5", size = 5) +
  geom_smooth() +
  geom_hline(aes(yintercept = 1), linetype = "dotted") +
  scale_y_continuous(breaks = c(-1, 0, 1, 3, 6)) +
  scale_x_continuous(breaks = c(0, 129), labels = c(80, 2000)) +
  xlab("Time in Days") + ylab("Simulated QMV Relative Hazard\n") +
  ggtitle("Simulated Relative Hazard for QMV from times 80-2000\n
          Based on Licht (2011) Fig. 2\n") +
  theme_bw(base_size = 15)

Licht points out that the relative hazard is substantively interesting for dummy variables (which have values 0 and 1), but we probably want other approaches for continuous variables. Only a few modifications to the code above are required to implement these.

Thanks to Jeff Chweiroth for pointing me to Amanda Licht's article and motivating me to put together this code.

Federal Register API/R Package Ideas?

2012-09-22T05:24:00.000+01:00

The other day Critical Juncture put up an API for the Federal Register. I thought it would be great if there was a package that could use this API to download data directly into R (much like the excellent WDI package).

This would make it easier to analyse things like:

the frequency of regulations issued on a particular issue over a given period of time,
the text of the actual regulations.

The nice people over at Critical Juncture tweeted me showing interest in the idea and wondering what would be useful.

I was thinking that in the package there could be commands such as getFedRegister and getMultiFedRegister that would do pretty much do what the API is set up to help now, except download the data into an R object rather than straight to JSON or CSV.

More Ideas?

Any other ideas for things that might be useful? Just leave them in the comments at my Tumblr site.

Create Beamer/knitr Lecture Slideshow with Bash, Explain the Script with knitr

2012-09-17T10:13:00.001+01:00

Setting up a beamer slideshow is tedious. Creating new slideshows with the same header/footer/style files every week for your course lectures is very very tedious.

To solve this problem I created a simple bash shell script. When you run the script in your terminal it asks whether you want to create a "Lecture" or "Seminar" and what number you want it to have. Then it does the rest.

You can find the script and all of the necessary files here.

To create the README file I used knitr version 0.8's new engine='bash' option. This allows you to knit bash code into your Markdown file the same what you would R code. It's pretty simple. See the R Markdown file for more details.

Please feel free to take and modify the files. Also, if you can help streamline them that would be great.

Oh kind of related tip: If you want a bash command to show up over more than one line in your knitted document place a backslash (\) at the end of the line.

The beamer theme I use is based on something I hammered together awhile ago. See this post for more details.