Contents lists available at ScienceDirect

International Journal of Medical Informatics

journal homepage: www.elsevier.com/locate/ijmedinf

Development and evaluation of an osteoarthritis risk model for integrationinto primary care health information technologyJason E. Blacka,*, Amanda L. Terryb, Daniel J. Lizotteca Graduate Program in Epidemiology & Biostatistics, Western University, 1151 Richmond Street, London, Ontario, N6A 5C1, Canadab Department of Family Medicine, Department of Epidemiology & Biostatistics, Schulich Interfaculty Program in Public Health, Western University, 1151 Richmond Street,London, Ontario, N6A 3K7, Canadac Department of Computer Science, Department of Epidemiology & Biostatistics, Schulich Interfaculty Program in Public Health, Department of Statistical and ActuarialSciences, 1151 Richmond Street, Western University, London, Ontario, N6A 3K7, Canada

A R T I C L E I N F O

Keywords:Prognostic prediction modelRiskElectronic medical recordOsteoarthritisCPCSSN

A B S T R A C T

Background: We developed and evaluated a prognostic prediction model that estimates osteoarthritis risk for useby patients and practitioners that is designed to be appropriate for integration into primary care health in-formation technology systems. Osteoarthritis, a joint disorder characterized by pain and stiffness, causes sig-nificant morbidity among older Canadians. Because our prognostic prediction model for osteoarthritis risk usesdata that are readily available in primary care settings, it supports targeting of interventions delivered as part ofclinical practice that are aimed at risk reduction.Methods: We used the CPCSSN (Canadian Primary Sentinel Surveillance Network) database, which containsaggregated electronic health information from a cohort of primary care practices, to develop and evaluate aprognostic prediction model to estimate 5-year osteoarthritis risk, addressing contextual challenges of dataavailability and missingness. We constructed a retrospective cohort of 383,117 eligible primary care patientswho were included in the cohort if they had an encounter with their primary care practitioner between 1January 2009 and 31 December 2010. Patients were excluded if they had a diagnosis of osteoarthritis prior totheir first visit in this time period. Incident cases of osteoarthritis were observed. The model was constructed topredict incident osteoarthritis based on age, sex, BMI, previous leg injury, and osteoporosis. Evaluation of themodel used internal 10-fold cross-validation; we argue that internal validation is particularly appropriate for amodel that is to be integrated into the same context from which the data were derived.Results: The resulting prediction model for 5-year risk of osteoarthritis diagnosis demonstrated state-of-the-artdiscrimination (estimated AUROC 0.84) and good calibration (assessed visually.) The model relies only on in-formation that is readily available in Canadian primary care settings, and hence is appropriate for integrationinto Canadian primary care health information technology.Conclusions: If the contextual challenges arising when using primary care electronic medical record data areappropriately addressed, highly discriminative models for osteoarthritis risk may be constructed using only datacommonly available in primary care. Because the models are constructed from data in the same setting where themodel is to be applied, internal validation provides strong evidence that the resulting model will perform well inits intended application.

1. Introduction

Prognostic prediction models (PPMs) estimate a patient’s risk ofdisease development [1,2] based on various predictors [3,4]. Predictorsmay include patient demographics (such as age and sex), family history,

lifestyle factors (such as smoking status or physical activity level), priormedical conditions, laboratory test results, radiographic imaging, orgenetic markers [5]. In turn, health care practitioners and patients canmake decisions informed by disease risk [6,7]. For example, a patientfound to be at high risk of lung cancer may be advised by their health

https://doi.org/10.1016/j.ijmedinf.2020.104160Received 30 December 2019; Received in revised form 28 February 2020; Accepted 24 April 2020

Abbreviations: PPM, prognostic prediction model; EMR, electronic medical record; CPCSSN, Canadian Primary Care Sentinel Surveillance Network; BMI, body massindex; ICD-9, International Classification of Disease, Ninth Revision; AUROC, Area Under the receiver operating characteristic Curve; TOARP, Tool for OsteoarthritisRisk Prediction; MOST, Multicenter Osteoarthritis Study

⁎ Corresponding author at: Graduate Program in Epidemiology & Biostatistics, Western University, 1151 Richmond Street, London, Ontario, N6A 5C1, Canada.E-mail addresses: [email protected] (J.E. Black), [email protected] (A.L. Terry), [email protected] (D.J. Lizotte).

International Journal of Medical Informatics 141 (2020) 104160

Available online 01 May 20201386-5056/ © 2020 Elsevier B.V. All rights reserved.

T

care provider to quit smoking. Several PPMs estimate a patient’s risk ofdeveloping osteoarthritis [8–11]; however, all existing PPMs for os-teoarthritis that we identified require information on predictors that arenot routinely collected in primary care, such as the Kellgren andLawrence grade (which requires radiographic imaging), and hence arenot suitable for integration into primary care health information sys-tems. Using existing data within primary care electronic medical re-cords (EMRs) instead would eliminate the need for collection of addi-tional, oftentimes burdensome, measures and would enable real-timerisk estimation at the point of care.

