LPA was conducted to identify statistically distinct patterns of service contact within a sample of individuals with a diagnosis of PD receiving treatment in secondary care in South London. In addition to identifying the optimal class solution in this setting, formal validation aimed to inform generalization estimation, reporting both average-level and individual-level classification errors. This paper further extends LPA research through consideration and evaluation of the structural model and extends such evaluation through protocolized resampling.

Validation-focused analyses assessed the structural model [39], with parameters tested for relevance when considering cluster membership. Out-of-sample approaches to prediction are well researched in wider structural equation modeling perspectives, for example, partial least squares structural equation modeling [40]. Importantly, while findings of mixture models may prove interesting, classes may not necessarily reside in the wider population [41], further motivating explicit consideration of generalization in study design.

LPA guidelines suggest using 3 sampling frames as good practice to validate findings [42,43]. To begin, the initial whole sample was “sample frame 1,” comprising the final analytical set. Validation samples randomly split the whole sample into two equal groups of 3941, resulting in “sample frame 2” [42,44] and 2 validation samples (1970/3941, 49.98%). Finally, cross-validation requires a holdout procedure [45], which partitions the dataset into test and training samples, to assess “out-of-sample” performance. This procedure split the whole sample into 5 parts, with 80% or k−4 parts designated as the training sample (3153/3941, 80%) and the remaining 20% or k−1 part designated as the test sample (788/3941, 20%). As shown in Figure 1, there were 3 sampling frames as a result.

In this study, posterior probabilities of cluster membership were identified using sample frames 1 and 2. The comparison between sampling frames (Tables S1 and S3 in Multimedia Appendix 1) describes the formal relevance of the model regarding validation. Posterior probabilities of cluster membership were calculated from the e-step of the expectation maximization algorithm [46] and describe how the person i belongs to the class k [47,48], also described as their modal class assignment, according to the selected model [42].

Patient data for the study were extracted from the Clinical Record Interactive Search platform at the South London and Maudsley NHS Foundation Trust. South London and Maudsley (SLaM) provides comprehensive, near-monopoly specialist mental health services to a geographic catchment area of 4 boroughs of southeast London (Croydon, Lambeth, Lewisham, and Southwark) with a population of approximately 1.3 million residents.

The Clinical Records Interactive Search (CRIS) system was set up in 2007-2008 through codevelopment of a data security and governance model with service users [49]. The CRIS system and related linked data have received successive research ethics approvals since 2008 as an anonymized data resource for secondary analyses (current approval: Oxford Research Ethics Committee C, approval number 23/SC/0257). The institutionally approved governance model has been described in detail elsewhere [50], but in summary, it includes automated deidentification of the record, standard researcher approval, a locally publicized opt-out facility, and an oversight committee with patient leadership and membership reviewing all data use.

The inclusion criteria comprised having an ICD-10 PD diagnosis [41] at any time, being aged between 20 and 101 years at initial presentation, and being in contact with SLaM due to referral, although 10.38% (409/3941) of the sample had no recorded team contacts in the analysis set. Given the population of interest and generalization intention, those aged <20 years (318/3941, 8.07%) and >101 years (8/3941, 2.02%) were excluded to reduce possible inclusion bias and to encourage greater sample representativeness.

The inclusion criteria were applied over an exposure period of 11 years, with care episodes sampled between 2007 and 2018. This study sought to identify a matched, complete set of service use data at service users’ first and final points of contact with services (time 1 and time 2, respectively). Time points were denoted by episode start and end dates, such that time 1 represented the start and end of the first episodes of care and time 2 represented the same at the final episodes of care. Dataset assembly resulted in a complete, matched set of data points, ascribed by time 1 and time 2 and organized by episode of care. An episode of care was calculated as the elapsed time between patients being referred to and in contact with services and subsequent discharge from secondary services. In the final matched pairs set, the average length of an episode of care was 131.86 (SD 309.69) days. Moreover, “discharge” is defined by SLaM Trust policy as “discharging a patient from an in-patient setting to a community team or from the Trust completely to primary care or another provider” [51]. Importantly, this follows when a patient is deemed to no longer benefit from either reduced risk or improved well-being [51] and suggests a shift in presentation and overall needs.

