Intraindividual pain variability and phenotypes of pain in sickle cell disease: a secondary analysis from the Pain in Sickle Cell Epidemiology Study.

ABSTRACT: Mean pain intensity alone is insufficient to describe pain phenotypes in sickle cell disease (SCD). The objective of this study was to determine impact of day-to-day intraindividual pain variability on patient outcomes in SCD. We calculated metrics of pain variability and pain intensity for 139 participants with <10% missing data in the first 28 days of the Pain in Sickle Cell Epidemiology Study. We performed Spearman rank correlations between measures of intraindividual pain variability and outcomes. We then used k-means clustering to identify phenotypes of pain in SCD. We found that pain variability was inversely correlated with health-related quality of life, except in those with daily or near-daily pain. Pain variability was positively correlated with affective coping, catastrophizing, somatic symptom burden, sickle cell stress, health care utilization, and opioid use. We found 3 subgroups or clusters of pain phenotypes in SCD. Cluster 1 included individuals with the lowest mean pain, lowest temporal instability and dependency, lowest proportion of days with pain and opioid use, and highest physical function. Cluster 2 included individuals with the highest mean pain, highest temporal dependency, highest proportion of days with pain and opioid use, and lowest physical function. Cluster 3 included individuals with high levels of mean pain, highest temporal instability, but with lower temporal dependency, proportion of days with pain and opioid use, and physical function compared with cluster 2. We conclude that intraindividual pain variability is associated with patient outcomes and psychological characteristics in SCD and is useful in delineating phenotypes of pain in SCD.

1. Introduction

Sickle cell disease (SCD), a chronic multisystem disorder, affects an estimated 100,000 individuals in the United States[15] and disproportionately impacts minorities and those of lower socioeconomic status.[28] Painful vaso-occlusive episodes are the hallmark of the disease[29] and contribute to morbidity[29] and poor health-related quality of life (HRQoL).[6,27] There is wide interindividual variability in pain episodes in SCD.[29] Many patients with SCD experience chronic pain,[34] which is associated with greater functional disability,[32] and worse pain-related patient-reported outcomes.[3] We have previously reported the inadequacy of mean pain intensity alone in defining pain phenotypes in SCD.[4] A recent review by Mun et al.[24] described approaches for the study of intraindividual variability in pain intensity, henceforth referred to as pain variability. Mun et al. reported that pain variability is associated with psychological well-being[31] and response to treatments[14,40] in arthritis and fibromyalgia, as well as postherpetic neuralgia and painful diabetic neuropathy.[11] Pain variability may be relevant in SCD as it may contribute to loss of control, a theme that has been encountered in SCD, and is related to the unremitting and unpredictable course of painful episodes[39] and the unpredictability of SCD pain.[1]

To understand the significance of intraindividual variability in pain intensity in SCD, we conducted a secondary analysis of diary data collected in the Pain in Sickle Cell Epidemiology Study (PiSCES), the largest epidemiological study of pain to date in SCD.[21,34] This study first demonstrated the high out-of-hospital burden of pain in SCD using daily pain diaries over a 6-month period. Of the 31,017 patient-days in this study, representing 232 patients submitting at least 30 diaries, pain was reported on 54.5% of days, and about 30% of this cohort reported pain on >95% of their diary days.[34] The objectives of this secondary analysis were to (1) determine whether intraindividual pain variability is associated with HRQoL, pain-related psychological characteristics, and health care utilization for pain and (2) determine whether an unsupervised learning statistical approach can identify subgroups or phenotypes of pain, based on measures of pain intensity and pain variability.

2. Methods

2.1. Pain in Sickle Cell Epidemiology Study

The design of the PiSCES has been previously published.[33,34] The Pain in Sickle Cell Epidemiology Study enrolled 308 individuals with SCD aged 16 years or older in Virginia between July 2002 and August 2004. Of the 308 individuals enrolled, 76 individuals were excluded from study analysis (23 did not send in any diaries, and 53 sent in <30 diaries), leaving 232 individuals to constitute the analysis sample.[34] Study participants completed daily pain diaries for a 6-month period, reporting daily pain, interference, distress from pain, presence of a self-reported crisis, health care utilization, and opioid use.[33,34] This study also collected baseline demographics, self-reported clinical comorbidities, and pain-related psychological characteristics.

2.2. Pain diaries

Paper pain diaries[33] were completed by patients daily for up to 6 months, and study participants received payment for each returned diary, with a higher payment in the latter 2 months of the study to encourage study completion. Study participants reported about the prior 24 hours. They rated their worst sickle cell pain intensity (0-9 numerical scale) and marked a body diagram to indicate where they hurt. In addition, they reported whether they were in a crisis and whether they had made a visit to the emergency department and hospital or either a scheduled or unscheduled clinic visit because of sickle cell pain. Patients also indicated what medication they had taken. Crisis days were self-defined by each patient using a check box on each daily diary.

2.3. Pain-related psychological characteristics and health-related quality of life

