Pain can’t be carved at the joints: defining function-based pain profiles and their relevance to chronic disease management in healthcare delivery design | BMC Medicine

Human population data source

The UK Biobank is a prospective epidemiological resource that provides a rich palette of information including high-quality brain imaging, genetics, and various biological and lifestyle measurements in a cohort of ∼500,000 participants recruited from across the UK, following ~ 9.2 million invitations sent to individuals registered with the UK’s National Health Service who were aged 40–69 years and lived within approximately 25 miles (40 km) of one of 22 assessment centers located throughout England, Wales, and Scotland [21]. Further information on the consent procedure can be referenced at biobank.ctsu.ox.ac.uk/crystal/field.cgi?id = 200. Analysis of UK Biobank non-identifiable data received Institutional Review Board (IRB) approval from MassGeneralBrigham and McGill IRBs as Not Human Subjects Research. Participants were aged 40–69 years when recruited (mean age 55, standard deviation 7.5 years). We focused on high-resolution T1-weighted structural brain scans the February 2020 data release of 40,681 (54% female) participants (see Supplementary Table 10 regarding demographic information of final sample, described below). UK Biobank data are available through a procedure described at The present analyses were conducted under UK Biobank application number 25163, “Structural genetic contributions brain function and behavior.” Before proceeding to mine the UKB for genetic contributions to pain condition (which, as stated above, have a high overlap with depression and anxiety), the present work is a first step towards more clearly defining population-level patterns of pain behavior that are anchored in biology. So defined, we expect that our pain profiles will facilitate later studies of the genetic architecture of pain conditions and, as focused on in this manuscript, guide clinical operations.

For the benefit of reproducibility and comparability, all our analyses of brain–behavior association were based on the precomputed and expert-vetted image-derived phenotypes [22]. For our analyses of gray matter structure, we relied on volume estimates in 100 cortical regions defined by the Schaefer-Yeo atlas (see Supplementary Table 5) [23]. All structural magnetic resonance imaging (MRI) data were preprocessed using the pipelines and quality-control workflows by the FMRIB team, Oxford University, UK [24]. The uniform preprocessing increases the comparability of our findings to other UK Biobank studies. We used common linear deconfounding to remove variation in all brain-imaging-derived phenotypes that could be explained by interindividual differences in head size and motion, following previously published UK Biobank research [25]. We applied the same normalization procedures to the pain experience variables, described below. As such, effects emerging from the subsequent modeling steps on the thus cleaned brain phenotypes cannot be explained by differences in brain size or motion.

Our primary goal was to define clinically relevant population-level modes of variability that as associated with pain conditions. As it is known that age [26], body mass index [27], and sex [28, 29] are associated with clinical pain in biologically and clinically important ways, we decided not to deconfound these variables from brain imaging-derived target measures with the logic that deconfounding would remove precisely the variability we were hoping to exploit to define clinically relevant population-level modes of co-variation. We further note that, consistent with our previous work, we have not regressed out sociodemographic variables with the goal of preserving rich insights into brain and behavior, by interlocking effects with other variables of interest [30,31,32,33].

Of the 40,681 UK Biobank participants with brain imaging data, 38,692 participants also had pain experience data (see Supplementary Table 1). We used two criteria to account for missing data: (1) we excluded participants with ≥ 90% missing values across all dataFields (i.e., all candidate pain variables) and (2) we excluded specific dataFields with ≥ 75% missing data across all participants (see Supplementary Table 2 and Supplementary Fig. 1 for further details). In addition to age, sex, and body mass index (BMI), a total of 154 pain-related phenotypes across 34,337 participants were available for analysis (cf. Supplementary Table 3). Demographic information for this final 34,337 sample is presented in Supplementary Table 10. We systematically normalized the pain-related variables by recoding each pain experience item from 0 to 10, with 0 being the “best” (reflective of health) and 10 being the “worst” (reflective of disease); see Supplementary Tables 3 and 4 for details and examples. Variable normalization was an important step as the UK Biobank codes used different variables to score similar responses. Following this variable normalization step, we imputed missing data by randomly sampling (with replacement) from existing, non-missing data.

