Flexible Estimation of Differences in Treatment-Specific Recurrent Event Means in the Presence of a Terminating Event

Summary .  In this article, we consider the setting where the event of interest can occur repeatedly for the same subject (i.e., a recurrent event; e.g., hospitalization) and may be stopped permanently by a terminating event (e.g., death). Among the different ways to model recurrent/terminal event data, the marginal mean (i.e., averaging over the survival distribution) is of primary interest from a public health or health economics perspective. Often, the difference between treatment-specific recurrent event means will not be constant over time, particularly when treatment-specific differences in survival exist. In such cases, it makes more sense to quantify treatment effect based on the cumulative difference in the recurrent event means, as opposed to the instantaneous difference in the rates. We propose a method that compares treatments by separately estimating the survival probabilities and recurrent event rates given survival, then integrating to get the mean number of events. The proposed method combines an additive model for the conditional recurrent event rate and a proportional hazards model for the terminating event hazard. The treatment effects on survival and on recurrent event rate among survivors are estimated in constructing our measure and explain the mechanism generating the difference under study. The example that motivates this research is the repeated occurrence of hospitalization among kidney transplant recipients, where the effect of expanded criteria donor (ECD) compared to non-ECD kidney transplantation on the mean number of hospitalizations is of interest.     

A Bayesian multi-risks survival (MRS) model in the presence of double censorings

Semi-competing risks data include the time to a nonterminating event and the time to a terminating event, while competing risks data include the time to more than one terminating event. Our work is motivated by a prostate cancer study, which has one nonterminating event and two terminating events with both semi-competing risks and competing risks present as well as two censoring times. In this paper, we propose a new multi-risks survival (MRS) model for this type of data. In addition, the proposed MRS model can accommodate noninformative right-censoring times for nonterminating and terminating events. Properties of the proposed MRS model are examined in detail. Theoretical and empirical results show that the estimates of the cumulative incidence function for a nonterminating event may be biased if the information on a terminating event is ignored. A Markov chain Monte Carlo sampling algorithm is also developed. Our methodology is further assessed using simulations and also an analysis of the real data from a prostate cancer study. As a result, a prostate-specific antigen velocity greater than 2.0 ng/mL per year and higher biopsy Gleason scores are positively associated with a shorter time to death due to prostate cancer.     

Exploring Migration Data Using Interval-Censored Time-to-Event Models

Abstract: Ecologists and wildlife biologists rely on periodic observation of radiocollared animals to study habitat use, survival, movement, and migration, resulting in response times (e.g., mortality and migration) known only to occur within an interval of time. We illustrate methods for analyzing interval-censored data using data on the timing of fall migration (from spring-summer-fall to winter ranges) for white-tailed deer (Odocoileus virginianus) in northern Minnesota, USA, during years 1991–1992 to 2005–2006. We compare both nonparametric and parametric methods for estimating the cumulative distribution function of migration times, and we suggest a parametric (cure rate) model that accounts for conditional (facultative) migrators as a potential alternative to traditional parametric models. Lastly, we illustrate methods for exploring the effect of environmental covariates on migration timing. Models with time-dependent covariates (snow depth, temp) were sensitive to the treatment of the data (as interval-censored or known event times), suggesting the need to account for interval-censoring when modeling the effect of these covariates.     

Dynamic prediction of cumulative incidence functions by direct binomial regression

In recent years there have been a series of advances in the field of dynamic prediction. Among those is the development of methods for dynamic prediction of the cumulative incidence function in a competing risk setting. These models enable the predictions to be updated as time progresses and more information becomes available, for example when a patient comes back for a follow‐up visit after completing a year of treatment, the risk of death, and adverse events may have changed since treatment initiation. One approach to model the cumulative incidence function in competing risks is by direct binomial regression, where right censoring of the event times is handled by inverse probability of censoring weights. We extend the approach by combining it with landmarking to enable dynamic prediction of the cumulative incidence function. The proposed models are very flexible, as they allow the covariates to have complex time‐varying effects, and we illustrate how to investigate possible time‐varying structures using Wald tests. The models are fitted using generalized estimating equations. The method is applied to bone marrow transplant data and the performance is investigated in a simulation study.     

