Abstract
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.
1. 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–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.
2. 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 . 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 iswhere , , 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 HRwith expected value .
The individual survival function is , where 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 and , the population survival function becomes
The population hazard rate is then derived as and the population HR can be expressed asa decreasing function of time, starting at when and declining to when , which is the long run population excess mortality. The decline towards is steeper for larger than for smaller .
The random HR in (2) is a constant function of time for each widowed person, with expected value . The population HR in (6) decreases towards over time . We have chosen a least squares’ (LS) approach to demonstrate estimation of , , and . The minimized sum of squares is , 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).
3. 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 , and 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 , where he/she is one of the less frail if or more frail if , such that the lowest possible individual mortality is . This limit is also the long run population HR. The size of decides how rapid the population HR declines towards . A large variance 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 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, is significantly positive, indicating a heterogeneous widowhood mortality. An 85+ years old widowed man is expected to have double risk of dying () compared to a married man of the same age, but according to the model it is in fact a mixture of low and high risk widowers.
Estimated standard deviations are between 15 and 92, and the declining factor in population HR one year after spousal death, , is between 0.014 and 0.466.
Survival times from model (4) are simulated for each sex and age group, with similar sample sizes as the widowhood data in Table 1. Simulations indicate that the LS estimator is unbiased and has a very good precision (Table 3). However, is overestimated for the smallest simulated value , and the confidence intervals (CIs) are quite wide in all cases. The bootstrap CIs for 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.
4. 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 , the long run population excess mortality for widowed compared to married, is accurate and has good precision. The frailty parameter estimate, , is less accurate and precise, but simulations show that the bootstrap CI is a good indicator of whether the frailty term is nonzero or not.
Data Availability
The data is not available to the public; however, monthly population HRs are provided in Supplementary Table 1.
Conflicts of Interest
The author declares that there are no conflicts of interest.
Acknowledgments
Thanks are due to Tormod Brenn for helpful suggestions.
Supplementary Materials
The marital status survival data is not available to the public; however, the population HRs for widowed compared to married individuals are provided in Supplementary Table 1. The time scale in the data is days since the loss of a spouse, here summarized as monthly changes in the population HR by age group and sex. Supplementary Table 2 contains an R-code for the LS estimation of in the parametric survival function . The sum of squares to be minimized is , where is the Nelson-Aalen survival function. (Supplementary Materials)