Knee osteoarthritis and time-to all-cause mortality in six community-based cohorts: an international meta-analysis of individual participant-level data

Background Osteoarthritis (OA) is a chronic joint disease, with increasing global burden of disability and healthcare utilisation. Recent meta-analyses have shown a range of effects of OA on mortality, reflecting different OA definitions and study methods. We seek to overcome limitations introduced when using aggregate results by gathering individual participant-level data (IPD) from international observational studies and standardising methods to determine the association of knee OA with mortality in the general population. Methods Seven community-based cohorts were identified containing knee OA-related pain, radiographs, and time-to-mortality, six of which were available for analysis. A two-stage IPD meta-analysis framework was applied: (1) Cox proportional hazard models assessed time-to-mortality of participants with radiographic OA (ROA), OA-related pain (POA), and a combination of pain and ROA (PROA) against pain and ROA-free participants; (2) hazard ratios (HR) were then pooled using the Hartung–Knapp modification for random-effects meta-analysis. Findings 10,723 participants in six cohorts from four countries were included in the analyses. Multivariable models (adjusting for age, sex, race, BMI, smoking, alcohol consumption, cardiovascular disease, and diabetes) showed a pooled HR, compared to pain and ROA-free participants, of 1.03 (0.83, 1.28) for ROA, 1.35 (1.12, 1.63) for POA, and 1.37 (1.22, 1.54) for PROA. Discussion Participants with POA or PROA had a 35–37% increased association with reduced time-to-mortality, independent of confounders. ROA showed no association with mortality, suggesting that OA-related knee pain may be driving the association with time-to-mortality. Funding Versus Arthritis Centre for Sport, Exercise and Osteoarthritis and Osteoarthritis Research Society International.


Introduction
The prevalence of musculoskeletal disorders (not including back pain) was ranked 19th for men and 20th for women in the 2017 Global Burden of Disease study. Knee OA made up 20% of this musculoskeletal burden. In terms of living with disability, musculoskeletal disorders ranked 10th and 11th for men and women, respectively [1]. The lifetime risk of knee osteoarthritis is estimated to be 45% [2], and the prevalence of knee OA is expected to rise in accordance with the increase in the ageing population and obesity epidemic in many parts of the world.
OA is a common debilitating joint disease, frequently associated with joint pain, functional limitation, and decreased quality of life [3]. It most commonly affects the knees, hips, hands, facet joints, and feet [4], with knee and hip OA causing the greatest burden to the population, as 1 3 pain and stiffness in these large weight-bearing joints often lead to significant physical dysfunction such as knee muscle weakness and limited flexion [5].
This variation in findings reflects differences in populations studied (clinical or general), the diagnostic methods used to define OA, statistical methodology used, and the use or inclusion of important confounders in each study. Traditional meta-analyses are valuable and efficient in terms of time and resources required, but do have several limitations, which have been widely recognised [17][18][19] including reliance by necessity on published data increasing the potential for publication bias as negative studies difficult to publish. Aggregate data are often not available, poorly reported, derived, and presented differently across studies (for example, odds ratio versus relative risk), and most studies vary in their definitions of exposures, confounders and outcomes [20].
Individual patient-level (IPD) meta-analysis utilises original raw data from cohorts and uses standardised statistical methods to analyse and produce pooled estimates [21]. IPD meta-analysis, although time-consuming and resource intensive, does not depend on previously published data, allows for a standardised definition of important variables and can be analysed using the same statistical approach. Within the current study, key measures of OA and relevant confounders are harmonised (based on expert consensus) [22], and consistent methods of analyses are used between cohorts to provide a more generalizable estimate of the association between OA and premature mortality in the general population.
This study seeks to overcome the limitations introduced when using aggregated results by gathering and analysing individual participant-level data from multiple international observational osteoarthritis cohort studies to describe the association between knee osteoarthritis and time-to all-cause-mortality.

