This example will illustrate how to fit the extended Cox Proportional hazards model with Gaussian frailty to continuous duration data (i.e. person-level data) and a discrete-time (longitudinal) data set. In this example, I will use the event of a child dying before age 5 in Haiti. The data for this example come from the Haitian Demographic and Health Survey for 2012 birth recode file. This file contains information for all live births to women sampled in the survey.

The longitudinal data example uses data from the ECLS-K. Specifically, we will examine the transition into poverty between kindergarten and third grade.

#Load required libraries
library(foreign)
library(survival)
library(car)
library(survey)

## Loading required package: grid
## 
## Attaching package: 'survey'
## 
## The following object is masked from 'package:graphics':
## 
##     dotchart

library(coxme)

## Loading required package: bdsmatrix
## 
## Attaching package: 'bdsmatrix'
## 
## The following object is masked from 'package:base':
## 
##     backsolve

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

#We form a subset of variables
sub<-data.frame(CASEID=haiti$caseid,v008=haiti$v008,bord=haiti$bidx,csex=haiti$b4,b2=haiti$b2, b3=haiti$b3, b5=haiti$b5, b7=haiti$b7, ibint=haiti$b11, rural=haiti$v025, educ=haiti$v106,age=haiti$v012,partneredu=haiti$v701,partnerage=haiti$v730, hhses=haiti$v190, weight=haiti$v005/1000000, psu=haiti$v021, strata=haiti$v022, region=haiti$v023)

sub$death.age<-ifelse(sub$b5==1,
                          ((((sub$v008))+1900)-(((sub$b3))+1900)) 
                          ,sub$b7)

#censoring indicator for death by age 5, in months (<=60 months)
sub$d.event<-ifelse(is.na(sub$b7)==T|sub$b7>60,0,1)
sub$d.eventfac<-factor(sub$d.event); levels(sub$d.eventfac)<-c("Alive at Age 5", "Dead by Age 5")
table(sub$d.eventfac)

## 
## Alive at Age 5  Dead by Age 5 
##          25981           3032

#recodes
sub$male<-ifelse(sub$csex==1,1,0)
sub$educ.high<-ifelse(sub$educ %in% c(2,3), 1, 0)
sub$age2<-sub$age^2
sub$partnerhiedu<-ifelse(sub$partneredu<3,0,ifelse(sub$partneredu%in%c(8,9),NA,1 ))
sub$region_name<-recode(sub$region, recodes = "1='Aire Metropolitaine';2='Reste-Ouest';3='Sud-Est';4='Nord'; 5='Nord-Est'; 6='Artibonite'; 7='Centre'; 8='Sud'; 9='Grand Anse';10='Nord-Ouest';11='Nippes';12='Camp' ", as.factor.result=T)
sub$region_num<-ifelse(sub$region==12, NA, sub$region)

Fit the ordinary Cox model

Here I fit the ordinary Cox model without frailty, just for comparison sake.

#using coxph in survival library
fit.cox2<-coxph(Surv(death.age, d.event)~bord+male+educ.high+I(age/5)+I(hhses>3), data=sub, weights=weight)
summary(fit.cox2)

## Call:
## coxph(formula = Surv(death.age, d.event) ~ bord + male + educ.high + 
##     I(age/5) + I(hhses > 3), data = sub, weights = weight)
## 
##   n= 29013, number of events= 3032 
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## bord              0.154211  1.166737  0.008420 18.315  < 2e-16 ***
## male              0.162654  1.176629  0.036940  4.403 1.07e-05 ***
## educ.high        -0.291035  0.747489  0.057436 -5.067 4.04e-07 ***
## I(age/5)         -0.034060  0.966514  0.013292 -2.562   0.0104 *  
## I(hhses > 3)TRUE -0.001588  0.998413  0.044693 -0.036   0.9716    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## bord                1.1667     0.8571    1.1476    1.1862
## male                1.1766     0.8499    1.0945    1.2650
## educ.high           0.7475     1.3378    0.6679    0.8366
## I(age/5)            0.9665     1.0346    0.9417    0.9920
## I(hhses > 3)TRUE    0.9984     1.0016    0.9147    1.0898
## 
## Concordance= 0.612  (se = 0.005 )
## Rsquare= 0.016   (max possible= 0.873 )
## Likelihood ratio test= 461.8  on 5 df,   p=0
## Wald test            = 503.2  on 5 df,   p=0
## Score (logrank) test = 516.8  on 5 df,   p=0