Working with the CRIS team, a custom data extract resulted in 82,265 episodes of care from 8510 service users. The inclusion of contacts from clinical teams data resulted in a reduced set of 6068 individuals, which was further processed to include complete “matched pairs” data for 3941 service users, at their first and final episodes of care (T1 and T2, respectively; Figure S1 in Multimedia Appendix 3). A complete dataset was preferred as there were no missing data, which strengthened the design using a large sample. T1 and T2 contained episode start and end dates, respectively, as discussed in the Statistical Analysis section. The authors acknowledge that the nature of this experimental design necessitated service users engaging with services on more than one occasion. Therefore, individuals with only a single episode of care (2127/3941, 53.97%) were omitted from the analysis. Notwithstanding, PD is typically a chronic condition [52,53], and in general, service users are likely to have ongoing needs and therefore present to services on more than one occasion.

Cross-validation uses the data to test “the predictive power” of the model [80]. For the procedure, cross-validation splits the dataset into k subsamples of equal size, repeating k times, each time using a different test set (Figure 1). The best-fitting parameters of the model are then fit using k−1 subsamples, with the model’s error then computed on the unused subsample or test set [81].

In line with the discriminant analysis procedure, cross-validation was conducted using sample frames 1 and 2 (in sample and out of sample, respectively), as shown in Figure 1. The Brier score [82,83] indicates the predictive usefulness of the model for individuals who, theoretically, may be unseen. Applied research in mental health has shown increased experimentation with the Brier score as a forecasting tool, such as for predicting suicidal ideation in individuals at near-term suicide risk [84] or identifying patients likely to drop out of psychotherapy [85]. In related psychotherapy research, the Brier score ranges from 0 (best prediction) to 1 (worst prediction), and it measures the accuracy of probabilistic predictions [85].

A detailed sampling strategy enabled both explanatory and predictive appraisal of the derived model solution, as shown in Figure 1. As well as the whole sample (frame 1), validation samples (frame 2) divided the dataset randomly, 2 × 50%. LPA was conducted on both frames 1 and 2 (total sample and validation samples, respectively), as shown in Figure 1, informing sample comparisons and confirming the profile structure [43]. Subsequently, LPA was rerun on the whole sample to aid interpretation of the findings.

Sample frame 3 for discriminant analyses split the total sample into 5 parts, with 80% (3153/3941) or k−4 parts designated as the training sample and the remaining 20% (788/3941) or k−1 part designated as the test sample. The holdout procedure was repeated over 10 folds, in line with best practices [45]. Discriminant analyses tests were then applied over sampling frames 1 and 3 (Figure 1), with cross-validation describing classification accuracy for both total and holdout sampling frames using the classification error and Brier scores. As discussed earlier, the Brier score is a proper score function that measures the accuracy of probabilistic predictions [85]. More technically, the Brier score comprises the mean square difference between the true classes and predicted probabilities [82,83].

Discussions on prediction-focused structural equation modeling consider discriminant validity in terms of in-sample versus out-of-sample perspectives [45]. Bootstrapping has gained prominence in sample validation; however, this remains a decidedly in-sample approach. In this study, larger sample sizes permitted a combination of k-fold assessments of discriminant validity, combined with cross-validation modeling, which ascertain both in-sample and out-of-sample classification accuracy.

LPA was then conducted on the “T1 and T2” dataset, while a holdout procedure (Figure 1) tested the predictive relevance of the profile solution across 10 partitions, or “folds,” of this procedure. Predictive performance was summarized by prediction metrics, namely, the MAE and MAPE [86,87]. The MAPE was calculated using the absolute error in each fold divided by the observed values that are evident for that fold. Averaging these fixed percentages yields the MAPE, which indicates how much error arises from the prediction model compared with the true value [72].

The flow of participants throughout the study is shown in Figure S1 in Multimedia Appendix 3. In brief, the initial analytical set comprised 82,265 episodes and 8510 patients. Data were prepared to include team contacts and ensure completeness, resulting in 170,365 episodes and 6421 patients. Further processing resulted in the final analytical set of 7882 episodes from 3941 patients to optimize the analysis, with this complete dataset having no missing values.