Study design
This study was designed to assess the relationship between knee osteoarthritis and time-to all-cause-mortality in multiple, prospective, longitudinal, community-based cohort studies from around the world. Subjects were stratified by the presence or absence of osteoarthritis at baseline, and time-to-mortality was compared between groups. Pooled estimates were produced using a two-stage individual participant-level meta-analysis framework consisting of two discrete steps: (1) analysing the individual cohorts separately; and (2) applying traditional meta-analysis methods to produce a pooled effect size [21].
A two-stage analysis can more easily handle cohort-specific characteristics such as heterogeneous populations, different risk relationships (such as direction and shape), and the effect of confounders, and can more overtly handle both sporadic and systematic missing data, unlike a one-stage analysis (i.e., pooling all data) [23]. A two-stage analysis allows for consistently defining the primary risk factors, outcome variables, adjusting for the same confounders, and using consistent statistical methods before producing a single pooled effect size. Unlike a traditional meta-analysis, it also allows for the inclusion of previously unpublished data.

Cohort and participant inclusion/exclusion criteria
Due to the type of data required (detailed pain and radiographic data), and the desire to use cohorts, including those which had not been previously published on the OA/mortality relationship, we identified cohorts using two sources: (1) published literature of cohort studies on knee osteoarthritis and mortality; and (2) contacting principal investigators of longitudinal osteoarthritis cohorts to see whether mortality data had been collected. We did not conduct a traditional systematic review, and as evidenced by the three recent systematic reviews and meta-analyses, several of the cohorts we have included in our study would not have been identified [7][8][9].
The inclusion criteria for cohorts were: (1) OA-related knee pain and knee radiographic data available at baseline for both OA and non-OA subjects; (2) time-to-mortality follow-up data for all participants; and (3) recruitment from the community (i.e., not identified through clinics, hospitals, or healthcare professionals). Exclusion criteria were: (1) cohorts where raw data could not be released for analysis; and (2) data not available for both OA and non-OA subjects. Cohorts were not selected with regard to previously published data on the relationship between OA and mortality.
We identified 40 cohorts via the two previously described sources as potentially having knee osteoarthritis data from the general population. Eighteen were excluded due to being a non-observational cohort or non-community-based or a case-control study. Thirteen lacked the appropriate knee X-ray or pain data at baseline after more detailed investigation, and two lacked available mortality or time-to-death data. Seven potentially eligible cohorts were identified, one of which had data access limitations, leaving six cohort studies available for analysis (see flow chart, Appendix 1). The six cohorts included were: three US community-based cohorts (Framingham and Johnston County Osteoarthritis Project) [24,25], one of which was enhanced for OA risk factors [Multicentre Osteoarthritis Study (MOST)] [26]; one community-based cohort from the United Kingdom (Chingford) [27]; one Chinese community-based cohort (Wuchuan) [28]; and one Australian community-based cohort [The Tasmanian Older Adult Cohort (TasOAC)] [29]. All cohorts provided data for all participants except Framingham which provided a random sample of 80%.
Key differences between cohorts (Appendix 2) are the year of baseline visit, length of follow-up, the baseline age of participants, and the lack of side-specific pain in a single cohort. Participants were included in the analysis if they were over 45 years of age, did not have evidence of rheumatoid arthritis, and had mortality data available. After initial data checks, subjects above the age of 80 were also excluded due to the extremely small numbers available (Appendix 2).

Data collection process
IPD was requested from the principle investigators of any identified cohort after submitting an analysis plan for their team to review. Principle investigators were also contacted directly in cases where data had never been previously released to outside research teams.
A subset of the full data containing only the pre-specified exposures, outcomes, and confounders was requested, transferred via encrypted online servers, and stored and managed centrally by the Oxford research team. A open email dialogue was maintained with principle investigators and key researchers from each cohort throughout the process of data acquisition, harmonisation, and analysis to ensure consistency between cohorts.

