SiZer Map to investigate significant features of body-weight profile changes in HIV infected patients in the IeDEA Collaboration

Objectives We extend the method of Significant Zero Crossings of Derivatives (SiZer) to address within-subject correlations of repeatedly collected longitudinal biomarker data and the computational aspects of the methodology when analyzing massive biomarker databases. SiZer is a powerful visualization tool for exploring structures in curves by mapping areas where the first derivative is increasing, decreasing or does not change (plateau) thus exploring changes and normalization of biomarkers in the presence of therapy. Methods We propose a penalized spline SiZer (PS-SiZer) which can be expressed as a linear mixed model of the longitudinal biomarker process to account for irregularly collected data and within-subject correlations. Through simulations we show how sensitive PS-SiZer is in detecting existing features in longitudinal data versus existing versions of SiZer. In a real-world data analysis PS-SiZer maps are used to map areas where the first derivative of weight change after antiretroviral therapy (ART) start is significantly increasing, decreasing or does not change, thus exploring the durability of weight increase after the start of therapy. We use weight data repeatedly collected from persons living with HIV initiating ART in five regions in the International Epidemiologic Databases to Evaluate AIDS (IeDEA) worldwide collaboration and compare the durability of weight gain between ART regimens containing and not containing the drug stavudine (d4T), which has been associated with shorter durability of weight gain. Results Through simulations we show that the PS-SiZer is more accurate in detecting relevant features in longitudinal data than existing SiZer variants such as the local linear smoother (LL) SiZer and the SiZer with smoothing splines (SS-SiZer). In the illustration we include data from 185,010 persons living with HIV who started ART with a d4T (53.1%) versus non-d4T (46.9%) containing regimen. The largest difference in durability of weight gain identified by the SiZer maps was observed in Southern Africa where weight gain in patients treated with d4T-containing regimens lasted 52.4 weeks compared to 94.4 weeks for those with non-d4T-containing regimens. In the other regions, persons receiving d4T-containing regimens experienced weight gains lasting 51-61 weeks versus 59-77 weeks in those receiving non-d4T-based regimens. Discussion PS-SiZer, a SiZer variant, can handle irregularly collected longitudinal data and within-subject correlations and is sensitive in detecting even subtle features in biomarker curves.


