Germán Rodríguez

Survival Analysis
Princeton University
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 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[2] :- 1 1 2 3 +----------------------------------------------+ 1 | -.2879554288 0 .3573350829 | +----------------------------------------------+

> p <- (1:3)/4 > q <- qgamma(p, shape = 1/v, scale = v); q [1] 0.6618945 0.9295689 1.2617365 > q / q[2] - 1 [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.

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 <- read.dta("https://data,princeton.edu/pop509/PebleyStupp2.dta") > 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[2] - 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.)

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[17] :+ b[1..5]) : w = (1, 5, 6, 12, 36) : H = runningsum(w :* h) : end -----------------------------------------------------------------------------------------------------

> h <- exp(b[1] + 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[6]*26 + b[7]*26^2 + b[8]*2 + b[13]

> xb <- b[6]*26 + b[7]*26^2 + b[8]*2 + b[13]

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[13] + b[10]) ) 1 2 3 4 5 +-----------------------------------------------------------------------+ 1 | .0763993311 .1028877058 .1409521526 .1846057143 .2072016182 | +-----------------------------------------------------------------------+

> 1 - exp( -H * exp(xb) ) [1] 0.03043380 0.04134353 0.05737145 0.07629658 0.08633646 > 1 - exp( -H * exp(xb - b[13] + b[10]) ) [1] 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[1] * exp(xb) ) 1 2 3 4 5 +-----------------------------------------------------------------------+ 1 | .0202495123 .0275605766 .0383528006 .0511760483 .058014883 | +-----------------------------------------------------------------------+ : 1 :- exp( -H * q[1] * exp(xb - b[13] + b[10]) ) 1 2 3 4 5 +-----------------------------------------------------------------------+ 1 | .0512447154 .0693431838 .0956710158 .1263555039 .1424560573 | +-----------------------------------------------------------------------+ : 1 :- exp( -H * q[3] * exp(xb) ) 1 2 3 4 5 +-----------------------------------------------------------------------+ 1 | .0382462283 .0518806776 .0718379424 .0952884784 .1076783111 | +-----------------------------------------------------------------------+ : 1 :- exp( -H * q[3] * exp(xb - b[13] + b[10]) ) 1 2 3 4 5 +-----------------------------------------------------------------------+ 1 | .0954132854 .1280228624 .174442573 .2270180579 .2539465035 | +-----------------------------------------------------------------------+ : end -----------------------------------------------------------------------------------------------------

> 1 - exp( -H * q[1] * exp(xb) ) [1] 0.02024903 0.02755993 0.03835190 0.05117485 0.05801353 > 1 - exp( -H * q[1] * exp(xb - b[13] + b[10] ) ) [1] 0.05124353 0.06934159 0.09566885 0.12635268 0.14245290 > 1 - exp( -H * q[3] * exp(xb) ) [1] 0.03824533 0.05187947 0.07183629 0.09528631 0.10767587 > 1 - exp( -H * q[3] * exp(xb - b[13] + b[10] ) ) [1] 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%.

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[13] + b[10]))) :^(1/v) 1 2 3 4 5 +-----------------------------------------------------------------------+ 1 | .0757812476 .1017714938 .1388706264 .181062986 .202757504 | +-----------------------------------------------------------------------+

> 1 - (1 / (1 + v * H * exp(xb)))^(1/v) [1] 0.03033502 0.04116154 0.05702187 0.07568015 0.08554838 > 1 - (1 / (1 + v * H * exp(xb - b[13] + b[10])))^(1/v) [1] 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.)

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[6]*29 + b[7]*29^2 + b[8]*3 + b[13] . 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[6]*29 + b[7]*29^2 + b[8]*3 + b[13] > 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[5] - S2[5] + S12[5] , S2[5] - S12[5] \ S1[5] - S12[5], S12[5]) . 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[5] - S2[5] + S12[5] , S2[5] - S12[5], S1[5] - S12[5], S12[5]), 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] 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.