Results of pain-related psychological characteristics and HRQoL measured at baseline in the PiSCES cohort have been previously reported.[8,19,22,23,33-35] Depression, anxiety, and somatic symptom burden were measured using the Patient Health Questionnaire (PHQ).[19] The 2 depression diagnoses generated by the PHQ (major depressive syndrome and other depressive syndrome) were combined into a single category of depression, and similarly, the 2 anxiety diagnoses (panic syndrome and other anxiety syndrome) were combined into a single category of anxiety.[19] Somatic symptom burden was measured by the PHQ, but the items related to pain in 4 sites (limb, back, stomach, and chest) were excluded as these were common sites of SCD pain, and the score was calculated based on 11 items.[35] Catastrophizing was measured using the catastrophizing subscale of the Coping Strategies Questionnaire.[8] Coping was measured by the Coping Strategies Questionnaire, that was adapted for SCD,[13,23] and scores from 3 subscales: active coping (ignoring pain sensations, reinterpreting pain sensations, calming self-statements, diverting attention, and increasing behavioral activities), affective or emotional coping (anger, fear, catastrophizing, isolation, and praying), and passive or behavioral adherence coping (taking fluids, resting, and heat/cold/massage) were calculated.[23] Health-related quality of life was measured using the SF-36 module,[42] and scores were reported for the 8 subscales: Bodily Pain (BP), General Health (GH), Mental Health (MH), Physical Function (PF), Role-Physical (RP), Role-Emotional (RE), Social Function (SF) and Vitality (VT),[22] as well as the physical component summary score and the mental component summary score.[23] Sickle cell disease stress was calculated using the Sickle Cell Stress Scale.[23]

2.4. Indices of pain variability

To describe affective instability seen in mood disorders, Jahng et al.[16] operationalized variability using several indices. The intraindividual standard deviation (iSD) measures the amplitude of fluctuations[24] and is a popular measure of variability[12] as it captures each participant's iSD over all of their observations and has been used to describe variability in SCD pain.[4] However, the iSD does not account for sequence or temporal dependency[16] and does not measure the frequency of fluctuations.[24] Temporal dependency, as measured by autocorrelation, indicates the extent to which current observations can be predicted from previous observations.[24] A positive autocorrelation at lag 1 (AR1) indicates that the pain experience is more likely to persist over time when it is above or below the mean pain level, whereas a negative value implies a back-and-forth pattern of pain ratings, and a value close to 0 indicates that pain level yesterday does not predict the pain level today.[24] Although autocorrelation takes into account the order of observations, and the persistency, it does not take into account the magnitude of fluctuations. Temporal instability is calculated based on measures of successive change, such as the Mean Square of Successive Differences (MSSD) and the Probability of Acute Change (PAC).[16] The MSSD is the average of the squared successive changes between 2 adjacent observations and takes into account both the amplitude of fluctuation and temporal dependency.[16] Mean square of successive difference is also a product of the population variance (iSD2) and (1-AR1),[24] and a large MSSD reflects a high iSD or a low AR1.[24] The PAC in pain intensity, which measures the number of acute changes in pain scores divided by the total number of successive changes,[16] can be calculated for an investigator-determined threshold of acute change. A change of 0.9 on the 0 to 10 Numerical Rating Scale has been found to be clinically meaningful in children with SCD,[25] and a change of 20% has been proposed as a minimally important change in chronic pain trials,[10] which would approximately be a change of 2 points on a 0 to 9 scale; therefore, in this study, it was calculated for both the probability of an acute change of pain intensity score of 1 (PAC1) and acute change of pain intensity score of 2 (PAC2). Both the MSSD and PAC can be used as measures of instability due to temporal fluctuations,[16] and the indices complement each other. Taken together, the above indices describe different facets of variability and may be calculated for each individual to measure intraindividual variability in pain intensity from diary data.

2.5. Data selection

For this analysis, we identified individuals in PiSCES who reported pain intensity scores for a minimum of 26 days of the first 28 days (up to 10% missing data) of diary data submission. The initial 28-day period was chosen apriori for analysis for several reasons: (1) it was likely that this period was associated with the least amount of missing data because missing data tend to increase over time,[4] and (2) this period was closest in time to baseline assessments of demographic, clinical, and psychological characteristics and HRQoL. Missing pain intensity scores were imputed with the simplest approach, last observation carried forward, for up to 2 missing data days. Imputation was performed to ensure a complete series and ability to calculate indices such as AR1, which cannot be calculated in the presence of missing values. If the first day was missing, the following day's pain intensity score was imputed. Imputation was not performed for other variables in the data set. As a check on sensitivity, indices of central tendency and pain variability were calculated with the original (unimputed) and imputed pain intensity variables and compared using the paired t-test.

2.6. Data analysis

