BivRec: an R package for the nonparametric and semiparametric analysis of bivariate alternating recurrent events.

BACKGROUND: Bivariate alternating recurrent event data can arise in longitudinal studies where patients with chronic diseases go through two states that occur repeatedly, e.g., care periods and break periods. However, there was no statistical software that provided tools for the analysis of such data. To meet this software need, we developed BivRec, a package for R that contains a set of tools for exploratory, nonparametric and semiparametric regression analysis of bivariate alternating recurrent events.
RESULTS: The BivRec package provides functions for nonparametric estimations for the joint distribution of bivariate gap times (bivrecNP) and semiparametric regression methods for evaluating covariate effects on the two types of gap times under the accelerated failure time model framework (bivrecReg). The package also provides exploratory data analysis tools such as a visualization of the gap times by groups. We utilize a subset of the South Verona Psychiatric Case Register (PCR) data to illustrate the use of the BivRec package for the reviewed methods.
CONCLUSIONS: We demonstrate BivRec's capability for data visualization, nonparametric and regression based analysis, as well as data simulation. The package has default methods with satisfactory performance despite the complexity of calculations and fills a gap in software for statistical analysis of bivariate alternating recurrent events. BivRec is accessible under the GPL-3 General Public License through CRAN, facilitating its installation.

Background

Bivariate alternating recurrent event data can arise in longitudinal studies where individuals go through two interchanging states that reoccur over time. For instance, patients with chronic conditions, such as addiction or mental illnesses, experience care periods and break periods repeatedly. The on and off periods of care or disease form a sequence of recurrent gap times. Consider a group of patients who enter a study at the time of a hospitalization due to a particular disease and during the study experience a series of recurrent hospitalizations caused by the same disease. In this case, patients experience alternating periods of care and breaks from care (referred to as Type I and Type II gap times, respectively, hereinafter) until a censoring event such as the end of study occurs.

First, we introduce some notation. Let and be random variables representing the length of Type I and Type II gap times of episode j experienced by subject i, respectively. Denote the collection of all episodes of subject i by . In addition, we assume that the bivariate recurrent event process of subject i is subject to right censoring C, which has survival function G(·) with maximum support τ= sup{t:G(t)>0}. Let m denote the number of episodes of bivariate alternating recurrence times for subject i, which satisfies the conditions: Note that, in the last bivariate pair , may or may not be censored, but is always censored, for which we use the hospitalization example to illustrate. When the censoring occurs during the last care period, , the break period afterwards, cannot be observed, and when the censoring occurs during the last break period, is only partially observed.

The observation of the gap times and is subject to the censoring time and , respectively, where . Due to the censoring, the observed data for subject i are , where and for j=1,…,m−1, while and . Figure 1 shows an illustration of this process where the Type I gap of the last episode is observed but the Type II gap following it is censored (i.e., ).

Fig. 1

Illustration of a bivariate alternating recurrent event process

The BivRec package was designed to analyze bivariate alternating recurrent data with the form as depicted in Fig. 1. It provides a consistent and user-friendly set of functions to: explore and visualize the data, estimate and plot the joint cumulative distribution function (cdf), the marginal survival and the conditional cdf using nonparametric methods [1] and fit semiparametric accelerated failure time models [2, 3] to estimate the effect of covariates on the two alternating gap times.

BivRec version 1.2.1 is available through CRAN at https://cran.r-project.org/package=BivRec. The reference manual can be found in both CRAN and Additional file 1 of this paper. The implementation in R was designed for users with some experience analyzing survival data and follows conventions used in similar R packages. Most functions in BivRec are S4 methods that produce S4 class objects. Where possible functions have been optimized using Fortran 90 to reduce running time.

Implementation

Nonparametric analysis

In 2005, Huang and Wang [1] developed a nonparametric method for estimating the joint distribution of the two types of alternating gap times, which is a useful data summary tool for bivariate alternating recurrent event data. In the BivRec package, we implemented Huang and Wang’s (2005) [1] nonparametric methods in the function bivrecNP() to estimate:

the joint cumulative distribution function (cdf) for the two types of gap times, , and its associated standard error (in an output data frame and a contour plot),

the marginal survival function for the Type I gap times, , and its associated standard error (in an output data frame and a survival plot), and

the conditional cdf for the Type II gap times given the Type I gap times, , and its associated standard error (in an output data frame and a conditional cdf plot).

For estimation for the cdf, assume that there exists a subject-level latent variable Ψ with an unspecified cdf, P(·) such that the bivariate gap times are identically and independently distributed (i.i.d.) given Ψ and that the censoring time C is independent of (N,Ψ). Define variables and , then their joint cdf is for z=x+y and =(x,y). The marginal survival function of is hence . Our interest lies in the estimation of the joint distribution , which is determined by through the identity . Let F(z,)=E[aI(Z≤z,≤,Δ=1)] and R(z)=E[aI(Z≥z)], where the weight a=a(C) is a non-negative function of C and satisfies , with a special case of a=1, i.e., no weights. Then, following Huang and Louis (1998) [4], it is shown that , where is the marginal survival function of [1]. Noticing that F(ds,) and R(s) can be replaced with their respective empirical estimators and that the survival function can be estimated by the estimator of [5] for univariate (i.e., single-type) recurrent gap times , one can estimate the joint distribution of interest for any (x,y) satisfying x+y≤τ as follows:

where are the distinct and uncensored recurrence times from with for m≥2 and for m=1, and .

Next, we focus on the estimation of the marginal survival function for Type I gap times, . Let denote the distinct and uncensored Type I gap times from . Following Wang and Chang (1999) [5], Huang and Wang in 2005 [1] proposed to estimate for x≤τ with: where and Finally, they argued that the marginal distribution of Type II gap times is not estimable due to induced dependent censoring, and hence proposed to estimate the conditional distribution of Type II gap times given Type I gap times, using the estimator

Similarly, one can estimate by for x1+y≤τ.

The standard errors of the estimators for the joint distribution and the marginal survival function of Type I gap times are estimated based on the large sample properties of these estimators proved in [1]. Briefly, for 0≤x+y≤L,L<τ<τ with τ being the maximal support of R(z), weakly converges to a Gaussian process with mean zero and variance-covariance function Σ=E[ϕ1(z1,1)ϕ1(z2,2)]. Similarly, converges weakly to a Gaussian process with mean zero and with variance-covariance function where x1,x2∈[0,L]. The definitions of ϕ1(z,) and ξ1(x) can be found in [1]. For the conditional distribution estimator, , the package provides the bootstrap standard error and confidence intervals.

Semiparametric regression

Researchers, especially in clinical settings, are often more interested in regression methods that allow them to understand the relationship between covariates and the recurrent event process. In this regard, Chang in 2004 [2] proposed a semiparametric accelerated failure time (AFT) model which allows the estimation of the covariate effects on the two types of alternating gap times to be done simultaneously (referred to as Chang’s method hereinafter). Recognizing that the estimation of the AFT model coefficients by Chang is based on a nonsmooth, rank-based estimating function, Lee et al. in 2018 [3] proposed a smooth, U-statistic-based estimating function whose solution is found to be more computationally tractable (referred to as Lee et al.’s method hereinafter). We now briefly review the AFT model for bivariate recurrent gap times and the estimation methods developed by [2] and [3] and implemented in the function bivrecReg().

Let denote a p×1 vector of baseline covariates of subject i and =(γ,γ)′, a subject-specific latent vector that carries information for within subject correlations among the recurrent gap times. The censoring time C is assumed to be independent of (N,,). The AFT model assumes that conditioning on and , the bivariate gap time pairs are i.i.d. within subject i. Furthermore, each (log) gap time is linearly related to the covariates as follows:

where 1,2 denote the regression coefficients for Type I and Type II gap times, respectively, and ε,k=1,2, are mutually independent random errors with mean zero. Both the errors and the latent vectors come from unspecified distributions.

In 2004, Chang [2] considered the transformed, complete gap times given by and , where . Their observed counterparts are:

The rank-based estimating functions [2] are:

where and for k=0,1 and for k=0; for k=1. Denote the solution to U1()=0 and U2()=0 by . Chang (2004) proposed to use the resampling method in [6] to estimate the covariance matrix for . Since both the point estimate and the resampling-based interval estimate rely on solving nonsmooth estimating functions, fitting the AFT model with Chang’s method can be computationally inefficient and even encounter a nonconvergence problem, in which case the R function will give the user an error message.

Motivated by a multi-state model [7], Lee et al. [3] defined the transformed gap times and with . Their observed counterparts are: The authors use O(·,·), a symmetric, continuous function on {(t,s):0≤t≤L,0≤s≤L} such that O(s,t) is monotonic in t if s is given and vice versa to derive the following U-statistic-based estimating equations: where and are Kaplan-Meier estimators of the survival function of the censoring time G(·) based on the data and , respectively, and L1<τ and L2<τ are limits that ensure respect for the support of G(·). The bivariate functions in the estimating functions, and are not necessarily the same. Current implementations used L1=L2=L and O(s,t)= log[ min{max(t,s),L}]− log(L). The regression coefficient estimators, and can be obtained by inductively solving and . Note that the smooth and monotonic nature of the estimating equations guarantees a unique solution, a property not possessed by Chang’s method. Moreover, Lee et al. proved the weak convergence of to a mean zero normal distribution with variance that can be consistently estimated by where the definitions of and can be found in [3].

Results

We use a subset from the South Verona Psychiatric Case Register (PCR) [8] to illustrate functions in BivRec for data exploration, visualization and analysis in R. Our PCR sample contains data on patients’ care and break periods and disease-related and socioeconomic factors such as age at disease onset, education level, and sex from 336 patients with schizophrenia or related disorders with their conditions first recorded between 1981 and 1995 in South Verona, Italy. We focus on two covariates, one categorical, EDU, and one continuous, Age10, that were previously studied in [3]. These correspond to the education level and age at onset (in 10 years), respectively. We also show a simulated data set and how to simulate data with a function in BivRec in a later section.

Data preparation

Even though bivariate alternating recurrent data may be displayed in a wide or long format, in line with various longitudinal methods, the BivRec package requires that the data is in a long format, with possibly multiple rows for each participant, reflecting the number of episodes that the participant experienced. In addition to a set of baseline covariates (which repeats for each row if a participant has more than one row in the data), the long format data should have six columns corresponding to , and which are defined in the previous section and can be specified in the statements of a data object function bivrecSurv() by “id=”, “episode=”, “xij=”, “yij=”, “d1=”, and “d2=”, respectively. See Table 1 for the detailed definitions of these arguments. All the functions for data exploration and analysis in BivRec use the data object created with this function. We used the PCR data to create a bivrecSurv() object as follows (see the package manual [9] for further details):

Table 1

Arguments and compatible standard functions for function bivrecSurv()

Argument	Description
id	Vector of subject’s unique identifier.
episode	Vector indicating the pair or episode number (j) for a subject (i); this will determine order of events for each subject.
xij	Vector with the lengths of time spent in event of Type I for individual i in episode j.
yij	Vector with the lengths of time spent in event of Type II for individual i in episode j.
d1	Vector of censoring indicator corresponding to Type I gap times (xij); d1 = 1 for uncensored, and = 0 for censored gap times.
d2	Vector of censoring indicator corresponding to Type II gap times (yij); d2 = 1 for uncensored, and = 0 for censored gap times. Note that in the last episode, yij is always censored (i.e., d2 = 0).
Compatible functions:	plot()

Note that in practice, the input data may have quality issues such as unequal length of variables, negative values for the gap time variables (xij and yij), gaps or non-integers in the episode variable within a subject (e.g., j=1,3,4 or j= a, b, c), unreasonable values in the censoring indicators (e.g., (d1, d2) = (0, 1) for any episode or (d1, d2) = (1, 1) for the last episode), or no subjects having any uncensored episodes observed (i.e., all d2 = 0). In these cases, the user will get an error message such as “Error: Data not cleaned” with possibly more details to help pinpoint the problem. However, missing values are allowed in the input data even though only subjects with complete data will be used in any subsequent analysis. It is worthwhile to mention that in applications, the time variables may not be observed continuously (e.g., in days); for example, a subject’s censoring event could occur on the same day as the last observed event, causing the situation of (d1, d2) = (1, 1) for the last episode of events (j=m). For these subjects, we suggest users to add a small quantity to the censoring time during the data cleaning process.