On restricted mean time in favor of treatment

The restricted mean time in favor (RMT-IF) of treatment is a nonparametric effect size for complex life history data. It is defined as the net average time the treated spend in a more favorable state than the untreated over a prespecified time window. It generalizes the familiar restricted mean survival time (RMST) from the two-state life–death model to account for intermediate stages in disease progression. The overall estimand can be additively decomposed into stage-wise effects, with the standard RMST as a component. Alternate expressions of the overall and stage-wise estimands as integrals of the marginal survival functions for a sequence of landmark transitioning events allow them to be easily estimated by plug-in Kaplan–Meier estimators. The dynamic profile of the estimated treatment effects as a function of follow-up time can be visualized using a multilayer, cone-shaped “bouquet plot.” Simulation studies under realistic settings show that the RMT-IF meaningfully and accurately quantifies the treatment effect and outperforms traditional tests on time to the first event in statistical efficiency thanks to its fuller utilization of patient data. The new methods are illustrated on a colon cancer trial with relapse and death as outcomes and a cardiovascular trial with recurrent hospitalizations and death as outcomes. The R-package rmt implements the proposed methodology and is publicly available from the Comprehensive R Archive Network (CRAN).     

Estimating cumulative treatment effects in the presence of nonproportional hazards

Summary .   Often in medical studies of time to an event, the treatment effect is not constant over time. In the context of Cox regression modeling, the most frequent solution is to apply a model that assumes the treatment effect is either piecewise constant or varies smoothly over time, i.e., the Cox nonproportional hazards model. This approach has at least two major limitations. First, it is generally difficult to assess whether the parametric form chosen for the treatment effect is correct. Second, in the presence of nonproportional hazards, investigators are usually more interested in the cumulative than the instantaneous treatment effect (e.g., determining if and when the survival functions cross). Therefore, we propose an estimator for the aggregate treatment effect in the presence of nonproportional hazards. Our estimator is based on the treatment-specific baseline cumulative hazards estimated under a stratified Cox model. No functional form for the nonproportionality need be assumed. Asymptotic properties of the proposed estimators are derived, and the finite-sample properties are assessed in simulation studies. Pointwise and simultaneous confidence bands of the estimator can be computed. The proposed method is applied to data from a national organ failure registry.     

Desiccation resistance of two semiterrestrial isopods, Ligia exotica and Ligia taiwanensis (Crustacea) in Taiwan

The ability to resist desiccation stress was examined in two semiterrestrial Ligia species, Ligia exotica Roux and L. taiwanensis Lee, in Taiwan, under a certain desiccation condition. L. exotica exhibited the longer survival time, lower weight-specific rates of water loss, and the slightly higher ability of tolerance to water loss, compared to L. taiwanensis. In each species, the animal size displays a positive correlation to the survival time and total ability to resist desiccation, yet this size effects on the weight-specific water loss rate is negative. Neither water content nor maximum tolerance to water loss shows the association with the animal size in both species. The path ways and magnitudes of the interactions between these traits of desiccation resistance are analyzed and diagrammed using a stepwise regression model. In this model, the body sizes of animal can explain the most part of the variations in the survival time. The body size has a direct effect and an indirect effect, through the effect on water loss rate, on the time that the experimental animals can survival under this desiccated condition. These results suggest that L. exotica attains larger size than does L. taiwanensis, a lower transpiration rate and, consequently, a greater ability in desiccation resistance. The performances of these interactions in the desiccated resistance are more advantageous for L. exotica to migrate and colonize in variable land habitats within a certain limit, and as a result that L. exotica shows a broader distribution pattern than did L. taiwanensis in Taiwan.     

Simulations of clonal species genotypic diversity – trembling aspen (Populus tremuloïdes) as a case study