All data analysis was conducted in CRAN R v.3.6.2 (R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria),[30] and statistical significance was assessed at the 0.05 threshold. We calculated mean, median (p50), 90th percentile of pain intensity reports (p90),[4,36] and indices of intraindividual variability (iSD, PAC1, PAC2, MSSD, and AR1) for every patient for the 28-day period. Autocorrelation at lag 1 could not be calculated in some cases (eg, division by 0). We also calculated the proportion of pain days (PPD), proportion of days with crisis, proportion of days with emergency department or hospital utilization, and proportion of days with opioid use. Spearman rank correlation coefficients were calculated to determine strength and direction of associations between (1) mean, median, and 90th percentile pain intensity, with indices of intraindividual pain variability, and (2) intraindividual pain variability with HRQoL and pain-related psychological factors and health care utilization for pain. Spearman rank correlations were also performed for 3 subgroups, those with <50% days with pain, ≥50 to <95% days with pain, and those with ≥95% days with pain during the 28-day period. The threshold of ≥50% of days was adapted to a 28-day period based on the frequency characteristic recently put forth by the AAPT guidelines for the diagnosis of chronic SCD pain, which defines the frequency criterion of chronic SCD pain as pain on more than half the days in the last 6 months[9]; however, we also separated out the subgroup with ≥95% days with pain, which represented those with daily or near-daily pain. Using centered and scaled data with Euclidean distance, we used k-means clustering procedures to identify unsupervised subgroups of the summarized pain measures using the factoextra package Version1.0.7 in R.[17] Variables entered into the algorithm included measures of pain intensity (mean, median, and p90), pain variability (iSD, PAC1, PAC2, MSSD, and AR1), and pain frequency (PPD). These unsupervised procedures allow for the discovery of patient groupings that may not have otherwise been considered.

Preliminary determination of the optimal number of clusters was based on a series of internal validation procedures, comprising the elbow method for within-cluster sum of squares, the silhouette method, and the Calinski–Harabasz criterion.[7] Robustness of the internal clustering solution was gauged via bootstrapping. Specifically, 1000 bootstrapped samples pulled participants from the original data set with replacement, approximating random samples drawn from the original population. The Calinski–Harabasz criterion was then calculated for each of the 1000 bootstrapped samples, and the number of selected clusters for each bootstrapped sample was plotted with a bar chart, with the majority rule compared with the original Calinski–Harabasz result. Following fit of the original (ie, internally valid) clustering solution, external validity was assessed using the Adjusted Rand Index (ARI), a metric that computes a similarity value between 2 cluster solutions (ie, predicted and true) and counts pairs that were assigned in the same or differing clusters in the predicted and true solutions. Namely, k-means clustering solutions were fit to each of the 1000 bootstrapped samples created for internal validation, and these bootstrapped solutions were compared with the original solution for similarity of cluster partitions. As a control, ARIs were also calculated for 1000 randomly generated clustering solutions relative to the original k-means solution. ARI results between the bootstrapped solutions and the randomly generated solutions were tested for a statistical difference using a paired t-test. Once the final clustering solution was robustly validated, both internally and externally, subsequent analysis compared differences between clusters on variables external to the algorithm, including demographic, clinical, and psychological variables using χ2 tests and Kruskal–Wallis tests, followed by an unadjusted post hoc multiple comparison procedure (Dunn test) as needed. The Institutional Review Board at Emory University provided a nonhuman subjects research waiver for secondary analysis of the deidentified PiSCES data set.

3. Results

3.1. Data selection and demographics

One hundred thirty-nine individuals of the total of 232 individuals had less than 10% missing data for pain intensity scores in the first 28 days of submission. A total of 112 pain intensity scores from 79 individuals were missing, of a possible total of 3892 pain intensity scores, representing 2.8% missing pain intensity data. Missing pain intensity scores were conservatively imputed using last value carried forward, as described. We did not find differences in the mean or medians between the imputed and unimputed series. However, we did find differences in the calculation of PAC1, PAC2, and MSSD (all P < 0.01), but the magnitude of the differences was moderately small based on effect size (ES range 0.25-0.36). We therefore conducted the analysis using the imputed variables, so the same series could be used to calculate indices of pain variability including AR1 because AR1 could not be calculated from the unimputed series due to the presence of missing data. Of the 139 individuals included, AR1 could not be calculated for 16 of them, as the denominator for calculation for AR1 was 0. Of the 139 participants included, we found that 55 has <50% days with pain (median PPD 0.11 [IQR 0-0.36]), 37 had ≥50 to <95% days with pain (median PPD 0.79 [IQR0.61-0.86]), and 47 had ≥95% days with pain (median PPD 1 [IQR 1-1]). Baseline demographic, clinical characteristics, scores on pain-related psychological measures, and HRQoL of the participants included in the analysis are reported in Table 1. With the exception of the RP subscale (P = 0.019) of the SF-36, there were no differences in demographic, clinical, psychological, or HRQoL characteristics between those included in this analysis and those who were excluded. Median scores on RP were lower in the individuals included as opposed to those who were excluded. We present descriptive statistics for the indices of pain intensity, variability, and pain frequency for the included participants in Table 2.

Table 1

Demographic, clinical, HRQoL, and psychological characteristics of the sample included in analysis and comparison with the excluded sample.

	Excluded, n= 93	Included, n=139	P