Primary risk factor: knee osteoarthritis
Due to the importance of using a consistent definition of osteoarthritis to avoid misclassification, we gained expert opinion on methods to harmonise knee osteoarthritis variables in prospective OA cohort studies, and all OA criteria used in this analysis were defined following a process of expert consultation, analysis, and agreement [22]. The key output of this meeting supported the use of both a binary self-reported pain question and the presence of radiographic OA to define knee OA in the general population. Thus, knee pain was defined using either an NHANES-type question (i.e., 'have you had pain for at least a month in the last month in your joint'), or a similar alternative pain question if an NHANES-type question had not been used to assess pain [30,31]. In cases where only WOMAC was available a threshold of 3 was used on the WOMAC pain subscale, this threshold was determined by the previous expert consensus and external validity study [22]. Radiographic OA was defined using the Kellgren and Lawrence (K/L) scoring method, grade 2, or above, and alternatively, an equivalent combination of radiographic features (osteophytes and joint space narrowing) from other validated scoring methods (such as the OARSI atlas) [32,33].
Subjects were divided into four categories: (1) no knee pain or radiographic OA (Pain-/ROA-); (2) radiographic OA with no pain (ROA); (3) knee pain with no radiographic OA (POA); (4) pain and radiographic OA (PROA). Person-level OA was calculated by assessing the OA status for each joint and using the 'highest' level of OA based on this system. For example, if a subject had no knee pain or radiographic OA (cat 1) in their right knee and radiographic OA with no pain (cat 2) in their left knee, their person-level knee OA status would be radiographic OA with no pain (cat 2).

Primary outcome: time-to-mortality
Each cohort contained a status variable (dead/alive) and a time-to-censoring variable for each participant. Three cohorts (Chingford, Johnston County, TasOAC) determined the date of death using nationally linked records, while the remaining cohorts used other methods to determine the date of death such as updates from Primary Care systems, death registries or municipal administration, family, medical records, and periodic examinations or contacts.
In cohorts where subjects were lost to follow-up at an unknown date, the previous visit when subjects had data was used as the last date where mortality status was known. Time-to-status was calculated from the baseline visit, determined by when knee X-rays and pain were assessed, to the last date that the subject's status was known. Survival was calculated using person-years attributing to the analysis.

Potential confounders
The potential confounders accounted for in this analysis were: age; sex; race; BMI; smoking; drinking; cardiovascular disease (CVD); and diabetes. These were based on clinical applicability and consistent availability across each cohort. To be modelled consistently between cohorts, variables were categorised into the broadest level of information available in any single cohort. For example, one cohort contained detailed data on the lifetime use of all tobacco products enabling the generation of a 'dose', while another cohort simply asked whether they were current, former, or never smokers. This second option was then generated for each cohort. Pain medication, such as NSAIDs, was not considered a potential confounder in this analysis, as it is on the causal pathway between painful OA and mortality, and a mediation analysis on this scale would not have been feasible due to both limitations in the data and in the methodology.
Age was defined as age at the time of baseline clinic visit when OA variables were assessed. Race was included as a potential confounder for any cohort which had more than one race category. Chingford, TasOAC, and Framingham have predominantly Caucasian participants; Johnston County and MOST have both Caucasian and African American subjects; and Wuchuan has predominantly Chinese subjects. BMI was calculated for each cohort using height and weight variables (weight/height in metres [2]). Extreme values were identified in several cohorts; however, due to the wide variety of subjects found in our dataset, we only excluded impossible (i.e. outside any known values) rather than improbable values. Smoking, Alcohol, Diabetes, and CVD were all generated as binary variables.
Smoking was calculated with current/ former smokers and never smokers. Alcohol was grouped by more than one drink per week versus none or one drink per week. Diabetes was based on the presence of self-reported clinically diagnosed diabetes, and CVD was calculated using self-reported responses to previous ischaemic heart disease, and general heart problems.