Multivariate pain-brain decomposition

Pattern-learning algorithms, such as partial least squares (PLS) canonical analysis, provide powerful approaches for extracting meaningful patterns from complex, heterogenous variable sets [34]. PLS is a natural choice for a deeply phenotyped dataset such as the UK Biobank for several reasons, including (1) we expected high autocorrelation between the input dimensions, (2) PLS allows for an effective use of data in the context of high-dimensional, phenome-wide data, as analyzed in our present study, and (3) PLS is highly data-efficient by estimating one coefficient corresponding to each input dimension. These 3 reasons also hold for linear discriminant analysis with the difference that linear discriminant analysis (LDA) is targeted at classification settings, whereas PLS is naturally targeted at explaining variation in continuous outcomes. PLS was implemented in Python package sklearn (version 0.21.3) to decompose the collection of > 150 pain-related participant responses (hereafter, pain features) with respect to cortex regional volumes (hereafter, brain features; cf. Supplementary Tables 1–4). PLS is particularly useful when trying to identify patterns in very large and strongly correlated data, such as in clinical populations where patients often report symptom domains or diseases that overlap [35, 36]. In our previous work, PLS was used to study how loneliness manifests itself in the brain’s anatomy [35] and, separately, how clinical risk factors for social isolation predispose people to dementias [36]. In our study, PLS was thus used to identify the most explanatory projections that yielded maximal covariance between sets of brain volumes in the context of participant reports of pain experience, aiming to yield neurobiologically grounded profiles of pain experience.

Intuitively, PLS has the ability to accommodate two variable sets and thus allows the identification of patterns that describe many-to-many relations. Although akin to principal components analysis (PCA), PLS maximizes the linear correspondence between two sets of variables. That is, a key feature of PLS is that it identifies the correspondence between two sets of variables, here capturing two different levels of observation—brain structure measurements (hereafter, brain features) and pain-related indicators (hereafter, pain features). Bringing multiple data modalities to the same table is a process that has previously been referred to as multi-modal fusion (see Yang et al., 2019 for a review [37]) The PLS algorithm, therefore, seeks dominant dimensions that describe shared variation across different sets of measures. In so doing, PLS can re-express a set of correlated variables in a smaller number of hidden factors of variation. These latent sources of variation are not always directly observable in the original measurements, but together explain dominant motifs of how the actual observations are intrinsically organized.

More technically, assuming X and Y are centered, the relationship between the two original variable sets (X and Y) and the derived canonical variates U and V can be understood as the best way to rotate the left variable set and the right variable set from their original measurement spaces to new spaces in a way that maximizes their linear correlation. The fitted parameters of PLS thus describe the rotation of the coordinate systems: the canonical vectors encapsulating how to get from the original measurement coordinate system to the new latent space, the canonical variates encoding the embedding of each data point in that new space. Similar to PCA, PLS re-expresses data in form of high-dimensional linear representations (i.e., the canonical variates). Each resulting canonical variate is computed from the weighted sum of the original variable as indicated by the canonical vector.

Concretely, to analyze the relationship between the brain features and the pain features, all data were systematically normalized (described above) prior to application of PLS to identify low-rank projections which link variable set $X$ corresponding to the brain features, to variable set $Y$ representing the participant-wise pain features. Two sets of linear combinations of the original variables are obtained as follows:

$$X\in \mathbbR^n\times p$$

$$Y\in \mathbbR^n\times q$$

$$\beginarrayccL_X=XV&L_Y=YU\endarray$$

$$\beginarrayccl_x,l=Xv_l&l_Yl=Yu_l\endarray$$

$$corr\left(l_X,l,l_Y,l\right)\propto l_X,l^Tl_Y,l=max$$