There is the potential for significant benefit to be derived by de-ploying such a risk engine in primary care. Affecting an estimated 13 %of adults over the age of 20, osteoarthritis causes significant morbidityin Canada [12]. This estimate increases to 29 % in adults 70 years ofage and older who receive primary care [13]. Symptoms of osteoar-thritis include joint pain and stiffness [14], commonly affecting thejoints of the hands, neck, lower back, hips, and knees. Osteoarthritistreatment largely consists of symptom management (e.g., non-steroidalanti-inflammatory medications for pain management), rather thantreatment of underlying disease mechanisms [15]. Total joint replace-ment is often required after significant degradation of the affected joint.To mitigate this burden, prevention strategies have shown potential inreducing the incidence of osteoarthritis. For example, a diet and ex-ercise program aimed at weight loss reduced the incidence of osteoar-thritis, though not statistically significantly [16]. Injury preventionprograms have been suggested as a potential strategy to prevent os-teoarthritis [17]. Interventions such as these may be improved by se-lectively targeting those at the greatest risk of osteoarthritis in order toreduce their risk. To perform this selective targeting, individualized riskestimates for osteoarthritis are required.

Designing and evaluating a PPM specifically in the primary carecontext leads to two important design decisions: 1) the PPM should useonly data that are easily available in the primary care context where itis to be deployed, for example those that exist in EMRs already, and 2)evaluation of the PPM should reflect the population where it is to bedeployed, that is, in primary care encounters. Researchers have begunto recognize the value of EMR data for research purposes more gen-erally [18]; however, data quality within these databases remains un-certain [19]. Issues such as implausible data and missing data arecommon in EMR data. When working with EMR data, researchers mustaddress these contextual challenges [20].

In this work, we developed and validated a PPM to estimate a pa-tient’s five-year risk of osteoarthritis development using primary careEMR data. Ultimately, we see this model being developed into a pur-pose-built tool to be used routinely by primary care practitioners duringpatient encounters to: 1) deliver a quantitative assessment of osteoar-thritis risk in patients where the patient and/or primary care practi-tioner is concerned about osteoarthritis risk; and 2) act as a passive riskscreening tool to identify high-risk patients who may have gone un-detected otherwise. We are confident that our learnings from this workwill be of use to others who are designing and evaluating PPMs usingEMR from and for primary care.

2. Methods

We developed and validated a prognostic prediction model to esti-mate the risk of osteoarthritis development within five years amongCanadian adults receiving primary care. Model development was in-formed by strategies of PPM development suggested by the TRIPODstatement [21], Steyerberg [22], Lee et al. [5], and Hendriksen et al.[3]. First, we compiled a list of risk indicators for osteoarthritis devel-opment based on the existing literature. Next, we identified a cohort ofpatients whose risk indicator status was known at baseline and assessedwhether they subsequently developed osteoarthritis within five years.

Based on this cohort, a multivariable model was constructed to enablethe estimation of individual patient risk. This model’s performance wasassessed in terms of its discrimination and calibration.

We began by examining existing literature to identify establishedpredictors of osteoarthritis development. We identified: BMI (bodymass index) [23–29], previous leg injury [23,25–27], leg length in-equality [30], older age [24–26,28,29], female sex [24–26,29], osteo-porosis [24], family history [29], occupation [29], and physicalworkload [28].

We used the CPCSSN (Canadian Primary Care Sentinel SurveillanceNetwork) database to develop our PPM. The CPCSSN database containsde-identified patient records from 12 regional networks across Canadaand includes more than 1.5 million patients [31]. Within this databaseare all structured patient records stretching back to 2008, includingpatient encounters, patient demographics, billing codes, laboratoryresults, prescriptions, referrals, risk indicator information, and medicalprocedures. The CPCSSN database contains only structured data; un-structured data (e.g., free-text clinician notes) are not available. Wefound CPCSSN contained data describing five of the nine risk indicatorsfor osteoarthritis: BMI, previous leg injury, older age, female sex, andosteoporosis.

2.1. Measures of risk indicators and outcomes

As is common when working with secondary data sources such asEMRs or health administrative data, risk indicators and outcomes maynot be found directly in EMR databases; we therefore developed astrategy to mitigate this contextual challenge. In consultation with ex-pert EMR users, we searched the EMR for data elements that we thoughtwere strongly associated with identified risk indicators. For example,we considered a patient with any of the following diagnostic codes tohave had a lower leg injury: ICD-9 (International Classification ofDiseases – Ninth Revision) 820-29 (fracture of lower limb), ICD-9 843(sprain or strain of hip and thigh), ICD-9 844 (sprain or strain of kneeand leg), or ICD-9 928 (crushing injury to lower limb). We compiledthese data elements into risk indicator definitions (Table 1).

We defined the risk indicator osteoporosis using evidence in the EMRfor any suspected bone disorder, despite this not referring specifically toosteoporosis. Our definition includes use of the ICD-9 code 733 (used tonote osteoporosis and other bone disorders) in the problem list, billingdata, or encounter diagnosis fields; the term “osteoporosis”; or pre-scription of a medication commonly used for the treatment or preven-tion of osteoporosis. Specific diagnoses of osteoarthritis (ICD-9 code:733.0) were unavailable.

