7 May 2008
Here's a link to a free, 18-hour mini-course on recent advances in econometrics and statistics from the National Bureau of Economic Research. It's co-taught by Guido Imbens and Jeffrey Wooldridge. The intended audience is obviously economists, but there are several topics (Bayesian inference, missing data, etc.) that are likely of interest to a wide range of social scientists. The course includes lecture videos, slides, as well as detailed notes on each topic.
One of the pesky things I've found in my (limited) experience with survival analysis is that it's almost impossible to plot several survival curves in the same space and include measures of uncertainty without the entire plot becoming incomprehensible. So, to build on the great R discussions Ellie and Andy have provided in recent blog posts, I'd like to offer an extension of my own. I've created a fairly flexible function that allows one to plot several survival curves along with estimation uncertainty from Zelig's Cox proportional hazards output (which was developed by Patrick Lam). Here are two examples of what my surv.plot() function can provide:
Hopefully this will be of some interest to a few readers. More details and example code below.
Here is the syntax for the command:
s.out: Simulated output from Zelig for each curve organized as a list()
duration: Surival time
censor: Censoring indicator
type: Display type for confidence bands. The default is "line" but "poly" is also supported (to create the shaded region in the right plot above).
plotcensor: Creates rug() plot indicating censoring times (Default is TRUE)
plottimes: plots a point for each survival time in the step function (Default is TRUE)
int: Desired uncertainty interval (Default is c(0.025,0.975) which corresponds to a 95% interval)
Here's the plot.surv() source code, and below I've copied the R code I used to create the plots above:
# Fit the Cox Model
z.out1 <- zelig(Surv(duration,ciep12)~invest+numst2+crisis,
# Set Low and High Quantities of Interest
low <- setx(z.out1,numst2=0)
high <- setx(z.out1,numst2=1)
# Simulate for Each
s.out1 <- sim(z.out1,x=low)
s.out2 <- sim(z.out1,x=high)
# Create list output that contains both simulations
out <- list(s.out1,s.out2)
# Plot the results