where n denotes the number of participants, p is the number of brain features (100), q is the number of pain features (154), $V$ and $U$ denote the respective contributions of X and Y to the derived latent “modes” of joint variation between the patterns derived from X and patterns derived from Y, $L_X$ and $L_Y$ denote the individual-level expression of the respective latent “modes,” $l_X,l$ is the lth column of $L_X$, and $l_Y,l$ is the lth column of $L_Y$. The goal of our pattern learning approach was to find pairs of latent vectors $l_X,l$ and $l_Y,l$ with maximal correlation in the derived latent embedding and quantify the strength of the relationship between the two component sub-groups. Since PLS was purely used as an exploratory analysis, uncertainty in effect sizes were not measured. This step was performed in Python, using SciKit-learn’s PLSCanonical [38]. Statistically significant pain-brain modes of co-variation were defined as p < 0.001, based on 1000 permutations (see Fig. 2C), as described below.

To assess the statistical robustness of the PLS-derived modes of pain-brain co-variation, we created an empirical null distribution following a previously published non-parametric permutation procedure [34, 39]. In each of the 1000 permutations, the brain feature matrix remained constant, while the pain feature matrix was randomly shuffled. This approach preserved the statistical structure of the brain features while selectively disrupting the structure of the pain features. The resulting distribution represented the null hypothesis of a random association between brain and pain features across participants. For each permutation iteration, four variables were recorded: the R score for each PLS mode (yielding 10 scores per iteration), pain feature loadings (one per pain feature per iteration), and brain feature loadings (one per brain feature per iteration). These values were used to create null distributions against which the statistical significance of our actual results could be assessed. For example, the Pearson correlation coefficient (rho) was recorded and contributed to the null PLS distribution. P-values were derived based on the 1000 rho values from the null PLS model. Multiple comparisons correction was not required, as the study focused exclusively on the leading structural brain variation mode estimated by PLS, which was statistically significant at p < 0.05.

Given the need to impute missing data (as described above), we conducted a sensitivity analysis to ensure that each pain profile’s characterization was not biased by imputed data. Following the same logic used to define the null PLS distribution, we generated an empirical null distribution for model weights for both pain and brain variables by applying the weights from the PLS permutation analysis (as shown in Fig. 2C). This distribution, shown in Fig. 3 (left column), allowed us to identify pain variables that were significantly loaded, defined as those in the top or bottom 5% of the empirical null distribution based on magnitude (shaded in grey in Fig. 3, left column). We then counted (Fig. 3, middle column) and summed (Fig. 3, right column) the weights of only the significant loadings within each domain, which collectively characterized each profile.

As a further validation of our pain-behavior profiles, we performed a split-half analysis (analogous to our previous UK Biobank studies (Kiesow 2020, Saltoun 2023)). We chose a split-half analysis given that our goal was to leverage the UKB’s highly unique choice (real-world) and extent (quantity) of pain variables to define clinically actionable pain profiles. We are unaware of any even remotely similarly rich pain neuroscience dataset, which makes an external validation challenging. As such, we performed a split-half analysis to further quantify the stability of our pain profiles. We randomly split the 34,336 participants into two halves across 1000 iterations, using one half to repeat our analysis pipeline and another for a second repetition of our analysis pipeline. To determine the stability of the 4-profile solution in each of these iterations, we then computed the correlation coefficient (rho) between the profiles as derived afresh from the first and second half of the cohort.

Pain profile clinical features

To investigate which clinical features were associated with each pain profile, we examined the potential associations between individual-level expressions of pain-brain modes and participants’ phenotypic, diagnostic, and medication histories. For this, we used an individual-level output from the PLS analysis, referred to as the “PLS Pain Score.” The PLS Pain Score is represented by an (n x q × 4) weight matrix, where n is the number of individuals, q is the number of pain features (154), and 4 is the number of significant loadings. Each weight value in this matrix indicates the relationship of an individual’s pain features to each of the four pain modes. For example, a higher loading (or PLS Pain Score) for a participant on Mode 1, which represents Pain Interference, suggests that this participant has a stronger expression of pain features associated with Pain Interference.