Tables 2 and 3 display the distribution of patient characteristics from the LPA and across sampling frames. Descriptive statistics for holdout samples were limited to the initial fold (of 10). Comparing latent profiles to the whole sample means and distribution provides an understanding of the characteristics of each group of patients across identified latent profiles (LPs)—LP1 and LP2—in the 2-profile solution. Service receipt, measured as the number of contacts from clinical teams, was markedly higher for LP2 than for LP1. The LPA best characterized service receipt as being either high or low, captured through LP2 and LP1, respectively. While LP1 counted no psychology contacts, LP2 on average had almost 6 per episode of care. This pattern was consistent for medical teams and, in particular, nursing teams.

In general, LP2 received treatment at a rate 10 times higher than LP1. The derived profiles were not differentiated by sociodemographic characteristics; both latent profiles had similar proportions of female and White individuals as well as comparable levels of neighborhood deprivation. In total, 11.3% (120/1062) of LP2 had internalizing diagnoses, with no such diagnoses present in LP1.

Analyses identified 2 profiles emerging from the contemporaneous (T1) dataset and 4 profiles from the repeated measures dataset (T1 and T2; Table 4). Profiles were enumerated according to the BIC using mclust [35] for T1 versus T1/T2 dataset, selecting 2 and 4 profiles, respectively. The comparison of classification and likelihood indices shown in Table 4 suggested that the 2-profile solution from T1 data was preferred, showing superior model fit. Stage-1 LRT tests compared 2 and 3 class solutions against the null hypothesis, that is, the 1-class model. Both solutions demonstrated superior likelihood than the null hypothesis; however, the 2-profile solution was preferred. Stage-2 LRT tests compared 3-profile (the alternative) versus 2-profile (the null) solutions and failed to demonstrate any advantage of a 3-profile solution over a 2-profile solution.

In model selection, comparative LRTs described how a 3-cluster solution significantly outperformed the 1-class model but not the 2-class model (Table 4); therefore, a 2-cluster solution was preferred.

Posterior probabilities were calculated as the proportion of correctly assigned patients to the respective clusters. These probabilities are “posterior,” in the sense that they are calculated after individuals have been assigned a cluster membership [91]. Table S2 in Multimedia Appendix 1 displays classification indicators and posterior probabilities, and the classification error of the model was also negligible, being near 0. The Brier scores further demonstrate near-perfect classification accuracy of the model in forecasting individuals’ likely engagement subtype, given the values of routinely collected indicators. High posterior probabilities are also considered clear indicators of classification accuracy, with values ideally exceeding 0.80, in line with established guidelines [66].

Findings from evaluation of the measurement model, taken with the near-perfect degree of classification accuracy (Table S3 in Multimedia Appendix 1), suggest the usefulness of combined parameters in predicting cluster membership.

Relating parameter values to clusters is expressed in the language of direct and controlled direct effects [92]. The direct effect describes the (unconditional) relationship between parameters and clusters, while the controlled direct effect describes this relationship based on the individual probabilities of cluster assignment. In the unconditional model, all parameters determined class assignment, apart from gender (Table 5). Figure 2 displays conditional model evaluation, where face-to-face contacts (t_3940,3939=5.73; P<.001; B=0.006; 95% CI 0.005-0.006) and the number of nursing contacts (t_3940,3939=5.52; P<.001; B=0.0005; 95% CI −0.003 to 0.0006) significantly determined class assignment.

Likelihood-based approaches tend to dominate the mixture modeling landscape; however, there can be inferential challenges when adhering to the orthodoxy [63,93]. A more Bayesian approach may consider the likely parameter value, given the cluster membership. Structural model evaluation using this approach followed the same 3-step procedure as advised by Nylund-Gibson et al [39] and described in the section Fitting the Structural Model. Results from these analyses are presented in Table 1 and further detailed in Figure 3.

In the unconditional model and of all care contacts, psychology demonstrated the strongest relationship to clusters (t_3940,3937=23.83; P<.001; B=0.02; 95% CI 0.02-0.02); however, this was not upheld in the conditional model (t_3940,3937=1.72; P=.09). Overall, key prognostic targets of service receipt given the knowledge of cluster membership were nurse contacts. In both unconditional (t₃₉₄₀=19.53; P<.001; B=7.27; 95% CI 6.54-8) and conditional (t₃₉₃₉=5.52; P<.001; B=−74.94; 95% CI −119.34 to −30.56) models, cluster membership was significantly related to the receipt of nursing contacts.

