FHtest: an R package for the comparison of survival curves with censored data

JSS Journal of Statistical Software November 2017, Volume 81, Issue 15. doi: 10.18637/jss.v081.i15 FHtest: An RPackage for the Comparison of Survival Curves with Censored Data Ramon Oller Universitat de Vic – UCC Klaus Langohr Universitat Politècnica de Catalunya Abstract The Fleming-Harrington class for right-censored data was first introduced by Harrington and Fleming (1982). This class is widely used in survival analysis studies and it is a subset of the so-called weighted logrank test statistics. Recently, Oller and Gómez (2012) proposed an extension of this class for interval-censored data. This paper introduces the Rpackage FHtest, which implements the Fleming-Harrington class for right-censored and interval-censored survival data. It provides an integrated approach for performing two-sample, k-sample and trend tests based on either counting process theory, likelihood theory, or permutation distributions. In this paper, we summarize the main aspects of the theory framework and present several examples with Rcodes to illustrate the usage of the main functions of FHtest. Keywords: interval-censored data, permutation tests, counting processes, R. 1. Introduction Survival times are encountered in many areas of research: for example, the time from cancer diagnosis until death; the time until first failure of a new machine; the sensory shelf life of yogurt, or the time an unemployed worker needs to find a new job. One common feature of all these examples is the possible presence of censored data, that is, incomplete information on the times of interest. Several types of censored data are distinguished: •Right-censored data: the (unknown) survival time lies above a certain value. •Left-censored data: the survival time lies below a certain value. •Interval-censored data: the survival time lies in an interval between two values. •Double-censored data: survival times are both left- and right-censored. 2FHtest: Comparison of Survival Curves with Censored Data in R As with other areas in statistics, common interests in survival analysis comprehend the estimation of distribution functions, hypothesis testing, or the fit of regression models. However, due to the presence of censored data, survival time data require specific methods that take into account that only partial information is available. In R(RCore Team 2017), there are plenty of packages that provide functions for the analysis of survival times including the survival package (Therneau 2017), which is one of the recommended packages of R. The CRAN task view on survival analysis (Allignol and Latouche 2017) gives an exhaustive overview on all these packages available on CRAN (https: //CRAN.R-project.org) or on Bioconductor (http://www.bioconductor.org/). The paragraph in the CRAN task view on hypothesis testing to compare survival functions lists several packages that offer functions for this purpose, such as function survdiff of the survival package. However, none of these packages included functions that implemented weighted logrank tests for interval-censored data based on the Fleming-Harrington class (Harrington and Fleming 1982). This lack of Rfunctions for this specific type of statistical tests to compare survival functions was the principal motivation for the development of the package FHtest (Oller and Langohr 2017). The package implements the Fleming-Harrington class for right-censored and interval-censored survival data and provides an integrated approach for performing twosample, k-sample and trend tests based on either counting process theory, likelihood theory, or permutation distributions. Several differences exist between functions provided by the FHtest package and those of other packages. For example, function survdiff of the survival package implements the Fleming-Harrington class with right-censored data and the counting process approach. Function survdiff, however, is limited to the one-parameter subclass instead of the whole class of two parameters. Recently, package survMisc (Dardis 2016) contains function comp, which implements the whole class of two parameters. However, neither survdiff nor comp consider one-sided alternatives, perform the permutation approach, or work with interval-censored data. Function ictest of the interval package (Fay and Shaw 2010) is programmed for interval-censored data, implements both the permutation and the score vector approach, and allows for one-sided alternatives. However, only two test statistics from the Fleming- Harrington class, the logrank and the Wilcoxon-Peto statistics, can be used. This paper introduces the FHtest package which is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=FHtest and it is organized around four sections: In Section 2, the theoretical background of the nonparametric tests of the Fleming-Harrington class for right and interval-censored data is set forth. Following, in Sections 3and 4, the four main functions of the FHtest package are presented and applied to real data examples, respectively. Section 5concludes the paper with some final remarks on the current state of the package and possible enhancements. 2. Comparison of survival curves with censored data We denote by Tthe time from a well-defined start point until the occurrence of an event of interest, E.Tis, hence, a positive random variable. One of the functions characterizing its distribution is the survival function, S, which corresponds to the probability of surviving a time t:S(t) = P(T > t)=1−F(t). Given the longitudinal nature of studies on survival times, a common feature of these are cen- Journal of Statistical Software 3 sored data. These arise whenever the occurrence of Eis not observed exactly, but only partial information is given. If Edoes not occur during the follow-up of a study, the corresponding time is right-censored. That is, it is only known that Tis larger than the study duration, but not by how much. Contrary to that, left censoring is present if Tis below a certain value, but it is not known by how much. If Eoccurs during a time interval, Tis said to be intervalcensored. The general convention is to consider the intervals left-open, that is T∈(L, R]. Both left and right-censored data can be considered particular cases of interval-censored data with L= 0 and R=∞, respectively. Alternatively, intervals can also be considered to be closed intervals, that is T∈[L, R]. Under this convention, exact observation can also be considered a particular case of interval-censored data, namely if L=R. In the current paper, we do not need to choose between both conventions because the FHtest package can handle both types of interval-censored data. For a more thorough introduction on survival times and censored data, see, e.g., Section 1.4 in Gómez, Calle, Oller, and Langohr (2009). In R, the nonparametric estimation of S(t)in the presence of right-censored data by means of the Kaplan-Meier estimator (Kaplan and Meier 1958) can be accomplished by using the function survfit of the survival package (Therneau 2017). For the estimation of the survival function with interval-censored data using the Turnbull estimator (Turnbull 1976), one may use the function icfit of the interval package (Fay and Shaw 2010) or different functions of the Icens package (Gentleman and Vandal 2017). When dealing with samples from two or more different populations, one is generally interested in comparing the corresponding distributions by testing one of the following hypotheses of interest: •Two-sample test: H0:S1(t) = S2(t)∀t > 0vs. H1:S1(t)6=S2(t)for some t > 0.(1) Alternatively, H1can also consider the case of a one-sided alternative. •k-sample test: H0:S1(t) = ··· =Sk(t)∀t > 0vs. H1:∃i, j, t such that Si(t)6=Sj(t).(2) •Increasing trend test: H0:S1(t) = ··· =Sk(t)∀t > 0vs. H1:S1(t)≤ ··· ≤ Sk(t)∀t > 0with, at least, one inequality. (3) Alternatively, H1can also consider the case of a decreasing alternative. The Rpackage FHtest provides functions for each of the three hypotheses above in the presence of right- or interval-censored data based on the Fleming-Harrington class. At this point, it is important to note that this family is geared to detect alternative hypotheses where the hazards between groups differ but do not cross. Under quite general conditions, the basic test statistic U=U1, . . . , Uk> 4FHtest: Comparison of Survival Curves with Censored Data in R is asymptotically normal with zero mean and covariance matrix V. The explicit form of U and Vis given below and depends on the type of censoring and the asymptotic approach. The test statistic for the hypothesis given in (1) is W1=U2 √v22 ,(4) where v22 is the second value in the diagonal of Vand the asymptotic distribution of W1is the standard normal distribution. The hypothesis in the k-sample problem (2) can be tested using the test statistic W2=U>V−U,(5) where V−is the generalized inverse of Vand the asymptotic distribution of W2is a χ2 k−1 distribution. The increasing trend alternative (3) can be tested with the statistic W3=a>U √a>V a,(6) where a= (a1, . . . , ak)>,a1< . . . < ak, are relevant covariate values for each group distribution. The asymptotic distribution of W3is the standard normal distribution. A one-tailed test in the left (or right) direction may show evidence of an increasing (or decreasing) trend. In the following, we briefly describe the test statistics of the Fleming-Harrington class for rightcensored (RC) and interval-censored (IC) survival data and derive the distribution approaches. 2.1. Right-censored data: Counting process approach The Fleming-Harrington class for right-censored data was first introduced by Harrington and Fleming (1982). This class of test statistics is a subset of the so-called weighted logrank test statistics. Fleming and Harrington (2005) give a theoretical framework for these statistics based on martingale and counting process methods. The Fleming-Harrington class has received much attention in the literature and many extensions have been proposed. For instance, Gillen and Emerson (2007) showed that these statistics are non-transitive under general alternatives and proposed an extension to avoid this drawback. More recently, Garès, Andrieu, Dupuy, and Savy (2015) proposed another extension which is designed to have good power against both late effects and proportional hazards alternatives. Consider a random sample of possibly right-censored data [Li, Ri], i = 1, . . . , n, with either Li=Ri(uncensored data) or Ri= +∞(right-censored data). Let t1< . . . < tmbe the m distinct uncensored times and suppose that there are n−mright-censored times. Then, the components of the weighted logrank test statistic URC are of the form URC j= m X r=1 w(tr)djr −njr nr dr, j = 1, . . . , k, (7) where drstands for the number of deaths at time trand nrgives the number of individuals at risk at time tr, that is, the number of individuals alive just prior to time tr. Both drj and nrj correspond to the same definitions but restricted to observations in the jth group. The weights suggested by Fleming-Harrington are w(tr) = ˆ S(tr−)ρ(1 −ˆ S(tr−))λ,(8) Journal of Statistical Software 5 where ˆ S(t)is the nonparametric estimation of the survival function obtained with the socalled Kaplan-Meier estimator. When (ρ, λ) = (0,0) and (ρ, λ) = (1,0), the statistic URC corresponds to the logrank and Prentice-Wilcoxon test statistics, respectively. The choice of the weights has to be made prior to the examination of the data taking into account that they should provide the greatest statistical power, which in turns depends on how the null is believed to be violated. For instance, in a given clinical trial, if one wanted to assess whether the effect of a treatment or therapy on the survival is stronger at the earlier phases of the therapy, λ= 0 should be chosen, with increasing values of ρemphasizing stronger early differences. If there were clinical reasons to believe that the effect of the therapy would be more pronounced towards the middle or the end of the follow-up period, it would make sense to choose ρ=λ > 0or ρ= 0, respectively, with increasing values of λemphasizing stronger middle or late differences. Statistical power may not be the only consideration in determining the weighting scheme, clinical aspects should also be contemplated. For instance, in cases where an intervention is thought to not alter the quality of the patient’s life, upweighting later differences in the survival curves might be important. The asymptotic null distribution of (7) is derived through a martingale central limit theorem. Under the null hypothesis, URC is asymptotically normal with zero mean and variancecovariance matrix VRCC given by VAR(Uj) = m X r=1 w(tr)2njr (nr−njr) n2 r (nr−dr)dr nr−1 and COV(Uj, Uj0) = − m X r=1 w(tr)2njr nj0r n2 r (nr−dr)dr nr−1, j, j0= 1, . . . , k. 2.2. Right-censored data: Permutation approach As shown in Abd-Elfattah and Butler (2007), the statistic (7) can also be written in a linear form URC = n X i=1 zici,(9) where zi=Z1i, . . . , ZKi>is a covariate vector of group indicators (Zji equals 1 or zero according to whether or not the ith individual is in the jth group) and ciis a score value associated with each individual: ci=            w(tr)− r X s=1 w(ts)ds ns if Li=Ri=tr, − r X s=1 w(ts)ds ns if Ri= +∞and tr≤Li< t(r+1). .(10) The permutation approach can be applied easily to this linear form. The idea is that if the null hypothesis is true and the underlying censoring process is identical across groups, the labels on the scores ciare exchangeable. The permutation distribution of URC is then obtained by permuting the labels and recomputing the test statistic for all the possible rearranged 6FHtest: Comparison of Survival Curves with Censored Data in R labels. The pvalues from the permutation distribution may be calculated either exactly using complete enumeration or a network algorithm, or by an approximation method like Monte Carlo resampling, or asymptotically using a version of the central limit theorem for exchangeable random variables (Sen 2006). This last method yields a normal approximation with zero permutation expectation and variance-covariance matrix VRCP =1 n−1n X i=1 c2 i n X i=1 (ziz> i−¯z¯z>), where ¯z =1 nPn i=1 ziis the sample mean. The permutation approach is, in fact, a conditional approach since the distribution of the test statistic is computed conditioned on the observed data. It is not obvious whether the permutation approach gives power properties similar to an unconditional approach. With right-censored data, Heimann and Neuhaus (1998) show that the permutation version of the logrank test and the unconditional version are asymptotically equivalent even under unequal censoring. A more recent contribution concerning this subject is given in Wang, Lagakos, and Gray (2010). 2.3. Interval-censored data: Permutation approach Oller and Gómez (2012) proposed a new class of test statistics for interval-censored data that extends the Fleming-Harrington class for right-censored data given in (7). The asymptotic behavior of the proposed tests was derived using a permutation approach and an observed Fisher information approach. Consider a random sample of interval-censored data (Li, Ri], let t1< . . . < tmdenote the unique ordered elements of {Li, Ri;i= 1, . . . , n}, and let ˆ S(t)be the nonparametric estimation of the survival function obtained with the so-called Turnbull estimator. Then, the components of UIC are of the form UIC j= m X r=1 wIC(tr)dIC jr −nIC jr nIC r dIC r, j = 1, . . . , k, (11) where dIC r=nˆ S(tr−)−ˆ S(tr)stands for the estimated number of deaths at time trand nIC r=nˆ S(tr−)gives the estimated number of individuals at risk at time tr. The computation of dIC rj and nIC rj is equivalent to dIC rand nIC rreplacing ˆ Sby ˆ Sj=1 NjPN i=1 Zji ˆ Siwhere Njis the sample size of the jth group and ˆ Siis an estimate of the individual survival function as defined by Fay and Shih (1998). The weights proposed by the authors are wIC(tr) = ˆ S(tr−)B(1 −ˆ S(tr); λ+ 1, ρ)−B(1 −ˆ S(tr−); λ+ 1, ρ) ˆ S(tr−)−ˆ S(tr), where B(t;a, b) = Rt 0xa−1(1 −x)b−1dx is an incomplete beta function. There are two reasons for these unusual weights. First, these weights follow from the score vector approach given in Section 2.4 for interval-censored data. Second, the equivalence between the weighted logrank form given in (11) and the linear form given in (12) is based on these weights. Further insights on these weights can be obtained from the following results. When (ρ, λ) = Journal of Statistical Software 7 (0,0) and (ρ, λ) = (1,0), the statistic UIC reduces to the logrank and Wilcoxon-Peto test statistics originally proposed in Peto and Peto (1972); for further details, see also Fay (1996, 1999). Asymptotically, as ˆ Sis nearly continuous at the point tr, the weight function wIC(tr) resembles [ˆ S(tr−)]ρ[1 −ˆ S(tr−)]λand, consequently, the interpretation of the parameters ρ and λreproduces the interpretation in the original right-censored family. This class UIC can also be written in a linear form (9) where cIC i=   ˆ S(Ri)B(1−ˆ S(Ri);λ+1,ρ)−ˆ S(Li)B(1−ˆ S(Li);λ+1,ρ) ˆ S(Li)−ˆ S(Ri)if ˆ S(Ri)6= 0, −B(1 −ˆ S(Li); λ+ 1, ρ)if ˆ S(Ri)=0. (12) Under the assumption that the underlying censoring process is identical across groups, this linear form of UIC allows the use of the permutation distribution (with zero expectation and variance-covariance matrix VICP). The approach is the same as explained in Section 2.2. Herein, it is considered that the observed intervals are half open intervals, however the definitions above are easily modifiable if we observe closed or open intervals. It is also worth mentioning that the methodology for interval-censored data can also be used for right-censored data as a particular case. In this case, both (10) and (12) would give similar but not identical results. 2.4. Interval-censored data: Score vector approach When we restrict the Fleming-Harrington class UIC to the subclass where λ= 0 (the so called Gρfamily), each test statistic (for a given value of ρ) can be seen as an efficient score test in a linear transformation model. We have to assume discrete or grouped continuous intervalcensored data: L, R ∈ {t0, t1, . . . , tm}where 0 = t0< t1< . . . < tm< tm+1 = +∞. We also assume a linear transformation model g(Ti) = −z> iβ+iwhere gis an unknown increasing function and ihas survival function S(t) = [1 + ρexp(t)]−1 ρ. Then, UIC =∂log(Lik(β, θ)) ∂ββ=0,θ=ˆ θ0 ,(13) where Lik(β, θ)is the likelihood function and θ= (θr)m r=1 are nuisance parameters such that θr=g(tr). Under discrete interval-censored data and the null hypothesis, the asymptotic behavior of the score vector UIC follows from maximum likelihood theory. This approach needs that none of the estimated θapproach the boundary of the parameter space, i.e., 1>ˆ S(t1)> . . . > ˆ S(tm)>0. The distribution is asymptotically normal with zero mean and variance-covariance matrix VICS. The explicit formula for VICS is not presented here. Instead, see Gómez and Oller (2008) for details. We note that URC and VRCC can also be derived by the score vector approach when λ= 0. In this case, however, a partial likelihood is used instead of a full likelihood. For this reason, though URC and VRCC or alternatively UIC and VICS may be used to analyze right-censored data, both approaches do not exactly coincide. 8FHtest: Comparison of Survival Curves with Censored Data in R 3. The FHtest package 3.1. Installation and dependencies The FHtest package can easily be installed from any CRAN repository by executing the R command R> install.packages("FHtest") To use FHtest, several Rpackages need to be available. As shown on the corresponding web page on CRAN (https://CRAN.R-project.org/package=FHtest), FHtest depends on the packages interval (Fay and Shaw 2010) and KMsurv (Yan, Klein, and Moeschberger 2012), and it imports functions from survival (Therneau 2017), perm (Fay and Shaw 2010), and MASS (Venables and Ripley 2002). The FHtest package is partially constructed over the interval package, which itself depends on the packages Icens (Gentleman and Vandal 2017) and MLEcens (Maathuis 2013): functions in both packages use similar structures in the input arguments and output values; some functions in FHtest are modifications of functions in interval and use the same internal procedures. The KMsurv package gives right-censored data sets that are used for illustration purposes. Both the survival and the MASS packages are recommended packages of Rand are, hence, included in all binary distributions of R. Finally, the package perm is from the same authors as interval and it is installed together with this package. Three of the packages – interval,MLEcens,KMsurv – will be automatically installed together with FHtest, if this has not been already done before. By contrast, package Icens is not available at the CRAN repositories, but at the Bioconductor web site (http://www. bioconductor.org/). For this reason, it has to be installed separately, for example, by executing the following instructions in R: R> source("http://bioconductor.org/biocLite.R") R> biocLite("Icens") Once all packages are installed, package FHtest can be loaded and used: R> library("FHtest") Loading required package: interval Loading required package: survival Loading required package: perm Loading required package: Icens Loading required package: MLEcens Loading required package: KMsurv 3.2. The main functions The four main functions of the FHtest package correspond to the censoring schemes and approaches exposed in Sections 2.1 through 2.4: Journal of Statistical Software 9 •FHtestrcc: The Fleming-Harrington test for right-censored data based on counting processes. •FHtestrcp: The Fleming-Harrington test for right-censored data based on permutations. •FHtesticp: The Fleming-Harrington test for interval-censored data based on permutations. •FHtestics: The Fleming-Harrington test for interval-censored data based on a score vector distribution. There are several features that are common to all four functions. First, the (censored) data and the group covariate can be specified in two different ways. One way is to assign two vectors, Land R, that define the censoring interval, and a group vector that gives the group covariate. The additional arguments Lin and Rin determine whether the endpoints are included or not. Another way is to use the common formula notation by including a ‘Surv’ object of the survival package; see the help page of function Surv for more details. Second, all functions allow for two-sample, k-sample and trend tests with hypotheses and alternatives as discussed in Section 2. The group vector determines the type of test: a twosample test when group has two unique values, otherwise a k-sample test when group is a character or factor vector, and a trend test when group is numeric. By default, the common argument alternative is set to "different", which corresponds to a two-sided alternative. For a one-sided alternative as in (3), the user can choose alternative = "increasing" for an increasing trend test or alternative = "decreasing" for a decreasing trend test. Third, the four functions have two arguments – rho and lambda – corresponding to the parameters ρand λof the weights of the Fleming-Harrington test statistic shown in (8). The only exception is function FHtestics, which does not allow values of lambda different from 0 as explained in Section 2.4. The default values are in all cases rho = 0 and lambda = 0. Finally, the outputs of all functions have the same structure indicating, among others, the type of test, the censoring scheme, the values of rho and lambda, and the alternative. In addition, the four functions, actually, return lists that can be stored as Robjects each with its own class according to the function in use. The arguments shared by the functions FHtestrcp and FHtesticp, but not available for FHtestrcc and FHtestics, have to do with the permutation approach. For example, argument method can be used to choose the way the pvalues are computed. There are four methods, the last two of which are only available for the two-sample test: "pclt" to use the permutational central limit theorem (Sen 2006); "exact.mc" to use an exact method via Monte Carlo; "exact.network" for an exact method using a network algorithm (see, e.g., Agresti, Mehta, and Patel 1990), and "exact.ce" for an exact method using complete enumeration, which is appropriate for very small sample sizes. The default value is "pclt". More details on all these aspects will be shown and explained in Section 4. 3.3. Data sets provided by package FHtest Two data sets with interval-censored data are provided by package FHtest:duser and illust3. 16 FHtest: Comparison of Survival Curves with Censored Data in R K-sample test for right-censored data Parameters: rho=1, lambda=1 Distribution: counting process approach Data: Surv(t2, d3) by group N Observed Expected O-E (O-E)^2/E (O-E)^2/V group=ALL 38 4.55 3.79 0.769 0.156 1.02 group=AML high risk 45 5.41 3.54 1.864 0.981 6.28 group=AML low risk 54 4.87 7.50 -2.633 0.924 8.99 Chisq= 9.9 on 2 degrees of freedom, p-value= 0.00697 Alternative hypothesis: survival functions not equal The use of function FHtestrcc for the k-sample test is apparently the same as before, but it is important to remark that the grouping variable must be either a factor or a character variable. Otherwise, a trend test is performed. Concerning the output of the k-sample test, it can be seen that the test statistic used, the one in (5), follows a χ2distribution with k−1=2degrees of freedom. It also shows a very small pvalue, that is, the null hypothesis can be rejected and we can assume that there are differences between the three risk groups with respect to disease-free survival after a bone marrow transplant. Moreover, we can observe that, contrary to the output of function survfit, the three rows of the output’s table are ordered alphabetically and not according to the internal coding of the factor group. Trend tests In the above situation or when dealing with an ordinal grouping variable, a trend test may be, actually, of greater relevance than the k-sample test. We reorder the factor group such that the risk group with best disease-free survival, low risk AML, is internally coded 1 and the one with worst disease-free survival, high risk AML, is coded 3: R> bmt$group <- relevel(bmt$group, ref = "AML low risk") R> levels(bmt$group) [1] "AML low risk" "ALL" "AML high risk" Hence, the null and alternative hypothesis of the trend test of interest are H0:S1(t) = S2(t) = S3(t)∀t > 0vs. H1:S1(t)≥S2(t)≥S3(t)∀t > 0with, at least, one inequality. (14) Contrary to the k-sample test before, in the case of the trend test, the grouping variable must be a numeric variable. For this reason, in the example below, function as.numeric is applied to factor group in order to use its internal numeric codes. It is also important to specify the value of argument alternative, the direction of the alternative hypothesis, which in (14) is decreasing: Journal of Statistical Software 17 R> FHtestrcp(Surv(t2, d3)~as.numeric(group), data = bmt, rho = 1, lambda = 1, + alternative = "decreasing", method = "exact.mc") Trend FH test for right-censored data Parameters: rho=1, lambda=1 Distribution: permutation approach (exact method using Monte Carlo with 999 replications) Data: Surv(t2, d3) by as.numeric(group) N O-E as.numeric(group)=1 54 -2.633 as.numeric(group)=2 38 0.769 as.numeric(group)=3 45 1.864 Statistic= 3, p-value= 0.001 99 percent confidence interval on p-value: 0.00000 0.00529 Alternative hypothesis: decreasing survival functions (higher group implies earlier event times) Among other aspects, the output shows in its first line that the trend test has been applied, and in its last line it indicates the direction of the alternative hypothesis. The function makes use of test statistic (6). The order of the rows of the output’s table follows the internal coding of the grouping variable. Therein, a negative sign in column O-E indicates less events observed than expected (under the null hypothesis), i.e., better survival than expected. Hence, the increase of the differences O-E from the first row to the last is in agreement with the alternative hypothesis of decreasing survival functions. This example, actually, not only illustrates the accomplishment of the trend test with rightcensored data, but also the use of an exact method via Monte Carlo to compute the pvalue. For this purpose, the value of FHtestrcp’s argument method was set to "exact.mc". As a result, in addition to the pvalue, the 99% confidence interval for the pvalue is also returned. Both, the number of the Monte Carlo replications (999) and the confidence level for the estimation of the pvalue, can be modified by changing the default values of argument permcontrol. To see which other characteristics of the permutation tests can be controlled through this argument, see the help page of (auxiliary) function permControl of the perm package, which is called by function FHtestrcp. Further examples for the functions FHtestrcc and FHtestrcp can be found and studied by calling example using the function names as topic argument: R> example("FHtestrcc") R> example("FHtestrcp") 4.2. Examples with interval-censored data For the illustration of functions FHtesticp and FHtestics, the data sets duser and illust3 of the FHtest package, which are presented in Section 3.3, will be used. As shown therein, 18 FHtest: Comparison of Survival Curves with Censored Data in R Years S ^(t) 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 Male Female 2.5 7.5 12.5 17.5 Years S ^(t) 0 1 2 3 4 5 6 7 0 0.2 0.4 0.6 0.8 1 Treatment Deferred therapy 500 mg/day dosage 1500 mg/day dosage Figure 2: Left panel: Survival functions of time from first injection drug use until HIV seroconversion among injection drug users. Right panel: Survival functions of time from treatment start until either AIDS onset or death among HIV-infected patients. Shaded rectangles represent the intervals within which the Turnbull estimate of S(t)is not identified. the right-censored observations in both data sets are identified by extremely large values for the right-endpoints of the censoring intervals: 9999 and 999, respectively. Even though not necessary for the application of FHtesticp and FHtestics, we first recodify these upper limits by NA values, in order to clearly distinguish rightfrom interval-censored data: R> is.na(duser$right) <- duser$right == 9999 R> is.na(illust3$right) <- illust3$right == 999 This recodification is also preferable for the nonparametric estimation of the survival functions by means of the Turnbull estimator, since no upper limit for the survival times is assumed. Two-sample tests For the two-sample test, we will use the data set duser, which contains the times from first injection drug use until HIV seroconversion among 940 injection drug users (IDUs). The left panel of Figure 2shows the estimation of the corresponding survival functions of both male and female IDUs. According to this graph, which is drawn with the plot method for ‘icfit’ objects from the interval package, the probability for seroconversion during the first three to four years is higher among female IDUs, after which the survival functions are nearly the same with the only difference that beyond 13 years, no more seroconversion is observed among male IDUs. In order to test the null hypothesis of equal survival functions among male and female IDUs, in the following two examples, we will use function FHtesticp, which implements the permutation approach. In the first example, we carry out the logrank test (ρ=λ= 0), use the permutational central limit theorem to compute the pvalue and choose the two-sided alternative hypothesis as shown in (1). In addition, we use arguments that are only available with functions FHtesticp and FHtestics: The arguments Lin and Rin are both set to TRUE to indicate that the censoring intervals are closed intervals, which is meaningful since time is measured in months. Moreover, the argument icontrol is used to increase the maximum Journal of Statistical Software 19 number of iterations used internally by the EM algorithm to estimate the survival functions. This is necessary here because convergence problems arise when using the default value of 10000 iterations. R> FHtesticp(Surv(left, right, type = "interval2") ~ zgen, data = duser, + Lin = TRUE, Rin = TRUE, icontrol = icfitControl(maxit = 50000)) Two-sample test for interval-censored data Parameters: rho=0, lambda=0 Distribution: permutation approach (asymptotic approximation using central limit theorem) Data: Surv(left, right, type = "interval2") by zgen N O-E zgen=0 759 -16.4 zgen=1 181 16.4 Statistic Z= 1.8, p-value= 0.0694 Alternative hypothesis: survival functions not equal The output is very similar to the one of function FHtestrcp with the title being the only difference. Here, we can see that the data set consists of 759 men (zgen = 0) and 181 women and that the pvalue is 0.0694. What the output does not show is that the computational cost for this example is high. The execution of the previous command, with 940 interval-censored data, takes more than 40 seconds under OS X Yosemite with an Intel processor of 2.5 GHz and the following session info: R> sessionInfo() R version 3.3.3 (2017-03-06) Platform: x86_64-apple-darwin13.4.0 (64-bit) Running under: OS X Yosemite 10.10.5 If we used function FHtestics instead, the computational cost would roughly be the same. That is, the high computation time is not due to the use of the permutational central limit theorem, but due to the use of interval-censored data, in particular due to the necessary nonparametric estimation of the common survival function. For this reason and in order to reduce computation times, if we are interested in carrying out further hypothesis tests with the same data, it is recommended to save the ‘FHtesticp’ object and use its element fit a posteriori as the value of FHtesticp’s argument icFIT. This is done in the following example, where we put more weight on early differences and use a one-sided alternative that could be motivated by some prior knowledge: 20 FHtest: Comparison of Survival Curves with Censored Data in R R> icpexam <- FHtesticp(Surv(left, right, type = "interval2") ~ zgen, + data = duser, Lin = TRUE, Rin = TRUE, + icontrol = icfitControl(maxit = 50000)) R> FHtesticp(Surv(left, right, type = "interval2") ~ zgen, data = duser, + rho = 1, Lin = TRUE, Rin = TRUE, alternative = "decreasing", + icFIT = icpexam$fit) Two-sample test for interval-censored data Parameters: rho=1, lambda=0 Distribution: permutation approach (asymptotic approximation using central limit theorem) Data: Surv(left, right, type = "interval2") by zgen N O-E zgen=0 759 -12.4 zgen=1 181 12.4 Statistic Z= 2.3, p-value= 0.0109 Alternative hypothesis: decreasing survival functions (group=1 has earlier event times) In this example, the use of the stored ‘FHtesticp’ object has greatly reduced the computational time to less than one second. k-sample tests The data set illust3, which contains the data from an AIDS clinical trial, will be used for the illustration of the k-sample test with interval-censored data. Since the group indicator is a numeric variable, we first put the corresponding labels to each of the three study groups, which are roughly of the same sample size: R> illust3$group <- factor(illust3$group, labels = c("Deferred therapy", + "500 mg/day dosage", "1500 mg/day dosage")) R> table(illust3$group) Deferred therapy 500 mg/day dosage 1500 mg/day dosage 541 538 528 The Turnbull estimates of the survival functions for time from treatment start until either AIDS onset or death are shown in the right panel of Figure 2. Therein, we can see that the estimated survival functions are nearly identical throughout the first year under treatment and then start to split: The survival function of the deferred therapy lies below the other two. These remain the same for nearly four years, after which the probabilities for disease-free survival are somewhat larger among the 1500 mg/day dosage group. To test, whether these difference allow us to assume that there are really differences between the three treatments with respect to disease-free survival, we now use function FHtestics, Journal of Statistical Software 21 which implements the score vector approach set forth in Section 2.4. The null and alternative hypotheses are the ones of the k-sample test as shown in (2) and for this reason, the grouping variable used in the following command syntax has to be a factor or a character vector. By setting rho = 3, we put more weight on early differences than to differences at long term, and since disease-free survival is given in months, arguments Lin and Rin are, again, set to TRUE. R> FHtestics(Surv(left, right, type = "interval2") ~ group, data = illust3, + rho = 3, Lin = TRUE, Rin = TRUE) K-sample test for interval-censored data Parameters: rho=3, lambda=0 Distribution: score vector approach Data: Surv(left, right, type = "interval2") by group N O-E group=1500 mg/day dosage 528 -17.27 group=500 mg/day dosage 538 -3.23 group=Deferred therapy 541 20.50 Chisq= 12.2 on 2 degrees of freedom, p-value= 0.00221 Alternative hypothesis: survival functions not equal Like before, when we applied the k-sample test to right-censored data, the value of test statistic (5) is computed, which asymptotically follows a χ2with 2 degrees of freedom. This value is fairly large and, consequently, the pvalue is quite small. Trend tests We also apply the trend test to this data set using the same null but the opposite alternative hypotheses as in (14): H0:S1(t) = S2(t) = S3(t)∀t > 0vs. H1:S1(t)≤S2(t)≤S3(t)∀t > 0with, at least, one inequality. (15) Here, S1is the survival function that corresponds to the deferred therapy, S3the one of the 1500 mg/day dosage treatment, and since H1states increasing disease-free survival functions, the argument alternative will now be set to "increasing". As with the trend test with right-censored data, the grouping variable must now be a numeric variable: R> FHtestics(Surv(left, right, type = "interval2") ~ as.numeric(group), + data = illust3, Lin = TRUE, Rin = TRUE, alternative = "increasing") Trend FH test for interval-censored data 22 FHtest: Comparison of Survival Curves with Censored Data in R Parameters: rho=0, lambda=0 Distribution: score vector approach Data: Surv(left, right, type = "interval2") by as.numeric(group) N O-E as.numeric(group)=1 541 51.4 as.numeric(group)=2 538 -8.0 as.numeric(group)=3 528 -43.4 Statistic Z= -4.2, p-value= 1.19e-05 Alternative hypothesis: increasing survival functions (higher group implies later event times) In both examples with function FHtestics, we either used equal weights for all differences between observed and expected events or put more weight on short-term differences. If, contrary to that, we wanted to put more weight on late differences, we would need to use function FHtesticp, because function FHtestics does not allow for values of lambda different from 0. Further examples for the functions FHtesticp and FHtestics can be found and studied by executing the following instructions: R> example("FHtesticp") R> example("FHtestics") 5. Final remarks In this paper we have presented the Rpackage FHtest, whose functions perform tests for right- and interval-censored survival data based on the Fleming-Harrington class. We have shown that the FHtest package is a flexible and user-friendly tool that integrates and extends methods in other Rpackages and we have given several detailed examples of the use of each function. The FHtest package has some limitations. For example, it covers left- and double-censored data as particular cases of interval-censored data, however, it does not cover neither truncated data nor doubly-censored (i.e., both the time origin and the event of interest are censored) data. Another intrinsic drawback is that weighted logrank tests are not appropriate for crossing hazards alternatives. In the future, we plan to perform a simulation study to compare the four FHtest methods. With right-censored data and for both small and large samples, we are interested in studying the differences between distinct approaches: (i) the RC counting process versus the IC score vector method; (ii) the RC permutation versus the IC permutation method; and (iii) the RC counting process versus the RC permutation method. With interval-censored data, we are also interested in studying the differences between both IC approaches: the IC permutation versus the IC score vector method. Journal of Statistical Software 23 Acknowledgments The work was partially supported by grants MTM2012-38067-C02-01, MTM2012-38067-C02- 02, and MTM2015-64465-C2-1-R of the Ministerio de Economía y Competividad (Spain) and by grants 2014 SGR 464 and 2014 SGR 598 from the Departament d’Economia i Coneixement de la Generalitat de Catalunya (Spain). The authors want to thank the research group Grup de Recerca en Anàlisi Estadística de la Supervivència (GRASS) for the fruitful discussions and are also thankful to the Data Analysis and Modeling Research Group for the valuable support. References Abd-Elfattah EF, Butler RW (2007). “The Weighted Log-Rank Class of Permutation Tests: p-Values and Confidence Intervals Using Saddlepoint Methods.” Biometrika,94(3), 543– 551. doi:10.1093/biomet/asm060. Agresti A, Mehta CR, Patel NR (1990). “Exact Inference for Contingency Tables with Ordered Categories.” Journal of the American Statistical Association,85(410), 453–458. doi:10. 1080/01621459.1990.10476220. Allignol A, Latouche A (2017). CRAN Task View: Survival Analysis. Version 2017-04-25, URL https://CRAN.R-project.org/view=Survival. Dardis C (2016). survMisc: Miscellaneous Functions for Survival Data.Rpackage version 0.5.4, URL https://CRAN.R-project.org/package=survMisc. Fay MP (1996). “Rank Invariant Tests for Interval Censored Data Under the Grouped Continuous Model.” Biometrics,52(3), 811–822. doi:10.2307/2533044. Fay MP (1999). “Comparing Several Score Tests for Interval Censored Data.” Statistics in Medicine,18(3), 273–285. doi:10.1002/(sici)1097-0258(19990215)18:3<273:: aid-sim19>3.0.co;2-7. Fay MP, Shaw PA (2010). “Exact and Asymptotic Weighted Logrank Tests for Interval Censored Data: The interval RPackage.” Journal of Statistical Software,36(2), 1–34. doi:10.18637/jss.v036.i02. Fay MP, Shih JH (1998). “Permutation Tests Using Estimated Distribution Functions.” Journal of the American Statistical Association,93(441), 387–396. doi:10.1080/01621459. 1998.10474120. Fleming TR, Harrington DP (2005). Counting Processes and Survival Analysis. John Wiley & Sons. doi:10.1002/9781118150672. Garès V, Andrieu S, Dupuy JF, Savy N (2015). “An Omnibus Test for Several Hazard Alternatives in Prevention Randomized Controlled Clinical Trials.” Statistics in Medicine, 34(4), 541–557. doi:10.1002/sim.6366. Gentleman R, Vandal A (2017). Icens: NPMLE for Censored and Truncated Data.doi: 10.18129/B9.bioc.Icens.Rpackage version 1.48.0. 24 FHtest: Comparison of Survival Curves with Censored Data in R Gillen DL, Emerson SS (2007). “Nontransitivity in a Class of Weighted Logrank Statistics Under Nonproportional Hazards.” Statistics & Probability Letters,77(2), 123–130. doi: 10.1016/j.spl.2006.06.001. Gómez G, Calle ML, Egea JM, Muga R (2000). “Risk of HIV Infection as a Function of the Duration of Intravenous Drug Use: A Non-Parametric Bayesian Approach.” Statistics in Medicine,19(19), 2641–2656. doi:10.1002/1097-0258(20001015)19:19<2641:: aid-sim527>3.0.co;2-p. Gómez G, Calle ML, Oller R, Langohr K (2009). “Tutorial on Methods for Interval-Censored Data and Their Implementation in R.” Statistical Modelling,9(4), 259–297. doi:10.1177/ 1471082x0900900402. Gómez G, Oller R (2008). “A New Class of Rank Tests for Interval-Censored Data.” Harvard University Biostatistics Working Paper Series. Working Paper 93, URL http://biostats. bepress.com/harvardbiostat/paper93/. Harrington DP, Fleming TR (1982). “A Class of Rank Test Procedures for Censored Survival Data.” Biometrika,69(3), 553–566. doi:10.1093/biomet/69.3.553. Heimann G, Neuhaus G (1998). “Permutational Distribution of the Log-Rank Statistic Under Random Censorship with Applications to Carcinogenicity Assays.” Biometrics,54(1), 168– 184. doi:10.2307/2534005. Kaplan EL, Meier P (1958). “Nonparametric Estimation From Incomplete Observations.” Journal of the American Statistical Association,53(282), 457–481. doi:10.2307/2281868. Maathuis M (2013). MLEcens: Computation of the MLE for Bivariate (Interval) Censored Data.Rpackage version 0.1-4, URL https://CRAN.R-project.org/package=MLEcens. Oller R, Gómez G (2012). “A Generalized Fleming and Harrington’s Class of Tests for Interval-Censored Data.” The Canadian Journal of Statistics,40(3), 501–516. doi:10. 1002/cjs.11139. Oller R, Gómez G, Calle ML (2004). “Interval Censoring: Model Characterizations for the Validity of the Simplified Likelihood.” The Canadian Journal of Statistics,32(3), 315–326. doi:10.2307/3315932. Oller R, Langohr K (2017). FHtest: Tests for Right and Interval-Censored Survival Data Based on the Fleming-Harrington Class.Rpackage version 1.4, URL https://CRAN. R-project.org/package=FHtest. Peto R, Peto J (1972). “Asymptotically Efficient Rank Invariant Test Procedures.” Journal of the Royal Statistical Society A,135(2), 185–207. doi:10.2307/2344317. RCore Team (2017). R: A Language and Environment for Statistical Computing.RFoundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Sen PK (2006). “Permutational Central Limit Theorems.” In Encyclopedia of Statistical Sciences. John Wiley & Sons. Journal of Statistical Software 25 Therneau TM (2017). survival: Survival Analysis.Rpackage version 2.41-3, URL https: //CRAN.R-project.org/package=survival. Turnbull BW (1976). “The Empirical Distribution Function with Arbitrarily Grouped, Censored and Truncated Data.” Journal of the Royal Statistical Society B,38(3), 290–295. Venables WN, Ripley BD (2002). Modern Applied Statistics with S. 4th edition. Springer- Verlag, New York. Volberding PA, Lagakos SW, Grimes JM, Stein DS, Rooney J, Meng T, Fischl MA, Collier AC, Phair JP, Hirsch MS, Hardy WD, Balfour HH, Reichman RC (1995). “A Comparison of Immediate with Deferred Zidovudine Therapy for Asymptomatic HIV-Infected Adults with CD4 Cell Counts of 500 or More per Cubic Millimeter.” New England Journal of Medicine,333(7), 401–407. doi:10.1056/nejm199508173330701. Wang R, Lagakos SW, Gray RJ (2010). “Testing and Interval Estimation for Two-Sample Survival Comparisons with Small Sample Sizes and Unequal Censoring.” Biostatistics, 11(4), 676–692. doi:10.1093/biostatistics/kxq021. Yan J, Klein JP, Moeschberger ML (2012). KMsurv: Data Sets From Klein and Moeschberger (1997), Survival Analysis.Rpackage version 0.1-5, URL https://CRAN.R-project.org/ package=KMsurv. Affiliation: Ramon Oller Departament d’Economia i Empresa Facultat d’Empresa i Comunicació Universitat de Vic – UCC Vic, Spain E-mail: [email protected] URL: https://www.uvic.cat/en/resultdir/ramon.oller Journal of Statistical Software http://www.jstatsoft.org/ published by the Foundation for Open Access Statistics http://www.foastat.org/ November 2017, Volume 81, Issue 15 Submitted: 2016-03-19 doi:10.18637/jss.v081.i15 Accepted: 2016-10-28