### Update to Graphing Non-Proportional Hazards in R

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:

### Showing results from Cox Proportional Hazard Models in R with simPH

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 th

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

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