Table 1Risk indicator definitions for CPCSSN database.

Risk Indicator Data Source Value

Age Patient Demographics NumericSex Patient Demographics Male or femaleBMI Patient Encounter NumericLeg Injury Billing ICD-9 Codes:

Health ConditionEncounter Diagnosis

• 820−29: fracture of lower limb• 843: sprain or strain of hip and thigh• 844: sprain or strain of knee and leg• 928: crushing injury to lower limbOsteoporosis Billing ICD-9 Code:

Health ConditionEncounter Diagnosis

• 733: Osteoporosis and other bonedisorders

Health ConditionEncounter DiagnosisRisk Factor

“osteoporosis”

Medications Alendronic acidRisedronic acidIbandronic acid

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

2

In an effort to increase the usability of CPCSSN data, Williamsonet al. previously constructed a case definition for osteoarthritis byconsulting with several expert EMR users [32]; this case definition wasvalidated by estimating the sensitivity and specificity of this definitioncompared to chart review by an expert EMR user as a “gold standard.”The case definition for osteoarthritis consisted of a combination ofbilling codes and problem list diagnoses (Table 2).

2.2. Cohort construction

We included any adult patient (known to be age 18 or older) whohad at least one interaction with a primary care practitioner between 1January 2009 and 31 December 2010 and had not been previouslydiagnosed with osteoarthritis. The data used for this work spanned2008–2016. The selection of interactions within this window allowedfor a wide window with sufficient time for follow-up. Risk indicatorswere assessed on the interaction within this window by examining re-cords up to and including the date of this interaction. For example, if apatient’s records indicated any diagnosis of osteoporosis prior to thisinteraction, they were considered to have the risk indicator osteo-porosis. Patients diagnosed with osteoarthritis in the five years fol-lowing this interaction date were considered cases of osteoarthritis.Thus, each patient had their own start date and follow-up period.Where there were multiple interactions between 1 January 2009 and 31December 2010, the earliest one was used.

2.3. Model building and evaluation

All available data were used in the development and validation ofthe model to maximize the predictive ability of the resulting model.Ten-fold cross validation was used to evaluate the model: all data wereused to create ten partitions; one partition was reserved for validation,the remaining nine were used to estimate the model parameters. Thisstep was repeated ten times, such that each partition was used for va-lidation once. Final estimates of validation measures were obtained byaveraging across each of the cross-validation models. For comparisonpurposes and to further investigate the question of evaluation in thiscontext, we also present abbreviated results of a single-split validation.

Missing data are a common contextual challenge of primary careEMR data analysis; we used a combination of complete case analysisand imputation to deal with missing data issues as follows. Some pa-tients’ birth years were zero; we treated these data as missing. Weeliminated any patient whose age was missing (n = 259), as we in-tended to estimate risk in confirmed adults. We used multiple im-putation to address missing data for BMI and sex [33,34]. For eachcross-validation step, separate imputation models were constructed forthe development and validation sets, as recommended by Wahl et al.[35]. We imputed five datasets for each partition using the MICEpackage in R [34]. In addition to the identified risk factors for os-teoarthritis, we included several covariates in the imputation model,including diagnoses of several chronic conditions and rurality to max-imize the accuracy of the imputations (Appendix Table A1). Diagnosisof osteoarthritis was also included in the imputation models, as

suggested by Moons et al. [36]. We then assessed the plausibility of theimputations by plotting the distribution of each variable with missingdata to ensure that the imputed data followed a similar distribution asthe original data. While the imputation process does smooth the cov-ariate distribution somewhat, the general location and scale of theimputed and observed data are similar.

Logistic regression was used to construct a prediction model for thedevelopment of osteoarthritis using the entire cohort; this model wasproduced by combining models from the five imputed datasets usingRubin’s rules. All risk indicators were included in the model in theiroriginal form (e.g., all continuous risk indicators were included ascontinuous variables; no binning was performed). Investigations per-formed on a similar dataset did not reveal non-linear associations be-tween risk indicators and the outcome: log-transformation of con-tinuous variables; the addition of polynomial transformations ofcontinuous variables; and the use of generalized additive models allresulted in models that performed similarly to logistic regression. Thesimpler logistic regression model was preferred.

Model discrimination and calibration were evaluated using cross-validation. As discussed, AUROC was used to assess discrimination; aprecision-recall curve was also used to investigate the precision andrecall of the model across various thresholds. Precision (or positivepredictive value) is the proportion of true positives amongst all pre-dicted positives. Recall (or sensitivity) is the proportion of predictedpositives amongst all true positives. The precision recall curve allows usto examine both precision and recall without setting an arbitrary riskthreshold [37]. To evaluate calibration, a calibration plot was usedbecause measures such as the Hosmer-Lemeshow goodness-of-fit testhave been shown to be oversensitive in large sample sizes [5]. A cali-bration plot displays the average predicted risk within each risk decileagainst the observed proportion of patients who develop osteoarthritiswithin that decile (i.e., observed risk). The calibration plot of a modelthat displays good calibration should closely follow a line with an in-tercept of zero and a slope of one, demonstrating strong agreementbetween the estimated and observed risks. Additionally, we presentcalibration in the large and calibration slope to further assess calibra-tion.

