Applied econometrics: Econ 508

Applied Econometrics 
Econ 508 - Fall 2014

Professor: Roger Koenker 

TA: Nicolas Bottan 

Welcome to a new issue of e-Tutorial. This e-TA will focus on on Duration Models (a.k.a. Survival Analysis) in the context of the PS5. 1


You can download the data set, called weco14.txt from the Econ 508 web site. Save it in your preferred directory.

Then you can load it in R. If you saved it in your hard drive in a folder called econ508 you can load it by typing:

weco<-read.table("C:/econ508/weco14.txt", header=T, sep="")

or you can retrieve it from the web

weco<-read.table("", header=T, sep="")
         y sex dex lex  kwit job_tenure status
1 13.72579   0  38  10 FALSE        277   TRUE
2 17.15373   1  55  11  TRUE        173   TRUE
3 13.62992   1  45  12 FALSE        410   TRUE
4 13.03998   1  41  11 FALSE        247   TRUE
5 13.19510   1  42  10 FALSE        340   TRUE
6 16.43010   0  38  12 FALSE        517   TRUE

Next generate the variables needed:


Survival Analysis


To do this kind of analysis in R we are going to use the package survival

Loading required package: survival
Loading required package: splines

Then we need to identify the “analysis time” variable, and the “failure” variable. The former indicates the duration of the process, while the latter indicates whether the data is censored. In the PS5 data set, “tenure” represents the “analysis-time” variable, i.e., the duration of the process, while “status” represents the “failure” variable, assuming values of 0 if it is censored, and 1 if it is failure.

Initially we need to generate the Kaplan-Meier estimator for men and women:

fit <- survfit(Surv(job_tenure,status)~sex, data=weco)
plot(fit, lty=c(1,3), xlab="Time", ylab="Survival Probability")
legend(2000,0.7, legend=c("Female", "Male"), lty=c(1,3) )

Next you may want to test formally for differences between this groups. To do so you can use the survdiff function:

survdiff(Surv(job_tenure,status)~sex, data=weco)
survdiff(formula = Surv(job_tenure, status) ~ sex, data = weco)

        N Observed Expected (O-E)^2/E (O-E)^2/V
sex=0 296      240      279      5.38      10.6
sex=1 387      332      293      5.12      10.6

 Chisq= 10.6  on 1 degrees of freedom, p= 0.00111 

Cox proportional hazard model

Next the PS asks for the estimation of a Cox proportional hazard model. You can estimate such model as follows:

coxph<-coxph(Surv(job_tenure,status)~sex+dex+lex+lex2, data=weco)
coxph(formula = Surv(job_tenure, status) ~ sex + dex + lex + 
    lex2, data = weco)

        coef exp(coef) se(coef)      z       p
sex   0.5444     1.724  0.08819   6.17 6.7e-10
dex  -0.0920     0.912  0.00668 -13.77 0.0e+00
lex  -1.1586     0.314  0.32682  -3.54 3.9e-04
lex2  0.0464     1.047  0.01317   3.52 4.2e-04

Likelihood ratio test=217  on 4 df, p=0  n= 683, number of events= 572 

And plot the results by calling:


  1. Please send comments to or