Age (median [IQR])	35 [28, 43]	33 [25, 41]	0.315
Female sex (n, %)	64 (68.8)	79 (56.8)	0.089
Education (n, %)
Less than high school	8 (8.6)	20 (14.4)	0.386
High school	38 (40.9)	50 (36)
More than high school	47 (50.5)	69 (49.6)
Income (n, %)
<$10,000	37 (41.6)	51 (37)	0.522
$10,000-$20,000	18 (20.2)	34 (24.6)
$20,000-$30,000	16 (18)	18 (13)
>$30,000	18 (20.2)	35 (25.4)
Married (n, %)	24 (25.8)	31 (22.3)	0.647
Genotype (n, %)
HbSS or HbS- β0 thalassemia	66 (71)	103 (74.6)	0.641
HbSC or S- β+ thalassemia	27 (29.0)	35 (25.4)
SF-36 (median [IQR])
General Health	40.00 [25, 55]	35.00 [23.75, 50]	0.269
Physical Function	65.00 [45, 85]	60.00 [45, 80]	0.602
Mental Health	76.00 [60, 88]	76.00 [64, 91]	0.430
Social Function	62.50 [50, 87.5]	62.50 [37.5, 87.5]	0.678
Bodily Pain	45.00 [32.5, 67.5]	45.00 [22.5, 67.5]	0.963
Vitality	40.00 [25, 50]	45.00 [22.5, 60]	0.712
Role-Physical	50.00 [0, 75]	25.00 [0, 75]	0.019
Role-Emotional	100.00 [33.33, 100]	66.67 [0, 100]	0.073
Physical component summary score	34.90 [31.09, 41.35]	33.93 [27.53, 41.05]	0.217
Mental component summary score	49.53 [41.44, 56.36]	49.01 [38.79, 57.26]	0.909
Depression (n, %)	27 (29)	37 (26.6)	0.800
Anxiety (n, %)	6 (6.5)	9 (6.5)	1.000
Catastrophizing (median [IQR])	14.00 [8.25, 18.00]	12.00 [7.00, 19.00]	0.934
Sickle cell stress (median [IQR])	20.00 [11.00, 26.00]	22.00 [12.22, 27.00]	0.174
Coping (median [IQR])
Affective/emotional focused	2.93 [2.20, 3.49]	2.70 [1.88, 3.60]	0.370
Passive/behavioral adherence	4.03 [3.40, 4.78]	4.19 [3.39, 4.83]	0.611
Active	2.67 [1.90, 3.33]	2.67 [1.80, 3.37]	0.808
Somatic symptom score* (median [IQR])	7.00 [4.00, 9.00]	7.00 [4.00, 10.00]	0.794
Avascular necrosis* (n,%)	18 (19.4)	30 (21.7)	0.785
Skin ulcers* (n, %)	11 (11.8)	15 (10.9)	0.989

For income n = 227, genotype n = 231, SF-36 subscales n = 226 to 232, catastrophizing n = 225, sickle cell stress n= 229, coping subscales n = 226 to 230, and somatic symptom score n = 230.

n = 138 (included), all available data presented.

HRQoL, health-related quality of life.

Table 2

Descriptive statistics for measures of pain intensity and pain variability.

	n	Median (IQR)
Mean	139	2.32 (0.5-4.59)
Median	139	2 (0-5)
Ninetieth percentile of pain intensity (p90)	139	5 (2-7)
Intraindividual standard deviation (iSD)	139	1.42 (0.78-2.13)
Probability of acute change of 1 point (PAC1)	139	0.22 (0.07-0.41)
Probability of acute change of 2 points (PAC2)	139	0.11 (0-0.22)
Mean square of successive differences (MSSD)	139	2.19 (0.83-4.96)
First-order autocorrelation (AR1)	123	0.35 (0.09-0.49)
Proportion of pain days (PPD)	139	0.68 (0.21-1)

3.2. Relationship of indices of pain intensity with pain variability

We examined Spearman rank correlations of mean pain intensity and measures of pain variability such as iSD, MSSD, PAC1, PAC2, and AR1. As reported in Table 3, there were statistically significant low to moderate correlations between measures of pain variability with mean pain for the overall sample. When stratified by pain frequency, we observed that the magnitude and significance of correlation between mean pain intensity and measures of pain variability was lower in those with more frequent pain.

Table 3

Spearman rank correlation coefficients between mean pain intensity and measures of pain variability for overall sample and stratified by pain frequency.

	Mean pain intensity
	All patients*	<50% pain days†	≥50-<95% pain days‡	≥95% pain days§
SD	0.53‖	0.97‖	0.80‖	−0.02
PAC1	0.54‖	0.92‖	0.34#	0.03
PAC2	0.37‖	0.86‖	0.47¶	0.02
MSSD	0.50‖	0.92‖	0.70‖	0.05
AR1	0.21#	0.62‖	0.11	−0.06

n = 139 (123 for AR1).

n = 55 (40 for AR1).

n= 37.

n = 47 (46 for AR1).

P < 0.001.

