{'fraction_non_alphanumeric': 0.03068124734340585, 'fraction_numerical': 0.016615788863557326, 'mean_word_length': 4.694169416941694, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Piecewise survival models: a change-point analysis on herpes zoster associated pain data revisited and extended

For many diseases it is reasonable to assume that the hazard rate is not constant across time, but also that it changes in different time intervals. To capture this, we work here with a piecewise survival model. One of the major problems in such piecewise models is to determine the time points of change of the hazard rate. From the practical point of view this can provide very important information as it may reflect changes in the progress of a disease. We present piecewise Weibull regression models with covariates. The time points where change occurs are assumed unknown and need to be estimated. The equality of hazard rates across the distinct phases is also examined to verify the exact number of phases. An example based on herpes zoster data has been used to demonstrate the usefulness of the developed methodology.


Introduction
Herpes Zoster is a viral disease characterized by a painful skin and vesicular rash. The occurrence of herpes zoster follows from recrudescence of latent varicella-zoster virus (VZV), which is reactivated and multiplied in the dorsal root of cranial nerve ganglia. After primary infection with varicella (i.e. chickenpox), the virus persists dormant in the central nervous system (CNS) of the infected patient causing damage to it. This neurological disease is a cause of considerable morbidity, especially in both elderly and immunosuppressed individuals [9,19]. The zoster-associated pain (ZAP) is the most problematic symptom of the disease as it can be highly debilitating, and can also limit the patient's productivity [11].
With regard to the classification of ZAP in 1994 by Dworkin and Portenoy [8], a division of pain into three phases, measured in days, has been proposed as: acute herpetic neuralgia [0,30), subacute herpetic neuralgia [30,120), and postherpetic neuralgia (PHN [120,∞). In addition, the same classification is verified in other research studies associated with herpes zoster pain data in order to model pain resolution [1,6]. Prompt antiviral regiments in patients who suffer from VZV, such as Acyclovir (ACV) and Valaciclovir (VACV), can accelerate rash healing, decrease viral shedding and pain duration, and prevent or reduce the incidence of postherpetic neuralgia [10,11].
In accordance to the aforementioned, there are various and slightly arbitrary definitions for these transition points for which, in most of the existing relevant literature, there is an assumption based on the clinical evolution of the disease rather than data supported evidence [2,18,20]. It is therefore of considerable interest to have a well-defined description and a valid separation for the distinct phases of zoster pain.
Consequently, it is very important in medical research to consider a model that evaluates the effect of antiviral therapies on each phase, while controlling for the effect of the other phases as well as other covariates of interest. The arising challenge is to detect the transition times where the disease passes to the next stage and characterize ZAP, and also to answer the clinical questions; "When acute becomes PHN pain?" and "What are the treatment effects for a specific phase?". This analysis of time-to-event data assumes different hazard rate over distinguished time intervals giving rise to piecewise survival models. From a practical point of view, such a change-point detection, where the curvature of the survival functions change, can provide significant information as it may reflect changes in the progress of the disease of interest.
To this direction, piecewise exponential models (PEM) which are able to model the hazard rate of ZAP throughout its distinct time periods have been proposed [1]. We proceed by considering the less developed piecewise Weibull models (PWM), as they provide a richer family having as special case the constant per period hazard rate, but we also examine the turning points of the disease based on such models. Such points are considered as unknowns that have to be estimated. We obtain change-points estimates optimizing the likelihood function of the piecewise survival models. We also provide methodology to evaluate an antiviral therapy and figure out how the treatment effects can be compared for a specific phase of ZAP. However, the fit of piecewise models is not straightforward, especially if the points of change of the hazard rate are not known. In such cases standard approaches may fail, since for example some discontinuities may arise and hence, derivative based maximization may fail. We avoid this obstacle working with the Simulated Annealing (SA) optimization method, which does not use derivatives. The structure of PEM relies on the approach of Arani et al. [1], which examined the existence of three phases.
We propose an improved approach in the change-point analysis of herpes zoster by using Weibull models. We aim at contributing towards two directions. Firstly, by reexamining the herpes zoster data using a piecewise Weibull model to show that the model improves with respect to the piecewise exponential used so far, providing data driven insights about the pain phases. Secondly, we estimate the change-points together with the parameters of the model via a Simulated Annealing method to allow for full inference.
The following sections are structured as follows. In Section 2 the statistical methods are described. The proposed piecewise regression models are extensively delineated for multiple and unknown change-points. In Section 3 we apply this methodology on Herpes Zoster associated pain data for both cases of known and unknown transition points. The available data set has already been analysed by Beutner et al. [2] as well as in a relevant analysis of ZAP conducted by Desmond et al. [6]. Section 4 summarizes the main statistical and clinical findings of the paper.

The Piecewise Weibull Model
Among the class of parametric statistical models, the Weibull model is one of the most widely used. This model has been proved proper enough to describe a wide range of time-to-event data [13]. The hazard function for Weibull distribution is given by where κ > 0 is the shape parameter, ρ > 0 is the scale parameter of the distribution, and λ = ρ −κ is a reparameterization to simplify the notation. Note that the hazard rate is constant over time if κ = 1, and the Weibull model reduces to an exponential model. If κ > 1, then the hazard increases as time increases, and decreases over time otherwise. The addition of this shape parameter offers to the Weibull model great flexibility. The mean, median and variance of a Weibull random variable T can be expressed as , respectively. Furthermore, the survival and probability densities are respectively defined as: The piecewise Weibull model (PWM) is a way to combine interpretability and more flexibility by lightening the strict assumption of the piecewise constant hazard model. In PWM case, the basic idea is the partition of the time scale into J intervals. In the bibliography, there is a variety of forms for such a piecewise Weibull baseline hazard rate model. According to Casellas [3], the first and simplest form of a piecewise Weibull baseline hazard rate model each phase shares the same shape and scale parameters. However, we give an extended approach that allows variety among Weibull distribution's parameters. Consequently, the piecewise Weibull baseline hazard function is given by where ). The number of change-points is K = J − 1, and τ = (τ 1 , τ 2 , ..., τ J−1 ) represents the change-point vector, while κ j and ρ j are the positive shape and scale parameters of j-th interval, respectively. The corresponding survival is Other formulations of piecewise Weibull models are analysed in [5,14], which are known as poly-Weibull models and the baseline hazard function arises as the sum of the J independent Weibull-type hazards.
The heterogeneity of the population is stated by the vector of p covariates x i for the i-th individual. In general, the form of piecewise Weibull instantaneous hazard of i-th individual, whose the observed survival time belongs to the j-th stage of ZAP is defined as where θ = (κ, τ, β) includes the unknown parameters of interest, which are involved into the model. These are the shape parameters κ, the regression coefficients β which constitute the scale parameters and the change-points τ, respectively. The log-likelihood function for PWM in any stage is given by where y i are the time to event observations, δ i is the censoring and x i is the vector with the covariate information for the i-th patient. The triples {(y 1 , Of course, the piecewise exponential model arises as special case when all κ's are equal to 1. PWMs have been used at the past. There are some applications with known change-points [see for example , 17]. A Bayesian approach to estimate the change-points was proposed in [3]. In [16] a likelihood-ratio test (LRT) statistic was proposed that can be used to test whether there is an abrupt change at an unknown point in the covariate coefficient vector in the Weibull hazard function. Similar theoretical results on the LRT are provided in [14]. For the PEM one can also see the work in [7]. Hitherto, the works mainly attempt to join time intervals with common hazard. In the present paper we attempt to estimate the change-points together with the parameters of the model.

Computational Details
Maximizing the log-likelihood can be demanding due to the change-points since the loglikelihood surface can have jumps. To avoid such problems we propose Simulated Annealing as the most appropriate derivative free method to optimize the log-likelihood function (θ), estimate the unknown transition times and compare the efficiency of the treatments, for both models. It is a generic probabilistic heuristic approach to global optimization problems, investigated by Kirkpatrick et al. [12]. It locates a good approximation to the global optimum in a large search space with reasonable probability. This means that this optimization algorithm allows moving to worse values of the objective function, utilizing a probability, which decreases exponentially as the time passes and as more iterations are performed. Accepting worse solutions is a fundamental property of metaheuristics because it allows a more extensive search for the optimal solution and improves the obstacles of Local search algorithms.
The SA Algorithm 1 includes a parameter which is called temperature. The SA is based on the Metropolis-Hastings algorithm and simulates the cooling process by gradually lowering the temperature of the system until it converges to a steady, frozen state [4].
• Select initial and final temperatures, let T 0 = 500 and T ∞ > 0. The temperature is the Local Variable which controls the probability of downward steps.
• Select the function T r+1 = Tr 1+0.01Tr in which the temperature decreases, T r denotes the temperature at iteration r. 1: At the r-th iteration, given the current value θ r , propose new points θ * = (κ * , τ * , β * ) by sampling from a N (θ r , Σ) where Σ is diagonal with the same variance σ 2 .
Note that (·) is given in (5). 3: Reduce the temperature T according to the selected function (cooling rate). 4: Stop when you cannot find any better solution after a large number of iterations.
Some important details are: • Simulated annealing injects just the right amount of randomness into things to avoid the trapping attraction of local optimum early in the process without getting off course late in the game, when a solution is nearby. This makes it pretty good at tracking down a decent answer, no matter its starting point.
• For better optimization, when initializing the temperature variable we should select a temperature that will initially allow for practically any move against the current solution. This gives the algorithm the ability to better explore the entire search space before cooling and settling in a more focused region.
• The algorithm converges asymptotically to globally optimum solutions after a series of iterations and it does not allow any move against the optimum.
• The initial temperature must be large enough to make the uphill and downhill transition probabilities about the same, which results in avoiding local optimums. On the other hand, the initial temperature must be such that not moving off the global optimum. Consequently, the definition of some control parameters (initial temperature, cooling rate, etc.) constitutes the main disadvantage, that is usual in local search algorithms, because it is subjective and must be based on empirical evidence [12].
Standard errors of the estimated parameters can be obtained with standard bootstrap approach.

About the Data
The data derived from a randomized clinical trial to investigate the therapeutic effects of Acyclovir and Valaciclovir used also in [2] and [6]. Overall, 1141 herpes zoster patients, aged 49 years or older, were randomized to orally receive three different dosing schemes of treatment (see Table 1). The three arms refered as ACV-7days (Acyclovir, 800mg, 5 per day for 7 days), VACV-7days (Valaciclovir, 1000mg, 3 per day for 7 days) and VACV-14days (Valaciclovir, 1000mg, 3 per day for 14 days). The primary clinical outcome is the time to complete cessation of ZAP (in days). Other demographic characteristics, such as age and sex, are also included in the dataset and can be seen in Table 1. Three observations were excluded since they reported no pain. All statistical hypothesis tests have been conducted in the typical level of significance 5% and the analyses were carried out with the R 4.0.5 software package. Table 1 provides some descriptive statistics about the data. The herpes zoster data consist of 779 (68.45%) events of completed pain cessation and 359 (31.55%) censored observations. The majority of study population claimed pain duration before the 75th day of observation.  In Figure 1, the patients administered Valaciclovir for 7 and 14 days seem to have similar behavior in their survival probabilities, which are considerably close to each other. This observation is further supported by the log-rank tests in Table A.1, which shows that K-M curves among the ACV, VACV7 and VACV14 treatment groups are statistically equivalent (p-value= 0.08), while the therapies with ACV and merged VACV have significantly different K-M survival curves (p-value= 0.03). This means that the two dosing schemes of VACV are not significantly different with respect to the resolution of herpes zoster-associated pain. Therefore, it is preferred to continue using the merged VACV treatment. Also, note that there is a sudden drop of survival probabilities close to the day 30, which gives the impression of a possible change in survival.
Applying the Cox proportional hazard (PH) models in both cases for i. ACV vs VACV7 vs VACV14 and ii. ACV vs VACV result that Valaciclovir for 7 and 14 days significantly accelerated the ZAP resolution (p-value= 0.0277 and p-value= 0.037) compared to ACV, respectively. Specifically, complete cessation of pain occurs 20.7% faster with Valaciclovir than Acyclovir p-value= 0.0147), as shown in Table A

Piecewise Weibull Model
The objective in this section is to explore the benefits of a parametric analysis of the ZAP data by using a piecewise Weibull distribution to model the times to complete cessation of ZAP. Two of these benefits are, firstly, to develop a data-driven approach to estimate the time-points where the hazard function changes, and, secondly, to evaluate the hazard rates via an exhaustive local search algorithm. Finally, bootstrap approach was employed in order to obtain standard errors and 95% confidence intervals of the estimates. Both approaches of known (τ 1 = 30, τ 2 = 120) and unknown change-points are examined.

Known time change-points
The log-likelihood for the PWM is used to test treatment effects within each phase of ZAP under the consideration of known transition times (τ 1 = 30, τ 2 = 120) [8] and estimated shape parameterκ = (0.90, 0.89, 0.93). We compare the effectiveness of Acyclovir at 800mg five times a day for 7 days and Valaciclovir at 1000mg three times a day for 7 and 14 days in a clinical trial. We set X i = 0 if the i-th patient was treated with Acyclovir, and X i = 1 if Valaciclovir is administered to the i-th patient. The hazard rate of the i-th patient in the j-th phase is given by In PWM case, the interest lies on the instantaneous hazard rate, which is computed individually for each patient due to its dependence on the observed time. The hazard rates for each treatment group and each phase have also been calculated at three indicative time points, t, which are tabulated in Table 2. The hazard rates of pain resolution over time for Acyclovir and Valaciclovir, respectively, have a downward trend of: (i) 0.0136 to 0.0114 and 0.0153 to 0.0128 in acute phase; (ii) 0.0099 to 0.0087 and 0.0134 to 0.0117 in subacute phase; (iii) 0.00463 to 0.00452 and 0.00546 to 0.00534 in chronic phase. There is an obvious downward trend in hazard rates over time for both treatments which tends to be eliminated by crossing to the next stages of the disease. Moreover, the decline of hazard rates from subacute to chronic pain (or phase of postherpetic neuralgia) is larger than the corresponding one from acute to subacute. Both treatment arms are again statistically significant in patients who came up against the middle stage of ZAP, (subacute: p ACV < 0.0001; p VACV = 0.0087), while it is observed that patients who received Acyclovir have slightly smaller hazard rate. Testing the hypothesis H 0 : β V ACV = 0 for each phase indicates statistically significant difference among the treatments for the subacute phase (see Table 2). The subacute phase of pain lasts longer until ZAP resolution as shown in Figure 2. Also note that the estimated shape parameters for the different Weibulls at each phase areκ = (0.90, 0.89, 0.93) and they are smaller than 1 as implied by an exponential survival model (constant hazards). A LRT statistic testing the piecewise exponential model against the piecewise Weibull (using the same known time points for change) has a value 14.62 with 3 degrees of freedom and a p-value equal to 0.0022 justifying the choice of the piecewise Weibull model.
The differences correspond to the mean time of each phase.

Unknown time change-points
Now the change-points, i.e. the times the hazard rates change need to be estimated. The SA algorithm was applied to the herpes zoster data. Here we assumed that the two change-points were unknown and we left the SA to estimate them as two extra parameters of the model. Ten independent chains of the exhaustive local search algorithm (SA) of 6000 iterations were executed, setting the initial values (κ 0 , τ 0 , β 0 ) on the parameters of interest. In addition, we set T 0 = 500 as an initial temperature, and the decreasing function was T new = T old 1+0.01T old [15]. Some details about the implementation of the SA algorithm follow. We used the 30th and 120th days as initial values for change-points, based on classification in Dworkin and Portenoy [8]. Nevertheless, these values were randomly varied with common variance σ 2 during the estimation procedure. Similarly, the initial values of the regression coefficients β (Equation (6)) were estimated through a simple log-likelihood optimization process given the change-points τ = (30, 120). Consequently, the selection of initial values plays substantial role in the estimation with SA algorithm. Table 3 summarizes the results of 100 bootstrap replications of the SA algorithm including the parameter estimates (i.e. transition times, hazard rates and shape parameters) across the phases. The first change in ZAP hazard function is estimated atτ 1 = 44 days on average (SE=9.48), while the second appears on average atτ 2 = 144 days (SE=3.06). There is a downward trend in hazard rate across the phases for both treatment arms. Moreover, there is still considerable difference in hazard rates among Acyclovir and Valaciclovir, but only over the first two phases (acute and subacute), while in the chronic phase there is no distinction between the two treatments (see Figure 2). In agreement with the literature, this study of zoster-associated pain shows that Valaciclovir is more effective than Acyclovir in healing acute lesions. Recall that the commonly used values for the change-points were 30 and 120, respectively. The value of 30 lies inside the 95% confidence interval for the first point (interval lies from (24,58)), while for the second one 120 is perhaps too early (interval lies from (137,149)). Also note that for the shape parameters the value 1 is included in the confidence interval corresponding to the first phase, marginally not inside the confidence interval of the second phase, while for the third phase the interval boundaries are far from 1 implying that the constant hazards assumption is not valid overall. A final point relates to the second change-point which appears too late for the data. After this time point we have not observed any other event.
It is also interesting to notice that the standard error of the first change-point ( τ 1 ) is larger than the standard error of the second ( τ 2 ), as given in Table 3. This might happen because most patients with herpes zoster declared pain duration before the second phase, as shown in Figure 1, and less events observed during the second and the third phases. Figure 3 illustrates the parametric piecewise Weibull cumulative hazard estimates for both treatments. There is considerable difference when comparing the fitted cumulative curve slopes from PWM between the proposed change-points τ p = (30, 120) [8] and the expected estimated change-points (44, 144) as well as the existence of the change-points is more evident in the right figure. Akaike's criterion (AIC) was implemented to compare the candidate models regarding the analysis of ZAP data. According to the Table 4, there is significant fit difference among the piecewise exponential and piecewise Weibull modelling. In fact, the 4th model (that is, PWM with change-points at 44 and 144 days) seems to be the selected one. Also note that from Figure 2 one can see that when estimating the change-points later than the 30 days, this has a large effect on the changes in the hazard rates.
Interest lies on the graphical presentation of the log-likelihood given each candidate pair of change-points and its behavior with respect to τ 1 and τ 2 , as shown in Figure 4. One can see that the log-likelihood is relatively flat around the maximum value indicating the increased uncertainty around the change-points as already noted.

SA Estimate
Mean Standard error 95% CI *

Equality of hazards across the phases
The equality of hazard rates across the distinct phases constitutes an important part of the change-point analysis because it verifies the number of phases. Thus, modelling can be improved concerning the goodness of fit and complexity of the model. Based on the hypothesis testing of Arani et al. [1], the comparison of hazards from the best fitted model; i.e. the PWM model with change-points at (44,144) can be tested as This means that under the null hypothesis the hazard rate of herpes zoster pain resolution in phase i is equal to the corresponding one of phase j. As shown in Table 5, the hazard rates of pain cessation are significantly different between all phases for both treatment groups.

Acyclovir
Parameter Ratio Estimate p-value *

Conclusions
In this paper we revisited a data set related to Herpes Zoster pain using the piecewise Weibull model that offers certain advantages and flexibility while keeping the idea of the three phases of the disease pain. The model has also covariates and it does not assume constant hazard for each stage and thus allows for different behaviors.
We have implemented an approach where the change-points of the piecewise model are considered as unknown parameters and we need to estimate them. This can also be the basis to confirm the theoretical part about the pain phases. To avoid numerical problems due to the change-points we implemented a Simulated Annealing algorithm to estimate simultaneously the change-points and the parameters at each stage.
An interesting question related to how many change-points exist is not treated here in detail. It seems that LRT approaches can be useful for that, especially since the SA method allows for fitting the model without problems. In our case, working with more than 2 changing points (3 phases alternatively) we did not find significant improvements in the log-likelihood.
The proposed piecewise Weibull model, accompanied by the Simulated Annealing estimation algorithm, presents an improved method to analyse ZAP data by incorporating in the modelling the covariate information and a variety of other model parameters in each phase. Based on the paper results, the piecewise Weibull model demonstrated its superiority over the commonly used piecewise exponential modelling. As far as the disease related findings, the hazard function is almost flat over the acute phase, one may think that applying a mixture of piecewise regression models such as a piecewise exponential model for the acute phase and and a piecewise Weibull model for the subacute and chronic phases could possibly lead to a better fit. Finally, given that Gnann and Whitley Gnann and Whitley [10] have introduced a prodromal phase lasting the 1-5 days before the rash onset, the extension of the proposed model to accommodate the number of change-points as an additional model parameter will increase the flexibility of the proposed model and it might improve further its accuracy.