Introduction
In the study of changes in longitudinal biomarkers in response to therapy or disease progression, it is useful to be able to identify the periods in time where changes occur. A key challenge arising from this effort is to isolate the underlying features of interest (say marker decreases or increases) in the presence of potentially large data variation. For example, in a data set of weight measurements in HIV-infected individuals initiating antiretroviral therapy (ART), which forms the core illustration of our methods in this paper, a scatterplot involving a mere 1% of the data ( (1)), which appear to capture the well-known features in such data involving rapid weight increase and subsequent plateau (2). However, it is unclear what the durability of weight gain is or whether there are decreases in weight after long-term exposure to therapy. In addition, each smoothing level produces a slightly different fit, particularly with respect to the timing of reaching the plateau in weight gain. As noted in Marron and Zhang (3), a hurdle in the application of smoothing methods is the selection of the smoothing parameter, because interesting features that are present in the data may be visible after applying some smoothing techniques or at some levels of smoothing but disappear in others, so choosing among the various smoothing techniques or the level of smoothing can be critical in extracting relevant features from the data; and of course, there is a tremendous computational burden associated with such data analyses, as the above conclusions were drawn from 2,000 patients and 12,000 observations, which, as mentioned above, constitute only 1% of the underlying database. The Significant Zero Crossings of Derivatives (SiZer) (4) was proposed to address many of the aforementioned issues (5). It is a useful Exploratory Data Analysis tool for understanding the significant features resulting from smoothed curves. SiZer simultaneously studies a family ofsmooth curves under a wide range of smoothing parameters (bandwidths) and produces inference on a smoothed version of the underlying curve viewed at varying levels of smoothing.
The standard implementation of SiZer (5) is based on the local linear smoother with a kerneltype smoothing method for a single predictor and a single outcome data (6). The SiZer map graphically explores structures in curves under study by mapping areas where the curve is significantly increasing, decreasing or does not change by studying its first derivative. Statistical inference is based on the derivatives of the smoothed curve by constructing confidence intervals at each location and also at each level of the smoothing parameters. The technique assembles these analyses at a wide range of smoothing parameters and synthesizes them in a single "map" where increase, decrease and plateau regions are identified by different colors, presenting an attractive global visualization of the data under many smoothing scenarios.
A number of extensions of the SiZer methodology have been proposed to increase inference precision (7,8). Of greater relevance to this paper is the extension proposed by Park and colleagues (9), which relaxes the assumption of independent errors, thus ignoring spurious features in the data which are caused by the existence of dependence in the data. These and other authors further extended the SiZer method into the area of time series (10,11). Another relevant extension is the SiZer for Smoothing Splines (SS-SiZer) (3), which uses splines to enhance detection of true features in the data.
Despite its attractiveness as a data visualization tool, the SiZer map has not been used widely in biological applications. This is unfortunate, since biological processes frequently involve changes in various measures (most notably biomarkers), which change over time, in response to disease progression or initiation and change of clinical therapies. One important reason for this is the fact that the SiZer and its extensions do not account for within-subject correlation. This frequently arises in longitudinal settings from measurements obtained repeatedly on the same individuals. A further technical complication is that the timing of these measurements becomes increasingly less regular with the passing of time (as subjects miss or reschedule clinical visits). It should be clear that this is a much different problem from time-series analysis, since longitudinal data are obtained from the same sample of study subjects repeatedly over time. Thus, neither the originally proposed method of SiZer nor its extensions in the area of time series fully address the challenges posed by longitudinal data.
It is a core aim of this paper to extend the method of SiZer maps to account for withinsubject correlation in the setting of irregularly collected longitudinal biomarker data, since SiZer offers an appealing global data visualization technique which can be tremendously useful in answering many important biological questions about the evolution of these data over time. We accomplish this by proposing a semiparametric extension of the SiZer methodology, named Penalized Spline SiZer (PS-SiZer), which combines a penalized spline regression model (10) with an embedded mixed-model representation of the marker evolution, coupled with methods which increase the computational efficiency of the standard SiZer. These computational advances are particularly attractive when analyzing massive databases with hundreds of thousands of patients and millions of observations. The paper is organized as follows: In the Methods we give a brief overview of the core ideas of the SiZer methodology as originally proposed (1,5) and the SS-SiZer methodology of Marron and Zhang (3) and present the proposed penalized spline PS-SiZer procedure for longitudinal data. Two simulation studies are presented in the Results where the proposed methodology is compared with the local linear smoother LL-SiZEr (1, 5) and the SS-SiZer (3) in, respectively, detecting changes and plateaus in longitudinal biomarker data. These are followed by the analysis of the IeDEA HIV data, where body weight measurements were collected repeatedly at each clinic visit. The clinical interest of this analysis is to describe the pattern of body weight changes after initiation of antiretroviral therapy as a surrogate of the treatment effectiveness and to determine the durability of weight gain by detecting a plateau in weight increases and the presence of possible decreases after long-term exposure to therapy. We conclude the paper with a brief discussion of our findings.
Insert Figure 1 about here

SiZer
More formally, for a given set of observed data and a smoothing function , we can consider a non-parametric regression model as follows: Here, is some "smooth" regression function that needs to be estimated from the data and The LL-SiZer applies the local linear regression method of Fan and Gijbels (12) to estimate and its derivative, . A common estimate of at each location of is given by the where is the minimum of the sum jointly over the regression coefficients, and . A line is fitted to the data for each using -weighted least-squares, where is a Gaussian ( ⋅ ) kernel. A SiZer map is then constructed by changing the value of the smoothing parameter .
The estimated regression function of and are obtained to construct a family of ( ), ' ( ) smooth functions at various levels of the smoothing parameter. Confidence limits for ' ( ) are obtained as where is a suitably defined Gaussian quantile (7). In the SiZer map, a pixel at and a specific smoothing level is colored blue if the confidence interval suggests that (implying that '( ) > 0 the curve at is increasing, red if the confidence interval suggests that (implying that ' ( ) < 0 the underlying curve is decreasing) and purple if the confidence interval contains zero (implying that no significant change in the curve can be detected).

SiZer for Smoothing Splines (SS-SiZer)
The SiZer for Smoothing Splines (SS-SiZer) (3) is an extension of a kernel-type estimation to the smoothing spline estimation. SS-SiZer incorporates the smoothing spline model and estimates the regression function by minimizing where is the smoothing spline parameter that determines the smoothness of the regression estimate and represents the roughness of the underlying function .
Here, the smoothing spline function is a natural cubic spline with knots at data locations ( ) . The smoothing parameter, acts similarly as the bandwidth in the LL-SiZer presented in 1 … the previous section.
SS-SiZer constructs point-wise confidence limits to produce the map. In our research, we apply the simultaneous confidence limit to the SS-SiZer model to address the multiplicity comparison issue (see next section). Otherwise, the interpretation remains the same as in the SS-SiZer maps (3). For other implementation details, such as the expression of first derivative estimate and its standard error, the reader is referred to the paper by Marron and Zhang (3).

