Predictors of outcome following psychological therapy for depression and anxiety in an urban primary care service: a naturalistic Bayesian prediction modeling approach

Abstract

Background

England's primary care service for psychological therapy (Improving Access to Psychological Therapies [IAPT]) treats anxiety and depression, with a target recovery rate of 50%. Identifying the characteristics of patients who achieve recovery may assist in optimizing future treatment. This naturalistic cohort study investigated pre-therapy characteristics as predictors of recovery and improvement after IAPT therapy.

Methods

In a cohort of patients attending an IAPT service in South London, we recruited 263 participants and conducted a baseline interview to gather extensive pre-therapy characteristics. Bayesian prediction models and variable selection were used to identify baseline variables prognostic of good clinical outcomes. Recovery (primary outcome) was defined using (IAPT) service-defined score thresholds for both depression (Patient Health Questionnaire [PHQ-9]) and anxiety (Generalized Anxiety Disorder [GAD-7]). Depression and anxiety outcomes were also evaluated as standalone (PHQ-9/GAD-7) scores after therapy. Prediction model performance metrics were estimated using cross-validation.

Results

Predictor variables explained 26% (recovery), 37% (depression), and 31% (anxiety) of the variance in outcomes, respectively. Variables prognostic of recovery were lower pre-treatment depression severity and not meeting criteria for obsessive compulsive disorder. Post-therapy depression and anxiety severity scores were predicted by lower symptom severity and higher ratings of health-related quality of life (EuroQol questionnaire [EQ5D]) at baseline.

Conclusion

Almost a third of the variance in clinical outcomes was explained by pre-treatment symptom severity scores. These constructs benefit from being rapidly accessible in healthcare services. If replicated in external samples, the early identification of patients who are less likely to recover may facilitate earlier triage to alternative interventions.

Introduction

Current first-line treatments for depression and anxiety disorders include both psychological and pharmacological interventions. Whilst these have roughly equivalent efficacy (Cleare et al., 2015), patients report a preference for psychological therapies (McHugh, Whitton, Peckham, Welge, & Otto, 2013), although these are more costly to deliver (Koeser, Donisi, Goldberg, & McCrone, 2015). In the UK, the National Institute for Health and Care Excellence (NICE) recommends a range of psychological therapies for depression (National Institute for Health and Clinical Excellence (NICE), 2022), including cognitive behavioral therapy (CBT) and interpersonal psychotherapy (IPT), and CBT for generalized anxiety disorders (NICE, 2011). Although people with depression and anxiety disorders may respond to first-line psychological and pharmacological treatments, this is not the case for all patients. Response rates are typically reported to be between 50% and 60% (Cuijpers, Stringaris, & Wolpert, 2020; Roy-Byrne, 2015).

The Improving Access to Psychological Therapies (IAPT) program was initially designed as a talking therapy service to implement evidence-based psychological interventions for people with mild-to-moderate symptoms of anxiety and depression. The service, based in England, has been relatively successful, with delivery scaled up and current recovery rates at 51.9% in August 2020 (NHS Digital, 2021). This is consistent with other evidence of psychological treatment outcomes for mood and anxiety disorders (Cuijpers et al., 2014), where rates are calculated from treatment ‘completers’; defined as people who get at least two therapy sessions, excluding those who are signposted elsewhere on assessment or those who do not attend any sessions.

Given this variability in treatment outcomes, one important question is whether it is possible to identify baseline variables prognostic of therapy working well (for depression and anxiety outcomes). Some demographic and clinical variables are routinely available and may help in tailoring effective treatment choices. Previous research has shown mixed results in different settings using differing methodologies. For example, systematic reviews of factors predicting treatment response for depression (Taylor, Marwood, Greer, Strawbridge, & Cleare, 2019) found inconsistent evidence for a wide range of predictors across studies. Severe symptoms seemed to be the most consistent predictor. Similar findings are apparent with anxiety disorders (Mululo, de Menezes, Vigne, & Fontenelle, 2012). The majority of these studies have taken place in settings which were not representative of real-world care (Doorn, Kamsteeg, Bate, & Aafjes, 2021). Where real-world care data have been used, Coley, Boggs, Beck, and Simon (2021) found little predictive information in baseline variables, although only a limited range of variables were assessed.

Given these mixed results and limited settings, our PROMPT study (Grant et al., 2014) purposefully combined a naturalistic approach (by recruiting participants within an IAPT service) with a collection of a rich set of baseline or pre-therapy variables. Initial findings suggested that people within IAPT commonly have a range of psychiatric comorbidities (Hepgul et al., 2016; Strawbridge, Alexander, Richardson, Young, & Cleare, 2023). The ultimate aim was to comprehensively characterize the predictive information in pre-treatment variables associated with post-treatment outcomes and take the first step toward development of a prediction model for the IAPT service outcomes (Kent, Cancelliere, Boyle, Cassidy, & Kongsted, 2020). The set of variables covered demographic, clinical, diagnostic, and psychosocial domains, using measures intended to be feasibly administered in clinical practice and alignment with the extant literature (Grant et al., 2014). We define this study as exploratory within the Prognosis Research Strategy (PROGRESS) framework (Hemingway et al., 2013; Riley et al., 2013).

