Frailty in Survival Analysis of Widowhood Mortality

Heterogeneity between individuals has attracted attention in the literature of survival analysis for several decades. Widowed individuals also differ; some are more frail than others and thereby have a higher risk of dying. The traditional hazard rate in a survival model is a measure of population risk and does not provide direct information on the unobservable individual risk. A frailty model is developed and applied on a large Norwegian data set of 452 788 widowed individuals. The model seemed to fit the data well, for both widowers and widows in all age groups. The random frailty term in the model is significant, meaning that widowed persons may have individual hazard rates.


Introduction
In survival analysis of time to death, or to any event, the hazard rate is a measure of the risk in the population at large.However, some individuals are more frail and thus have a higher risk than others.Individuals with a high risk are more likely to die early and the remaining individuals are mainly those with lower risk.This will give a "pulled down" population hazard rate as time passes [1, chapter 6].It is therefore important to keep in mind that the observable population risk is not a measure of the individual risk.Some of the individual variation can be expressed by covariates and then included in the model.The unobservable differences between individuals are usually considered as random variation, but in survival analysis, observing over time, this individual variation cannot be ignored.
Several studies of mortality after the death of a spouse, summarized in a meta-analysis [2], demonstrated that the population hazard rate is substantially elevated immediately after the death of the spouse, before gradually declining.This is also the case after controlling for age, sex, and other covariates.Thus, the hazard rate for the widowed population is observed as being "pulled down" by time in widowhood.Many of the studies find this "pulled down" pattern to fade out at a substantial higher level than that of a married population [3][4][5][6][7].
The unobservable individual hazard rate can be modeled as a random variable from which the population hazard is derived.Such models are called frailty models [8], first suggested by Vaupel, Manton and Stellard [9] in mortality studies, and by Lancaster [10] in a study of unemployment duration.The mathematical properties of frailty models are thoroughly discussed in [11].
The aim of this study is to develop a frailty model to fit the characteristic widowhood mortality.Properties are discussed, and the model is applied on a large Norwegian data set.

A Frailty Model
A commonly used and simple model for the unobserved heterogeneity between individuals is the individual proportional hazard rate: (), where () is a basic rate and  a gamma distributed random variable with expectation () = 1.It is shown by [11] that this frailty effect will "pull down" the population hazard rate over time.Widowed mortality is previously shown to be elevated immediately after the death of the spouse and thereafter decrease, but still remaining higher than the married mortality several years into widowhood [7].Assuming constant individual hazard ratio (HR) compared to married, but with the characteristic "pulled down" shaped population HR, our proposed individual hazard rate is where  0 > 0,  1 ≥ 0, and   () is a known mortality rate for a married population of same sex and age as the widowed.This will give us a constant individual HR with expected value is the known cumulative mortality rate for married population.
The population survival function is found by integrating over the distribution of , after replacing the individual married rates   () by an average rate () for the group at risk The expectation on the right hand side is the Laplace transform of , and for a gamma distributed  with E() = 1 and var(Z) =  2 , the population survival function becomes The population hazard rate () is then derived as and the population HR can be expressed as a decreasing function of time, starting at  0 +  1 when  = 0 and declining to  0 when  → ∞, which is the long run population excess mortality.The decline towards  0 is steeper for larger  2 than for smaller  2 .The random HR in (2) is a constant function of time for each widowed person, with expected value  0 +  1 .The population HR in (6) decreases towards  0 over time .We have chosen a least squares' (LS) approach to demonstrate estimation of  0 ,  1 , and  2 .The minimized sum of squares is ∑  [ Ŝ() − ()] 2 , with the Nelson-Aalen survival estimator Ŝ() and the survival function () in (4).Bootstrapping and simulation results are also included.The population HR in ( 6) is then estimated using these LS estimates.Comparing other estimation methods is beyond the scope of this paper.
Many of the common frailty models in the literature have a multiplicative random element in the hazard rate, and a variety of estimation procedures are available in standard statistical software [12].However, to the author's knowledge, there is no software available for our model (1).

Data Example
The population of married individuals in Norway on January 1, 1975, were followed until the end of 2006 for changes in marital status and date of death.Data were provided by Statistics Norway and are described in detail elsewhere [7].Individuals dying on the same day as their spouse are not recorded as a widower or widow.During the follow-up period, 452 807 individuals became widowed at the age of 55 to 94 years.We further excluded 19 widowed individuals dying within three weeks after their spouse and of a similar external cause, assuming their death was caused by the same accident as their spouse.The remaining sample containing 452 788 widowers and widows is used here for estimation of the frailty model.
The parameters  0 ,  1 , and  2 in model ( 4) are estimated separately for widowers and widows in four age groups.The total numbers in each group are presented in Table 1.Time  is days after spousal death.The LS estimations are performed in the statistical computing language R [13], and the R-code is available in Supplementary Table 2.The cumulative mortality rates for married () are calculated from life tables given in [7].
The individual HR in (2) is random but constant over time.A particular widowed individual  is thus assumed to have a constant HR  =  0 +    1 , where he/she is one of the less frail if 0 <   < 1 or more frail if   > 1, such that the lowest possible individual mortality is HR  =  0 .This limit is also the long run population HR.The size of  2 decides how rapid the population HR declines towards  0 .A large variance  2 indicates a high amount of very frail widowed persons, who are expected to die early, resulting in a rapid decline in population HR.
The observed population HR is presented in Figure 1 as monthly measurements Δ Â()/Δ(), where Δ Â() are the increments of the Nelson-Aalen estimator and Δ() the increments of the cumulative mortality rate from life tables for a married population.The Norwegian data has a "pulled down" shaped population HR, and the frailty model ( 1) seems to fit this shape quite well.
Estimates show that  0 is significantly greater than one in all age groups and both sexes, meaning that even the least frail widowed individual has higher mortality than a married of same age and sex (Table 2).Furthermore,  1 is significantly   1. Simulations indicate that the LS estimator β0 is unbiased and has a very good precision (Table 3).However,  1 is overestimated for the smallest simulated value 0.4, and the confidence intervals (CIs) are quite wide in all cases.The bootstrap CIs for β1 are also wide and of the same magnitude as the simulation CIs.
Our frailty model with time constant individual HR seems to fit the widowhood data well.However, there may be time dependent bereavement mechanisms causing the elevated and decreasing population HR, such as loss of income, social support, grief, and emotional stress [14].It is not possible to conclude on what is causing the "pulled down" population HR, without more information from each widowed individual.

Conclusion
The characteristic decreasing population HR over time since spousal death can be explained as heterogeneity in individual mortality.A gamma frailty model is fitted to a large Norwegian data set, and it is demonstrated how this model fits the observed "pulled down" population HR.According to the particular frailty model we use here, each person has a constant risk of dying at any time in widowhood, except for the effect of aging.However, some widowed have higher risk than others.High risk persons will eventually die sooner than those who are less frail, causing the population risk to decrease.
The simulation study shows that, for our particular frailty model in (1), the LS estimator for  0 , the long run population excess mortality for widowed compared to married, is accurate and has good precision.The frailty parameter estimate, β1 , is less accurate and precise, but simulations show that the bootstrap CI is a good indicator of whether the frailty term  1 is nonzero or not.

Table 1 :
Number of widowed men and women.Norway 1975-2006.

Table 2 :
Estimated parameters with Bootstrap confidence intervals (CI), 1000 repetitions.(365) is the cumulative mortality rate for the married at time  = 365.

Table 3 :
Least square estimates based on simulated survival times from () in (4).Averaged estimates and confidence intervals (CI) from 1000 repetitions in two groups: men age 55-64 and women age 75-84.