Other modelling techniques were considered, including Cox pro-portional hazards; however, our goal is to predict and convey risk ra-ther than to establish the relationship between outcomes and covari-ates, and using a Cox model for this purpose would require theadditional step of estimating baseline risk. Furthermore, although asurvival analysis method could take into account censoring, the EMRdata lack censoring information in the way it would be conceived in atraditional longitudinal cohort study because it is not clear when aprimary care patient is “lost to follow-up”; we illustrate this issue usinga sensitivity analysis that demonstrates that any such censoring wouldhave limited impact on our model. Therefore, we chose to forego sur-vival analysis for this work and opt for the simpler logistic regressionmodel.

3. Results

The final cohort was composed of 383,117 patients (Fig. 1). Patientcharacteristics were typical of a primary care population, as they wereslightly older and more likely to be female [38–40] (Table 3). After fiveyears of follow-up, 12,803 (3.3 %) patients developed osteoarthritis.

Data were commonly missing for BMI, while sex was almost nevermissing (Table 4). Multiple imputation was used to address missingdata for BMI and sex.

Kernel density estimates of the distribution over the five imputeddatasets (Fig. 2) demonstrate that the distributions of the imputed va-lues are similar to the distributions of the original (unimputed) values.

Table 2Validated case definition for osteoarthritis.

Billing Problem List

Any occurrence of the following codes: Any occurrence of the following codes:

• 715, Osteoarthritis and allieddisorders

• 721, Spondylosis and allieddisorders

• 715, Osteoarthritis and allieddisorders

• 721, Spondylosis and allied disorders

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

3

The final model is given in Table 5. All risk indicators were sig-nificantly associated with incident osteoarthritis diagnosis.

Internal validation of the model based on cross-validation yielded ahigh AUROC (0.84, 95 % bootstrapped CI 0.83 to 0.85). Calibration inthe large (−0.00002) and the calibration slope (0.999) indicated ex-cellent calibration, as would be expected due to our use of validationdata derived from the same population as for model development.Visual inspection of the calibration plot shows that predicted risk clo-sely followed the observed risk (Fig. 3). Risk was slightly under-estimated among lower-risk patients, while the risk of moderate tohigh-risk patients was slightly overestimated. We also performed asingle-split evaluation of our model for comparison purposes using238,875 for training and 144,242 for validation to ensure a preciseestimate of AUROC. This resulted in an identical estimated AUROC(0.84) with a 95 % bootstrapped CI of 0.83 to 0.84.

As seen in the precision-recall curve shown in Fig. 4, the model hasmoderate precision-recall characteristics but is able to identify thosewith higher risk.

4. Discussion

We produced a prognostic prediction model for the diagnosis ofosteoarthritis using EMR data that are commonly available in primarycare. We see this model ultimately being used in two ways in Canadianprimary care settings. First, the model can provide estimates of os-teoarthritis risk on demand when requested by a provider or patientwho is interested in osteoarthritis risk. Second, the model can operatein the background of the provider’s EMR system during all patient visitsand automatically flag any patients whose osteoarthritis risk surpassessome prespecified threshold, or to order patients from highest-risk tolowest-risk.

Our model compares extremely favourably with existing work.Existing models for osteoarthritis risk estimation include the Tool forOsteoarthritis Risk Prediction (TOARP) [8]; the Nottingham knee os-teoarthritis risk prediction models [9]; and models derived from datafrom the Rotterdam Study-1 [10] and the Multicenter OsteoarthritisStudy (MOST) [11]. These models estimate risk of knee osteoarthritis,whereas our model estimates risk of osteoarthritis in any joint. All ex-isting models were constructed using population-based cohorts rangingin size from 400 to 3000 people who were assessed using interviews,physical examinations, and laboratory tests, including radiographicimaging. In contrast, our model was constructed using existing EMRdata from a primary care population of over 380,000 patients. Allmodels, including ours, used multivariable logistic regression to con-struct the prediction model. Internal validation of the existing modelsreported AUROC ranging from 0.70 to 0.79. Our model had the highestdiscriminative ability (0.84) and hence we consider it state-of-the-artfor its purpose of a primary care tool.

Our model has several limitations. First, not all risk indicators forosteoarthritis were available within the CPCSSN database: data de-scribing family history, occupation, leg length inequality, and physicalworkload were not available. We were unable to link to additionaldatabases to obtain these data. Thus, estimated risk for those whopossess these uncaptured risk indicators will be an underestimation oftheir true risk. Second, there are elements of a patient’s treatment thatwe may not be able to observe. As such, we cannot adjust for the impactthis treatment potentially has on patient risk. Third, the model weproduced estimates a patient’s risk of osteoarthritis without specifyingthe joint affected. However, interventions to reduce osteoarthritis riskare typically not joint specific [16], thus knowing which joint is likelyto develop osteoarthritis does not inform prevention.