This paper addresses the primary aim of the PROMPT study in two parts. Firstly, we describe and summarize key baseline variables and use these as univariate predictors of therapy outcome status (recovered or not/improved or not). Secondly, we use a robust Bayesian prediction model strategy to assess multivariate predictors of outcome. Defining a set of useful predictors is a variable selection problem to which the Bayesian projective prediction approach is particularly well suited (Piironen, Paasiniemi, & Vehtari, 2020). Firstly a ‘reference’ model is developed which predicts the data (and associated uncertainties) as well as possible. Secondly, information from the reference model is ‘projected’ to sub-models to find the most parsimonious or simplest model which gives similar predictive performance to the reference model. This method is attractive both clinically and statistically. Clinically, it provides simple, interpretable models. Statistically, projective prediction models incorporate (rather than ignore) the uncertainty associated with the variable selection process. In previously used methods such as frequentist forward or backward selection, this selection uncertainty is ignored leading to over-optimistic models which do not generalize well to other samples. Other methods such as the LASSO contain little information about the uncertainty associated with model parameters. In short, this projective prediction method compares well with other methods of variable selection, both Bayesian (Piironen & Vehtari, 2017a) and non-Bayesian (Bartonicek, Wickham, Pat, & Conner, 2021) with the advantage of allowing valid estimates of the uncertainty of model parameters.

Method

Participants

The study was a naturalistic cohort study of patients referred to the Southwark IAPT Psychological Therapies Service in South London. Participants were recruited by PROMPT study researchers accessing electronic patient records. Lists of individuals who had agreed to take part in research when registered for therapy were obtained. Patients were contacted from this list and asked to take part in PROMPT. If willing, they were invited to a baseline interview with a trained postgraduate researcher (prior to therapy initiation) where assessments selected to measure variables predictive of subsequent therapy outcomes were collected. Subsequently, on starting therapy, patients' progress through therapy sessions was monitored using the treatment outcome data collected routinely within IAPT electronic patient records. The study team accessed and used this follow-up data as outcomes for PROMPT. Retention of blindness in relation to the collection of outcome data was not possible. All research visits took place at the National Institute for Health Research King's Wellcome Clinical Research Facility. The first patient was enrolled on 5 February 2014 and recruitment continued until 29 July 2016.

Participants were eligible to participate if they were aged 18 years or over, had provided informed consent, and were on a waiting list to receive psychological therapy. Exclusions from the research study (but not the IAPT service) were made if participants had already begun therapy or if they could not speak English well enough to carry out the baseline interview.

Treatments

Southwark IAPT talking therapies service implements a stepped care provision in which most individuals are initially offered low-intensity therapy. If individuals do not respond to low-intensity treatment, they are ‘stepped up’ to high-intensity therapy. Step 1 involves registering the individual with the service and a pre-treatment assessment. Individuals are then assigned to either a step 2 or 3 level treatment. Step 2 (low-intensity therapy) entails interventions such as four to six sessions of guided self-help based on CBT principles, and/or workshops such as CBT for sleep with a trained psychological well-being practitioner. Step 3 (high-intensity therapy) offers a course of, e.g. individual CBT, counseling, or IPT with a trained psychological therapist. Individuals can be stepped up from step 2 to 3, and/or stepped down (National Collaborating Centre for Mental Health, 2018).

Outcomes

We utilized four outcomes in this analysis: the IAPT recovery index (primary outcome), IAPT reliable change, and total severity scores from the depression and anxiety self-report scales. In fact, outcomes 3 and 4 were the self-report scales from which composite IAPT outcomes are derived; the Patient Health Questionnaire for depression (PHQ-9; Kroenke, Spitzer, & Williams, 2001) and the Generalized Anxiety Disorder assessment for anxiety (GAD-7; Spitzer, Kroenke, Williams, & Löwe, 2006). The PHQ-9 has nine Likert scale items and the GAD-7 has seven items. Patients rate each symptom from 0 ‘not at all’ to 3 ‘nearly every day’ with total scores (maximum 27 on the PHQ-9 and 21 on the GAD-7) used as indicators of symptom severity. These scales have recently been shown to be unidimensional with temporal measurement invariance (Stochl et al., 2020).

The recovery index is a joint consideration of the clinical thresholds of these two scales. Patients are considered a ‘clinical case’ if they score 10 or more on the PHQ-9 or score 8 or more on GAD-7. Recovery means moving below these clinical cut-offs on both questionnaires at discharge from the IAPT service. The recovered YES/NO classification reflects individuals who have recovered from case-ness at the end of treatment v. those who have not. Here we relax the condition of pre-treatment case-ness as including non-cases, in alignment with the naturalistic approach of our study. Additionally, patients who are symptomatic but not cases may still have a clinically meaningful change in outcomes. However, given the IAPT definition of recovery, we include a sensitivity analysis in which recovery is predicted for cases only.

The second outcome, IAPT reliable improvement, considers whether individuals undergo a positive ‘reliable change’ (Jacobson & Truax, 1991) in either PHQ-9 or GAD-7 scores. To count as reliable improvement for one outcome there must be no negative change in the other. This outcome putatively identifies patients whose symptoms improve sufficiently beyond that expected by measurement error, even if they do not meet IAPT definition of recovery. Note that the reliable recovery index also exists, a further composite of both the recovery index and change score (Gyani, Shafran, Layard, & Clark, 2013). As a composite of a composite, we consider reliable recovery to be particularly difficult to interpret, and accordingly, we only report this in supplementary information also as preliminary analysis showed this outcome to have little association with predictors.