Further testing modeled an interaction term between nursing contacts and cluster membership in determining nonattendance, as shown in Figure 3. The interaction term was significant both for face-to-face contacts (t_3940,3937=5.81; P<.001; B=0.74; 95% CI 0.49-0.99) and did-not-attends (t_3940,3937=3.62; P<.001; B=0.12; 95% CI 0.06-0.20), with confirmatory testing indicating that the rate of did-not-attends did not significantly differ from a Poisson distribution (N=3941, χ²₈=10.3; P=.20).

Figure 3 shows this relationship visually and indicates how did-not-attend rates vary between clusters as a function of nursing contacts. The assignment of profiles is shown probabilistically between either latent profile 1 and latent profile 2, which allows for assignment likelihood or uncertainty in clinical practice. For LP2, the number of nursing contacts made less of a difference in fitted rates of did-not-attends, compared to LP1. As it represented the majority of the sample, LP1 was more sensitive to the effect of nursing contacts regarding fitted rates of attended appointments.

Discriminant analysis modeling indicated high classification accuracy of the selected model (Table S3 in Multimedia Appendix 1). Both the MAE and MAPE were taken as summative indicators regarding the overall prediction error of the model. Over 10 folds of the holdout procedure, or train × test runs, the “grand” MAE was 0.004 (Table S3 in Multimedia Appendix 1), meaning that the model predicted correct class assignment with near-perfect accuracy. Conversely, out-of-sample tests on the replicative utility of engagement measures were less favorable. For face-to-face contacts, the MAE was 11.27, whereas for the number of did-not-attends, the MAE was 1.97. The MAPE may aid interpretation of these results, specifying that the number of face-to-face contacts erroneously predicts cluster membership in 11% (86/788) of cases (of 10 “folds”; Table S2 in Multimedia Appendix 1).

In the explanation-focused in-sample approach, the Brier score was 0.005. From an out-of-sample perspective, the average Brier score was 0.19 (Table S1 in Multimedia Appendix 1). The average classification error over 10 folds was 0.005, indicating substantially high classification accuracy. In addition, the average (out-of-sample) Brier score was 0.19 (SE 0.005), suggesting a nearly 80% rate of forecasted classification accuracy from an out-of-sample perspective.

This study aimed to characterize the patterns of service receipt in a large sample of patients with a diagnosis of PD in an urban, secondary care setting. The main findings of this study were identification of 2 LPs that best characterize service receipt in specialist PD treatments. LP2, at 26% (1062/3941) of the sample, received almost 5-fold as many nursing contacts as LP1 (9.43 vs 3.98) and 6-fold as many medical contacts (5.90 vs 1). In addition, general service use was almost 10 times greater in LP2 than LP1. Sociodemographic characteristics were not associated with LPs, but diagnostic factors were associated. Specifically, LP1 had no internalizing diagnoses. These findings suggest important distinctions of service receipt as a function of routinely collected data and suggest the discriminatory relevance of such descriptors in determining likely profile membership.

LPA identified statistically reliable categories of service use using routine observational data for those with a diagnosis of PD. Profiles distinguish service receipt for those with a diagnosis of PD and were replicated across a range of sampling frames. Discriminant analyses and cross-validation modeling suggest substantial discriminatory performance in secondary settings using these types of data. This study described 2 LPs, primarily distinguished by diagnostic information. Of the approximately 3000 patients in LP1 (n=2879, 73.05%), none had internalizing PD diagnoses. In inferential testing, patients in LP1 were twice as likely to have less service contact.

The literature discussing PD and service use describes greater service use in the general presence of PD [94,95]; however, this has yet to be robustly applied between different PD diagnoses. One study delineates the interplay between personality (disorder) and health [96], and there is rigorous work, such as comparisons of service use based on ethnicity in PD treatments [97]. Despite economies of scale offered by population-based interventions, these programs should be customized for subgroups on the basis of personality [96]. This study details foundational work in this area using a large sample in secondary care and attempts to export this comparative focus to diagnostic factors within PD treatments.