The precision of our model was somewhat limited for larger recallthresholds (Fig. 4). However, this is not of great concern for this par-ticular application because the cost associated with the treatment offalse positives in this context is minimal; most treatments for osteoar-thritis consist of lifestyle modifications with little to no risk of harm[16]. Thus, the lower precision of our model is not alarming. Ad-ditionally, the burden on primary care practitioners and patients usingour model is low; no additional measures beyond those routinely col-lected in the EMR are required. The risk information gained from ourmodel comes at minimal cost to the practitioner and may be used toinform further targeted screening that gathers additional tailored in-formation.

Our model did not account for censoring of patients. This followedfrom the assumption that if a primary care patient did not seek carefrom their primary care practitioner and receive an osteoarthritis di-agnosis, they did not develop osteoarthritis. Given the population (pa-tients who visit their doctors) and given that osteoarthritis diagnosis is

Fig. 1. Patient flowchart.

Table 3Descriptive statistics for the dataset.

CPCSSN(n = 383,117)

Age, median (IQR) 44 (30–59)BMI, median (IQR) 26.6 (23.3–30.7)Female, n (%) 221,021 (57.6 %)Prior leg injury, n (%) 10,893 (2.8 %)Prior diagnosis of osteoporosis, n (%) 11,647 (3.0 %)

Table 4Predictors with missing data.

Predictor CPCSSN(n = 383,117)

BMI, n (%) 256,413 (66.9)Sex, n (%) 59 (0.02)

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

4

often instigated by patients reporting their symptoms, we argue thatthis is a reasonable assumption. We nonetheless recognize that patientsmay have been censored through “loss to follow-up” (perhaps becausethey changed providers) or died during the study; however, unlike insome prospective studies, we are not able to directly observe this cen-soring. To investigate the impact of censoring, we conducted a sensi-tivity analysis whereby we excluded patients who did not have at leastone interaction with their primary care practitioner after the end oftheir follow-up window. This ensured that patients in the remainingcohort were unlikely to have been lost to follow-up. This sensitivityanalysis revealed that the potential impact of censoring on the modelwas minimal: model estimates based on this restricted cohort were si-milar to those of the original cohort (Appendix Table A2); model per-formance remained strong as well (AUROC: 0.83, 95 % CI 0.82 to 0.84;see Appendix Fig. A1 for calibration plot).

One consideration that may or may not be a limitation depending onone’s goal is the potential lack of generalizability of our model to thegeneral population, since it was derived from primary care data. If thegoal of the model was to be deployed in a population-level public healthcampaign aimed at reducing osteoarthritis incidence, for example, itmay not be an appropriate tool. However, the use of primary care EMRdata positions our model strongly for deployment in the primary caresetting. The data used for risk estimation are those already currentlycollected within EMR systems; no additional measures need to be col-lected in order to use our prediction model. This aspect is unique to ourmodel: all other existing models require some radiographic imagingthat is not routinely collected. This consideration also supports the useof internal validation as an approach for evaluation; given that the

model is anticipated to be applied to the same population from which itwas derived, internal validation is an appropriate way to evaluateperformance. Specifically, our model is most appropriate for use in the

Fig. 2. Kernel density estimates for the marginal distribution of the five imputed datasets (red, shorter peaks) and the original data (blue, taller peaks) (left:development sets; right: validation sets) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).

Table 5Model estimates.

Reference category/units β coefficient (95% CI) Odds ratio (95 % CI) P value

(Intercept) −8.13 (−8.23 to −8.02) — <0.001Age Years 0.059 (0.058 to 0.060) 1.061 (1.059–1.062) <0.001BMI kg/m2 0.042 (0.039 to 0.044) 1.043 (1.040–1.045) <0.001Sex Male (Reference) (Reference)

Female 0.21 (0.17 to 0.25) 1.24 (1.19–1.29) <0.001Prior leg injury No (Reference) (Reference)

Yes 1.61 (1.54–1.67) 5.00 (4.68–5.34) <0.001Prior diagnosis of osteoporosis No (Reference) (Reference)

Yes 0.92 (0.86 to 0.99) 2.52 (2.37–2.68) <0.001

CI: confidence interval; BMI: body mass index.

Fig. 3. Calibration plot based on validation sets.

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

5

Canadian primary care setting, as data were derived from Canadianprimary care EMRs. Use in new regions would ideally be preceded bymodel ‘updating’ using data from the new region [41].

A related issue arises given that we have used what might be termed“surrogate” measures of known risk factors. For example, consider thecollection of diagnosis codes we have used to indicate a leg injury.Because we do not have gold standard information for leg injuries, wedo not know the accuracy of this surrogate; however, we have shownthat the surrogate is predictive of osteoarthritis in its own right, and solong as it is available at the time of prediction (which it would be if themodel were deployed in a similar primary care context as that whichgenerated the data) then it is a valid and useful risk indicator.