We built two models to follow clonal species genotypic diversity (G/N) over long periods of time at the stand and landscape levels. The models were then validated with empirical data from trembling aspen (Populus tremuloides) populations in Quebec’s boreal forest. Data was collected using a chronosequence approach in seven sites that burned in 1717, 1760, 1797, 1823, 1847, 1944, and 1916. Genetic identification was done by using four microsatellite loci. At the stand scale, simulations were repeated for a genet size of 5, 25, 50 and 100 ramets each. At the landscape level, we simulated the cumulative genet survival rate under different fire cycles (5–500 years) for 500 years after fire. Stand simulations indicated that ramet mortality within genets rather than genet mortality accounts for the increase in G/N with time since fire. Both the initial genet size and the recurrent suckering of some genets (or ramet recruitment) play an important role in maintaining high G/N levels for long periods of time. In general, the larger the number of ramets per genet, the longer the genet survives under a gap disturbance regime and a minimum of 100 ramets per genet is required to maintain aspen genet survival for 500 years. At the landscape level, genet loss increases as the fire cycle gets longer. In Quebec’s boreal forest, short rotation even-aged management practices seem to maintain a genet survival rate similar to that produced by the natural succession regime.     

Cause‐Specific Cumulative Incidence Estimation and the Fine and Gray Model Under Both Left Truncation and Right Censoring

Summary The standard estimator for the cause‐specific cumulative incidence function in a competing risks setting with left truncated and/or right censored data can be written in two alternative forms. One is a weighted empirical cumulative distribution function and the other a product‐limit estimator. This equivalence suggests an alternative view of the analysis of time‐to‐event data with left truncation and right censoring: individuals who are still at risk or experienced an earlier competing event receive weights from the censoring and truncation mechanisms. As a consequence, inference on the cumulative scale can be performed using weighted versions of standard procedures. This holds for estimation of the cause‐specific cumulative incidence function as well as for estimation of the regression parameters in the Fine and Gray proportional subdistribution hazards model. We show that, with the appropriate filtration, a martingale property holds that allows deriving asymptotic results for the proportional subdistribution hazards model in the same way as for the standard Cox proportional hazards model. Estimation of the cause‐specific cumulative incidence function and regression on the subdistribution hazard can be performed using standard software for survival analysis if the software allows for inclusion of time‐dependent weights. We show the implementation in the R statistical package. The proportional subdistribution hazards model is used to investigate the effect of calendar period as a deterministic external time varying covariate, which can be seen as a special case of left truncation, on AIDS related and non‐AIDS related cumulative mortality.     

Phenology simulation model of Scotinophara lurida (Hemiptera: Pentatomidae)

A phenology simulation model was developed for Scotinophara lurida (Burmeister). The components for the model were a degree-day immigration flight model of overwintered adults, temperature-dependent developmental models of each stage, survival rates of each stage, and an adult oviposition model. A degree-day model for immigration flight of overwintered adults was developed with blacklight trap catch data by a Weibull function. Laboratory experiments using seven constant temperature regimens were conducted to determine the effect of temperature on the development of immature stages. Developmental rates of each immature stage fit well to a linear model. Distribution of developmental time for each immature stage was successfully modeled against physiological age by a Weibull function. To determine the temperature effect on longevity, fecundity, and survival of female adults, laboratory and greenhouse experiments were conducted. The adult developmental rate (1/median longevity) was described by a linear model. The oviposition model was developed incorporating the three components of average total fecundity, cumulative oviposition rate function, and survival rate function. The simulation model predicted the time of peak occurrences of life stages of S. lurida well.     

A mechanistic model based analysis of cotton aphid population dynamics data

1 A 2‐year field study was conducted to generate data on seasonal abundance patterns of cotton aphids Aphis gossypii Glover and to develop a mechanistic model based on cumulative population size. The treatments consisted of three irrigation levels (Low, Medium and High) with 65%, 75% and 85% evapotranspiration replacement and three nitrogen fertility treatments (blanket‐rate‐N, variable‐rate‐N and no nitrogen). 2 A nonlinear regression equation, the analytical solution of a cumulative size mechanistic model, was fitted to each of the 27 individual data sets collected in 2003 and in 2004. The size and time of the peak, the cumulative aphid density, and the birth and death rates were estimated for each population, and each of these five variables was analyzed as a response variable in the analysis of variance. 3 For 2003 (a dry year), the Water (irrigation) main effect was found to be significant for the time of peak, the death rate and the cumulative density. The lower aphid death rate at low water levels might be due to the water stress in plants. 4 For 2004 (a year with moderate precipitation), the Nitrogen main effect was significant for both the birth and death rates. As nitrogen applications were increased, the decrease in both the aphid birth and death rates translates into a decrease in crowding and an increase in aphid survival. 5 The fact that treatment effects may be manifested through birth and death rate parameters in the new mechanistic model opens up new avenues for analyzing population size data of this kind.     