plot(survfit(fit.cox2),  ylim=c(.8,1), xlim=c(0,60), ylab="S(t)", xlab="Age in Months")
title(main="Survival Function for Child Mortality")

Fit the Cox model with frailty at the regional level

The coxme() function in the coxme library Link will fit the Cox model with shared frailty, assuming a Gaussian frailty term. The model would look like:

$h_j (t) = h_{0j} e^{(x'\beta + u_j)}$

where

$u_j \sim N(0, \sigma^2)$

is a Normally distributed random effect, identical for each person in the jth group. This term raises or lowers the average hazard function the same way for each person within each group, but allows the overall risk for people in different groups to be different. This would be considered to be a random intercept model, if we were considering an ordinary linear or genearlized linear model.

fit.cox.f<-coxme(Surv(death.age, d.event)~bord+male+educ.high+I(age/5)+I(hhses>3)+(1|region_name), data=sub, weights=weight)
summary(fit.cox.f)

## Cox mixed-effects model fit by maximum likelihood
##   Data: sub
##   events, n = 3032, 29013
##   Iterations= 7 38 
##                     NULL Integrated    Fitted
## Log-likelihood -29900.32  -29640.46 -29624.44
## 
##                    Chisq    df p    AIC    BIC
## Integrated loglik 519.71  6.00 0 507.71 471.61
##  Penalized loglik 551.76 14.39 0 522.98 436.40
## 
## Model:  Surv(death.age, d.event) ~ bord + male + educ.high + I(age/5) +      I(hhses > 3) + (1 | region_name) 
## Fixed coefficients
##                         coef exp(coef)    se(coef)     z       p
## bord              0.15495566 1.1676062 0.008438343 18.36 0.0e+00
## male              0.16259318 1.1765579 0.036945296  4.40 1.1e-05
## educ.high        -0.29376289 0.7454532 0.057496194 -5.11 3.2e-07
## I(age/5)         -0.03182241 0.9686786 0.013336771 -2.39 1.7e-02
## I(hhses > 3)TRUE -0.17325594 0.8409224 0.052800758 -3.28 1.0e-03
## 
## Random effects
##  Group       Variable  Std Dev    Variance  
##  region_name Intercept 0.18601623 0.03460204

This gives us the variance in child mortality by region, which is honestly pretty substantial. We can use a likelihood ratio test to see if the frailty model is significantly better than the ordinary Cox model:

anova(fit.cox.f, fit.cox2)

## Analysis of Deviance Table
##  Cox model: response is  Surv(death.age, d.event)
##  Model 1: ~bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 | region_name)
##  Model 2: ~bord + male + educ.high + I(age/5) + I(hhses > 3)
##   loglik  Chisq Df P(>|Chi|)    
## 1 -29640                        
## 2 -29669 57.948  1 2.687e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(fit.cox.f)

## [1] 59277.65

AIC(fit.cox2)

## [1] 59348.87

Which it is, and this is supported by the AIC difference between the two models of 71.22 points.

So, what are the frailties in this case? We can get those from the frail portion of the model structure:

fit.cox.f$frail

## $region_name
## Aire Metropolitaine          Artibonite                Camp 
##         0.369173163         0.136504990         0.171973528 
##              Centre          Grand Anse              Nippes 
##         0.032839617        -0.122132751        -0.003100147 
##                Nord            Nord-Est          Nord-Ouest 
##        -0.042252525        -0.092178225        -0.314393906 
##         Reste-Ouest                 Sud             Sud-Est 
##         0.025396603        -0.090208004        -0.071622344

Which shows the region region_name.Aire Metropolitaine has the highest frailty, which means that the average level of childhood mortality is highest in that region, while the lowest frailty is in region_name.Nord-Ouest . This is nice, but wouldn’t it look nice to map these?

library(maptools)

## Loading required package: sp
## Checking rgeos availability: TRUE

htshp<-readShapePoly("~/Google Drive/dem7223/data/sdr_subnational_data_2015-03-18/shps/sdr_subnational_data_dhs_2012.shp")
plot(htshp)
text(getSpPPolygonsLabptSlots(htshp), labels=as.character(htshp$DHSREGFR), cex=0.6)

## Warning: use coordinates method

#Since the region variable in the shapefile doesn't include two of the areas, I need to make them an ID number:
htshp$reg_merge<-ifelse(htshp$DHSREGFR=="Aire Métropolitaine", 1, ifelse(htshp$DHSREGFR=="Reste-Ouest",2, htshp$REGCODE+1))
htshp$reg_merge

##  [1]  3  4  5  6  7  8  9 10 11  2  1

#make a data frame to hold the frailties, and the accompanying ID numbers
frail<-data.frame(frail=exp(unlist(fit.cox.f$frail)), region=c(1,6, 12, 7, 9, 11, 4, 5, 10, 2, 8, 3), regname=gsub(names(unlist(fit.cox.f$frail)), pattern = "region_name.", replacement = ""))
frail

##                                     frail region             regname
## region_name.Aire Metropolitaine 1.4465381      1 Aire Metropolitaine
## region_name.Artibonite          1.1462606      6          Artibonite
## region_name.Camp                1.1876464     12                Camp
## region_name.Centre              1.0333848      7              Centre
## region_name.Grand Anse          0.8850309      9          Grand Anse
## region_name.Nippes              0.9969047     11              Nippes
## region_name.Nord                0.9586277      4                Nord
## region_name.Nord-Est            0.9119426      5            Nord-Est
## region_name.Nord-Ouest          0.7302313     10          Nord-Ouest
## region_name.Reste-Ouest         1.0257218      2         Reste-Ouest
## region_name.Sud                 0.9137411      8                 Sud
## region_name.Sud-Est             0.9308824      3             Sud-Est

#just check the names are the same
htshp@data[, c("DHSREGFR", "reg_merge")]

##               DHSREGFR reg_merge
## 0              Sud-Est         3
## 1                 Nord         4
## 2             Nord-Est         5
## 3           Artibonite         6
## 4               Centre         7
## 5                  Sud         8
## 6          Grande-Anse         9
## 7           Nord-Ouest        10
## 8               Nippes        11
## 9          Reste-Ouest         2
## 10 Aire Métropolitaine         1

#merge the shapefile's data frame to the frailties
mdat<-htshp@data
mdat<-merge(mdat, frail, by.x="reg_merge", by.y="region", all.x=T, sort=F)

#stick it back on the shapefile
htshp@data<-mdat

library(sp)
library(RColorBrewer)
brks<-quantile(htshp$frail); brks[1]<-brks[1]-.1; brks[5]<-brks[5]+.1
rownames(htshp@data)<-as.character(htshp$DHSREGFR)


sp.label <- function(x, label) {
    list("sp.text", coordinates(x), label)
}

ISO.sp.label <- function(x) {
    sp.label(x, row.names(x@data["frail"]))
}

make.ISO.sp.label <- function(x) {
    do.call("list", ISO.sp.label(x))
}

spplot(htshp,"frail",at=brks, col.regions=brewer.pal(n=4,"RdBu" ), sp.layout = make.ISO.sp.label(htshp), main="Frailty Values for Haitian Regions")

Which shows the highest and lowest frailty areas, including the area of Port-au-Prince, in dark blue to the south, called “Aire Métropolitaine”.

Random slopes

If we were interested in whether a predictor variable had heterogeneous effects across the various groups within our data, we could include that in our model as well, and we would have effectively a random slope model:

$h_j (t) = h_{0j} e^{(x'\beta + u_j+\gamma_j 'x)}$

where $\gamma_j$ is a group-specific effect of a particular predictor variable, and these two random effects will be distributed as:

$\left[\begin{array} {rrr} u_j \\ \gamma_j \end{array}\right] \sim MVN(0, \Sigma)$

#See if higher birth order children face equal disadvantage in all regions
fit.cox.f2<-coxme(Surv(death.age, d.event)~bord+male+educ.high+I(age/5)+I(hhses>3)+(1+bord|region_name), data=sub, weights=weight)
summary(fit.cox.f2)

## Cox mixed-effects model fit by maximum likelihood
##   Data: sub
##   events, n = 3032, 29013
##   Iterations= 16 83 
##                     NULL Integrated    Fitted
## Log-likelihood -29900.32  -29639.86 -29619.97
## 
##                    Chisq    df p    AIC    BIC
## Integrated loglik 520.92  8.00 0 504.92 456.79
##  Penalized loglik 560.69 17.46 0 525.77 420.71
## 
## Model:  Surv(death.age, d.event) ~ bord + male + educ.high + I(age/5) +      I(hhses > 3) + (1 + bord | region_name) 
## Fixed coefficients
##                         coef exp(coef)   se(coef)     z       p
## bord              0.15503765 1.1677019 0.01013823 15.29 0.0e+00
## male              0.16247804 1.1764225 0.03694877  4.40 1.1e-05
## educ.high        -0.28940245 0.7487108 0.05758529 -5.03 5.0e-07
## I(age/5)         -0.03176316 0.9687360 0.01333913 -2.38 1.7e-02
## I(hhses > 3)TRUE -0.17257940 0.8414915 0.05276513 -3.27 1.1e-03
## 
## Random effects
##  Group       Variable  Std Dev       Variance      Corr         
##  region_name Intercept  0.1979265059  0.0391749017 -0.3304455077
##              bord       0.0172420452  0.0002972881

anova(fit.cox.f, fit.cox.f2)

## Analysis of Deviance Table
##  Cox model: response is  Surv(death.age, d.event)
##  Model 1: ~bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 | region_name)
##  Model 2: ~bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 + bord |      region_name)
##   loglik  Chisq Df P(>|Chi|)
## 1 -29640                    
## 2 -29640 1.2116  2    0.5456

And it looks like there is significant variation in this effect, because the model with the additional term fits the data significantly better than the model with only the random “intercept”. What does this effect look like?

spplot(htshp,"frail",at=brks, col.regions=brewer.pal(n=4,"Reds" ), sp.layout = make.ISO.sp.label(htshp), main="Random Slopes for Birth Order for Haitian Regions")

And we basically see that the birth order vaiable has stronger negative effects (higher hazards), generally in the southern portions of the country.

Using Longitudinal Data

As in the other examples, I illustrate fitting these models to data that are longitudinal, instead of person-duration. In this example, we will examine how to fit the Cox model to a longitudinally collected data set.

First we load our data

load("~/Google Drive/dem7903_App_Hier/data/eclsk.Rdata")
names(eclsk)<-tolower(names(eclsk))
#get out only the variables I'm going to use for this example
myvars<-c( "childid","gender", "race", "r1_kage","r4age", "r5age", "r6age", "r7age","c1r4mtsc", "c4r4mtsc", "c5r4mtsc", "c6r4mtsc", "c7r4mtsc", "w1povrty","w1povrty","w3povrty", "w5povrty", "w8povrty","wkmomed", "s2_id", "c1_5fp0", "c15fpstr", "c15fppsu")
eclsk<-eclsk[,myvars]


eclsk$age1<-ifelse(eclsk$r1_kage==-9, NA, eclsk$r1_kage/12)
eclsk$age2<-ifelse(eclsk$r4age==-9, NA, eclsk$r4age/12)
#for the later waves, the NCES group the ages into ranges of months, so 1= <105 months, 2=105 to 108 months. So, I fix the age at the midpoint of the interval they give, and make it into years by dividing by 12
eclsk$age3<-recode(eclsk$r5age,recodes="1=105; 2=107; 3=109; 4=112; 5=115; 6=117; -9=NA")/12

eclsk$pov1<-ifelse(eclsk$w1povrty==1,1,0)
eclsk$pov2<-ifelse(eclsk$w3povrty==1,1,0)
eclsk$pov3<-ifelse(eclsk$w5povrty==1,1,0)

#Recode race with white, non Hispanic as reference using dummy vars
eclsk$hisp<-recode (eclsk$race, recodes="3:4=1;-9=NA; else=0")
eclsk$black<-recode (eclsk$race, recodes="2=1;-9=NA; else=0")
eclsk$asian<-recode (eclsk$race, recodes="5=1;-9=NA; else=0")
eclsk$nahn<-recode (eclsk$race, recodes="6:7=1;-9=NA; else=0")
eclsk$other<-recode (eclsk$race, recodes="8=1;-9=NA; else=0")
eclsk$male<-recode(eclsk$gender, recodes="1=1; 2=0; -9=NA")
eclsk$mlths<-recode(eclsk$wkmomed, recodes = "1:2=1; 3:9=0; else = NA")
eclsk$mgths<-recode(eclsk$wkmomed, recodes = "1:3=0; 4:9=1; else =NA")

Now, I need to form the transition variable, this is my event variable, and in this case it will be 1 if a child enters poverty between the first wave of the data and the third grade wave, and 0 otherwise. NOTE I need to remove any children who are already in poverty age wave 1, because they are not at risk of experiencing this particular transition.

eclsk<-subset(eclsk, is.na(pov1)==F&is.na(pov2)==F&is.na(pov3)==F&is.na(age1)==F&is.na(age2)==F&is.na(age3)==F&pov1!=1&is.na(eclsk$c15fpstr)==F)
eclsk$povtran1<-ifelse(eclsk$pov1==0&eclsk$pov2==0, 0,1)
eclsk$povtran2<-ifelse(eclsk$povtran1==1, NA,ifelse(eclsk$pov2==0&eclsk$pov3==0,0,1))

Now we do the entire data set. To analyze data longitudinally, we need to reshape the data from the current “wide” format (repeated measures in columns) to a “long” format (repeated observations in rows). The reshape() function allows us to do this easily. It allows us to specify our repeated measures, time varying covariates as well as time-constant covariates.

#make an id that is the combination of state and strata 
eclsk$sampleid<-paste(eclsk$c15fpstr, eclsk$c15fppsu)
#within each sampling unit, sum the weights
wts<-tapply(eclsk$c1_5fp0,eclsk$sampleid,sum)
#make a data frame from this
wts<-data.frame(id=names(unlist(wts)), wt=unlist(wts))
#get the unique sampling location ids'
t1<-as.data.frame(table(eclsk$sampleid))
#put all of this into a data set
wts2<-data.frame(ids=wts$id, sumwt=wts$wt, jn=t1$Freq)
#merge all of this back to the original data file
eclsk2<-merge(eclsk, wts2, by.x="sampleid", by.y="ids", all.x=T)
#In the new data set, multiply the original weight by the fraction of the
#sampling unit total population each person represents
eclsk2$swts<-eclsk2$c1_5fp0*(eclsk2$jn/eclsk2$sumwt)



e.long<-reshape(eclsk2, idvar="childid", varying=list(age=c("age1","age2"), age2=c("age2", "age3"), povtran=c("povtran1", "povtran2")), times=1:2, direction="long" , drop = names(eclsk)[4:20])
e.long<-e.long[order(e.long$childid, e.long$time),]

#find which kids failed in the first time period and remove them from the second risk period risk set
failed1<-which(is.na(e.long$povtran1)==T)
e.long<-e.long[-failed1,]
e.long$age1r<-round(e.long$age1, 0)
e.long$age2r<-round(e.long$age2, 0)
head(e.long, n=10)

##            sampleid  childid gender race c1_5fp0 c15fpstr c15fppsu pov1
## 0001002C.1     63 1 0001002C      2    1  352.35       63        1    0
## 0001002C.2     63 1 0001002C      2    1  352.35       63        1    0
## 0001007C.1     63 1 0001007C      1    1  358.99       63        1    0
## 0001007C.2     63 1 0001007C      1    1  358.99       63        1    0
## 0001010C.1     63 1 0001010C      2    1  352.35       63        1    0
## 0001010C.2     63 1 0001010C      2    1  352.35       63        1    0
## 0002004C.1     63 1 0002004C      2    1  228.05       63        1    0
## 0002004C.2     63 1 0002004C      2    1  228.05       63        1    0
## 0002006C.1     63 1 0002006C      2    1  287.25       63        1    0
## 0002008C.1     63 1 0002008C      2    1  259.19       63        1    0
##            pov2 pov3 hisp black asian nahn other male mlths mgths    sumwt
## 0001002C.1    0    0    0     0     0    0     0    0     0     0 21844.73
## 0001002C.2    0    0    0     0     0    0     0    0     0     0 21844.73
## 0001007C.1    0    0    0     0     0    0     0    1     0     1 21844.73
## 0001007C.2    0    0    0     0     0    0     0    1     0     1 21844.73
## 0001010C.1    0    0    0     0     0    0     0    0     0     1 21844.73
## 0001010C.2    0    0    0     0     0    0     0    0     0     1 21844.73
## 0002004C.1    0    0    0     0     0    0     0    0     0     1 21844.73
## 0002004C.2    0    0    0     0     0    0     0    0     0     1 21844.73
## 0002006C.1    1    1    0     0     0    0     0    0    NA    NA 21844.73
## 0002008C.1    0    0    0     0     0    0     0    0     0     1 21844.73
##            jn      swts time     age1     age2 povtran1 age1r age2r
## 0001002C.1 72 1.1613419    1 6.433333 7.894167        0     6     8
## 0001002C.2 72 1.1613419    2 7.894167 9.750000        0     8    10
## 0001007C.1 72 1.1832273    1 5.300000 6.825000        0     5     7
## 0001007C.2 72 1.1832273    2 6.825000 8.916667        0     7     9
## 0001010C.1 72 1.1613419    1 5.885833 7.385833        0     6     7
## 0001010C.2 72 1.1613419    2 7.385833 9.333333        0     7     9
## 0002004C.1 72 0.7516504    1 5.450000 6.977500        0     5     7
## 0002004C.2 72 0.7516504    2 6.977500 9.083333        0     7     9
## 0002006C.1 72 0.9467730    1 5.158333 6.685833        1     5     7
## 0002008C.1 72 0.8542875    1 5.905833 7.433333        0     6     7

Now we fit the Cox model and doing fraily by the sample strata (to mimic the sample design). I use the weights calculated above, standardized to the within cluster sample size.

#Fit the model
fitl1<-coxme(Surv(time = age1r, time2 = age2r, event = povtran1)~mlths+mgths+black+hisp+other+nahn+(1|sampleid), e.long)

## Warning in Surv(time = age1r, time2 = age2r, event = povtran1): Stop time
## must be > start time, NA created

summary(fitl1)

## Cox mixed-effects model fit by maximum likelihood
##   Data: e.long
##   events, n = 616, 12978 (159 observations deleted due to missingness)
##   Iterations= 21 130 
##                     NULL Integrated    Fitted
## Log-likelihood -5318.838  -5041.363 -4939.557
## 
##                    Chisq    df p    AIC    BIC
## Integrated loglik 554.95  7.00 0 540.95 509.99
##  Penalized loglik 758.56 96.93 0 564.70 135.94
## 
## Model:  Surv(time = age1r, time2 = age2r, event = povtran1) ~ mlths +      mgths + black + hisp + other + nahn + (1 | sampleid) 
## Fixed coefficients
##             coef exp(coef)   se(coef)     z       p
## mlths  0.7983549 2.2218827 0.11413223  6.99 2.7e-12
## mgths -0.9214421 0.3979448 0.09388894 -9.81 0.0e+00
## black  1.3401604 3.8196561 0.12782313 10.48 0.0e+00
## hisp   0.9333661 2.5430549 0.11046930  8.45 0.0e+00
## other  0.3262753 1.3857968 0.26887849  1.21 2.2e-01
## nahn   0.9889591 2.6884347 0.23186045  4.27 2.0e-05
## 
## Random effects
##  Group    Variable  Std Dev   Variance 
##  sampleid Intercept 0.5085613 0.2586346

Event History Analysis - Example 7 Frailty in the Cox model

coreysparks

February 4, 2015

Fit the ordinary Cox model

Fit the Cox model with frailty at the regional level

Random slopes

Using Longitudinal Data