This example will illustrate how to fit parametric hazard models to continuous duration data (i.e. person-level data). In this example, I use the time between the first and second birth for women in Haiti. The data for this example come from the Haitian Demographic and Health Survey for 2012 individual recode file. This file contains information for all women sampled in the survey.

#Load required libraries
library(foreign)
library(survival)
## Loading required package: splines
library(car)
library(survey)
## Loading required package: grid
## 
## Attaching package: 'survey'
## 
## The following object is masked from 'package:graphics':
## 
##     dotchart
library(muhaz)
library(eha)

#load the data
haiti<-read.dta("/Users/ozd504/Google Drive/dem7223/data//HTIR61FL.DTA", convert.factors = F)

head(test, n=6) In the DHS individual recode file, information on every live birth is collected using a retrospective birth history survey mechanism. Since our outcome is time between first and second birth, we must select as our risk set, only women who have had a first birth. The bidx variable indexes the birth history and if bidx_01 is not missing, then the woman should be at risk of having a second birth (i.e. she has had a first birth, i.e. bidx_01==1). I also select only non-twin births (b0 == 0). The DHS provides the dates of when each child was born in Century Month Codes. To get the interval for women who acutally had a second birth, that is the difference between the CMC for the first birth b3_01 and the second birth b3_02, but for women who had not had a second birth by the time of the interview, the censored time between births is the difference between b3_01 and v008, the date of the interview.

We have 8671 women who are at risk of a second birth.

table(is.na(haiti$bidx_01))
## 
## FALSE  TRUE 
##  8671  5616
#now we extract those women
sub<-subset(haiti, haiti$bidx_01==1&haiti$b0_01==0)

#Here I keep only a few of the variables for the dates, and some characteristics of the women, and details of the survey
sub2<-data.frame(CASEID=sub$caseid, 
                 int.cmc=sub$v008,
                 fbir.cmc=sub$b3_01,
                 sbir.cmc=sub$b3_02,
                 marr.cmc=sub$v509,
                 rural=sub$v025,
                 educ=sub$v106,
                 age=sub$v012,
                 partneredu=sub$v701,
                 partnerage=sub$v730,
                 weight=sub$v005/1000000,
                 psu=sub$v021, strata=sub$v022)

Now I need to calculate the birth intervals, both observed and censored, and the event indicator (i.e. did the women have the second birth?)

sub2$secbi<-ifelse(is.na(sub2$sbir.cmc)==T, ((sub2$int.cmc))-((sub2$fbir.cmc)), (sub2$fbir.cmc-sub2$sbir.cmc))
sub2$b2event<-ifelse(is.na(sub2$sbir.cmc)==T,0,1) 
plot(survfit(Surv(secbi, b2event)~1, sub2), conf.int=T, ylab="S(t)", xlab="Months")
title(main="Survival Function for Second Birth Interval, Haiti", sub="All Women")

Estimating Parametric Hazard Models

While parametric models are not so common in demographic research, fundamental understanding of what they are and how they are constructed is of importance. Some outcomes lend themselves very readily to the parametric approach, but as many demographic duration times are non-unique (tied), the parametric models are not statistically efficient for estimating the survival/hazard functions, as they assume the survival times are continuous random variables. In this section, we first estimate the empirical hazard function and then fit a variety of parametric models to it (Exponential, Weibull, Log-normal and Piecewise exponential). Ideally, a parametric model’s hazard function should approximate the observed empirical hazard function, if the model fits the data.

#since these functions don't work with durations of 0, we add a very small amount to the intervals
fit.haz.km<-kphaz.fit(sub2$secbi[sub2$secbi>0], sub2$b2event[sub2$secbi>0] , method = "product-limit")
#this is a version of the hazard that is smoothed using a kernel-density method
fit.haz.sm<-muhaz(sub2$secbi[sub2$secbi>0], sub2$b2event[sub2$secbi>0] )

#Empirical hazard function (product-limit estimate) plot
kphaz.plot(fit.haz.km)
#overlay the smoothed version
lines(fit.haz.sm, col=2, lwd=3)

So now we see what the empirical hazard function looks like.

#Create some predictor variables: Woman's education, Woman's age^2, Partner's education (> secondary school)
sub2$educ.high<-ifelse(sub2$educ %in% c(2,3), 1, 0)
sub2$age2<-sub2$age^2
sub2$partnerhiedu<-ifelse(sub2$partneredu<3,0,ifelse(sub2$partneredu%in%c(8,9),NA,1 ))