In a large study using administrative data in Quebec, emergency department visits were compared between “externalizing” or Diagnostic and Statistical Manual cluster B PD and schizophrenia, those diagnosed with both conditions, and the general population [98]. General practitioner visits were highest in the cluster B group, while the use of psychiatric (vs nonpsychiatric specialist or generalist) services differed in the presence of a PD diagnosis, with 1 mean annual visit versus 4, for schizophrenia.

Similar findings were captured in this study, which demonstrates lower service use for noninternalizing diagnostic groups. Greater service use was shown for the smaller latent profile, LP2. In practice, greater patterns of service use from the minority group arguably have predictive relevance in profile determination and tailoring treatment pathways accordingly. “Hotspotting,” as it is known in general practice [99], is a data-driven practice that aims to identify patterns of “extreme” use [100]. Complex needs and PD populations are similarly characterized by high use [101]. This responsive concept is well known in general practice but has yet to be specifically applied to mental health care. In the case of PD and complex needs, general practice is often a key contact [102]. As demonstrated in this study, understanding individual (stratified) patterns of service use offers an example for expanding hotspotting practices, commonly used by general practitioners, to address historical challenges in generalist service provision for complex needs. By applying advanced techniques such as LPA, this approach can help create more tailored treatment pathways.

Therefore, this study argues that more data-enabled, proactive service provision, engendered by hotspotting practices, may assist modernization of mental health provision and encourage more individually intelligible and tailored treatment pathways. The assertion is that the determination of cluster membership, given the routinely measured parameters, may help inform similar practices in secondary care with complex needs. Related work in improving access to psychological therapies highlights 8 LPs [43], which contrasts with the solution presented in this study.

Interestingly, the rates of nonattendance were stable across profiles, averaging around 1 in 10, which is broadly in line with the literature in general practice [103]. Service use was poorly distinguished as a function of the treating team. Structural evaluation revealed nursing contacts and did-not-attends of key interest; modeling the interaction term showed that the fitted rate of did-not-attends differs between clusters as a function of nursing contacts. The implication is that more intensive treatment pathways with more nursing contacts result in greater fitted nonattendance for LP1. However, for LP2, greater nursing contacts had much less of an effect on did-not-attends.

The use of a large sample, the multiplicity of analyses, and the development of a protocolized resampling scheme for validation comprise key strengths of this study. Discriminant validity and classification accuracy tests showed near-perfect classification for individuals with a diagnosis of PD into LPs (Table S2 in Multimedia Appendix 1). Generalizability was explicitly considered pursuant to resampling and validation, which may have implications regarding wider (national) policy implementation. From an out-of-sample perspective, for example, the model accurately predicted cluster assignment in nearly 80% of cases, across 10 “folds” of the holdout procedure, or train × test runs, compared with a Brier score of 0.005 for in-sample model performance (Table S1 in Multimedia Appendix 1). At least, for explanatory purposes, the model was able to classify all cases correctly.

A limitation of this study is that both exposure and outcome were measured over the same periods. This contemporaneous design was indicated in model selection; however, more frequent measurement may have resulted in more precise model solutions. Comparative modeling between both T1 and T1/T2 data suggested that a 2-cluster solution using T1 measurement showed superior model performance. In addition, the analytical design only included those who had more than 1 episode of care, which may exclude more sporadic or crisis presentations. Covariates were also drawn over time, rather than at presentation, rendering this study more suitable for validation than for strict prediction.

In addition, clinical outcomes were not assessed in this study, although there is scope to include these in future research. Finally, the direction of causality remains unclear—specifically, whether service receipt determines latent profile, or vice versa—along with the relationship between service receipt, profile, and clinical outcomes.

Truly “predictive” models of care are currently uncommon in mental health [5]. However, this study plays an important role in moving toward more predictive models of care in the medium term, following the formal validation of LPs and deployment of a rigorous resampling procedure in testing classification uncertainty at both average and individual levels. Through the identification and description of LPs, services are envisioned to be equipped to include such information in more personalized treatment approaches. More reflexively, such approaches can be seen as facilitating preventative models of care rather than reactive ones [104,105], representing a radical and, if not timely, necessary shift in care delivery.