Data exploration

We begin the exploration of this data set by obtaining a visualization of the care and break periods (Fig. 2) using the plot() function on a bivrecSurv object in the following way:

Fig. 2

Care and break periods in the South-Verona Psychiatric Case Register (PCR) data sorted by the overall follow-up time of each individual

The resulting plot shows that during the study, the majority of the observation time was spent out of the hospital or care facility, as accentuated by the larger portions of blue in the graph.

The data can also be viewed in subgroups defined by a categorical covariate of six or less levels if a user specifies such a covariate using the by argument of plot(). An example of this feature is shown by looking at times in and out of care based on education levels in the PCR data set. Figure 3 shows the resulting subgroup plots for the levels of the variable EDU, an indicator with values of one for participants with secondary or higher education levels and zero otherwise. Note that in this example, the covariate Age10 is automatically dropped from the by statement, since the function detects that Age10 is possibly a continuous variable. In addition, 10 subjects had missing values for the EDU variable, which is reflected in a message letting the user know only 326 of the subjects were used instead of the full 336 sample.

Fig. 3

Care and break periods in the PCR data stratified by education (1 = secondary education or higher; 0 = less than secondary education)

The “bluer” shade of the left panel compared to the one on the right in Fig. 3 indicates that patients with higher education level might have longer break periods than those with less education. Thus, considering education as a factor that affects the re-hospitalization process is reasonable. Note that more than one categorical covariate can be specified in the by argument, but subgroup plots are created for each covariate separately. If subgroups defined by the combinations of multiple categorical variables (e.g., sex and race) are desired, users can define a new categorical variable (with ≤6 levels) before applying the plot() function.

We use the function bivrecNP() to estimate the joint distribution for all combinations of values given by the options u1 and u2. The non-negative weight in the nonparametric estimator in Equation (1), a=a(C), is specified by the option ai. For the data example, we use a simple unit weight a(C)=1 by setting ai=1 (the default). If one sets ai=2, the weight will be the censoring time of each subject, a(C)=C. The function will automatically estimate the marginal survival probability for Type I gap times for all distinct and uncensored recurrence times. We set conditional=TRUE to request the conditional distribution estimate for Type II gap times given that Type I gap times fall into a certain interval specified by the given.interval argument. In the following example, we set given.interval=c(100, 500) to estimate , and the confidence level of the point-wise confidence intervals is set by level=0.99 as an illustration that the confidence level can be changed to values different than 0.95. Table 2 and the package manual [9] provide further details of the function bivrecNP().

Table 2

Arguments and compatible standard functions for function bivrecNP()

Argument	Description
response	A response object of the bivrecSurv class.
level	The confidence level for the point-wise confidence interval; must be between 0.50 and 0.99; the default value is 0.95.
ai	Value 1 or 2 to indicate which weight function to use in the nonparametric estimator; 1 indicates that the weights are 1 for all subjects, a(Ci)=1 (default); 2 indicates that the weight is the subject’s censoring time, a(Ci)=Ci.
u1	A vector (or single number) of time values to be used for the estimation of the joint cdf, Pr(X0≤u1,Y0≤u2).
u2	A vector (or single number) of time values to be used for the estimation of the joint cdf, Pr(X0≤u1,Y0≤u2).
conditional	A logical value. If TRUE, this function will calculate the conditional cdf for the Type II gap time given an interval of the Type I gap time and the bootstrap standard error and confidence interval at the specified confidence level; the default is FALSE.
given.interval	A vector c(v1, v2) that must be specified if conditional=TRUE. The vector indicates an interval for the Type I gap time to be used for the estimation of the cdf of the Type II gap time given this interval.
	If given.interval=c(v1, v2), the function calculates Pr(Y0≤y|v1≤X0≤v2). The given values v1 and v2 must be in the range of gap times in the estimated marginal survival.
