**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: