Germán Rodríguez
Survival Analysis Princeton University

Child Mortality in Guatemala: A Shared Frailty Model

We will show how to fit a shared frailty model using data on child mortality in Guatemala. The data come from a retrospective survey conducted in 1974-76 by the Instituto de Nutrición de Centroamérica y Panamá (INCAP) and RAND, and include data on multiple children per family.

The data were first analyzed ignoring the clustering by Pebley and Stupp (1987), Reproductive Patterns and Child Mortality in Guatemala. Demography,24:43-60. (This link requires access to JSTOR.)

The data were then reanalized by Guo and Rodríguez (1992), Estimating a Multivariate Proportional Hazards Model for Clustered Data Using the EM Algorithm, with an Application to Child Survival in Guatemala, Journal of the American Statistical Association,87:969-976, using what we now call a shared frailty model.

The dataset is available as a Stata file in long format (one record per child):

. use http://data.princeton.edu/pop509/PebleyStupp
(Child mortality in Guatemala)
 
. desc
 
Contains data from http://data.princeton.edu/pop509/PebleyStupp.dta
  obs:         3,120                          Child mortality in Guatemala
 vars:            19                          8 Mar 2005 20:48
 size:       137,280 (99.6% of memory free)   (_dta has notes)
-------------------------------------------------------------------------------
              storage  display     value
variable name   type   format      label      variable label
-------------------------------------------------------------------------------
momid           float  %9.0g                  
kid             byte   %9.0g                  
kids            byte   %8.0g                  
time            float  %9.0g                  
income          float  %9.0g                  
death           byte   %8.0g                  
sex             byte   %8.0g                  
mage            float  %9.0g                  
borde           float  %9.0g                  
home            byte   %8.0g                  
pdead           byte   %8.0g                  
edyears         float  %9.0g                  
p0014           byte   %8.0g                  
p1523           byte   %8.0g                  
p2435           byte   %8.0g                  
p36up           byte   %8.0g                  
f0011           byte   %8.0g                  
f1223           byte   %8.0g                  
kidid           float  %9.0g                  
-------------------------------------------------------------------------------
Sorted by:  momid  kid

A Piece-wise exponential model

The first step is to reproduce exactly the results in Pebley and Stupp. We start by stseting the data noting that time is age at death or at interview, death is an indicator of whether the child died, and kidid uniquely identifies each kid. (This is the 10 times the momid in the original dataset plus the kid number.)

. stset time, fail(death) id(kidid)
 
                id:  kidid
     failure event:  death != 0 & death < .
Obs. Time interval:  (time[_n-1], time]
 exit on or before:  failure
 
------------------------------------------------------------------------------
     3120  total obs.
        0  exclusions
------------------------------------------------------------------------------
     3120  obs. Remaining, representing
     3120  subjects
      403  failures in single failure-per-subject data
 131512.5  total analysis time at risk, at risk from t =         0
                             earliest observed entry t =         0
                                  last observed exit t =        60

Pebley and Stupp used a piece-wise exponential model where the hazard is constant in intervals with boundaries at 1, 6, 12 and 24 months. We can fit that kind of model in Stata using the standard Poisson trick or using streg to fit an exponential model with dummy variables for the duration categories. We will do the latter because it allows us to introduce gamma frailty later. So we stsplit the dataset and generate dummies for all categories. (Note that the original paper shows the hazard for each age, so we generate dummies for all categories and omit the constant.)

. stsplit dur, at(0 1 6 12 24 60)
(10474 observations (episodes) created)
 
. gen a0      = dur==0
 
. gen a1to5   = dur==1
 
. gen a6to11  = dur==6
 
. gen a12to23 = dur==12
 
. gen a24up   = dur==24
 
. local age "a0 a1to5 a6to11 a12to23 a24up"

The model includes mother's age and age squared, a linear term on birth order, and an indicator of whether the previous kid in the family died. There are also indicators for the length of the previous birth interval with categories 0-14, 15-23, 24-35 and 36 or more, with first births serving as the reference cell. We need to compute only one of these variables. It is useful to create macros for all:

. gen mage2 = mage^2
 
. local fixed "mage mage2 borde pdead"
 
. local previous "p0014 p1523 p2435 p36up"

The model also includes a time-varying covariate with time-varying effects. The variable in question is the length of the following birth interval, with indicators f0011 for intervals under a year and f1223 for intervals of one to two years. These dummies capture very short and short intervals where a child would face competition from a subsequent sibling. Pebley and Stupp treat these as time-varying covariates because a sibling is assumed to affect the life of the index child only after it is born. So we consider very short intervals (< 12) only at ages 12 months and higher, and short intervals (12-23 months) only at ages 24 months and higher. This is the time-varying part. But the effect of very short intervals is also allowed to vary, with intervals < 12 having different effects at ages 12-23 and 24 or more. This means that we need to create three variables that interact the length of the following interval with the age of the child:

