Germán Rodríguez
Survival Analysis Princeton University

## Interpretation of Frailty Models

We continue our discussion of shared frailty models with a focus on interpretation of the results using various calculations based on the parameter estimates.

### The Distribution of Frailty

The estimated variance of frailty at the family level is 0.2142. In terms of Clayton's model this means that the hazard given that a sibling died at age a is 21.4% higher than the hazard given that the sibling survived age to a.

Using the result of Oakes, we can also say that the Tau correlation between sibling survival times is estimated to be 0.097, so the probability of a concordant pair is about 10 percentage points higher than the probability of a discordant pair. (The "pair" here refers to a pair of mothers, say A and B, each with two children, say 1 and 2. The pair is concordant if both children of mother A die at younger (or older) ages than the children of mother B. The pair is discordant if child 1 of A lives longer than child 1 of B but child 2 of A dies younger than child 2 of B. The interpretation is not terribly useful in this application because most observations are censored and thus one can't establish concordance or discordance.)

A more direct approach is to look at the actual distribution of frailty. Stata has a function `gammaden(a, b, g, x)` to compute the density of the gamma distribution with shape `a`, scale `b` (which is 1/ß in our notation) and location `g` (here 0). `b = 1/a`. R has a function `dgamma(x, shape, rate = 1, scale = 1/rate)` to compute the density of the gamma distribution with given shape and scale (or its reciprocal the rate). Recall that the estimated variance is the reciprocal of the shape parameter.

```. scalar v = exp( _b[ln_the:_cons] )

. twoway function y=gammaden(1/v,v,0,x), range(0 3) ///
>         title("Gamma density with mean 1 and variance 0.214")

. graph export gfr.png, width(500) replace
(file gfr.png written in PNG format)
``` ```> v <- exp( -1.5407144 )

> theta <- seq(0, 3, 0.02)

> density <- dgamma(theta, shape = 1/v, scale = v)

> gd <- data.frame(theta, density)

> ggplot(gd, aes(theta, density)) + geom_line() +
+   ggtitle("Gamma Density (mean 1, var 0.214)")

> ggsave("gfrr.png", width = 500/72, height = 400/72, dpi = 72)
``` We can also compute quantiles. Stata has a function `invgammap(a, p)` to compute quantiles of the standard gamma distribution with shape `a`, which has scale 1 and location 0. A gamma with shape `a` and scale `b` is just `b` times a standard gamma with shape `a`. In R we use `qgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE)` to compute quantiles. So here are the quartiles of the distribution depicted above.

```. mata v = st_numscalar("v")

. mata q = invgammap(1/v, (0.25,.50,.75) )*v

. mata q
1             2             3
+-------------------------------------------+
1 |  .6618944858   .9295688957   1.261736474  |
+-------------------------------------------+

. mata q :/ q :- 1
1              2              3
+----------------------------------------------+
1 |  -.2879554288              0    .3573350829  |
+----------------------------------------------+
```
```> p <- (1:3)/4

> q <- qgamma(p, shape = 1/v, scale = v); q
 0.6618945 0.9295689 1.2617365

> q / q - 1
 -0.2879554  0.0000000  0.3573351
```

So the quartiles are 0.66, 0.93 and 1.26. Computing the ratios Q1/Q2 and Q3/Q2 we see that families with frailty at Q1 have 29% lower risk, and families with frailty at Q3 have 36% higher risk, than families with median frailty. Clearly, unobserved family characteristics have a very substantial effect on child survival.

### Distribution of Observed Risks

It may be interesting to contrast the above results with the risks that can be attributed to observed characteristics.

At this point we need the results from the piecewise exponential model with gamma frailty. I saved the file generated in the previous log as `PebleyStupp2.dta`. We read the file and fit the model again. and then read the parameter estimates from `PebleyStupp2.dat`.

```. use https://data.princeton.edu/pop509/PebleyStupp2, clear
(Child mortality in Guatemala)

. streg i.dur mage mage2 borde pdead p0014 p1523 p2435 p36up i011a1223 i011a24p i1223a24p , ///
>         dist(exponential) frailty(gamma) shared(momid)

failure _d:  death
analysis time _t:  time
id:  kidid

Fitting exponential model:

Fitting constant-only model:

Iteration 0:   log likelihood = -2471.3265  (not concave)
Iteration 1:   log likelihood = -2335.4575
Iteration 2:   log likelihood = -2283.1245
Iteration 3:   log likelihood = -2274.2706
Iteration 4:   log likelihood = -2273.9186
Iteration 5:   log likelihood = -2273.9183
Iteration 6:   log likelihood = -2273.9183

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     =     13,594
Gamma shared frailty                   Number of groups  =        851
Group variable: momid
Obs per group:
No. of subjects =         3120                                min =          1
No. of failures =          403                                avg =   15.97415
Time at risk    =     131512.5                                max =         40

LR chi2(15)       =    853.56
Log likelihood  =   -1847.1382                  Prob > chi2       =    0.0000

------------------------------------------------------------------------------
_t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
dur |
1  |    .073227   .0115575   -16.56   0.000     .0537434     .099774
6  |    .090922   .0127802   -17.06   0.000     .0690275    .1197612
12  |   .0546844   .0075736   -20.98   0.000     .0416845    .0717385
24  |   .0098223   .0021615   -21.01   0.000     .0063811    .0151194
|
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
_cons |   .3710262   .2874212    -1.28   0.201     .0812846    1.693562
-------------+----------------------------------------------------------------
/ln_the |  -1.540714    .632364    -2.44   0.015    -2.780125   -.3013038
-------------+----------------------------------------------------------------
theta |    .214228   .1354701                      .0620307     .739853
------------------------------------------------------------------------------
LR test of theta=0: chibar2(01) = 3.29                 Prob >= chibar2 = 0.035
```
```> library(foreign)

> gux\$death <- gux[,"_d"]

> b <- scan("https://data.princeton.edu/pop509/PebleyStupp2.dat")
```

Next we simply predict the hazard holding family frailty constant at one. Remember, however, that each child contributes one or more pseudo observations, but should be counted once only, so we pick first segments.

```. predict haz, haz
(option alpha1 assumed)

. sum haz if _t0==0, detail

Predicted hazard
-------------------------------------------------------------
Percentiles      Smallest
1%     .0309073       .0280617
5%      .033402       .0297215
10%     .0361795        .029956       Obs               3,120
25%     .0413038        .029956       Sum of Wgt.       3,120

50%     .0470627                      Mean           .0516725
Largest       Std. Dev.      .0182111
75%     .0552813       .1704806
90%     .0726199       .1706315       Variance       .0003316
95%     .0934537        .171702       Skewness       2.365153
99%     .1179419        .187827       Kurtosis       10.73898

. di r(p25), r(p50), r(p75)
.04130384 .04706271 .05528133

. di  r(p25)/r(p50) - 1, r(p75)/r(p50) - 1
-.12236582 .17463127
```
```> mf <-  death ~ factor(dur) + mage + mage2 + borde + pdead +
+   p0014 + p1523 + p2435 + p36up + i011a1223 + i011a24p + i1223a24p

> X <- model.matrix(mf, data = gux)

> haz <- exp(X %*% b)

> first <- gux\$dur == 0

> r <- quantile(haz[first], probs = 1:3/4); r
25%        50%        75%
0.04130319 0.04706190 0.05528053

> r/r - 1
25%        50%        75%
-0.1223645  0.0000000  0.1746346
```

So children who are at Q1 in terms of the observed risk factors have 12% lower risk, and children at Q3 of observed risk factors have 17% higher risk, than children at the median. Looks like characteristics that are observed at birth have a smaller impact on survival than unobserved family characteristics. (We have not considered time-varying covariates, such as the birth of another child, but one can calculate their impact for any given child by specifying a trajectory.)

### Subject-Specific Probabilities

Another approach to presenting results is to calculate survival probabilities. I will calculate probabilities of surviving to ages one and five, the complements of the infant and child mortality "rates". First let us get the baseline hazard, the width of the intervals, and thus the cumulative baseline hazard. (The current run includes a constant, so we need to take that into account.)

```. mata
------------------------------------------------- mata (type end to exit) ---------------------------
:         b = st_matrix("e(b)")

:         h = exp(b :+ b[1..5])

:         w = (1, 5, 6, 12, 36)

:         H = runningsum(w :* h)

: end
-----------------------------------------------------------------------------------------------------
```
```> h <- exp(b + c(0, b[2:5]))

> w <- c(1, 5, 6, 12, 36)

> H <- cumsum(w * h)
```

Next I will focus on a mother who is 26 years old, pretty close to the mean age, is having a second child (so we can include the length of the first interval as a predictor), and has not experienced a child death before giving birth to the second child. We consider first the case where the preceding birth interval is 3 years or longer.

The relevant parameters for age, age-squared, birth order and previous birth interval are in positions 6, 7, 8, and 13 of the vector of parameters.

```. mata xb = b*26 + b*26^2 + b*2 + b
```
```> xb <- b*26 + b*26^2 + b*2 + b
```

We can now compute the survival function for the average family with these characteristics. We will also change the previous birth interval to one year, by removing the coefficient in position 13 and adding the one in position 10:

```. mata 1 :- exp( -H * exp(xb) )
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0304345123   .0413444931   .0573727829   .0762983366   .0863384417  |
+-----------------------------------------------------------------------+

. mata 1 :- exp( -H * exp(xb - b + b) )
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0763993311   .1028877058   .1409521526   .1846057143   .2072016182  |
+-----------------------------------------------------------------------+
```
```> 1 - exp( -H * exp(xb) )
 0.03043380 0.04134353 0.05737145 0.07629658 0.08633646

> 1 - exp( -H * exp(xb - b + b) )
 0.07639758 0.10288538 0.14094905 0.18460173 0.20719721
```

So the probabilities of infant and child death for the average 26-year old mother having a second child are 5.7 and 8.6% with a three-year interval and 14.1% and 20.7% with a one-year interval.

If we consider instead mothers at the first and third quartiles of frailty we obtain

```. mata
------------------------------------------------- mata (type end to exit) ---------------------------
:         1 :- exp( -H * q * exp(xb) )
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0202495123   .0275605766   .0383528006   .0511760483    .058014883  |
+-----------------------------------------------------------------------+

:         1 :- exp( -H * q * exp(xb - b + b) )
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0512447154   .0693431838   .0956710158   .1263555039   .1424560573  |
+-----------------------------------------------------------------------+

:         1 :- exp( -H * q * exp(xb) )
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0382462283   .0518806776   .0718379424   .0952884784   .1076783111  |
+-----------------------------------------------------------------------+

:         1 :- exp( -H * q * exp(xb - b + b) )
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0954132854   .1280228624    .174442573   .2270180579   .2539465035  |
+-----------------------------------------------------------------------+

: end
-----------------------------------------------------------------------------------------------------
```
```> 1 - exp( -H * q * exp(xb) )
 0.02024903 0.02755993 0.03835190 0.05117485 0.05801353

>       1 - exp( -H * q * exp(xb - b + b ) )
 0.05124353 0.06934159 0.09566885 0.12635268 0.14245290

>       1 - exp( -H * q * exp(xb) )
 0.03824533 0.05187947 0.07183629 0.09528631 0.10767587

>       1 - exp( -H * q * exp(xb - b + b ) )
 0.09541112 0.12802001 0.17443881 0.22701329 0.25394127
```

So for a low-risk (Q1) family the probability of child death goes from 5.8 to 14.2 when we compare long and short intervals. For a high-risk (Q3) family the corresponding probability goes from 10.8 to 25.4%.

## Population-Average Probabilities

We can also calculate the average probabilities in the population of 26-year old mothers having a second birth. From the results in the notes, the (complements of) the survival probabilities under gamma frailty are

```. mata  1 :- (1 :/ (1 :+ v * H :* exp(xb))) :^(1/v)
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0303357356   .0411625024   .0570231867   .0756818763   .0855503244  |
+-----------------------------------------------------------------------+

. mata  1 :- (1 :/ (1 :+ v * H :* exp(xb - b + b))) :^(1/v)
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .0757812476   .1017714938   .1388706264    .181062986    .202757504  |
+-----------------------------------------------------------------------+
```
```> 1 - (1 / (1 + v * H * exp(xb)))^(1/v)
 0.03033502 0.04116154 0.05702187 0.07568015 0.08554838

> 1 - (1 / (1 + v * H * exp(xb - b + b)))^(1/v)
 0.07577953 0.10176922 0.13886761 0.18105915 0.20275328
```

So on average the probabilities of infant and child deaths are 5.7 and 8.6% with three-year intervals and 13.9 and 20.3% with one-year birth intervals. These are a bit lower than the corresponding probabilities for the average family with the given characteristics. (The difference is modest because by ages one and five there hasn't been a lot of time for selection to operate.)

### Marginal, Joint and Conditional Probabilities

One last calculation we can do explicitly under gamma heterogeneity involves the marginal and joint probabilities of infant and child death for two children in the same family.

We could do these calculations for 26 year-old mothers having a second birth three or more years after the first, but unless they have twins the calculation of bivariate probabilities doesn't make a lot of sense. So I will do the calculations for a second birth at age 26 and a third at age 29, everything else being the same. (A simpler approach would be to do the calculations for mothers whose observed risk factors put them at the median.)

Applying the results in the notes

```. mata xb2 = b*29 + b*29^2 + b*3 + b

. mata S1 = (1 :/ (1 :+ v * H :* exp(xb))) :^(1/v)

. mata S2 = (1 :/ (1 :+ v * H :* exp(xb2))) :^(1/v)

. mata S12 = (1 :/ (1 :+ v * H :* (exp(xb) + exp(xb2)))) :^(1/v)

. mata S1 :* S2 \ S12
1             2             3             4             5
+-----------------------------------------------------------------------+
1 |  .9392544348   .9180445945   .8874193438    .852069807    .833669972  |
2 |  .9394506446   .9184008084   .8880888478   .8532200758    .835120237  |
+-----------------------------------------------------------------------+
```
```> xb2 <- b*29 + b*29^2 + b*3 + b

> S1  <- (1 / (1 + v * H * exp(xb)))^(1/v)

> S2  <- (1 / (1 + v * H * exp(xb2)))^(1/v)

> S12 <- (1 / (1 + v * H * (exp(xb) + exp(xb2))))^(1/v)

> matrix( c(S1 * S2 , S12), 2, 5, byrow=TRUE)
[,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.9392560 0.9180467 0.8874222 0.8520734 0.8336740
[2,] 0.9394522 0.9184029 0.8880916 0.8532237 0.8351242
```

So the probability of two children surviving to age five is slightly higher than the product of the two marginals. A better way to see the correlation is to calculate a two by two table with survival to age five for two children in the same family.

```. mata M = ( 1 - S1 - S2 + S12 , S2 - S12 \ S1 - S12, S12)

. mata M
1             2
+-----------------------------+
1 |  .0090075205   .0765428038  |
2 |  .0793294387    .835120237  |
+-----------------------------+

. mata (M[1,1]/M[1,2]) / (M[2,1]/M[2,2])
1.238840858
```
```> M <- matrix( c(1 - S1 - S2 + S12 , S2 - S12, S1 - S12, S12), 2, 2)

> M
[,1]       [,2]
[1,] 0.009007067 0.07932741
[2,] 0.076541313 0.83512421

> (M[1,1]/M[1,2]) / (M[2,1]/M[2,2])
 1.23884
```

The last calculation is an odds ratio: the odds of one child dying by age five are 23.9% higher if the other child died by age five than if the other didn't!

A separate note discusses the use of log-normal frailty, which has the important advantage that it extends easily to multi-level models. Calculation of unconditional survival, however, requires numerical integration. The results are very similar to those presented here.