Compatible functions:	plot(), head(), print()

Using the function head() we show a snapshot of all the output elements. To see the entire data frame for all the output elements, use the R function print(). The individual output elements (joint_cdf, marginal_survival, and conditional_cdf) can also be retrieved using the $ operator following the output object.

We use the plot() function on the resulting object to simultaneously generate plots for all the estimated distributions. To individually plot the joint cdf, the marginal survival probability and the conditional cdf, one can use the plotJoint, plotMarg and plotCond functions, respectively.

Figure 4 is the contour plot of the joint cdf in the half-plane where estimations meet the condition x+y≤τ, for τ=5697. The joint cdf shows that, for instance, the probability that participants were in care for 2000 days or less and out of care for 3000 days or less is over 0.70. Note that the choices of u1 and u2 in the code will determine if the plot of the joint distribution is drawn on a half plane or not. For instance, if we had coded u1=seq(100, 2000, 50) and u2=seq(100, 3000, 50), then all the (x,y) pairs defined by u1 and u2 would be estimable because max(u1)+ max(u2)≤τ. In addition, the resulting object (npresult in this example) will show values of NA for the joint probability for any combination of u1 and u2 that does not meet the x+y≤τ condition.

Fig. 4

Contour plot of the joint cdf of the care and break periods for combinations that meet the condition x+y≤τ, where τ=5697

The line plots of the estimated marginal survival probability and conditional cdf (with corresponding 99% confidence intervals in this example) are shown in Fig. 5. The left panel shows that care times shorter than a year and a half are the most likely as the survival probability drops sharply before 500 days. With this in mind, we look at the conditional cdf plot for patients who received between 100 and 500 days of care and conclude that the probability that those patients spent 3000 or fewer days out of care was close to 80%. The bivrecNP() function with conditional=TRUE specification for the PCR data produced results in 248.39 seconds (4.14 minutes, using an AMD Ryzen 5 3500U with Radeon Vega Mobile Gfx 2.10 GHz processor with 16 GB RAM for this and other analyses in this paper). Note that the computing time would be shorter if a less fine grid was specified for u1 and u2. The visualizations for the results, shown in Figs. 4 and 5 ran in less than 2 seconds each.

Fig. 5

Plots of the marginal survival probability of the care period and the conditional cdf of the break period given the care period

The function bivrecReg() provides point estimates, standard error estimates, an estimated variance-covariance matrix and confidence intervals at a specified confidence level for the effects of covariate(s) on the two types of gap times, based on the Chang’s method and Lee et al.’s method. Using the PCR data, we show how to estimate the effects of the education and age of disease onset on the two types of gap times (care time) and (break time) by fitting the AFT model using the function bivrecReg(). We again need to specify a bivrecSurv object as the response, but now we add the covariates to the right hand side of the formula (for further details on the syntax of this function see Table 3 and the reference manual of this package [9]). Note that only baseline (i.e., time-invariant) covariates can be fit in the AFT model in Eq. (2), otherwise, user will get the message, ~Error: Time-varying covariates not allowed.~ We use the function’s default estimation method, Lee et al.’s method, to obtain estimates and Wald confidence intervals based on the asymptotic standard error estimate. This method can also be explicitly specified by the option method=~Lee.et.al~. The messages shown while the function is running and the model fitting output from summary() are as follows.

Table 3

Arguments and compatible standard functions for function bivrecReg()

Argument	Description
formula	A formula with a bivrecSurv object on the left of a ’ ∼’ operator as response, and the covariate(s) on the right.
data	A data frame that includes all the variables listed in the formula.
method	A string indicating which method (“Lee.et.al” or “Chang”) to estimate the effects of covariates; the default is “Lee.et.al”.
Compatible functions:	summary(), vcov(), coef(), confint(), print()

We show the confidence interval for education and the full variance-covariance matrix using the functions confint() and vcov(), respectively. Similar to the standard R function confint(), one can specify the confidence level and an individual parameter of interest. If needed, one could also see coefficient results using the function coef().