Competing Risks and Time‐Dependent Covariates

Time‐dependent covariates are frequently encountered in regression analysis for event history data and competing risks. They are often essential predictors, which cannot be substituted by time‐fixed covariates. This study briefly recalls the different types of time‐dependent covariates, as classified by Kalbfleisch and Prentice [The Statistical Analysis of Failure Time Data, Wiley, New York, 2002] with the intent of clarifying their role and emphasizing the limitations in standard survival models and in the competing risks setting. If random (internal) time‐dependent covariates are to be included in the modeling process, then it is still possible to estimate cause‐specific hazards but prediction of the cumulative incidences and survival probabilities based on these is no longer feasible. This article aims at providing some possible strategies for dealing with these prediction problems. In a multi‐state framework, a first approach uses internal covariates to define additional (intermediate) transient states in the competing risks model. Another approach is to apply the landmark analysis as described by van Houwelingen [Scandinavian Journal of Statistics 2007, 34 , 70–85] in order to study cumulative incidences at different subintervals of the entire study period. The final strategy is to extend the competing risks model by considering all the possible combinations between internal covariate levels and cause‐specific events as final states. In all of those proposals, it is possible to estimate the changes/differences of the cumulative risks associated with simple internal covariates. An illustrative example based on bone marrow transplant data is presented in order to compare the different methods.     

Estimation of separable direct and indirect effects in continuous time

Many research questions involve time-to-event outcomes that can be prevented from occurring due to competing events. In these settings, we must be careful about the causal interpretation of classical statistical estimands. In particular, estimands on the hazard scale, such as ratios of cause-specific or subdistribution hazards, are fundamentally hard to interpret causally. Estimands on the risk scale, such as contrasts of cumulative incidence functions, do have a clear causal interpretation, but they only capture the total effect of the treatment on the event of interest; that is, effects both through and outside of the competing event. To disentangle causal treatment effects on the event of interest and competing events, the separable direct and indirect effects were recently introduced. Here we provide new results on the estimation of direct and indirect separable effects in continuous time. In particular, we derive the nonparametric influence function in continuous time and use it to construct an estimator that has certain robustness properties. We also propose a simple estimator based on semiparametric models for the two cause-specific hazard functions. We describe the asymptotic properties of these estimators and present results from simulation studies, suggesting that the estimators behave satisfactorily in finite samples. Finally, we reanalyze the prostate cancer trial from Stensrud et al. (2020).     

Using multiple traits to estimate the effects of heat shock on the fitness of Aphidius colemani

Extreme temperature events are a great challenge for most ectotherms, particularly for the immature stages of parasitoids, as they do not possess the ability to behaviourally thermoregulate. In this study, we measured the effect of an acute heat shock, combined with desiccation stress (34 °C and 35% r.h. during 10 h) during the mummy stage on several fitness‐related traits of emerging adults of the aphid parasitoid Aphidius colemani Viereck (Hymenoptera: Braconidae: Aphidiinae). Our results showed that the emergence rate was strongly affected by the heat shock (63 ± 2.3 vs. 96.1 ± 0.7% for stress vs. control conditions, respectively), and the resulting population was male biased. Heat stress reduced the lifespan of emerged parasitoids but had no effect on fat reserves and female size. Egg load at emergence and parasitism rate were also reduced by heat treatment. The antennal symmetry was disturbed by the heat treatment, and stressed individuals had reduced mating success compared to control parasitoids. Specifically, time before wing fanning, a typical male courtship behaviour, was significantly longer in parasitoid pairs involving females that had emerged from heat‐treated mummies. Additionally, olfactometry tests on experimental females indicated that their response to host odour was similar to that of control individuals, but they required a longer time for perception. These results highlight that exposure to high temperatures during the mummy stage has cumulative detrimental effects that may strongly impact parasitoid populations under natural conditions and hinder the success of biological control programs.     