It is important to note that our osteoarthritis prediction modelquantifies patient risk to guide decision making and identify high riskpatients. As with all prediction models, it would be inappropriate toinfer causality between risk indicators and osteoarthritis based on ourmodel. The recommendation of any intervention to reduce risk of os-teoarthritis should be based on evidence demonstrating the effective-ness of that intervention in reducing osteoarthritis risk.

5. Conclusions

Primary care EMRs are a rich, yet underutilized, source of long-itudinal health data that can support the development of novel tools forintegration into primary care health information systems. Our workdemonstrates the utility of these data for constructing PPMs despitecontextual challenges such as missing data, using an osteoarthritis riskmodel as a success story, and provides a strategy and rationale for in-ternal validation. Two key future directions for this work will be to 1)design and evaluate strategies for incorporating the model into clinical

workflows as appropriate and 2) to consider and evaluate more com-plex PPMs both for osteoarthritis and other diseases, for example thosederived using machine learning techniques; this second direction willrequire a careful trade-off between potential performance improve-ments and interpretability, and will thus go hand-in-hand with activ-ities that evaluate how best to enhance primary care health informationtechnology.

Summary Points

– Developing prognostic prediction models (PPMs) presentscontextual challenges that depend simultaneously on dataprovenance and on the target population.

– We present a new, state-of-the-art PPM for osteoarthritis de-rived from primary care data that can integrate into primarycare health information technology.

– We describe how using primary care EMR data and deployingto primary care influences model design and evaluationchoices.

Author statement

Lead author was JB. All authors initiated the research idea. JBdrafted the research idea. JB wrote the article that is being submitted.AT and DL revised the article. JB extracted and analyzed the CPCSSNdata. DL and AT supported the methodology development and con-tributed to the writing of the article. Each author has read and approvedthe final version of this article.

Ethics

Ethics approval was obtained from the Western University ResearchEthics Board #107572.

Consent for publication

Not applicable as no individual person’s data were presented.

Availability of data and materials

The Canadian Primary Care Sentinel Surveillance Network(CPCSSN) database is not publicly accessible, in keeping with the intentof the agreement made with primary health care practitioners con-tributing to the CPCSSN database. Researchers can request data fromCPCSSN directly.

Funding

Funding for this research was provided by the Natural Sciences andEngineering Research Council of Canada. The funding source played norole in this research.

Declaration of Competing Interest

The authors have no conflicts of interest to declare.

Fig. 4. Precision-recall curve based on validation sets.

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

6

Appendix A

Table A1Covariates included in imputation model.

Covariates for imputation model

• Previous diagnosis of:o Diabeteso Hypertensiono Depressiono Alcohol use disordero Epilepsyo Schizophreniao Anxiety disordero Cancero Cardiovascular diseaseo Chronic obstructive pulmonary diseaseo Rheumatoid arthritiso Chronic kidney disease

• Rurality• Income

Table A2Model estimates based on sensitivity analysis.

Reference category/units β coefficient (95% CI) Odds ratio (95 % CI) P value

(Intercept) −8.22 (−8.35 to −8.10) — <0.001Age Years 0.064 (0.062 to 0.065) 1.066 (1.064–1.067) <0.001BMI kg/m2 0.040 (0.037 to 0.043) 1.041 (1.038–1.043) <0.001Sex Male (Reference) (Reference)

Female 0.25 (0.20 to 0.29) 1.28 (1.23–1.34) <0.001Prior leg injury No (Reference) (Reference)

Yes 1.59 (1.52–1.66) 4.90 (4.56–5.28) <0.001Prior diagnosis of osteoporosis No (Reference) (Reference)

Yes 0.82 (0.75 to 0.89) 2.27 (2.12–2.44) <0.001

CI: confidence interval; BMI: body mass index.

Fig. A1. Calibration plot based on sensitivity analysis.

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

7

Appendix B. Supplementary data

Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.ijmedinf.2020.104160.

References

[1] Public Health Agency of Canada, Canadian Chronic Disease Indicators, Quick Stats,(2017) Edition. Ottawa, ON; 2017.

[2] R. Birtwhistle, R. Morkem, G. Peat, T. Williamson, M.E. Green, S. Khan, et al.,Prevalence and management of osteoarthritis in primary care: an epidemiologiccohort study from the Canadian Primary Care Sentinel Surveillance Network, COpen 3 (3) (2015) E270–5.

[3] M. Doherty, J.W.J. Bijlsma, N. Arden, D. Hunter, N. Dalbeth, Osteoarthritis andCrystal Arthropathy, 3rd ed., Oxford University Press, 2016 507 p..

[4] M.C. Hochberg, R.D. Altman, K.T. April, M. Benkhalti, G. Guyatt, J. Mcgowan,et al., American College of Rheumatology 2012 recommendations for the use ofnonpharmacologic and pharmacologic therapies in osteoarthritis of the hand, hip,and knee, Arthritis Care Res. (Hoboken) (2012).