To improve prediction performance and stratify for differences in severity at the start of treatment, PHQ-9 and GAD-7 scores at the start of treatment (first therapy session) were included in present analyses.

Baseline predictors

Measures included as predictors were collected at the baseline visit and were divided into three domains: demographic, clinical/diagnostic, and psychosocial. Demographic variables were age, gender, ethnicity, education, employment status, income, and relationship status. Clinical status was indicated by patient age (in years) from the first onset of psychiatric symptoms, the duration of the current episode in months, and current anti-depressant use. Relevant co-occurring diagnoses selected from the MINI diagnostic interview (Sheehan et al., 1998) were bipolarity, diagnoses of social phobia, panic disorder, agoraphobia, obsessive compulsive disorder (OCD), post-traumatic stress disorder (PTSD), alcohol abuse, substance abuse/dependence, anorexia, and bulimia. Other clinical assessments included were personality disorder traits (Standardized Assessment of Personality – Abbreviated Scale [SAPAS]; Moran et al., 2003) and physical health (Cumulative Illness Rating Scale [CIRS]; Linn, Linn, & Gurel, 1968).

Psychosocial predictors were: current beliefs about mental illness (Brief Illness Perceptions questionnaire [B-IPQ]; Broadbent, Petrie, Main, & Weinman, 2006); self-criticism by negative cognitions and self-reassurance (Forms of Self Reassurance and Self Criticism scale [FSRC]; Gilbert, Clarke, Hempel, Miles, & Irons, 2004), self-efficacy (General self-efficacy scale; Schwarzer & Jerusalem, 1995), quality of life (EuroQol questionnaire [EQ5D]; Herdman et al., 2011), social support (Oslo Social Support scale [OSS]; Dalgard et al., 2006), stressful life events (List of Threatening Events Questionnaire [LTE]; Brugha, Bebbington, Tennant, & Hurry, 1985), and adverse events during childhood (Childhood Trauma Questionnaire [CTQ]; Bernstein et al., 2003). Further description of the measures and variables are detailed in the supplementary information.

Statistical analysis

Descriptive statistics