Based on the above results, we conclude that patients in the secondary education or higher (EDU=1) group tended to have 1.76 (=e0.56, p-value =0.007) times longer care periods and 1.21 (=e0.19, p-value >0.05) times longer break periods than patients with less education (EDU=0) after adjusting for covariate Age10. In addition, a ten year delay in age of disease onset was associated with a 1.27-fold (=e0.24, p-value =0.001) increase in the length of break period and a minimal change in the care period of less than 1% (=e0.007, p-value >0.05), holding EDU constant.

We also fit the AFT model using Chang’s method for point estimation along with Parzen’s resampling algorithm to obtain standard error estimates and construct Wald confidence intervals. The results are as follows.

It is important to point out that due to the resampling algorithm needed to obtain variances, applying Chang’s method to the PCR data (n=326 after omitting missing values for EDU) led to substantially longer computing times than applying Lee et al.’s method (56 minutes vs. 67 seconds). Similarly, applying these methods to a simulated data set with a smaller sample size, n=150, led to a much shorter computing time (602 and 16 seconds for the two methods, respectively). In addition, as discussed earlier, the convergence of Chang’s method is not guaranteed for any dataset, and when the convergence is not achieved, users get the following error message, ~Error: Max Iterations reached. Did not converge.~ For these reasons, we recommend to use the default ~Lee.et.al~ method for fitting the AFT model. For a comparison of the estimation bias, standard error, and coverage probability between these two methods with more extensive simulation studies, see [3]. Nonetheless, we show how to use function simBivRec() provided in the BivRec package to simulate a bivariate alternating recurrent event data. This provides a way for users to perform their own simulation studies.

Following [3], in the example below, we show how to simulate a dataset with a categorical covariate a1 from a binomial distribution with success probability 0.5 and a continuous covariate a2 from a uniform (0,1) distribution using the function simBivRec(). This function has options to simulate data based on the various scenarios outlined in [3]. As an example, we set the parameters in the same way as for the scenario presented in the top panel of Table 3 in [3]: the sample size is n=150 (nsize=150); the regression coefficients for the effect of covariate a1 on the Type I and Type II gap times are set as beta1=c(0.5, 0.5); those for covariate a2 are set as beta2=c(0, -0.5); and the support of the uniform distribution (0,τ) for the censoring time is set with tau_c=63, which yields censoring rate of 15% for the first bivariate gap time pairs, on average. Additional parameters for the within-subject correlation structure of the gap times are set through the option, set=1.1, meaning the first parameter setting of simulation scenario 1 in [3].

The first few lines of the simulated data are shown as follows, where id is the subject ID (which can take both numerical and string values in the package), epi is the episode number of a gap time pair, xij and yij are the jth episode of the observed Type I and Type II gap times of subject i, respectively, ci is the overall censoring time of subject i, d1 and d2 are their corresponding censoring indicators, and a1 and a2 are two baseline covariates.

Conclusions

Despite the growing need and interest in the study of recurrent event data along with the development of R packages for its analysis such as survrec [10], reda [11] and reReg [12], software for analyzing bivariate alternating recurrent events has been lacking. In this paper, we reviewed nonparametric and semiparametric regression methods for gap times between alternating recurrent events and demonstrated how to use the BivRec package that we developed in R to perform these analyses. We also demonstrated BivRec’s capabilities for data visualization and simulation.

There is still a need for further package development such as additions of univariate recurrent gap times (i.e., gap times of recurrent events of the same type) methods as the degenerated case of the bivariate alternating gap times data such as the seminal nonparametric work by Wang and Chang (1999) [5] and the subsequent semiparametric regression methods by various authors [2, 13–16]. Finally, some competing, intensity-based models for bivariate alternating recurrent event data, which are not included in our package, can be found in [17] and [18].

Availability and requirements

Project name: Bivariate Alternating Recurrent Event Data Analysis

Project home page: https://cran.r-project.org/web/packages/BivRec/index.html

Operating system(s): Platform independent

Programming language: R

Other requirements: R 3.5.0 or higher

License: GPL-3

Any restrictions to use by non-academics: None