P < 0.01.

P < 0.05.

AR1, first-order autocorrelation; MSSD, mean square of successive differences; PAC1, probability of acute change of 1 point in pain intensity score; PAC2, probability of acute change of 2 points in pain intensity score.

In Supplementary Table 1 (available at http://links.lww.com/PAIN/B493), we report the relationship of pain intensity and pain variability with each other. For the overall sample, we also found that mean, p50, and p90 were highly correlated with each other. Among measures of pain variability, iSD was highly correlated with MSSD, PAC1, and PAC2. Low correlations were observed between AR1 and other pain variability indices.

3.3. Relationship of intraindividual pain variability with health-related quality of life

In Table 4, we report the associations of the subscales of the SF-36 with mean and pain variability, for the overall sample, and stratified by pain frequency (<50% days with pain, ≥50 to <95% days with pain, and ≥95% days with pain). Measures of pain variability were negatively correlated with HRQoL, and the strength of the associations was low, as shown in Table 4 and Figure 1, and further in Supplemental Table 2 (available at http://links.lww.com/PAIN/B493). In the subgroup with <50% days with pain, similar to mean pain intensity, measures of variability were negatively correlated with HRQoL, particularly General Health, RP, and BP subscales and the physical component summary score. In those with ≥50 to <95% days with pain, scores on the BP subscale were negatively correlated with several measures of variability. In the subgroup with ≥95% days with pain, mean pain intensity was negatively correlated with physical function and SF, but not with other subscales. However, the GH subscale was positively correlated with PAC1, PAC2, and MSSD and negatively correlated with AR1. With increasing PAC1, PAC2, and MSSD, overall GH scores increased, suggesting that some fluctuation in pain in those with daily or near-daily pain is associated with better GH, and increasing AR1 (or persistence) was associated with lower GH. First-order autocorrelation, however, was negatively associated with scores on the RE subscale, and the mental component summary score in those with ≥95% days with pain, but none of the other measures of variability nor correlated with MH scores in any subgroup.

Table 4

Spearman rank correlation coefficients between mean pain and measures of pain variability with HRQoL, for all patients, and stratified by pain frequency.

All	Mean	iSD	PAC1	PAC2	MSSD	AR1
General Health	−0.42§	−0.24‖	−0.21¶	−0.14#	−0.21¶	−0.15#
Physical Function	−0.53§	−0.19¶	−0.22‖	−0.13	−0.19¶	−0.14
Role-Physical	−0.46§	−0.23‖	−0.20¶	−0.15#	−0.20¶	−0.22¶
Role-Emotional	−0.44§	−0.19¶	−0.19¶	−0.11	−0.17#	−0.17#
Vitality	−0.45§	−0.24‖	−0.22‖	−0.11	−0.18¶	−0.21¶
Mental Health	−0.31§	−0.22‖	−0.26‖	−0.18¶	−0.20¶	−0.05
Social Function	−0.48§	−0.22¶	−0.15#	−0.10	−0.14#	−0.28‖
Bodily Pain	−0.58§	−0.38§	−0.29§	−0.20¶	−0.31§	−0.33§
Physical component summary score	−0.56§	−0.29§	−0.26‖	−0.20¶	−0.27‖	−0.24‖
Mental component summary score	−0.37§	−0.22‖	−0.21¶	−0.13	−0.17#	−0.1

n = 136-139 (120-123 for AR1).

n = 55 (40 for AR1).

n = 37.

n = 47 (46 for AR1).

P < 0.001.

P < 0.01.

P < 0.05.

0.05 > P < 0.1.

AR1, first-order autocorrelation; iSD, intraindividual standard deviation; MSSD, mean square of successive differences; PAC1, probability of acute change of 1 point in pain intensity score; PAC2, probability of acute change of 2 points in pain intensity score.

Figure 1.

Correlation between pain variability and outcomes. Heatmap representation of Spearman rank correlation coefficients between mean pain and measures of pain variability with HRQoL, psychological characteristics, health care utilization, and opioid use for all patients (statistically significant [P < 0.05] correlations are represented in color). AR1, first-order autocorrelation; HRQoL, health-related quality of life; MCS, mental component summary score; MSSD, mean square of successive differences; PAC1, probability of acute change of 1 point in pain intensity score; PAC2, probability of acute change of 2 points in pain intensity score; PCS, physical component summary score.

3.4. Relationship of intraindividual pain variability with psychological factors

In Table 5, we report association of pain variability with psychological factors, including anxiety, depression, coping, somatic symptom burden, and sickle cell stress for the overall sample. Pain variability was most consistently associated with affective coping, catastrophizing, somatic symptom burden, and sickle cell stress; however, the magnitude of the association was low. In Supplementary Table 3 (available at http://links.lww.com/PAIN/B493), we explored these relationships, stratified by pain frequency, and found that in those with ≥95% days with pain, AR1 was positively correlated with affective coping, catastrophizing, and sickle cell stress.

Table 5

Spearman rank correlation coefficients between mean pain intensity and pain variability with psychological factors, health care utilization, and opioid use.

	Mean	iSD	PAC1	PAC2	MSSD	AR1
Psychological factors
Anxiety	0.27†	0.18‡	0.08	0.14	0.18†	0.04
Depression	0.28†	0.11	0.10	0.06	0.09	0.03
Affective coping	0.31*	0.25†	0.23†	0.16§	0.20‡	0.18§
Passive coping	0.12	0.15§	0.09	0.08	0.13	0.17§
Active coping	0.13	0.10	0.07	0.08	0.11	0.02
Catastrophizing	0.22†	0.24†	0.22‡	0.16§	0.19‡	0.14
Somatic symptom burden	0.44*	0.28*	0.30*	0.24†	0.30*	−0.01
Sickle cell stress	0.44*	0.27†	0.23†	0.14§	0.19‡	0.27†
Health care utilization and opioid use
Proportion of days with ED or hospital utilization	0.30*	0.39*	0.25†	0.26†	0.29*	0.28†
Proportion of days with opioid use	0.81*	0.44*	0.40*	0.26†	0.39*	0.27†

n = 135 to 139 (119-123 for AR1).

P < 0.001.

P < 0.01.

P < 0.05.

0.05 > P < 0.1.

AR1, first-order autocorrelation; ED, emergency department; iSD, intraindividual standard deviation; MSSD, mean square of successive differences; PAC1, probability of acute change of 1 point in pain intensity score; PAC2, probability of acute change of 2 points in pain intensity score.

3.5. Relationship of pain variability with health care utilization and proportion of days with opioid use

In Table 5 and Figure 1, we also report that there were low to moderate correlations between pain variability, opioid use, and health care utilization in the overall sample, but mean pain intensity retained the strongest association with opioid use. In general, the strength and significance of the correlation coefficients, however, declined with increasing pain frequency, as shown in Supplemental Table 3 (available at http://links.lww.com/PAIN/B493).

3.6. Phenotypes of pain in sickle cell disease, based on measures of pain intensity, pain variability, and pain frequency

Of the 139 individuals, 16 individuals had missing scores for AR1, and they were removed from this analysis. Clusterability was first evaluated via the Hopkins statistic (values range from 0 to 1, with higher values indicating higher clusterability), and clustering was determined to be adequate (0.77).[17] Euclidean distance was calculated. Internal validation methods (ie, (1) elbow method for within-cluster sum of squares; (2) silhouette method; and (3) Calinski–Harabasz criterion) unanimously returned a 3-cluster solution for the original data set (Fig. 2). This result was internally confirmed by bootstrapping the Calinski–Harabasz criterion and finding 3 clusters chosen most commonly (60%), followed by 4 clusters (22%) and 2 clusters (14%), shown in Figure 2. External validity was demonstrated by the ARI, comparing bootstrapped k-means solutions vs randomly selected solutions (Fig. 2). ARI values near 1 indicate high cluster consistency, and the mean ARI for the bootstrapped solutions was 0.87 (SD = 0.12), relative to the mean ARI for the random solutions (0.00, SD = 0.01) (P < 0.001). The results of the internal and external bootstrapping simulations indicate that the 3-cluster solution was robust to changes in the data and suggested a replicable cluster structure in the broader population.

Figure 2.

Internal and external validation of cluster solution. Internal validation for a 3-cluster solution by Elbow (A) and Silhouette (B) methods, by the Calinski–Harabasz Index (C), and C-H index calculated by bootstrapping 1000 samples. (D) External validation, as demonstrated by the Adjusted Rand Index, comparing randomly selected solutions (E) vs bootstrapped k-means solutions (F).

As shown in Table 6, we analyzed individual cluster characteristics in variables internal to the algorithm. Cluster 1 included individuals with the lowest mean pain and lowest temporal dependency and instability. Cluster 2 included individuals with the highest levels of mean pain and highest temporal dependency (AR1), whereas cluster 3 included individuals with high levels of mean pain and the highest temporal instability (MSSD and PAC). Clusters 2 and 3 had individuals with the high proportion of days with pain (median of 100% and 83.9%, respectively), whereas, cluster 1 had individuals with the lowest proportion of days with pain (median of 21.4%). The cluster solution is represented in Figure 3.

Table 6

Individual cluster characteristics, in variables internal and external to the clustering algorithm.

Metrics of pain intensity (median [IQR])	Cluster 1 (n=44)Lowest mean pain, lowest temporal instability/dependency	Cluster 2 (n=37)Highest mean pain, highest temporal dependency	Cluster 3 (n=42)High mean pain, highest temporal instability	P *	Post hoc pairwise comparison significant†
Mean	0.59 [0.25, 1.26]	5.04 [4.18, 6.11]	3.27 [2.39, 5.12]	<0.001	‡ § ‖
50th percentile (p50)	0.00 [0.00, 0.12]	5.00 [4.00, 6.00]	3.00 [2.00, 5.38]	<0.001	‡ § ‖
90th percentile (p90)	2.00 [0.60, 3.00]	7.00 [6.00, 8.00]	6.15 [5.08, 8.00]	<0.001	‡ §
Intraindividual standard deviation (iSD)	1.12 [0.73, 1.50]	1.33 [0.94, 1.75]	2.39 [2.07, 2.73]	<0.001	§ ‖
Probability of acute change of 1 point (PAC1)	0.11 [0.07, 0.22]	0.26 [0.15, 0.30]	0.44 [0.37, 0.55]	<0.001	‡ § ‖
Probability of acute change of 2 points (PAC2)	0.07 [0.04, 0.12]	0.07 [0.00, 0.11]	0.33 [0.22, 0.40]	<0.001	§ ‖
Mean square of successive differences (MSSD)	1.50 [0.69, 2.44]	1.89 [1.15, 2.52]	7.57 [5.15, 8.99]	<0.001	§ ‖
First-order autocorrelation (AR1)	0.30 [0.00, 0.44]	0.45 [0.24, 0.65]	0.36 [0.03, 0.46]	0.018	‡ ‖
Proportion of pain days (PPD)	0.21 [0.11, 0.47]	1.00 [1.00, 1.00]	0.84 [0.61, 0.92]	<0.001	‡ § ‖

χ2 (categorical variables) or Kruskal–Wallis (continuous variables) test used as appropriate.

Post hoc testing–Dunn test, 2 sided. Dunn test without adjustment for multiple comparisons.

Pairwise comparison cluster 1 vs 2.

Pairwise comparison cluster 1 vs 3.

Pairwise comparison cluster 2 vs 3.

ED, emergency department.

Figure 3.

Pain phenotypes (clusters) in PiSCES. Visual representation of the cluster solution with description of individual cluster characteristics. PiSCES, Pain in Sickle Cell Epidemiology Study.

We analyzed differences between clusters in variables external to the algorithm (Table 6). We found that the proportion of married individuals differed across clusters. Individuals in cluster 1 tended to be younger (median age 31.5) and have the best HRQoL scores, along with lower levels of catastrophizing, stress and use of affective coping strategies. Cluster 1 had the highest PF scores, whereas cluster 2 had the lowest PF scores. Cluster 1 most consistently differentiated itself from clusters 2 and 3 on subscales of HRQoL and psychological factors. Differences between clusters 2 and 3 were limited to PF and SF, with cluster 3 having higher median scores as compared with cluster 2 on these subscales, but individuals in cluster 3 were younger, with a median age of 31.5 years in cluster 3 as compared with 39 years in cluster 2. Psychological characteristics were similar between clusters 2 and 3, but different from cluster 1. Clusters 1 and 3 were different in the proportion of days with health care utilization. The proportion of days with opioid use was different across all 3 clusters, with cluster 2 having the highest proportion of days with opioid use, despite having lower proportion of days with crisis as compared to cluster 3. Last, in cluster 1, most individuals had <50% days with pain, whereas cluster 2 had almost all individuals with ≥95% days with pain. In cluster 3, most individuals had ≥50% days with pain, but only a small proportion of this cluster had ≥95% days with pain. There was no difference in proportion of patients with avascular necrosis or skin ulcers or median number of SCD comorbidities between the 3 clusters.

4. Discussion

In this study, we present results of a systematic investigation of the role of intraindividual variability in pain intensity, using data from PiSCES, the most contemporary and largest pain diary data set of individuals with SCD. Our findings provide empirical evidence for the association of pain variability, particularly first-order autocorrelation (AR1), with HRQoL and psychological well-being in SCD pain. Our exploratory analysis also shows that there are meaningfully distinct subgroups or pain phenotypes in SCD, which could be identified on the basis of pain intensity, pain variability, and pain frequency.

Previously published results from the PiSCES study demonstrated that mean pain intensity was highly predictive of HRQoL, controlling for other sociodemographic variables such as age, sex, and years of education.[22] In this study, we show that intraindividual pain variability is associated with HRQoL, psychological characteristics, and clinical outcomes. Although mean pain intensity was associated with HRQoL overall, it was less relevant in the subgroup of those with ≥95% days with pain, suggesting that mean pain intensity alone is inadequate in describing outcomes with those with daily or near-daily pain. We found that the relationships between variability and outcomes were different in the subgroup of individuals with more frequent pain. In those with ≥95% days with pain, AR1 seems to be negatively associated with some subscales of HRQoL and suggests that higher AR1, which indicates persistence of pain, may be associated with poorer outcomes. In this subgroup, some fluctuation in pain intensity seems to be associated with improved scores in the GH domain of HRQoL. This is somewhat different from what Schneider et al.[31] postulated, which was that patients may find pain more manageable if it varied less from day to day. However, our study suggests that the relationship of variability and outcomes may be nuanced and that pain variability may have a differential impact based on pain frequency.

We found associations between pain variability and psychological factors in this study. Schneider et al.[31] reported that in rheumatological illnesses and osteoarthritis, increased variability was associated with depression. Zakoscielna et al.[43] reported that in older adults, pain variability was correlated with higher depression. The results of this study contrast with these 2 previously published studies, as we did not find an association between pain variability and depression, although a quarter of individuals in the sample were depressed. We also found that increased pain variability correlated with greater use of affective coping strategies, somatic symptom burden, and sickle cell stress. In addition, even in the subgroup with ≥95% days with pain, higher autocorrelation continued to be associated with sickle cell stress, use of affective coping, and catastrophizing, whereas other measures of variability were not significantly associated with psychological characteristics. This suggests that the tendency for pain to persist, as measured by the first-order autocorrelation, may be a determinant of outcomes in SCD. Our findings are different from the findings of Mun et al.,[24] where higher autocorrelation was associated with higher predictability, and lower autocorrelation, which reflected more fluctuation and unpredictability, was proposed to be associated with lower psychological well-being. We, however, report positive associations of autocorrelation with psychological well-being, namely maladaptive coping strategies, catastrophizing, and stress. Given the nature of our study, it is possible to tell whether autocorrelation is causal of maladaptive coping, catastrophizing, and stress or whether these psychological characteristics lead to higher autocorrelation of pain. Our results provide the rationale for further examination of the relationship of pain-associated psychological characteristics and AR1, especially in those with daily or near-daily SCD pain, and must be closely examined in future studies in both SCD and non-SCD pain. Future studies may also examine whether pain variability indices remain stable over longer periods of time, particularly in those with chronic daily or near-daily pain. This would provide some insight into whether pain variability is a trait or a state and how much it changes over time.

Using unsupervised clustering methods, we identified underlying subgroups or clusters among individuals in this study. These clusters differed based on variables that were both internal as well as external to the clustering algorithm. Cluster analysis methods have been used to identify subgroups among painful conditions based on measures of pain sensitivity and psychological distress in the OPPERA study,[2] based on quantitative sensory testing profiles in neuropathic pain,[5] based on symptoms in fibromyalgia,[41] and based on pain and psychological symptoms in older adults with chronic pain.[18] In this analysis, we found 3 distinct clusters or phenotypes of pain, each of which had different characteristics. It seems that the cluster 2 phenotype with the highest mean pain, highest proportion of days with pain, highest temporal dependency, and highest proportion of days with opioid use tends to fare the worst in terms of HRQoL, particularly PF. It is also interesting that cluster 3 phenotype had the highest temporal instability, but did not have appreciably different scores on measures of psychological function related to pain, and had only differences physical and SF subscales of HRQoL as compared with cluster 2. Using the threshold of ≥50% days with pain, we find that most individuals with ≥50% days with pain were found in clusters 2 and 3, whereas cluster 1 had the most number of individuals with <50% days with pain. Our results may be confounded by age because cluster 2 had a higher median age as compared with cluster 3, but did not have any difference in number of individuals with comorbidities such as avascular necrosis or skin ulcers, which may be associated with persistent pain in SCD. Although our analysis is exploratory and a separate validation cohort is needed for confirmation, our results seem robust and point to the presence of pain phenotypes in SCD and the potential role of pain variability in characterizing these clusters. Despite similar pain frequency, there were notable differences in PAC and MSSD between clusters 2 and 3. Thus, the role of temporal dependency and instability on quality of life in SCD chronic pain bears further study. Future research should also examine the biological characteristics of cluster membership and whether it predicts pain trajectories or response to treatments.

Although pain diary data were collected prospectively in PiSCES, they were collected using paper diaries. Compared with ecological momentary assessments (EMA) through electronic pain diaries, paper pain diaries may be subject to errors and omissions,[26] poor compliance,[26] inflation in retrospective reports,[20] and inaccuracies due to backfilling of entries.[38] Thus, the results of this study should be confirmed with EMA-based studies of SCD pain. Individuals with chronic painful conditions with higher intraindividual pain variability tend to have higher recalled retrospective ratings of pain,[37] and a similar phenomenon is seen among those who experience acute pain.[12] Thus, it is possible that paper pain diary data itself may not reflect the true underlying pain extent of pain intensity or variability. Although most psychological characteristics associated with persistent pain, and measures of HRQoL were included in this study, we were limited by the availability of psychological measures and HRQoL assessments collected at baseline in PiSCES. Being an exploratory, secondary analysis, we presented unadjusted P values, so statistical inferences should be made with caution and confirmed with future studies. Future studies should also investigate the relationship of pain variability with patient-reported outcome measures specific to SCD. We also calculated pain variability after imputation of missing values using the last value carry forward method, the simplest and most conservative option, which may have artificially modified the variability. Using EMA, future studies should also examine the effect of other variables on pain variability, such as time of day of reporting, methods of imputation of missing data, concurrent effects of mood and sleep on pain variability, and the impact of pain variability on pain interference.

5. Conclusion

Intraindividual pain variability is associated with patient outcomes and psychological characteristics in SCD and is useful in delineating phenotypes of pain in SCD.

Conflict of interest statement

The authors have no conflicts of interest to declare.

Supplemental digital content associated with this article can be found online at http://links.lww.com/PAIN/B493.