Now we fit the models. I use the eha package to do this, since it fits parametric proportional hazard models, not accellerated failure time models. I prefer the interpretation of regression models on the hazard scale vs. the survival time scale.

Exponential Model

#exponential distribution for hazard, here we hard code it to be
#a weibull dist with shape ==1 
fit.1<-phreg(Surv(secbi, b2event)~educ.high+partnerhiedu+I(age/5), data=sub2[sub2$secbi>0,], dist="weibull", shape=1)
summary(fit.1)
## Call:
## phreg(formula = Surv(secbi, b2event) ~ educ.high + partnerhiedu + 
##     I(age/5), data = sub2[sub2$secbi > 0, ], dist = "weibull", 
##     shape = 1)
## 
## Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
## educ.high           0.332    -0.391     0.676     0.032     0.000 
## partnerhiedu        0.072    -0.397     0.673     0.070     0.000 
## I(age/5)            7.053     0.004     1.004     0.008     0.651 
## 
## log(scale)                    4.067               0.062     0.000 
## 
##  Shape is fixed at  1 
## 
## Events                    5856 
## Total time at risk        380253 
## Max. log. likelihood      -30158 
## LR test statistic         10.35 
## Degrees of freedom        3 
## Overall p-value           0.0158268
plot(fit.1)

Which shows us what the constant hazard model looks like, it assumed the hazard is constant with respect to time, which after seeing the plots above, we know is false. We see the effects of both woman’s and partner’s education are negative, which makes sense. More educated people have lower risks of having more children. We also see the age effect is insignificant, which doesn’t make sense.

Weibull Model

#weibull distribution for hazard
fit.2<-phreg(Surv(secbi, b2event)~educ.high+partnerhiedu+I(age/5), data=sub2[sub2$secbi>0,], dist="weibull")
summary(fit.2)
## Call:
## phreg(formula = Surv(secbi, b2event) ~ educ.high + partnerhiedu + 
##     I(age/5), data = sub2[sub2$secbi > 0, ], dist = "weibull")
## 
## Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
## educ.high           0.332    -0.459     0.632     0.032     0.000 
## partnerhiedu        0.072    -0.492     0.612     0.070     0.000 
## I(age/5)            7.053    -0.070     0.932     0.009     0.000 
## 
## log(scale)                    3.747               0.045     0.000 
## log(shape)                    0.367               0.009     0.000 
## 
## Events                    5856 
## Total time at risk        380253 
## Max. log. likelihood      -29547 
## LR test statistic         399.01 
## Degrees of freedom        3 
## Overall p-value           0
plot(fit.2)

plot(fit.2, fn="haz")
lines(fit.haz.sm, col=2)

Here, we see a more realistic situation, where the hazard fucntion changes over time (Weibull allows this), but compared to the empirical hazard, the model is a very poor fit, as empirically, the hazard goes up, but then goes down. The Weibull hazard just goes up, as the model does not allow the hazard to change direction, only rate of increase (i.e. it can incraese at a slower or faster rate, but not change direction). We also see the effect of mom’s age is significant and negative (older women have lower risk of having a second birth)

The Log-normal distribution is more flexible and allows the hazard to change direction. Log-Normal Model

#log-normal distribution for hazard
fit.3<-phreg(Surv(secbi, b2event)~educ.high+partnerhiedu+I(age/5), data=sub2[sub2$secbi>0,], dist="lognormal", center=T)
summary(fit.3)
## Call:
## phreg(formula = Surv(secbi, b2event) ~ educ.high + partnerhiedu + 
##     I(age/5), data = sub2[sub2$secbi > 0, ], dist = "lognormal", 
##     center = T)
## 
## Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
## (Intercept)                  -2.734               0.186     0.000 
## educ.high           0.332    -0.442     0.643     0.033     0.000 
## partnerhiedu        0.072    -0.408     0.665     0.070     0.000 
## I(age/5)            7.053    -0.034     0.967     0.009     0.000 
## 
## log(scale)                    2.691               0.032     0.000 
## log(shape)                    1.540               0.072     0.000 
## 
## Events                    5856 
## Total time at risk        380253 
## Max. log. likelihood      -28412 
## LR test statistic         323.21 
## Degrees of freedom        3 
## Overall p-value           0
plot(fit.3)