Dynamic model for predicting survival of mature larvae of Tribolium confusum during facility heat treatments

Structural heat treatment, a viable alternative to methyl bromide fumigation, involves raising the ambient temperature of food-processing facilities between 50 and 60 degrees C by using gas, electric, or steam heaters, and holding these elevated temperatures for 24 h or longer to kill stored-product insects. A dynamic model was developed to predict survival of mature larvae, which is the most heat-tolerant stage of the confused flour beetle, Tribolium confusum (Jacquelin du Val), at elevated temperatures between 46 and 60 degrees C. The model is based on two nonlinear relationships: 1) logarithmic survival of T. confusum mature larvae as a function of time, and 2) logarithmic reduction in larval survival as a function of temperature. The dynamic model was validated with nine independent data sets collected during actual facility heat treatments conducted on two separate occasions at the Kansas State University pilot flour and feed mills. The rate of increase of temperature over time varied among the nine locations where mature larvae of T. confusum were exposed, and the approximate heating rates during the entire heat treatment ranged from 1.1 to 13.2 degrees C/h. The absolute deviation in the predicted number of larvae surviving the heat treatment was within 3-7% of the actual observed data. Comparison of the absolute deviation in the time taken for equivalent larval survival showed that the model predictions were within 2-6% of the observed data. The dynamic model can be used to predict survival of mature larvae of T. confusum during heat treatments of food-processing facilities based on time-dependent temperature profiles obtained at any given location.     

Relationship between survival rates and numbers of natural teeth in an elderly Japanese population

Objectives: The objective of this study was to assess whether elderly people with 20 or more natural teeth were more likely to live longer than a cohort with less than 20 teeth. Materials and methods: Groups of elderly people over 80 years of age (24 males and 35 females) with 20 or more teeth (≥20 group) were compared with elderly people (24 males and 35 females) with less than 20 teeth (<20 group). Follow‐up studies were conducted at regular intervals for 10 years from July 1992 to July 2002. The cumulative survival rate of the ≥20 group (average ± SE tooth number of teeth – males, 23.9 ± 0.6; females, 23.8 ± 0.4) was compared with the <20 group (average number of teeth – males, 3.8 ± 1.1; females, 2.6 ± 0.8). The multivariate Cox proportional hazard models with the number of teeth in a group (≥20 group or <20 group). Smoking status and alcohol intake as covariates were used to adjust the cumulative survival rate. Results: The male participants in the ≥20 group had a significantly higher cumulative survival rates (p < 0.05) than the <20 group at 18 and 21 months from baseline. There were no significant differences in survival rates between the female groups. Adjusted cumulative survival rate was significantly different at 72, 75 and 78 months between the ≥20 group and <20 group for males but not for females. Conclusion: Having 20 or more natural teeth was associated with increased survival rate in elderly males, but not among the elderly females.     

Bayesian approaches to joint cure-rate and longitudinal models with applications to cancer vaccine trials

Complex issues arise when investigating the association between longitudinal immunologic measures and time to an event, such as time to relapse, in cancer vaccine trials. Unlike many clinical trials, we may encounter patients who are cured and no longer susceptible to the time-to-event endpoint. If there are cured patients in the population, there is a plateau in the survival function, S(t), after sufficient follow-up. If we want to determine the association between the longitudinal measure and the time-to-event in the presence of cure, existing methods for jointly modeling longitudinal and survival data would be inappropriate, since they do not account for the plateau in the survival function. The nature of the longitudinal data in cancer vaccine trials is also unique, as many patients may not exhibit an immune response to vaccination at varying time points throughout the trial. We present a new joint model for longitudinal and survival data that accounts both for the possibility that a subject is cured and for the unique nature of the longitudinal data. An example is presented from a cancer vaccine clinical trial.     