None of the clinical features (phenotypic, diagnostic, or medication histories) were included in the original PLS analysis. Each association with individual-level clinical histories therefore represents an independent profile of each pain mode, and thus an independent validation of the meaningfulness of the derived pain profiles.

Phenome-wide association study (PheWAS) of participants’ derived pain scores were defined by reference to a wide variety of almost 1000 lifestyle factors, demographic indicators, and laboratory values. Saltoun et al. recently published this tool to characterize brain asymmetries [40]. Feature extraction was carried out using two utilities designed to obtain, clean, and normalize UK Biobank phenotype data according to predefined rules. Initially, a raw collection of ~ 15,000 phenotypes from 11 major categories was available for analysis. We used the FMRIB UK Biobank Normalisation, Parsing And Cleaning Kit (FUNPACK version 2.5.0; to harmonize and clean these phenotypes using predefined, standard arguments available within the FUNPACK workflow. Cleaning steps included excluding any brain-imaging-derived or mental health-related variables (as these were both related to the PLS analysis), and deploying built-in tools to remove “do not know” responses and unasked dependent data. FUNPACK therefore yielded ~ 3300 phenotypes, which was then input into PHEnome Scan ANalysis Tool (PHESANT), a tool designed specifically for curating UK Biobank phenotypes [22]. The PHESANT toolkit combined phenotypes across visits, normalized, cleaned, and categorized the data as belonging to one of four datatypes: categorical ordered, categorical unordered, binary, and numerical. All categorical unordered columns were converted into binary columns to encode a single response. For example, the employment status phenotype was originally encoded as a set of values representing different conditions (e.g., retired, employed, on disability). The output of categorical one-hot encoding on unordered phenotypes was then combined with all measures classified by PHESANT as binary, numerical, or categorical ordered. By default, FUNPACK assessed pairwise correlation between phenotypes and discarded all but one phenotype of a set of highly correlated (> 0.99 Pearson’s correlation rho) phenotypes. For example, left and right leg fat percentages were highly correlated (Pearson’s rho 0.992). Hence, only right leg fat percentage was included in the final set of phenotypes. Also by default, both FUNPACK and PHESANT exclude variables with fewer than 500 participants; we kept this default. There were 11 categories available for analysis, but because some of the mental health phenotypes may overlap with pain-related responses, we removed all 220 phenotypes belonging to the “mental health” category. Overall, 757 phenotypes from 10 major categories were available for PheWAS.

Each PheWAS examined a battery of associations with subject-wise pain scores with the 757 extracted phenotypes. By computing Pearson’s correlation for each phenotype and pain score pair, we were able to extract both the association strength and accompanying statistical significance of the given phenotype-pain profile pairing. For each profile, two standard tests were used to adjust for the multitude of associations being assessed. First, we used Bonferroni’s correction for multiple comparisons, adjusting for the number of tested phenotypes (0.05/757 = 6.6e-5). Second, we further evaluated the significance our correlation strength using the false discovery rate (FDR), another popular method of multiple comparison correction which is less stringent than Bonferroni’s correction. The false discovery rate [41] was set as 5% and computed for each asymmetry pattern in accordance with standard practice [22, 42]. For visualization purposes, phenotypes in Manhattan plots were colored and grouped according to the category membership defined by FUNPACK.

Diagnosis-wide Association Study (DiWAS) of participants’ derived pain scores were defined by reference to a wide variety of almost 1,500 disease phecodes extracted from participants’ electronic health records. Feature extraction for both ICD-9 and ICD-10 diagnoses was carried out using the FMRIB UK Biobank Normalisation, Parsing And Cleaning Kit (FUNPACK version 2.5.0; to extract ICD codes. Once extracted, ICD-9 and ICD-10 codes were grouped into 1477 hierarchical phenotypes which combine related ICD-9 and ICD-10 codes into a single “phecode,” using previously established definitions [43]. In addition to combining related ICD codes, we use built-in exclusion criteria to reduce case contamination of control groups [44]. In our UKB subset of people with both pain variables and brain imaging (n = 34,337), 22 diseases were not represented in our subset. We excluded these 22 ICD codes from the final analysis. The final set of 1425 diagnostic phecodes was later compared to discovered pain profile expression to probe for relations between patient-reported pain profiles and clinical diagnoses.

More concretely, we used FUNPACK on the UK Biobank sample to extract phenotype information related to ICD-9 diagnosis codes and ICD-10 diagnosis codes. The FUNPACK utility contained built-in rules tailored to the UK Biobank which resulted in the extraction of a total of 542 ICD-9 codes and 4105 ICD-10 codes. ICD-10 codes contain more granular information than ICD-9 codes, and separate ICD-9 and −10 codes may reflect common etiologies [45]. To consolidate related medical diagnoses across and within ICD coding paradigms, we used the Phecode Map 1.2 with ICD-10 codes made available from https://phewascatalog.org/phecodes_icd10 [43]. The phecode system additionally established clinically guided exclusion criteria relevant to creating a health comparison cohort associated with each phecode. As a specific example, the exclusion criteria associated with “Type 2 diabetes” (phecodes 250.1) includes individuals with related diseases, including Type 1 diabetes and secondary diabetes mellitus, as well as patients with symptoms commonly associated with Type 2 diabetes, such as abnormal glucose [43]. Individuals with related diseases, but not the disease of study, are therefore excluded from the control group based on the built-in exclusion criteria, which may otherwise decrease the statistical power for finding genotype–phenotype associations [44]. The final set of 1425 clinical phecodes spanned 17 disease classes as defined by Wu and colleagues [43], ranging from congenital anomalies and neoplasms, to mental disorders and infectious diseases.

We charted robust cross-links between the subject-wise expression of a given pain profile and the 1425 extracted phecodes, with appropriate correction for multiple comparisons. For each extracted phecode, the Pearson correlation between the given phecode and the inter-individual variation in pain archetype revealed both the association strength and accompanying statistical significance of the given phenotype-pain archetype association. Individuals fulfilling the exclusion criteria defined by Wu et al. (2019) were removed from the control group, separately for each phecode. For each uncovered pain profile, two standard tests were used to adjust for the multitude of associations being assessed. First, we used Bonferroni’s correction for multiple comparisons, adjusting for the number of tested phenotypes (0.05/1425 = 3.5e − 5). Second, we further evaluated the significance of our correlation strength using the false discovery rate (FDR), another popular method of multiple comparison correction which is less stringent than Bonferroni’s correction. The false discovery rate [41] was set as 5% [22, 40] and computed for each pain archetype in accordance with standard practice [42]. For visualization purposes, phecodes in Manhattan plots were colored and grouped according to the category membership defined by the phecode map developed by Wu and colleagues (2019).

Medication-wide Association Study (MedWAS) of participants’ derived pain scores were defined by reference to the Anatomical Therapeutic Chemical (ATC) classification published by the World Health Organization (www.who.int/tools/atc-ddd-toolkit/atc-classification). The ATC classification groups active substances according to the organ or system on which they act and to their therapeutic, pharmacological, and chemical properties. Drugs are classified in groups at five different levels, based on target organ (level 1), therapeutic subgroup (2), pharmacological subgroup (3), chemical subgroup (4), and chemical substance (5). We chose to use level 3 in our analyses as we were interested in drug category. The UKB includes 6745 different types of drug or supplement entries across 374,891 participants (listed under Data-Field 20,003). Using in-house code, each of 6745 medications was searched for its corresponding ATC codes at each ATC level. This search was performed solely on the medication name and not on the medication dose (see Supplementary Methods, Supplementary Table 7). Because medication dose was not collected by the UKB team (see 12.4, all medications were collapsed into medication name independent of dose and formulation. Individual medications were subsequently collapsed into ATC codes thereby providing another step which summarized medications on the individual level. Medication frequency (i.e., how many times per day a medication was taken), treatment duration (how many days a medication was taken), or whether a medication was taken at the time of interview (i.e., that day) was not available for analysis.

Because some medications fall into more than one ATC category depending on which dose or route of administration is used (e.g., finasteride, see whocc.ni/atc/structure_and_principles/), we selected the first ATC category in alpha-numeric order. Medications listed in the UKB that did not have corresponding ATC codes were excluded from analysis (e.g., many dietary supplements, free text entries). Medications which were combinations of two drugs (e.g., budesonide and formoterol are found together as an inhalant for COPD) were excluded from analysis. Given that our analysis revealed an association between opioids and pain interference, we noted that the ATC classification of codeine was unclear in some instances. Codeine is an opioid that (in combination with paracetamol or aspirin or ibuprophen) can be used as an analgesic; however, because it is a weak opioid, it is most often used as a cough suppressant (R05D). To not bias our results, we insured that codeine was not included under the larger category of N02A Opioids. We did note that insulin was listed in the UKB as “insulin product” (coded as 1,140,883,066); because this did not specify which type of insulin, we coded this ATC Level 3 Category “A10A Insulins and Analogues.” Overall, the UKB medication data represents 14 ATC Level 1 Domains (all of them), 73 Level 2 categories, and 137 Level 3 categories. At the individual participant level, medication profiles were coded as follows: the total number of medications a participant was prescribed within each ATC category was represented as a whole number. For example, a participant who reported taking escitalopram, bupropion, and oxycodone would correspond with an ATC Level 3 profile would be a 137-place vector with a “N06A Antidepressant” value of 2 and an “N02A Opioid” value of 1.

Each MedWAS examined a battery of associations with subject-wise pain scores with the 137 ATC Level 3 codes. By computing Pearson’s correlation for each ATC code and pain score pair, we were able to extract both the association strength and accompanying statistical significance of the given medication-pain profile pairing. For each profile, two standard tests were used to adjust for the multitude of associations being assessed. First, we used Bonferroni’s correction for multiple comparisons, adjusting for the number of tested phenotypes (0.05/137 = 3.64e − 4). Second, we further evaluated the significance our correlation strength using the false discovery rate (FDR), another popular method of multiple comparison correction which is less stringent than Bonferroni’s correction. The false discovery rate [41] was set as 5% and computed for each asymmetry pattern in accordance with standard practice [22, 42].

Regarding opioid and antidepressant versus pain interference and mood disorder,we wanted to better understand the observation that opioids and antidepressants were associated with different profiles of pain interference and mood disorder. Across the 502,507 participants in the UK Biobank, we restricted our analysis to those who had data about medication, pain interference (the Brief Pain Inventory, or BPI), and mood (the Patient Health Questionnaire 9, or PHQ-9). Similar to our original criteria, we excluded participants with either ≥ 90% missing data and further excluded dataFields with ≥ 75% missing data. This resulted in 167,203 participants. We then imputed missing data by randomly sampling without replacement using the same procedure reported above. To facilitate visualization, for each individual, we then used established methods to sum, round, and bin the total score: the BPI into mild (1–4), moderate (5–6), and severe (7–10) categories; the PHQ-9 into minimal (1–4), mild (5–9), moderate (10–14), moderately severe (15–19), and severe (20–27) categories. These data were visualized as a matrix of pie charts.

Comparison of pain profiles with classic body part framework

We next sought to provide evidence that the four pain profiles add information beyond the classic body part framework. We used LDA to train and test models that predict individuals’ dominant pain profile or painful body part (class definition described below) based on behavioural and lifestyle factors used in the PheWAS. LDA was a natural choice of method as it allows the prediction of multiple classes from a large set of features which, in this case, comprised brain structure and real-world behavioural and lifestyle measures. Below, we describe class definition and performance measurement.

Class definition: painful body part

We converted the UKB’s pain location question into body part labels as described by the International Association for the Study of Pain (IASP) Axis 1: Regions. Given that our intent was to associate each participant’s phenotypes with the body part of worst pain, whenever an entry in Data-Field 120,037 (“area most bothered by pain in last 3 months”) was present, this was used as default. When 120,037 was not present, we mapped UKB’s Data-Field 120,021 (“Pain or discomfort all over the body in the last 3 months”) to IASP’s 900 (“More than three major sites”). We further mapped all participants who had more than a single body part identified in the Data-Field 6159 as IASP’s 900. We did this because, except for the 69,917 participants who responded to question 120,037, the UKB did not identify which body part was most painful, but instead simply listed them in numerical order. 157,042 people identified > 1 painful body part. For those with only two body parts identified, we selected as the primary location UKB response 6159–0.1 (or the first response). Mappings of these codes are presented in Supplementary Table 8 [46]. Counts shown in Supplementary Table 8 are based on the UKB’s Data Showcase website, accessed 3/20/2023: https://biobank.ndph.ox.ac.uk/showcase/field.cgi?id=6159.

Class definition: dominant pain profile

Individuals were classified according to the pain profile which most strongly reflects their experience of pain. As not all individuals are expected to be in clinically relevant amounts of pain, individuals which do not strongly express any pain profiles were categorized into a separate no pain class. The proportion of individuals in the no pain class was informed by the number of people who experience no pain in the location-based approach.

By construction, the developed pain profiles rely on pain-location information (cf. above). As we were interested in comparing pain experience to pain location, we elected to use modified pain scores which exclude pain location information in the LDA analysis. To do so, we zeroed-out all data fields which recorded pain location before applying the pain profile definitions obtained through the PLS analysis ($U$ in the equations above). The resultant “pain-agnostic” individual-level pain scores maintain the correlation structure between data-fields established from the full dataset, but do not contain information about individual-level pain location.

To facilitate the interpretations of individual level pain scores, we flipped the signs of pain profiles such that the resultant pain-agnostic pain scores are positively correlated with increasing amounts of location-based pain. Through this procedure, we are ensured that increasing individual level scores are associated with increasing bodily pain, though the pain-agnostic pain scores themselves do not contain any body-location pain information. Note that flipping the sign does not change the magnitude of covariance between pain variables, eigenvalues of the pain profiles, or interindividual variation in pain score. Flipping the sign is purely to aid interpretation.

All individuals express all pain profiles to various extents. To assign individuals to particular pain profile which most closely reflects their experience of pain, all individuals were categorized according to their largest pain score. For example, an individual that expresses the pain interference pain profile more strongly than any of the other three pain profiles would belong to the pain interference class through this procedure. Through this approach, everyone belongs to one of four pain profile classes. However, we know that not all individuals are in clinically relevant amounts of pain. We therefore leveraged information from the pain-location based approach, wherein 34.0% of individuals reported no pain in any location, to establish an additional class of individuals with no pain. The 34.0% of individuals with the smallest pain scores across all pain profiles were grouped into a new “no pain” class instead of their respective “strongest” pain profile. These resultant 5 classes therefore reflect the pain/no-pain distribution seen in the location-based approach and categorize individuals according to their experience of pain rather than where it is located.

Cross-validation and model performance

To improve the generalizability of the LDA models, we used a fivefold out-of-sample cross validation wherein the available data were divided into 5 equal samples, trained on 4/5, and tested on the left out 1/5. Class distribution was maintained across all splits, such that all classes were represented in the held out test set for all five splits. Given substantial class imbalance observed in dominant pain profile and painful body part membership (see above), data was resampled to balance the classes prior to fitting the LDA, and was tested on held-out data with and without resampling. For each training iteration, 10 separate resampling protocols were performed for a total of 50 iterations across splits and samples. Performance was measured as accuracy (computed as correctly predicted class labels divided by the total number of labels). We further report Area Under the Curve (AUC) using a micro-averaging approach (AUC(micro)). AUC(micro) is a common method for evaluating multi-class classification performance. AUC(micro) calculates AUC by pooling together true positive rate and false positive rates across all classes using class size as weights [47,48,49].

link