. gen i011a1223 = f0011 * (dur == 12)
 
. gen i011a24p  = f0011 * (dur == 24)
 
. gen i1223a24p = f1223 * (dur == 24)
 
. local follow "i011a1223 i011a24p i1223a24p"
We are now ready to fit the model. You should verify that this reproduces exactly the results in Pebley and Stupp, which is quite remarkable.

. streg `age' `fixed' `previous' `follow', dist(exponential) noconstant
 
         failure _d:  death
   analysis time _t:  time
                 id:  kidid
 
Iteration 0:   log likelihood =  -2819.899  
Iteration 1:   log likelihood = -2115.6375  
Iteration 2:   log likelihood = -1864.5861  
Iteration 3:   log likelihood = -1849.9913  
Iteration 4:   log likelihood = -1848.7872  
Iteration 5:   log likelihood = -1848.7838  
Iteration 6:   log likelihood = -1848.7838  
 
Exponential regression -- log relative-hazard form 
 
No. of subjects =         3120                     Number of obs   =     13594
No. of failures =          403
Time at risk    =     131512.5
                                                   Wald chi2(16)   =   9905.61
Log likelihood  =   -1848.7838                     Prob > chi2     =    0.0000
 
------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      Z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          a0 |   .3375566   .2532262    -1.45   0.148     .0775884    1.468575
       a1to5 |   .0245332   .0186095    -4.89   0.000     .0055473    .1084995
      a6to11 |   .0303346   .0229222    -4.63   0.000     .0068981    .1333968
     a12to23 |   .0180935   .0136752    -5.31   0.000     .0041132    .0795915
       a24up |   .0032229   .0025034    -7.39   0.000     .0007032    .0147712
        mage |   .8612339   .0501274    -2.57   0.010     .7683828    .9653051
       mage2 |   1.002581    .001035     2.50   0.013     1.000555    1.004612
       borde |   1.063634   .0355569     1.85   0.065     .9961777    1.135658
       pdead |   1.103248   .1649263     0.66   0.511     .8230493    1.478839
       p0014 |    1.71365   .3640741     2.54   0.011     1.130004     2.59875
       p1523 |   .8845005   .1648386    -0.66   0.510     .6138541    1.274474
       p2435 |   .7718101   .1424211    -1.40   0.160     .5375754    1.108107
       p36up |   .6764905    .141174    -1.87   0.061     .4493947    1.018346
   i011a1223 |   2.246906   1.609712     1.13   0.258     .5517892    9.149482
    i011a24p |   4.933551   3.633516     2.17   0.030     1.164816    20.89594
   i1223a24p |    1.07556   .4053648     0.19   0.847     .5138396    2.251343
------------------------------------------------------------------------------
note: no constant term was estimated in the main equation
 
. estimates store pwe

A Shared-Frailty Model

We now introduce a shared frailty term at the mother's level to allow for intra-family correlation in child survival. Following standard practice we assume that the frailty term has a gamma distribution independent of observed covariates. Guo and Rodriguez derive an EM algorithm to fit this model (and Guo's dissertation includes the Fortran program used to fit the model). I now show how to reproduce their results using streg. All we need to do is specify the frailty distribution and the variable that identifies the clusters:

. streg `age' `fixed' `previous' `follow', dist(exponential) noconstant ///
>         frailty(gamma) shared(momid)
 
         failure _d:  death
   analysis time _t:  time
                 id:  kidid
 
Fitting exponential model:
 
Fitting full model:
 
Iteration 0:   log likelihood = -1864.1367  
Iteration 1:   log likelihood = -1848.4841  
Iteration 2:   log likelihood = -1847.2191  
Iteration 3:   log likelihood =  -1847.139  
Iteration 4:   log likelihood = -1847.1382  
Iteration 5:   log likelihood = -1847.1382  
 
Exponential regression --
         log-relative hazard form               Number of obs      =     13594
         Gamma shared frailty                   Number of groups   =       851
Group variable: momid
 
No. of subjects =         3120                  Obs per group: min =         1
No. of failures =          403                                 avg =  15.97415
Time at risk    =     131512.5                                 max =        40
 
                                                Wald chi2(16)      =   8796.54
Log likelihood  =   -1847.1382                  Prob > chi2        =    0.0000
 