Merits of Modelling Multivariate Survival Data Using Random Effects Proportional Odds Model

Recently, there has been a great deal of interest in the analysis of multivariate survival data. In most epidemiological studies, survival times of the same cluster are related because of some unobserved risk factors such as the environmental or genetic factors. Therefore, modelling of dependence between events of correlated individuals is required to ensure a correct inference on the effects of treatments or covariates on the survival times. In the past decades, extension of proportional hazards model has been widely considered for modelling multivariate survival data by incorporating a random effect which acts multiplicatively on the hazard function. In this article, we consider the proportional odds model, which is an alternative to the proportional hazards model at which the hazard ratio between individuals converges to unity eventually. This is a reasonable property particularly when the treatment effect fades out gradually and the homogeneity of the population increases over time. The objective of this paper is to assess the influence of the random effect on the within‐subject correlation and the population heterogeneity. We are particularly interested in the properties of the proportional odds model with univariate random effect and correlated random effect. The correlations between survival times are derived explicitly for both choices of mixing distributions and are shown to be independent of the covariates. The time path of the odds function among the survivors are also examined to study the effect of the choice of mixing distribution. Modelling multivariate survival data using a univariate mixing distribution may be inadequate as the random effect not only characterises the dependence of the survival times, but also the conditional heterogeneity among the survivors. A robust estimate for the correlation of the logarithm of the survival times within a cluster is obtained disregarding the choice of the mixing distributions. The sensitivity of the estimate of the regression parameter under a misspecification of the mixing distribution is studied through simulation. (© 2004 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

A Model for Germination Rate during Dormancy Loss in Hordeum vulgare

A quantitative model for change in germination rate of barley(Hordeum vulgare L.) during and after loss of primary dormancyis presented. Change in mean germination time on a logarithmicscale is normally distributed within the period of after-ripeningand the standard deviation of this distribution is shown tobe a quantitative function of after-ripening temperature. Therate of change of mean germination time is in inverse proportionto the product of the standard deviation and a parameter whichis characteristic of the seed population. The latter parameteris the rate constant for change in probit cumulative germinationas a negative linear function of the logarithm of mean germinationtime. A model based solely on dormancy loss is combined withan existing model of change in probit viability as a functionof mean germination time to produce a model which predicts thetime to and optimum value of mean germination rate of a populationas it after-ripens. The model provides a quantitative link betweenthe effect of pre-germination and germination environments ontotal and rate of germination of an initially dormant population.Experimental data from dormant barley (cv. Triumph) stored at27, 38, 45, 50 and 60 °C, and germinated at 18 °C wereused to validate the model. The data show that germination ratecontinues to increase after primary dormancy is lost until itreaches an upper limit determined by the intrinsic germinativevigour of the seed lot. Rate of loss of primary dormancy andincrease in germination rate thus appear to be quantitativelylinked as a function of after-ripening temperature and factorswhich may be specific to the mode of induction of dormancy withina seed lot prior to harvest.Copyright 1995, 1999 Academic Press Hordeum vulgare L. barley, germination rate model, dormancy, vigour, after-ripening temperature     

Corridors maintain species richness in the fragmented landscapes of a microecosystem

Theory predicts that species richness or single-species populations can be maintained, or at least extinctions minimized, by boosting rates of immigration. One possible way of achieving this is by establishing corridors of suitable habitat between reserves. Using moss patches as model microecosystems, we provide here probably the first field experimental test of the idea that corridors can reduce the rate of loss of species, and therefore help to maintain species richness. Connecting patches of habitat with corridors did indeed slow the rate of extinction of species, preserving species richness for longer periods of time than in disconnected habitat patches. The pattern of γ-diversity, the cumulative species richness of entire connected systems, is similarly higher than that of fragmented systems, despite the homogenizing effects of movement. Predators are predicted to be more susceptible to fragmentation because of their greater mobility and smaller population sizes. Our data are consistent with this prediction: the proportion of predator species declined significantly in disconnected as compared with connected treatments.