[5] J. Runhaar, M. van Middelkoop, M. Reijman, S. Willemsen, E.H. Oei,D. Vroegindeweij, et al., Prevention of knee osteoarthritis in overweight females:the first preventive randomized controlled trial in osteoarthritis, Am. J. Med. 128(8) (2015) 888–895.e4.

[6] C.A. Emery, E.M. Roos, E. Verhagen, C.F. Finch, K.L. Bennell, B. Story, et al., OARSIClinical Trials Recommendations: design and conduct of clinical trials for primaryprevention of osteoarthritis by joint injury prevention in sport and recreation,Osteoarthr. Cartil. 23 (5) (2015) 815–825.

[7] D.M. Lloyd-Jones, P.W.F. Wilson, M.G. Larson, A. Beiser, E.P. Leip, R.B. D’Agostino,et al., Framingham risk score and prediction of lifetime risk for coronary heartdisease, Am. J. Cardiol. 94 (1) (2004) 20–24.

[8] S.A.M. Nashef, F. Roques, P. Michel, E. Gauducheau, S. Lemeshow, R. Salamon,European system for cardiac operative risk evaluation (EuroSCORE), Eur. J. Cardio-Thorac. Surg. 16 (1) (1999) 9–13.

[9] J.M.T. Hendriksen, G.J. Geersing, K.G.M. Moons, J.A.H. de Groot, Diagnostic andprognostic prediction models, J. Thromb. Haemost. 11 (Suppl 1) (2013) 129–141.

[10] E.W. Steyerberg, Y. Vergouwe, Towards better clinical prediction models: sevensteps for development and an ABCD for validation, Eur. Heart J. 35 (29) (2014)1925–1931.

[11] Y.H. Lee, H. Bang, D.J. Kim, How to establish clinical prediction models,Endocrinol. Metab. (Seoul, Korea). 31 (1) (2016) 38–44.

[12] D.T. Felson, Y. Zhang, J.M. Anthony, A. Naimark, J.J. Anderson, Weight loss re-duces the risk for symptomatic knee osteoarthritis in women. The Framinghamstudy, Ann. Intern. Med. 116 (7) (1992) 535–539.

[13] D.T. Felson, Weight and osteoarthritis, Am. J. Clin. Nutr. 63 (3 Suppl) (1996)430S–432S.

[14] G.B. Joseph, C.E. McCulloch, M.C. Nevitt, J. Neumann, A.S. Gersing,M. Kretzschmar, et al., Tool for osteoarthritis risk prediction (TOARP) over 8 yearsusing baseline clinical data, X-ray, and MRI: data from the osteoarthritis initiative,J. Magn. Reson. Imaging 47 (6) (2018) 1517–1526.

[15] W. Zhang, D.F. McWilliams, S.L. Ingham, S.A. Doherty, S. Muthuri, K.R. Muir, et al.,Nottingham knee osteoarthritis risk prediction models, Ann. Rheum. Dis. 70 (9)(2011) 1599–1604.

[16] H.J.M. Kerkhof, S.M.A. Bierma-Zeinstra, N.K. Arden, S. Metrustry, M. Castano-Betancourt, D.J. Hart, et al., Prediction model for knee osteoarthritis incidence,including clinical, genetic and biochemical risk factors, Ann. Rheum. Dis. 73 (12)(2014) 2116–2121.

[17] D.L. Riddle, P.W. Stratford, R.A. Perera, The incident tibiofemoral osteoarthritiswith rapid progression phenotype: development and validation of a prognosticprediction rule, Osteoarthr. Cartil. 24 (12) (2016) 2100–2107.

[18] H. Carr, S. de Lusignan, H. Liyanage, S.-T. Liaw, A. Terry, I. Rafi, Defining di-mensions of research readiness: a conceptual model for primary care researchnetworks, BMC Fam. Pract. 15 (2014) 169.

[19] S. de Lusignan, S.-T. Liaw, P. Krause, V. Curcin, M.T. Vicente, G. Michalakidis, et al.,Key concepts to assess the readiness of data for international research: data quality,lineage and provenance, extraction and processing errors, traceability, and cura-tion. Contribution of the IMIA Primary Health Care Informatics Working Group,Yearb. Med. Inform. 6 (2011) 112–120.

[20] A.L. Terry, M. Stewart, S. Cejic, J.N. Marshall, S. de Lusignan, B.M. Chesworth,

et al., A basic model for assessing primary health care electronic medical recorddata quality, BMC Med. Inform. Decis. Mak. 19 (1) (2019) 30.

[21] K.G.M. Moons, D.G. Altman, J.B. Reitsma, J.P.A. Ioannidis, P. Macaskill,E.W. Steyerberg, et al., Transparent reporting of a multivariable prediction modelfor individual prognosis or diagnosis (TRIPOD): explanation and elaboration, Ann.Intern. Med. 162 (1) (2015) W1.

[22] E.W. Steyerberg, Clinical Prediction Models, Springer New York, New York, NY,2009 (Statistics for Biology and Health).