The Penalized SiZer (PS-SiZer)
In this section, we present our extension of the SiZer map to handle data that arise in the longitudinal setting. In the proposed model, we consider subject-specific correlation arising from data obtained repeatedly on the same individuals. We utilize an approach similar to the standard SiZer in which a family of smooth functions is used at various levels of smoothing parameters .
We enhance the underlying model through the use of a computationally efficient smoothing model (presented below). In the PS-SiZer map, we also apply simultaneous confidence limits to resolve the issues related to multiple comparisons. Our proposed methodology extends SiZer in the following areas:

1)
Adding a random intercept component to summarize subject-specific correlation 2) Applying a P-spline (13) as the underlying smoothing function 3) Constructing a simultaneous 95% confidence limit addressing multiple-comparison issues

Model specification
Let denote measurements on subject at time , . We model the = 1, 2, … = 1,2. .., responses as, where is a smooth function indexed by a smoothing parameter and is a vector of ( ) random normal error terms with mean and variance . The model in (2) where is a vector of coefficients, and is the B-spline basis function of degree .
For the penalty term, the P-spline model of Eilers & Marx uses a base penalty on higher-order finite differences, (13). Consequently, the difference penalty matrix with order can be ∆ ∆ written as, . Here, is a matrix such that constructs the vector of difference of ∆ ∆ ∆ ∆ ℎ the coefficients i.e., ; and so on.
The second component of the model is the addition of a subject-specific random effect ~ . This results in the penalized least square objective function minimizing where and

Mixed model representation
The minimization problem discussed in the previous section can be handled using the mixedmodel framework (14). Equation ( penalty matrix is singular and has rank . A singular value decomposition (0, 2 The model in (4)

Estimate and variability bands of the derivatives
The derivatives of the smooth function are obtained by defining Here, ( ) and . The first derivative estimate of is: The estimated standard error of is .

Confidence bands
The where is the corresponding empirical best linear unbiased predictor (EBLUP) obtained ( ) from the mixed model framework. Here, is the quantile of the random variable at a smoothing level , i.e., which is the supremum on the set . The quantile was approximated The simulated values were sorted from smallest to largest, and = 10,000 the one with rank was used as . For a PS-SiZer map, we obtained the 95% quantile of equation (4) based on a simulation of size at each level of The confidence limits .
for were obtained as follows:

Simulation studies
In practice, the fundamental function of a SiZer map is to detect the underlying features in the data. For this reason, it is natural to compare the SiZer maps according to which ones detect the correct number of underlying features.
We have conducted Monte Carlo simulation studies to evaluate the performance of PS-SiZer map under various scenarios. The key objective of this simulation study was to compare the PS-SiZer with the LL-SiZer and SS-SiZer. Our simulation studies were designed to mimic the HIV data analyzed later as part of the illustration of the methodology.
Using these simulated data, the relative performance of the SiZer maps was evaluated by the following approaches: 1) By making the SiZer maps comparable at a similar level of "Effective Degrees of SS-SiZer) were generated with the same range of EDFs.

2)
By comparing the performance of the three SiZer maps at various levels of EDF according to which flags more features of a curve when a curve changes its status (increasing, decreasing or stable) from one to another.

3)
By comparing the performance of the three SiZer maps which are most sensitive to detect plateaus that are really present in the data.
In this research, performance of the PS-SiZer maps are presented in two different simulation studies: "Simulation Study 1" which addresses item 2 above and "Simulation Study 2" which addresses item 3.
Simulation study 1 Longitudinal data were simulated as chosen to be equally spaced in the interval [0, 1] with ( ) = 65 + 25 -2.0 * * sin(5 ( + 5)) + + where is random noise, is the subject-specific random intercept and denote the time of measurement. The function is a periodic function which has sin (5 ( + 5)) five features. By this term, we mean changes in the curve from increasing to decreasing or vice versa. The quantity, is a function to control the spread of the periodic sine function, 25 -2.0 which has the effect of diminishing the size of the features at higher time intervals.
Insert Table 1