------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      Z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          a0 |   .3710262   .2874212    -1.28   0.201     .0812846    1.693562
       a1to5 |   .0271691   .0212782    -4.60   0.000     .0058538    .1261004
      a6to11 |   .0337345   .0263341    -4.34   0.000     .0073047    .1557928
     a12to23 |   .0202894   .0158516    -4.99   0.000     .0043878    .0938193
       a24up |   .0036443   .0029231    -7.00   0.000     .0007566    .0175539
        mage |   .8559554   .0514063    -2.59   0.010     .7609048    .9628796
       mage2 |   1.002689   .0010689     2.52   0.012     1.000596    1.004786
       borde |   1.059148   .0376965     1.61   0.106     .9877825     1.13567
       pdead |   .9272842   .1655893    -0.42   0.672     .6534463    1.315878
       p0014 |     1.7738   .3859226     2.63   0.008     1.158004    2.717059
       p1523 |   .9075818   .1718846    -0.51   0.609     .6261508    1.315505
       p2435 |   .7960679    .149776    -1.21   0.225     .5505554    1.151063
       p36up |   .6898139   .1463948    -1.75   0.080     .4550791    1.045628
   i011a1223 |   2.210407   1.593665     1.10   0.271     .5379858    9.081837
    i011a24p |   4.959532   3.677954     2.16   0.031       1.1593    21.21708
   i1223a24p |   1.076996   .4067595     0.20   0.844     .5137269    2.257852
-------------+----------------------------------------------------------------
     /ln_the |  -1.540714    .632364    -2.44   0.015    -2.780125   -.3013038
-------------+----------------------------------------------------------------
       theta |    .214228   .1354701                      .0620307     .739853
------------------------------------------------------------------------------
Likelihood-ratio test of theta=0: chibar2(01) =     3.29 Prob>=chibar2 = 0.035
note: no constant term was estimated in the main equation
 
. estimates store gamma

You should verify that the results agree quite closely with the paper. (The labels for the complex interaction terms involving the following birth interval are out of order in the published paper, but OK in Guo's dissertation.)

The estimated variance of frailty of 0.21 (or 0.22 in the paper) implies modest association between the lifetimes of siblings. Following Oakes it corresponds to a rank correlation of 0.09. Even this small association, however, translates into substantial hazard ratios in the sense of Clayton. The conditional hazard for a child knowing that a sibling died at age a is 22% higher than if the sibling survived to age a, holding all covariates constant. Knowing that two siblings died raises the hazard by 44%. Note that to test the significance of the variance we have to be careful because the test is on a boundary of the parameter space. In the JASA paper we obtain a t-ratio of 1.6 for the variance, which agrees with Stata (0.214/0.135=1.6, although working in the log-scale may be better). Stata's likelihood ratio test gives a chi-squared of 3.3 with an estimated p-value of 0.035, essentially the same result.

. estimates table pwe gamma, t  eform
 
----------------------------------------
    Variable |    pwe         gamma     
-------------+--------------------------
_t           |                          
          a0 |  .33755656    .37102625  
             |      -1.45        -1.28  
       a1to5 |  .02453322    .02716914  
             |      -4.89        -4.60  
      a6to11 |  .03033465    .03373445  
             |      -4.63        -4.34  
     a12to23 |  .01809347    .02028935  
             |      -5.31        -4.99  
       a24up |  .00322289    .00364434  
             |      -7.39        -7.00  
        mage |  .86123393    .85595543  
             |      -2.57        -2.59  
       mage2 |  1.0025814    1.0026886  
             |       2.50         2.52  
       borde |  1.0636339    1.0591481  
             |       1.85         1.61  
       pdead |  1.1032484    .92728419  
             |       0.66        -0.42  
       p0014 |  1.7136501    1.7737999  
             |       2.54         2.63  
       p1523 |  .88450054    .90758181  
             |      -0.66        -0.51  
       p2435 |  .77181013    .79606791  
             |      -1.40        -1.21  
       p36up |  .67649048    .68981395  
             |      -1.87        -1.75  
   i011a1223 |  2.2469057    2.2104071  
             |       1.13         1.10  
    i011a24p |  4.9335514    4.9595321  
             |       2.17         2.16  
   i1223a24p |  1.0755599    1.0769956  
             |       0.19         0.20  
-------------+--------------------------
ln_the       |                          
       _cons |                 .214228  
             |                   -2.44  
----------------------------------------
                             legend: b/t

The other aspect of interest in the results is that the estimates of the covariate effects are remarkably stable. The one change worth mentioning is the coefficient for 'previous child died', which changes sign. This variable was clearly acting as a proxy for unobserved family effects. Note also that the t-ratios for observed covariate effects have generally declined in magnitude, which is consistent with the notion that correlated observations have less information than independent observations.

Guo and Rodríguez also fit a non-parametric frailty model using a two-point distribution. This distribution is not built into Stata, but the model can be fit in Stata with only modest effort by following the method described in Section 3.2 of the paper. Details are left as an exercise for the reader.

We continue our analysis of this dataset with a discussion of model interpretation, including calculation of predicted risks and probabilities.