To report our findings we use the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD) guidelines (Collins, Reitsma, Altman, & Moons, 2015). Table 1 summarizes patient characteristics by recovery condition. Median and interquartile range (IQR) are reported for continuous measures, and frequency and proportion for categorical. Group differences in continuous measures were assessed with the Wilcoxon rank-sum test and for categorical variables using the c² test (or Fisher's exact test if expected counts were less than 5 in table cells).

Table 1. Potential predictor variables by recovery status at the end of treatment

PHQ, Patient Health Questionnaire; GAD, Generalized Anxiety Disorder scale; SAPAS, Standard assessment of Personality – Abbreviated Scale; CIRS, Cumulative Illness Rating Scale; IPQ, Illness Perceptions Questionnaire; FSCRS, Forms of Self-Criticizing/Attacking and Self-Reassuring Scale.

^a Median (IQR); n (%).

^b Wilcoxon rank-sum test; Pearson's c² test; Fisher's exact test.

Pre-processing of predictors and univariate associations

Prior to developing the prognostic models and to reduce the opportunity for overfitting, a set of 30 candidate variables were identified from those collected at baseline (Grant et al., 2014) by authors A. J. C., R. S., A. H. Y., and J. H. These were intended to cover as wide a range of variable domains as possible whilst excluding potentially redundant information; for example, relationship support would be captured better by the social support scale than categorical relationship status. For the prognostic modeling, pre-processing of predictors was made to ease both computation and interpretation, i.e. categorical variables simplified (cf. SI methods). Due to the low frequency of occurrence, eating disorder diagnoses were not included, leaving 28 variables. Following this, independent associations between the selected predictor variables and binary outcomes (recovery and reliable improvement) were assessed using Firth penalized logistic regression, adjusted for age, gender, and pre-therapy PHQ-9 and GAD-7 scores.

Prediction modeling

To build a robust and interpretable prognostic model we used Bayesian projective prediction for variable selection. This is a recently developed method with two distinct phases (Piironen et al., 2020). Firstly, a reference model is constructed including all candidate variables as the best model of the data (Pavone, Piironen, Bürkner, & Vehtari, 2022). Best predictive performance is evaluated as the predicted out-of-sample performance (how well the model would predict new data), estimated via leave-one-out cross-validation (LOO-CV). Prediction performance was measured by the expected log predicted density (ELPD; Gelman, Hill, & Vehtari, 2020). In constructing the reference model, we addressed potential overfitting by using regularized horseshoe priors, which have been shown to perform well in comparison to other informative priors (Piironen & Vehtari, 2017a). These function by weighting model coefficients toward zero but allowing larger coefficients with more information to escape this weighting. In setting the prior, the number of expected non-zero variables is defined. Here, an expectation of five non-zero variables seemed reasonable with the aim of producing a parsimonious model.

Having defined a reference model, the second phase variable selection was undertaken via projective prediction (Piironen et al., 2020). The latter involves both a search and evaluation component. In the search, variables are ranked for inclusion by the degree to which they reduce the Kullback–Leiber divergence (a measure of information) between the reference and null model. Sub-models (with variables input by rank) are compared to the reference models using the ELPD. Non-inferiority in predictive performance was set at least 1 standard error of the ELPD difference and determined the number of variables in the model. Cross-validation was used for both search and evaluation procedures to avoid overfitting, but this also meant it was possible to evaluate model stability. Model stability was estimated by how often variables were included in cross-validation sets for each different size sub-model, i.e. the number of times a particular variable was included in a sub-model of size 2, 3, 4, and so on up to 30. The coefficients are interpretable as associations between the predictor and outcome conditional on the model selection process (and as is clear for prediction models, do not represent estimates of causal associations).

Bayesian models were specified with the R package rstanarm (Goodrich, Gabry, Ali, & Brilleman, 2022) which fits models using Hamilton Monte Carlo Markov chains (HMC). HMC chains work by repeatedly sampling from a combination of prior information and data (posterior distribution) to reach a best (and stable) solution. Models were run with four chains of 1000 warm up and 2000 sampling iterations and checked for adequacy with diagnostic tests for convergence: visual, Gelman and Rubin statistics (Rhat) and effective sample size (ESS). All models reported below passed these checks with credible results (Rhat near 1 and ESS > 2000). Model coefficients are reported as median and 95% credible intervals of the posterior distributions.

Model performance

To assess logistic regression model performance for recovery and reliable improvement outcomes, we generated predicted probabilities adjusted for overfitting by cross-validation. Using these predicted probabilities, metrics for model validation and calibration were calculated. Firstly, the C statistic or area under the curve (AUC) indicates how well the model can discriminate between recovered and non-recovered cases. The AUC is the proportion of pair-wise comparisons between participants with an event (recovery) with to those without. To address calibration (how well outcome probability matches the true underlying probability), the calibration intercept and slope are estimated by a logistic regression of the outcomes by predicted probabilities. For the calibration slope, deviations from 1 indicate the degree to which predicted probabilities do not agree with those expected across the range of risk of the outcome. The calibration intercept gives a measure of overall calibration with 0 good, and less or greater than 0 reflecting over or underestimation of risk, respectively. For continuous outcomes (PHQ-9 depression and GAD-7 anxiety severity scores), the root mean squared error of the fitted linear models was reported. Finally, for all models, the Bayesian R ² (Gelman, Goodrich, Gabry, & Vehtari, 2019) was calculated as a measure of explained variance.

Missing data

There were missing data in both the outcome and predictor variables. Participants with missing outcomes were removed from further analysis (45/308, 14.6%) as they had not had therapy (cf section ‘Sample size and candidate variables’). In total, there were 446 missing predictor values (6.3%) which varied according to predictor. For 20 predictors, the total number of missing values was low (<6%) with only BME status (22%), current episode duration (23%), age of onset (19%), and CIRS (19%) above this. In total, 132 (50%) of participants had a full set of predictor variable data and 248 (94%) with three or less predictor values missing. We generated a single imputation of missing data prior to modeling using a non-parametric random forest algorithm available in the MICE package. This algorithm can handle datasets with both continuous and categorical variables and makes few assumptions about the multivariate structure of the data (Waljee et al., 2013).

Sample size and candidate variables

This was a naturalistic, observational study. A priori sample size estimates were based on patient throughput, and it was estimated that 600 patients could be recruited into the baseline PROMPT interview across the recruitment period. This was an overestimation, and 371 were ultimately recruited. It was anticipated that 28% of patients would drop-out of treatment (Grant et al., 2014). Ultimately, drop-out and enrolment of patients who had already started therapy led to the exclusion of 108 (29%) individuals (see Fig. 1) which gave an event-to-predictor ratio of 4.1 (124 non-recovered/263). In fact, this reduced sample size supported the adoption of a Bayesian prediction modeling approach with informative priors. Regularized horseshoe priors were used as these have been shown to perform well (i.e. produce stable and reliable models) when the ratio of events to predictors is quite low (Piironen & Vehtari, 2017b).

Figure 1. Flowchart of participants through the study.

Results

Participant flow

Of 1806 patients identified and 922 assessed for eligibility, 371 (40.2% of those deemed eligible) were enrolled into the study. In total, 263 patients completed therapy, defined by attendance at least two therapy sessions, Fig. 1. The median number of therapy sessions was 8 (IQR: 6, 14; minimum 2, maximum 32).

Outcomes

Of the 263 eligible patients who completed the study, 139 (52.9%) were classed as recovered. Median PHQ-9 scores were 4 (IQR: 2, 6) in the recovered group and 15 (IQR: 11, 19) in the non-recovered group. Median GAD-7 scores were 3 (IQR: 2, 5) in the recovered group and 13 (IQR: 9, 16) in the non-recovered group. In total, 53.6% of patients (141/263) were classified as having made a reliable improvement. In total, 99/141 (70.2%, or 37.6% of the whole sample) met both response definitions of recovery and reliable improvement.

Therapy allocation

Of the 263 patients enrolled in the longitudinal study, 99 (37.6%) were allocated to step 2 therapy (commonly computerized CBT or guided self-help) and the remaining 164 (62.4%) to step 3 therapy (largely face-to-face counseling or CBT). Patients had up to three separate treatment courses with the majority receiving one course of treatment. Recovery was more likely after one treatment (123/220, 57% recovered) relative to two (12/30, 40% recovered) or 3 treatments (3/12, 25% recovered). Patients who had recovered/responded after step 2 treatment were discharged whereas those who did not respond were likely to move to step 3.

Sample characteristics for recovered and non-recovered groups

Table 1 shows differences in baseline characteristics between recovered and non-recovered patients. These group differences were observed for several variables. In terms of demographics there were differences in employment and income with the non-recovered group having larger proportions of people in the unemployed or long-term sick categories relative to recovered. The recovered group had higher household income than non-recovered. Clinically, for recovered v. non-recovered patients, their lifetime age of onset was higher (median 18 v. 15) and current episode duration was shorter (median 7 v. 18 months). Further, recovered patients were less likely to have diagnoses of panic disorder (40% v. 56%), agoraphobia (33% v. 57%), social phobia (27% v. 45%), OCD (11% v. 30%), PTSD (6.5% v. 20%), and reported substance abuse (3.6% v. 11%). Recovered patients also had fewer personality disorder traits indicated (SAPAS score 3 v. 4), and better scores on psychosocial measures: they had more positive illness perception scores, higher self-efficacy ratings, were less self-critical (lower negative cognition and higher self-reassurance), had higher QoL scores, better social support, and lower reports of childhood trauma.

At the start of therapy, recovered patients had lower PHQ-9 and GAD-7 scores. For the PHQ-9, the median score was 10 in the recovered group v. 18 in the non-recovered group (p < 0.001). Similarly, median GAD-7 scores were 9 in the recovered v. 15 in the non-recovered group (p < 0.001). Not all patients were cases at pre-treatment (treatment session 1): 211 of the 263 (80.2%) patients were IAPT cases. Of the 52 non-cases, 48 (92.3%) were recovered at the end of therapy. Four patients who were non-cases at the start of treatment were cases at outcome.

Univariate predictors of recovery and reliable improvement

Univariate logistic regressions with recovery and reliable improvement as outcomes showed associations with several variables, but a more limited range than the group differences described in section ‘Sample characteristics for recovered and non-recovered groups’ when adjusted for age, gender, and pre-treatment depression (PHQ-9) and anxiety (GAD-7); cf. Table 2. For recovery, unemployment through long-term sickness, higher personality disorder trait scores, and diagnoses of agoraphobia and OCD and lower QoL were associated with lower odds of recovery. Only agoraphobia and OCD diagnoses remained statistically significant after FDR adjustment for multiple comparisons. Patients were less likely to indicate reliable improvement if they had lower educational qualifications, higher personality disorder scores, were taking antidepressant medication, met criteria for agoraphobia, social phobia, and OCD, and had lower self-efficacy and QoL ratings. Following adjustment for multiple comparisons, associations with personality disorder, OCD diagnosis, self-efficacy, and QoL remained statistically significant.

Table 2. Univariate predictors of Recovery and Reliable Improvement outcomes, adjusted for age, gender, and pretreatment depression (PHQ 9) and anxiety (GAD 7) scores

PHQ, Patient Health Questionnaire; GAD, Generalized Anxiety Disorder scale; SAPAS, Standard assessment of Personality – Abbreviated Scale; CIRS, Cumulative Illness Rating Scale; IPQ, Illness Perceptions Questionnaire; FSCRS, Forms of Self-Criticizing/Attacking and Self-Reassuring Scale; EQ-5D, EuroQol Quality of Life Instrument 5D (3 level) Utility Index; OSS, Oslo Social Support scale; LTE, List of Threatening Events; CTQ, Childhood Trauma Questionnaire.

^a False discovery rate correction for multiple testing.

Bayesian prediction modeling

Recovery

For the recovery outcome, the reference model had an AUC of 0.80 and Bayesian R ² of 0.26 (95% Bayesian credible intervals [bCI] 0.16–0.36). Internal calibration was good with slope and intercept at 0 and 1 (Table 3). Regularization via horseshoe priors meant most variables had an OR of 1 with only pre-treatment PHQ-9, OCD, and agoraphobia diagnoses showing a discrepancy from 0 of greater than 0.01 (cf. SI Table 2.2.1). Projective inference resulted in a parsimonious model with two predictors for performance (Fig. 2). Recovery was negatively associated with higher PHQ-9 scores (OR 0.30, 95% bCI 0.20–0.41) and OCD diagnosis (OR 0.75, 95% bCI 0.52–0.98). Using the divide-by-4 rule for the log-odds-ratio (Gelman et al., 2020), a 1 s.d. increase in baseline PHQ-9 (6.5 points in this sample) decreases the probability of recovery by 32%. A diagnosis of OCD reduces the probability of recovery by a maximum of 7.2%. In terms of model stability, both variables were included in the two-variable model for 100% of the cross-validation sets (SI Table 2).

Table 3. Performance statistics for reference and projection models for recovery, reliable improvement, depression (PHQ 9), and anxiety (GAD 7) scores

C (ROC), C statistic (area under the curve); ELPD (s.e.), expected log predictive density (standard error); PHQ, Patient Health Questionnaire; GAD, Generalized Anxiety Disorder.

Figure 2. Odds ratios (OR) for posterior distributions from projection models for recovery and reliable improvement outcomes, visualizing the posterior distribution with the black dot and bar representing the median and 95% Bayesian credible intervals (CI). The pre-treatment predictors are: PHQ, Patient Health Questionnaire; GAD, Generalized Anxiety Disorder scale; SAPAS, Standard assessment of Personality; QoL: EQ-5D, EuroQol Quality of Life Instrument 5D (3 level) Utility Index; OCD, diagnosis of Obsessive Compulsive disorder.

The model showed similar results when including only cases at baseline; only PHQ-9 and OCD diagnosis were included in the projection model. However, there was some concomitant reduction in model performance (AUC = 0.72, Bayesian R ²: 0.14, bCI 0.06–0.23). The association with pre-treatment PHQ-9 was weaker (OR 0.45, 95% bCI 0.30–0.62) but the OCD association was near identical (OR 0.73, 95% bCI 0.48–0.99). PHQ-9 and OCD were in 100% of the cross-validation sets (SI Table 4).

Reliable improvement

The reference model for reliable improvement had a Bayesian R ² of 0.12 (95% bCI 0.04–0.19) and AUC of 0.69; predictive performance was notably lower than that for recovery. The projection model included four terms (Fig. 2). The odds of reliable improvement increased for higher baseline GAD-7 scores (OR 2.27, 95% bCI 1.59–3.14), and for higher baseline QoL (OR 1.21, 95% bCI 0.99–1.69). The odds of reliable improvement decreased with an OCD diagnosis (OR 0.86, 95% bCI 0.59–1.02) and increased SAPAS scores (OR 0.85, 95% bCI 0.59–1.02). In terms of variable stability, these four variables were included in 100% of cross-validation folds (SI Table 6).

Continuous depression outcome (PHQ-9)

For PHQ-9, the reference model (SI Table 7) had a Bayesian R ² of 0.37 (95% bCI 0.25–0.46). Projective prediction resulted in a sub-model with two variables, baseline depression (PHQ-9) and QoL scores; see Fig. 3a. More severe depressive symptoms (higher pre-treatment PHQ-9) were associated with higher PHQ-9 scores at outcome (b = 3.45; 95% bCI 2.67–4.18), as were lower QoL scores (b = −1.33; 95% CI −2.12 to −0.39). Note that a sub-model additionally including OCD resulted in equivalent performance to the reference model (b = 0.77; 95% CI 0.04–1.51). The selected variables were included in the 100% of cross-validation folds (SI Table 8).

Figure 3. Coefficient posterior distributions from projection models for depression (PHQ-9) and anxiety (GAD-7) continuous outcomes, visualizing the posterior distribution with the black dot and bar representing the median and 95% Bayesian credible intervals (CI). The pre-treatment predictors are: PHQ, Patient Health Questionnaire; GAD, Generalized Anxiety Disorder scale; QoL: EQ-5D, EuroQol Quality of Life Instrument 5D (3 level) Utility Index; Agoraphobia diagnosis.

Continuous anxiety outcome (GAD-7)

For the anxiety reference model (SI Table 9) Bayesian R ² was 0.31 (95% CI 0.21–0.40). Using projective prediction, a sub-model with four variables gave performance equivalent to the reference model. More severe anxiety symptoms (higher GAD-7 scores) at outcome were associated with higher pre-treatment GAD-7 (b = 1.39, 95% bCI 0.21–2.42) and PHQ-9 scores (b = 1.42, 95% bCI 0.14–2.45), lower QoL (b = −0.76, 95% bCI −1.49 to −0.10), and meeting criteria for agoraphobia (b = 0.83, 95% bCI 0.18–1.53). Again, model stability was good with variables appearing in 100% of cross-validation sets (SI Table 10).

Discussion

Summary of findings

This study set out to assess whether clinically useful baseline predictors of psychotherapy outcomes could be identified in a naturalistic setting (an English psychotherapy ‘IAPT’ service). Including a rich set of baseline predictor variables, we showed in a descriptive analysis that many of these predictors carry information predictive of outcomes (cf. Tables 1 and 2). However, combining these predictors in a multivariable prediction model showed that baseline predictor information adds little beyond that contained in pre-treatment depression and anxiety severity scores. This held for all versions of the clinical outcome. Not only were potentially useful baseline predictors few in number (projection models contained a maximum of 3 baseline variables in addition to pre-treatment depression and anxiety), but the strength of the associations between these baseline predictors and outcomes was relatively small (Figs 2 and 3). The usefulness of these baseline predictors in developing prognostic models for treatment outcome may be only modest.

Selection of predictor variables

There was some variation in the variables included in the prediction models by outcome. Considering baseline symptoms, pre-treatment depression severity (PHQ-9) predicted both recovery and PHQ-9 severity outcomes; both pre-treatment anxiety and depression severity predicted GAD-7 outcome and only pre-treatment anxiety for reliable improvement. These differences are due to the correlation structure between PHQ-9 and GAD-7 scores at pre-treatment and outcomes. For example, the partial correlations between post-treatment PHQ-9 and pre-treatment GAD-7 score are independent given (conditioned on) pre-treatment PHQ-9 scores. In the case of reliable improvement, since it is a change score, the effect of pre-treatment scores may be due to residual confounding. In fact, higher GAD-7 pre-treatment scores predicted greater improvement consistent with regression to the mean.

Although contributions to prediction performance were relatively modest, there were some variables providing additional information, although this varied by outcome. Not meeting the criteria for OCD was associated with better odds of recovery and reliable improvement. Secondly, higher QoL was associated with higher odds of reliable improvement, and with less severe depression and anxiety scores at the end of treatment. Thirdly, fewer personality disorder traits (lower SAPAS scores) were associated with higher odds of reliable improvement. Finally, not meeting criteria for agoraphobia was associated with lower endpoint anxiety (GAD-7). Whilst there were some discrepancies between the variables included per outcome, it should be noted that these four variables were identified for all models in the top 5 or 6 ranked items in the search algorithm and all these variables could contribute positively to the prediction of outcomes. For example, for the PHQ-9 severity outcome, if QoL was not included in variable selection, a combination of OCD, agoraphobia, and SAPAS would give match reference model performance (Bayesian R ² = 0.37, SI Table 12). These other variables, although their predictive information is relatively modest, may suggest targets for treatment. For example, Simkin, Hodsoll, and Veale (2022) showed the interdependence of change in OCD and depression symptoms through therapy (in an observational cohort of OCD patients).

As indicated in the descriptive analysis in Tables 1 and 2, this does not mean that variables not selected in the projection prediction model cannot be informative for outcomes. There were several univariate associations for both recovery and reliable improvement as shown in Table 2. All these variables contain predictive information about the outcome which overlaps with the pre-treatment variables in the projection models. To illustrate this point, we re-ran the Bayesian variable selection process excluding the PHQ-9 and GAD-7 pre-treatment scores from the projection pathway for the PHQ 9 outcome. In estimating the predictive performance of the remaining variables (relative to the reference model), whilst not able to match the reference model performance for the recovery outcome, a three- or four-variable model with QoL (EQ5D), illness perceptions, OCD diagnosis, and income reached within 1 s.e. of maximum performance at Bayesian R ² = 0.24 (SI Table 13). Similarly for GAD-7, a four- or five-variable model with QoL (EQ5D), illness perceptions, OCD, agoraphobia diagnoses, and income had an R ² of 0.24 (SI Table 14).

Comparisons of outcomes

Overall, model performance differed between the outcomes, with recovery and dimensional continuous scores performing better than reliable improvement (as Bayesian R ² can be calculated for all models, Table 3). One point is that reliable improvement is based on change scores which in terms of predictive performance are less statistically efficient than endpoint scores as they contain two sources of variability (pre and post treatment). A further loss of power results from the dichotomization of the continuous measure of change and the threshold for change being conservative. Indeed, reliable change assumes homogeneity of measurement error and is liable to recording false positives when measurement error is high, missing true smaller change when measurement error is low (McAleavey, 2024). Whilst recovery showed better performance than reliable improvement it should be noted that excluding non-cases (as per the service definition of recovery) reduced performance considerably (R ² = 14% for cases only v. 26% including non-cases). The best predictive performance was for continuous depression and anxiety outcomes. Arguably these outcomes are the best choice for predicting patient outcomes. Whilst the recovery index is typically preferred by stakeholders for purposes of interpretation, classification of patient response could be made using predicted scores from the projective prediction models for PHQ-9 and GAD-7. Note that the mean absolute difference between observed and predicted scores from the projection models here was smaller than the criteria for reliable improvement; 4.5 for the PHQ-9 and 3.8 for the GAD-7 relative to the range of the scales (0–27 PHQ-9 and 0–21 GAD-7).

Relation to literature

The finding that a simple model based primarily on pre-treatment depression or anxiety severity predicts outcomes after psychological therapy mirrors other recent attempts to identify predictors of anxiety and depression treatment outcomes. Several studies are also relevant from the psychiatric prediction modeling literature. In a naturalistic psychiatric hospital setting, Webb et al. (2020) predicted 2-week outcomes (PHQ-9 scores) for 484 inpatient depression patients undergoing various multimodal therapies, predominantly intensive group CBT over a 2-week period with adjunct individual CBT therapy and medication if needed. As here, the most important predictor for the trained model was PHQ-9 baseline scores but the inclusion of 13 other demographic and clinical variables (out of a possible 51) significantly improved prediction performance (increasing R ² by 8%). Interestingly, there were similarities and differences to our results with some reflecting the different set of included variables (e.g. QoL). OCD was included in Webb et al.'s elastic net model, along with other anxiety diagnoses (social anxiety, PTSD, and GAD scores); whilst our reference model for PHQ-9 indicated predictive information in the social anxiety, agoraphobia, and PTSD variables (SI Table 8), the projection model showed that these were not needed for equivalent prediction performance. One advantage of the Bayesian projection approach here is that a minimal set of predictive variables with associated uncertainty can be found in addition to standard regularization methods.

Conversely, other studies have been more negative about the information value of prognostic variables other than baseline outcome. The NESDA study (Dinga et al., 2018) predicted MDD diagnosis (2 years post baseline, n = 804) for unipolar depression patients recruited from the general population, undergoing psychotherapy, pharmacotherapy or no treatments. Prediction performance was modest overall (AUC = 0.66) with other variables from psychological, clinical, and biological domains adding only 0.01 to that prognostic value. Similarly, two anti-depressant treatment studies found that baseline depression severity items were the strongest predictor of remission, namely the GENDEP study (Iniesta et al., 2016), and Chekroud et al. (2016), the latter using data from the STAR*D trial. Whilst Iniesta et al. showed clinically useful prediction performance (AUC = 0.72), the prediction model for the Chekroud study showed poor performance with external data from other trials (max AUC = 59.7%). The pre-dominance of baseline severity symptom score appears to generalize over both therapy and anti-depressant treatments.

Finally, two recent studies focused on psychotherapy outcomes using routine data are particularly relevant. Coley et al. (2021) developed and evaluated prediction models for psychotherapy (typically CBT) outcomes in depression. The strongest predictor for PHQ-9 outcomes was baseline PHQ-9 scores with only anxiety and medication further included in the model. In contrast to the results here, model performance both for this outcome and dichotomized 50% reduction was poor. The model only explained about 12% of PHQ-9 outcome variance whilst the AUC for 50% reduction was 57% and calibration poor. Similarly, an examination of outcomes to outpatient CBT (Hilbert et al., 2020) found that performance of the prediction models was robust, but not at the level of clinical utility (for remission AUC < 0.6 and dimensional scores R ² < 0.05). Both these studies used routine care databases with separate validation sets rather than estimating out-of-sample performance using cross-validation. Given the contrast in these findings to the relatively good performance of our prognostic models identified, the translation of findings from smaller research-led studies such as this one to routine care will need well-designed validation.

Strengths

The PROMPT naturalistic study was designed as a minimal adjunct to the IAPT service delivery and as such reflects care and outcomes in real-world service provision. Consequently, our prognostic factors and model identified have potential applicability to clinical practice. Secondly, the use of modern Bayesian methods and the use of a reference model allowed an evaluation of variable selection which was statistically valid. In estimating of out-of-sample performance through cross-validation we were able to evaluate the stability of variable selection and uncertainty of the final model parameters. Thirdly, we evaluated a wide range of outcomes, adding dimensional scores in addition to the composite IAPT recovery and reliable improvement outcomes to the analysis. As continuous measures (as well as being the outcomes the IAPT recovery index was derived from), they hold more information and likely support better prediction of patient outcomes.

Limitations

Arguably the most important limitation of this study is the fact there was no control group, meaning it is not possible to move beyond prognosis, that is, distinguish predictors of recovery due to treatment from those associated with recovery but not specifically related to treatment. A randomized design, with repeated baseline measures of predictors to account for measurement error, would help to resolve this. This study would also have been improved with a new test dataset to externally validate the prediction model. Whilst the out-of-sample performance can be estimated using LOO-CV, a more robust test of performance would be using data collected as part of a different study. Thirdly, the PROMPT sample may not be representative in terms of all IAPT service-user settings and population. Recruitment was only undertaken in South London, whilst IAPT is a nationwide service. There may also be differences in the recruited sample (who both engaged in the IAPT process and consented to take part in additional research interviews) compared with non-engagers and non-consenters. Including non-cases at baseline, in our primary analysis, brings challenges although the exclusion of non-cases did not strongly affect our results. Another issue is measurement error, a pervasive problem in clinical research and it is unclear how measurement bias in predictor variables may affect variable selection. Finally, recruitment did not reach the intended target meaning the sample size was smaller than planned. Because of this we implemented a clinician-led pre-analysis selection of candidate variables from all those measured. Whilst this was a change to the planned protocol, recent studies have shown selection by clinician to perform as well as machine learning or data-driven methods (e.g. Fusar-Poli et al., 2019). On the other hand, we may have excluded relevant predictor variables in this analysis. It may further be the case that symptom-level predictors are important, and these questions will be the subject of further work. There may also be factors unmeasured in our study, e.g. biomarkers which may have had predictive utility.

Implications and conclusion

The present study showed that Bayesian prediction models performed well in predicting depression and anxiety (and to a lesser extent recovery) at outcome for a cohort from a South London IAPT service. Variables prognostic of outcome were predominantly pre-treatment depression and anxiety symptoms, but with some more modest contribution from other predictors (poor outcomes predicted by lower QoL, OCD, and agoraphobia diagnoses and higher personality disorder traits). These findings have potential implications for IAPT treatment. Specifically, patients with more severe symptoms may require earlier stepping up of care, longer treatments or adjunctive/combination therapy approaches. Patients with particular comorbidities (such as OCD) may be less likely to respond to more generic therapy for depression and anxiety symptoms, raising the possibility of screening for those and offering targeted therapy instead of, or in addition to, standard IAPT approaches. The poorer outcomes in those with agoraphobia diagnoses may be indicative of a need to address avoidance behaviors. Similarly, higher scores on SAPAS screening could inform treatment choices. Clearly this work requires replication, and potential interventions based on screening would need further evaluation. Nevertheless, if sufficiently replicated in further development and validation studies, such an approach could have implications for improving future clinical outcomes by tailoring the stepped care program to better serve patients identified by those factors prognostic of worse outcome.