Statistical methods: descriptive statistics
Descriptive statistics [percentages, means (standard deviations), and medians (inter-quartile ranges)] were calculated for baseline characteristics of all cohorts using all available data. The difference between baseline characteristics in subjects with and without complete data (OA and confounders) was calculated using t tests (or Wilcoxon-Mann-Whitney) for continuous variables and Chi-square tests (or Fisher's exact) for binary and categorical variables. Descriptive statistics for baseline characteristics and time-to-mortality data were stratified by the categories, no pain/no ROA, POA, ROA and PROA.

Statistical methods: missing data
There were three potential types of missing data to consider within our analyses. To identify data that was missing at random (MAR) and missing completely at random (MCAR), we tested patterns and predictors of missingness for all exposures and potential confounders. We identified several MAR variables and ensured to include any required predictors in the imputation model. All other variables were assumed to be MCAR, a non-testable assumption. There were also systematically missing variables, which were missing in their entirety in a single cohort. Appendix 2 shows the systematically missing and MAR/MCAR variables for each cohort.
Multiple imputation with chained equations (MICE) was used to impute any missing data for both the primary risk factor and for confounders [34,35]. Systematically missing variables (i.e., variables which were missing in their entirety) were excluded from all models and analyses.
Participants with missing mortality data were excluded from all analysis (cohorts had no more than three percent missing mortality data). The Nelsen-Aalen estimator was used to approximate the baseline hazards in the imputation models [36]. Variables used for the imputation models were congruent with the analysis model described in the next section. Missing PROA and race were modelled using multinomial logistic regression; BMI by linear regression; sex, smoking, alcohol, CVD, and diabetes by logistic regression. Age was modelled by predictive mean matching due to non-normality from being restricted between ages 45 and 80 [37].

Statistical methods: survival analysis
Cox proportional hazard regression models were used to estimate hazard ratios (HRs) and 95% confidence intervals (95% CIs) between three OA categories (POA, ROA, PROA) and the time-to all-cause-mortality using no pain/no ROA as the comparator group for each analysis. Three models were run: (1) univariable models assessed OA alone; (2) adjusted for age, sex and race; (3) adjusted for age, sex, race, BMI, smoking, alcohol, CVD, and diabetes. Models run in the Johnston County cohort also included a variable for recruitment wave. Several cohorts were systematically missing key potential confounders (primarily smoking and alcohol) (Appendix 2).
To satisfy the assumptions of the Cox proportional hazard model linearity was assessed between continuous variables (age and BMI) and time to death using fractional polynomials and kernel. The proportional hazard assumption of the primary risk factor (OA) was tested using Schoenfeld residuals. Due to the violation of this proportionality assumption, Johnston County was truncated to the 13-year follow-up post hoc, which was the maximum follow-up time of one of the recruitment waves. This corrected the violation of proportionality for the PROA variable, although reduced the power of this cohort. A priori interactions of OA and age, and OA and BMI were tested in all cohorts.

Statistical methods: individual participant data analysis
Individual participant-level meta-analysis methods were utilised, using a two-staged approach [21,38]. In the first stage, hazard ratios (HR) and 95% confidence intervals (CI) were first produced for each individual cohort. Data were pooled in the second stage using random-effects analysis, using the Hartung-Knapp estimation to account for uncertainty around the tau statistic [39,40].

3
The Stata admetan command was used to produce the pooled estimates in addition to forest plots which graphically demonstrate the results [41]. All analyses were conducted using Stata version 13·0 statistical software (StataCorp, College Station, Texas, USA).

Role of funding source
Versus Arthritis UK (formally Arthritis Research) had no role in study design, data collection, data analysis, data interpretation, or writing of the report. Members of the PCCOA steering committee from Osteoarthritis Research Society International (a non-profit scientific organisation) had roles in study development and interpretation as outlined in the author contribution section with all contributors named in the writing group. The corresponding author had full access to all the data in the study and had final responsibility for the decision to submit for publication ( Fig. 1).

Results
10,723 participants in six cohorts from four countries were included in the analyses. All cohorts had less than 3% missing mortality data and less than 12% missing risk factor or confounder data. Participants with missing mortality data were excluded, whilst those missing risk and confounder data were included in imputed analyses. In several cohorts, there was a statistically significant difference in OA, age, BMI, diabetes, and CVD in subjects with and without missing data (Appendix 3). Table 1 shows the baseline demographics for all cohorts stratified by baseline OA. Median follow-up for this analysis ranged from 5·6 to 20·0 years after baseline. There was substantial variability in the baseline age (54.3-62.7 years), BMI (22.5-30.7 kg/m 2 ), prevalence of PROA (6.7-33.3%), and the duration of follow-up in each cohort, such that the percentage of subjects that died in each cohort ranged from 2.9 to 22.3% ( Table 2).

Key results
This individual participant-level meta-analysis of over 100,000 people from four countries revealed that participants with knee pain only, or a combination of knee pain and radiographic OA, had an increased association with reduced time-to-mortality, independent of age, sex, and race (HRs of 1.36-1.44). To explore whether the association could be explained by co-morbid conditions, the models were further adjusted for BMI, smoking, alcohol, CVD, and diabetes. The results remained consistent with HRs of 1.35 for those with POA and 1.37 for those with PROA, compared to participants without knee pain and ROA (pain-/ROA-). Interestingly, we did not observe an association with time-to-mortality in participants with radiographic changes alone (ROA), suggesting that it is pain or some functional consequence of pain such as walking disability or reduced physical activity, rather than the structural aspect of knee OA, that may be driving the increased association with premature mortality [42,43]. While many studies have found an association between OA-related pain and premature mortality, the potential pathways that explain this association is still unclear. A study using large population-based data sets to investigate the effect of pain phenotype on the association between pain and mortality found that the impact of pain in daily life was more important than the presence or extent of pain in the relationship between pain and mortality [44]. Findings from one of the same cohorts examining the potential mechanisms between OA and all-cause mortality highlighted frequent walking as a potential target to reduce all-cause mortality. While anxiety, depression, and unrefreshed sleep had statistically significant effects, the extent of their mediation effect had low clinical significance [45].

Results in the context of other studies
Three recent meta-analyses found no association between OA and mortality, with pooled effect sizes of 0.91 (0.68, 1.23), 1.06 (0.88, 1.28), and 1.21 (0.82, 1.78) in a knee only analysis [7][8][9]. All three articles combined results of studies which used multiple forms of OA diagnosis, including clinician diagnosed OA, self-reported clinical diagnosis, pain, and radiographic OA, increasing the measurement error of the OA variable. Two of the meta-analyses combined studies with knee, hip, and hand data into a single effect size [7,9]. Individual studies, which assessed knee pain and radiographic OA with mortality, tend to report higher effect sizes more consistent without our results. For instance, Liu et al. [15] [11,47], supporting the concept that the association of knee OA with reduced time-to-mortality may be driven by pain rather than by structural changes identified by radiographs. In this analysis, we treated a number of comorbidities as potential confounders; however, the relationship between these comorbidities and OA is poorly understood and may ultimately be part of the causal pathway; therefore, the associations we found here may not represent a causal association between OA and mortality. We could be underestimating this association if some of the potential confounders are actually mediators on the causal pathway, and we may be overestimating the association depending on how well our adjusted models are accounting for confounding. Patients with OA have on average 2·6 moderate-tosevere co-morbidities [48] and 31% of patients have five or more other chronic conditions [49]. Our fully adjusted models, which included lifestyle factors and cardiovascular conditions, did not change substantially from the models adjusted for age, sex, and race. This may indicate that the additional potential confounders we have adjusted for do not have a substantial confounding effect on the association between OA and reduced time-to-mortality. This suggests that either the association is driven by OA or is due to residual confounding caused by measurement error in self-reported variables and/or by the lack of potential confounders such as physical activity and occupation. An additional potential source of unmeasured confounding is the pain sensitization, which may effect the relation of painful knee OA and mortality, and should be pursued in future research. The current study focuses on the knee; however, it is known that limitations in activities of daily living and mobility vary according to hip or knee site [50]; previous studies have also found an increased risk of mortality in individuals with hip symptoms [51].

Strengths and limitations
A limitation of this study is that the included cohorts were designed as independent studies and were not originally Fig. 2 a-c Forest plots of fully adjusted models: a ROA; b POA; c PROA compared to Pain-/ROA-designed to be directly compared to one another. Therefore, osteoarthritis was assessed differently between cohorts. It is known that even small variations in the way a pain question is worded, or X-rays are graded, can result in differences in OA prevalence [22,52]. To minimise this variation, we made every effort to harmonise pain and ROA variables between cohorts by conducting an international expert consensus study [22].
One of the strengths of our study, unlike traditional metaanalyses, is that we actively sought cohorts that had not previously published on the association between OA and mortality, to avoid publication bias. To also capture people without the symptomatic aspects of OA, we restricted our studies to those that included the general population and one enhanced risk factor cohort. These people would not be included in clinical OA cohorts, which is a known issue in the accurate reporting of the true burden of OA [53,54].
The MOST cohort included additional focussed recruitment to include a larger proportion of participants that were older, female, overweight, or had knee surgery/injury, all factors associated with an increased risk of OA. Therefore, the reference group (without pain or ROA) is likely to have a higher prevalence of OA risk factors than the pain and ROA-free group in other cohorts, which may have biased our results toward the null in this cohort.
The follow-up of the included studies ranged from 5.6 to 20.0 years; however, only baseline knee OA and confounders were included in the analysis, meaning that participants may have changed OA categories after the baseline visit, resulting in possible misclassification bias. A further potential limitation is that the age of our participants at baseline ranged between 45 and 80; however, the mean age between cohorts was relatively similar with lowest having a mean age of 56.0 and the highest with a mean age of 64.4 (Table 1).
Both a strength and limitation of the current study is that we included cohorts from different countries, with different races, cultures and health care systems. Confounders were harmonised using the least detailed information available in any single cohort at the baseline visit only, which likely increased our risk of residual confounding in our models. However, by harmonising the individual confounders and adjusting for them consistently between studies, we have reduced unnecessary heterogeneity between studies. Therefore, remaining differences between cohorts are more likely to reflect racial, country, and/or cultural variations rather than how variables were defined, or which statistical models were used.
Previous individual cohort or meta-analysis studies have suggested that a large proportion of the increased risk of mortality is due to cardiovascular mortality [9,14]. Cause-specific mortality was not available in the majority of our cohorts and justifies further investigation. Likewise, medical detail was not available across all cohorts to consider the effect of painrelieving medications. Pain medication is on the casual pathway between painful OA and mortality, and by not including it in our model, our associations are ultimately combining both the direct effect of OA on mortality and the indirect effect of OA through pain medication on mortality. Future research using mediation analysis will help to clarify this pathway.
IPD meta-analyses are time-consuming and resource intensive compared with traditional meta-analyses; however, they allow for standardising exposures, outcomes, and statistical methods, and, more importantly, avoid publication bias by not being limited to the inclusion of previously published studies, which is rarely done in the traditional meta-analysis.

Conclusion
This study is the first individual participant-level data metaanalysis of knee osteoarthritis and premature mortality. It demonstrates that participants with knee pain only or a combination of knee pain and radiographic OA had a 35-37% increased association with reduced time-to all-cause-mortality independent of age, sex, race, BMI, smoking, alcohol, CVD, or diabetes.
With the increasing prevalence of knee OA, it is essential that clinicians and public health bodies are aware of the potential that people with OA may have an increased burden of premature mortality compared to people without OA. This finding highlights that osteoarthritis is a serious disease and supports the need for further research to identify whether OA-related mechanisms are causally associated with premature mortality.    Funding This study was funded by the Centre for Sport Exercise and Osteoarthritis Research Versus Arthritis (grant 21595).