[23] C. Cooper, H. Inskip, P. Croft, L. Campbell, G. Smith, M. McLaren, et al., Individualrisk factors for hip osteoarthritis: obesity, hip injury and physical activity, Am. J.Epidemiol. 147 (6) (1998) 516–522.

[24] K.M. Lee, C.Y. Chung, K.H. Sung, S.Y. Lee, S.H. Won, T.G. Kim, et al., Risk factorsfor osteoarthritis and contributing factors to current arthritic pain in South Koreanolder adults, Yonsei Med. J. 56 (1) (2015) 124.

[25] T. Neogi, Y. Zhang, Osteoarthritis prevention, Curr. Opin. Rheumatol. 23 (2) (2011)185–191.

[26] V. Silverwood, M. Blagojevic-Bucknall, C. Jinks, J.L. Jordan, J. Protheroe,K.P. Jordan, Current evidence on risk factors for knee osteoarthritis in older adults:a systematic review and meta-analysis, Osteoarthr. Cartil. 23 (4) (2014) 507–515.

[27] E. Vignon, J.-P. Valat, M. Rossignol, B. Avouac, S. Rozenberg, P. Thoumie, et al.,Osteoarthritis of the knee and hip and activity: a systematic international reviewand synthesis (OASIS), Joint Bone Spine 73 (4) (2006) 442–455.

[28] I. Vrezas, G. Elsner, U. Bolm-Audorff, N. Abolmaali, A. Seidler, Case-control study ofknee osteoarthritis and lifestyle factors considering their interaction with physicalworkload, Int. Arch. Occup. Environ. Health 83 (3) (2010) 291–300.

[29] G.J. Leung, K.D. Rainsford, W.F. Kean, Osteoarthritis of the hand I: aetiology andpathogenesis, risk factors, investigation and diagnosis, J. Pharm. Pharmacol. 66 (3)(2014) 339–346.

[30] W.F. Harvey, M. Yang, T.D.V. Cooke, N.A. Segal, N. Lane, C.E. Lewis, et al.,Association of leg-length inequality with knee osteoarthritis a cohort study, Ann.Intern. Med. 152 (5) (2010) 287–W92.

[31] R. Birtwhistle, K. Keshavjee, A. Lambert-Lanning, M. Godwin, M. Greiver,D. Manca, et al., Building a pan-Canadian primary care sentinel surveillance net-work: initial development and moving forward, J. Am. Board Fam. Med. 22 (4)(2009) 412–422.

[32] T. Williamson, M.E. Green, R. Birtwhistle, S. Khan, S. Garies, S.T. Wong, et al.,Validating the 8 CPCSSN case definitions for chronic disease surveillance in a pri-mary care database of electronic health records, Ann. Fam. Med. 12 (4) (2014)367–372.

[33] M.J. Azur, E.A. Stuart, C. Frangakis, P.J. Leaf, Multiple imputation by chainedequations: what is it and how does it work? Int. J. Methods Psychiatr. Res. 20 (1)(2011) 40–49.

[34] S. van Buuren, K. Groothuis-Oudshoorn, Mice: multivariate imputation by chainedequations in R, J. Stat. Softw. 45 (3) (2011) 1–67.

[35] S. Wahl, A.-L. Boulesteix, A. Zierer, B. Thorand, M.A. van de Wiel, Assessment ofpredictive performance in incomplete data by combining internal validation andmultiple imputation, BMC Med. Res. Methodol. 16 (1) (2016) 144.

[36] K. Moons, R. Donders, T. Stijnen, F. Harrell, Using the outcome for imputation ofmissing predictor values was preferred, J. Clin. Epidemiol. 59 (2006) 1092–1101.

[37] T. Saito, M. Rehmsmeier, The precision-recall plot is more informative than theROC plot when evaluating binary classifiers on imbalanced datasets, PLoS One 10(3) (2015) e0118432.

[38] C.A. Mustard, P. Kaufert, A. Kozyrskyj, T. Mayer, Sex differences in the use of healthcare services, N. Engl. J. Med. 338 (23) (1998) 1678–1683.

[39] K.D. Bertakis, R. Azari, L.J. Helms, E.J. Callahan, J.A. Robbins, Gender differencesin the utilization of health care services, J. Fam. Pract. 49 (2) (2000) 147–152.

[40] J.X. Nie, L. Wang, C.S. Tracy, R. Moineddin, R.E. Upshur, Health care service uti-lization among the elderly: findings from the Study to Understand the ChronicCondition Experience of the Elderly and the Disabled (SUCCEED project), J. Eval.Clin. Pract. 14 (6) (2008) 1044–1049.

[41] K.J.M. Janssen, Y. Vergouwe, C.J. Kalkman, D.E. Grobbee, K.G.M. Moons, A simplemethod to adjust clinical prediction models to local circumstances, Can. J. Anesth.56 (3) (2009) 194–201.

J.E. Black, et al. International Journal of Medical Informatics 141 (2020) 104160

8