about here
We compared three SiZer maps through the various levels of combination of error variance and subject-specific variance, ( respectively (Table 1)  most of the time. We were mainly interested to find how sensitive PS-SiZer is to detect the fourth and the fifth features of the true curve compared to LL-SiZer and SS-SiZer, as these were significantly diminished by the addition of the phasing-out component in the data as described above. When subject-specific variation is small (i.e., ), PS-SiZer detected all five features 2 = 2 51% of the time, whereas the SS-SiZer and LL-SiZer were able to detect all features 2% and 14% of the time respectively. At the same variability level, four features were detected by PS-SiZer 88% of the time, compared to 47% and 34% by SS-SiZer and LL-SiZer respectively.
When the subject-specific variation is high (i.e., ), PS-SiZer was still able to detect four 2 = 15 features in the data almost 68% of the time compared to 40% and 24% by SS-SiZer and LL-SiZer respectively. Interestingly, the fifth feature was not detected by LL-SiZer at all in this variability level, compared to 8% by PS-SiZer and 5% by SS-SiZer (Table 1).
Results from the above table are illustrated in Figure 2 randomly selected simulated trial are presented in Figure 4. Three maps were able to plot the pattern of the curve by moving from the blue region to the purple region at all levels of EDF.
Combined, Simulation Studies 1 and 2 demonstrate that the PS-SiZer map not only detects the significant changes of the true curve, but is also sensitive enough to detect the true time point where the curve reaches its plateau. Even though all three SiZer maps were able to detect the dominant features of the underlying curves, (that is, the trajectory of the curve from significantly increasing -blue area -to decreasing -red area -to non-significant -purple area), they were not able to detect less pronounced features at almost all levels of the EDF and were less sensitive than PS-SiZer in locating the true plateau of the curve.
Insert Figure 4 about here

Illustration
As an illustration of the proposed methodology of the PS-SiZer, we analyze data on weight changes in people living with HIV who initiate ART. In addition to detecting features in the data corresponding to body weight increases after the start of therapy, an important clinical question pertains to the durability of weight gains under different treatment regimens. More specifically, we explore possible differences in the durability of weight gain between stavudine (d4T) containing ART regimens versus non-d4T-containing regimens. Previous literature suggests that d4T is associated with lipodystrohy, a problem with the way the body produces and stores fat (19) and long-term weight loss compared to other regimens such as, for example, those containing Insert Table 2

about here
SiZer Maps for the Southern IeDEA region were generated for each ART group, i.e., one map each for the groups of patients initiating ART with a regimen containing or not containing d4T. To address the issue of durability of weight gain, we need to determine the first time point (in weeks from start of ART) at which weight gain stops, i.e., the time when either weight stops increasing or starts to decrease. The SiZer map provides an overall visual representation of the longitudinal weight change after the start of ART. However, to reach a final conclusion on the durability of weight increases after ART start, we need to decide on a single optimum level of smoothing. In this paper, the optimum smoothing parameter was estimated using the 'Rule of thumb' approach for local linear regression modelling described in Section 4.2 of Fan and Gijbels (12). Statistical analyses were performed using SAS Software version 9.3 and R software (2.13.2). SAS was used to create the analysis data sets for each of the five IeDEA regions. The R package SiZer (22) was used to fit local linear kernel model and corresponding SiZer maps using the analysis data sets.
The R package locpol (23) was used to obtain the best possible smoothing parameter for each of the treatment regimens in all five different regions so that estimated body weight changes over time could be generated along with their associated confidence interval.
Insert Figure 5

about here
The SiZer map for the Southern-Africa region is presented in Figure 5 The PS-SiZer map of body-weight changes among individuals initiating ART with a non-d4Tcontaining regimen shows that, at lower smoothing levels, there are some blue and purple areas suggesting an intermittent weight increase. Otherwise, the map consists of mostly blue areas (indicating weight increases) for medium and higher levels of smoothing for up to about 100 weeks after ART initiation. This indicates that patients starting ART with a non-d4T-containing regimen experience sustained body-weight increases for a period possibly double that of patients treated with d4T-containing regimens.
Insert Table 3

about here
To reach a conclusion about the comparison of the durability of weight changes in the d4Tcontaining versus not-containing ART regimens, we choose the SiZer analysis (row of the SiZer map) at the optimum level of smoothing for the Southern Africa IeDEA region (Figure 2 panel   b1). This analysis shows that weight in patients treated with d4T-containing ART regimens increased rapidly after ART initiation and plateaued afterwards. Consulting the first derivative ( Figure 2 panel b2), we observed that the 95% CI of the curve includes zero after 52.4 weeks in the group of patients who received a d4T-containing regimen compared to 94.4 weeks for patients treated with non-d4T-containing regimens (panel b2). A numerical summary of these results is also shown in the first row of Table 3. Supplementary Figures 1-4 for the East-Africa, West Africa, Central-Africa, and Asia-Pacific IeDEA regions respectively. The SiZer maps corresponding to the East and West Africa regions are very similar. For d4T-containing regimens (panel a1 in Supplementary Figures 1 and 2), blue areas are followed by purple areas after about 50-60 weeks for most levels of smoothing, indicating significantly increasing weight during this period. After this point, weight gain diminishes. By contrast, the blue areas in the SiZer maps corresponding to the non-d4T-containing regimens (panels a2 in Supplementary Figures 1 and 2) extend past week 60, indicating that weight continues to increase past 60 weeks after initiation of ART. Analyses at the optimum smoothing level produced the estimated curves of weight measurements shown in panels b1 and b2 of Supplementary Figures 1 and 2 and rows 2 and 3 in Table 3. For East Africa, results at the optimal smoothing levels showed that the weight in patients treated with d4Tcontaining regimens did not significantly increase after 54.5 weeks compared to 76.5 weeks for patients treated with non-d4T-containing regimens. For West Africa, the results are similar, with d4T-containing regimens estimated to weight gain for 51.4 weeks versus 65.8 weeks for the non-  Table 3. The estimated duration of weight gain in d4T-containing regimens was 51.3 weeks versus 58.7 weeks in non-d4T-containign regimens.

Discussion
This paper presents a significant extension of the SiZer methodology, the penalized SiZer or PS-SiZer. Current SiZer methods, such as the standard LL-SiZer (5) and SS-SiZer (3) do not account for the correlation induced by repeated measurements obtained on the same patient, which invariably arise in longitudinal settings with particular frequency in biomarker data. In addition, to developing a SiZer variant which can take into account within-subject correlation, our efforts were also centered on developing computationally efficient methods to address analyses involving massive databases from tens of thousands of subjects and millions of individual measurements.
The fundamental motivation of the originally proposed SiZer map is to detect the underlying features in the data and present a global visualization of changes in quantitative data for a spectrum of smoothing levels. The key goal of this research is to show propose a SiZer variant which can detect more real features in data in the context of data collected repeatedly from the same subjects at irregular time points longitudinally. From the simulation results, it was evident that both the standard LL-SiZer formulation and the SS-SiZer method, while able to detect large dominant features in the data, missed more subtle features, because neither method appropriately addresses within-subject variability resulting in wider confidence intervals and a diminished sensitivity when features in the data become attenuated (i.e., smaller changes from increases to decreases or vice versa). (3) have also attempted to compare these two maps by carrying out a number of simulations studies. The authors concluded that the original local linear version (here, LL-SiZer) of the SiZer and the smoothing spline SiZer (here, SS-SiZer) often performed similarly, without one method dominating the other in all cases. Similar findings were observed in our own simulation studies. By contrast, the PS-SiZer maps identified more underlying features in the simulation data than the other two SiZer map methods at a wide range of smoothing levels. In addition, both LL-SiZer and SS-SiZer detected a plateau in the data later compared to the PS-SiZer map, which detected the plateau almost exactly at the true time that it occurred in the simulated data. The simulation studies thus clearly demonstrate that, at a wide range of smoothing levels, PS-SiZer was more sensitive to small features in the data than the other two methods, presumably due to its improved ability to account for the presence of correlation in the data.

Marron & Zhang
The main idea of SiZer maps is to detect significant changes in the data by mapping areas where the 95% confidence intervals of the first derivative is significantly different from zero. The combination of the penalized spline regression model with random intercepts in the PS-SiZer map results in narrower confidence intervals, which, in turn, lead to more sensitive detection of even less prominent features present in the data compared to standard SiZer maps. In summary, PS-SiZer is a reasonably accurate addition to the family of SiZer map methods particularly when analyzing data from longitudinal settings.
In the application of the PS-SiZer methodology, we analyzed a database involving more than 185,000 adult HIV-infected patients and well over two million longitudinal weight measurements. Our ability to handle such a large data set, underlines the computational advantages of the proposed methodology. In addition to a global visualization of the data, the PS-SiZer analysis produced meaningful clinical results by showing that the durability of weight gain experienced by after starting ART with regimens containing d4T is likely significantly shorter than among persons who start ART with regimens which do not contain d4T.
Specifically, within the Southern Africa region, weight increases in the former regimens were observed to end after about 50 weeks from initiating of ART compared to almost 100 weeks among patients who started ART with regimens not containing d4T. While the clinical importance of this finding is less pronounced, given the almost universal phasing out of stavudine as a first-line regimen, weight gain among people living with HIV is a relevant topic, particularly with the wide adoption of integrase inhibitors and dolutegravir in particular, as main line antiretroviral therapies, all of whom are known to result in significant weight gain in these patients (24)(25)(26).        (a1) (a2) (b1) (b2) and its first derivative (bottom row), for d4T-containing and non-d4T-contiainging ART regimens (left and right column respectively).