We followed the PROBAST (Prediction Model Risk of Bias Assessment Tool) [33] and TRIPOD+AI (Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis+AI) guidelines [34] when conducting analyses and reporting findings. The study was preregistered on the Open Science Framework [35]. Minor changes made to the preregistered protocol are listed in Multimedia Appendix 1.

This study was approved by the Stockholm Regional Ethical Review Board (REPN 2009/1089-31/2 and 2014/1897‐31) and adhered to the Declaration of Helsinki. All participants provided written informed consent before entering the study. Data were pseudonymized prior to being used for research purposes in accordance with Swedish regulations.

This study used data from the MULTI-PSYCH cohort comprising 2668 patients treated with ICBT for depression (n=1300), panic disorder (PD; n=727), or social anxiety disorder (SAD; n=641) at the Internet Psychiatry Unit of the Psychiatric Clinic Southwest at Karolinska University Hospital, Huddinge, Sweden, between 2008 and 2020 [36]. All treatments were provided over 12 weeks with therapist guidance via asynchronous messaging and consisted of similarly structured psychoeducation and exercise modules, with disorder-specific differences in intervention content (eg, behavioral activation for depression and exposure therapy for anxiety disorders). During screening, patients completed multiple questionnaires and answered interview questions about their symptomatology, medical history, and socioeconomic status, as well as provided blood or saliva samples for DNA analysis. Each patient’s unique Swedish personal identity number was used for record linkage across multiple nationwide registers: the National Patient Register [37], the Stockholm regional health care database VAL [38], the National Prescribed Drug Register [39], and the LISA (Longitudinal Integrated Database for Health Insurance and Labor Market Studies), which is managed by Statistics Sweden and combines data from the Register-based Employment Statistics, the Total Population Register, the Education Register, the Register of Income and Taxation, the Swedish Public Employment Service, and the Swedish Social Insurance Agency [40].

Given the substantial etiological and pathological overlap among internalizing disorders, treatment similarities, and a growing interest in transdiagnostic clinical prediction tools, the 3 disorder-specific subsamples were combined into a single dataset.

We constructed 2 sets of predictors—phenotypic and genetic—and combined their predicted probabilities via ensembling to produce the final model. The phenotypic set included variables collected during pretreatment screening and from national registers. To enable baseline prediction, we restricted variables to those available before treatment initiation. First, to reduce the initial predictor space, we performed partial a priori selection of candidate predictors informed by empirical evidence and domain expertise [29,41]. Upon manually screening over 1000 variables, we retained 56 phenotypic predictors (Multimedia Appendix 2), which were subsequently refined using elastic net variable selection. For genetic predictors, we employed the approach suggested by Albiñana et al [42], who leveraged high pleiotropy of complex traits and overlapping underlying mechanisms to demonstrate that incorporating hundreds of PGSs for psychiatric, neurological, and behavioral phenotypes can increase polygenic predictions compared to single-phenotype PGS. Following this methodology, we constructed a library of 615 traits, comprising 15 PGSs from the latest disorder-specific psychiatric GWAS (Multimedia Appendix 3) and 600 PGSs from the UK Biobank GWAS Atlas that covers sociodemographic, lifestyle, and health-related traits [16]. Genotyping, imputation, and quality-control procedures are described elsewhere [43]. PGSs were constructed using PRS-CS, a Bayesian polygenic prediction method that places continuous shrinkage priors on SNP effect sizes and improves performance over traditional approaches [44].

Symptom severity was estimated with the Montgomery-Åsberg Depression Rating Scale, Self-Rating Version (MADRS-S) for depression [45], Panic Disorder Severity Scale, Self-Report (PDSS-SR) for PD [46], and Liebowitz Social Anxiety Scale, Self-Report (LSAS-SR) for SAD [47]. The primary outcome was a binary Clinically meaningful improvement, coded as 1 if either of the following criteria was met: (1) symptom reduction from pretreatment to posttreatment (≥40% for depression and PD and ≥30% for social anxiety) or (2) remission (MADRS-S≤10, PDSS-SR≤7, and LSAS-SR≤35). The composite outcome was chosen to capture 2 clinically relevant pathways to a favorable treatment outcome: symptom reduction of a magnitude deemed clinically relevant or low posttreatment symptom burden. A remission-only criterion may classify patients as nonimproved despite substantial symptom reduction if they remain above the cutoff level, whereas a large relative change may result in patients being classified as improved despite residual symptoms. Moreover, the remission component of the composite outcome is inherently influenced by baseline severity, as patients with lower pretreatment scores require a smaller absolute reduction to meet a remission cutoff, whereas higher baseline severity may make remission harder despite substantial relative improvement. The composite definition was therefore intended to partially mitigate limitations of both end points and reflect a clinically meaningful positive outcome either by reaching an acceptable end-state symptom level or by improving enough that the treatment outcome can be judged clinically beneficial. The primary binary outcome was constructed from original disorder-specific scales: raw pretreatment and posttreatment scores for relative symptom reduction and raw-score cutoffs for remission. Disorder-specific cutoffs for symptom reduction and remission were applied to account for differences in clinically significant symptom reduction across disorders, following the SibeR (Swedish National Quality Register for internet-based psychological treatment) guidelines and existing literature [48-50]. A binary outcome was chosen because it is more directly interpretable and actionable for clinicians than a continuous posttreatment score.

As a secondary outcome, we analyzed the continuous posttreatment score to reduce the information loss inherent in dichotomization. To enable pooling of the 3 disorder-specific symptom scales, MADRS-S, PDSS-SR, and LSAS-SR scores were linearly rescaled to a 0 to 100 range and combined into a single harmonized symptom severity score.

Data preprocessing was performed in R (version 4.3.1; R Foundation for Statistical Computing) and Python (version 3.10.12). Model training and validation were performed in Python, primarily using scikit-learn [51] (1.5.2), NumPy [52] (2.0.0), and pandas [53] (2.2.3) libraries.

To mitigate the risk of overfitting, expert knowledge was applied to filter the original GWAS Atlas of 600 PGSs and preselect those with higher relevance for the ICBT outcome. While based on face validity and theory, this selection was intentionally liberal to allow unexpected associations to emerge. To address collinearity, predictor pairs with a variance inflation factor >5 were excluded. The resulting final set contained a total of 293 PGSs. Following the methodology in Albiñana et al [42], upon standardizing the PGSs, we trained a linear (LASSO‐penalized regression) and a nonlinear model (XGBoost). The hyperparameter tuning procedure mirrored the phenotypic models. Sex, age, and the first 20 principal components were included to control for population stratification, and models with and without covariates were compared by the increment in AUC (ΔAUC). Finally, the model with the best performance across the 10 outer folds of the nested CV was retrained on the entire training set and evaluated on the holdout test set.

Phenotypic and genetic predictor spaces differ substantially in their structure, dimensionality, and signal-to-noise ratio, requiring separate optimization pipelines. Therefore, we used a late fusion strategy to combine predictions from the 2 data modalities. Specifically, out-of-fold predicted probabilities from the best-performing phenotypic and genetic models were combined in an LR stacking meta-learner. This approach was selected because it allows the meta-learner to downweight a constituent model that contributes little independent predictive signal and provides a parsimonious way to estimate the relative contribution of each predictor type. It also enables a direct evaluation of whether genetic predictions add incremental value beyond phenotypic predictions.

Key clinical and demographic sample characteristics are summarized in Table 1, with detailed descriptive statistics and comparisons between the training and test samples provided in Multimedia Appendix 6. A total of 1115 (62.3%) patients achieved clinically meaningful improvement, of whom 753 (67.5%) met both the symptom reduction and remission criteria, 258 (23.1%) achieved symptom reduction only, and 104 (9.3%) achieved remission only, illustrating that the 2 criteria capture partially distinct clinically relevant outcomes.

The training set included patients treated from 2008 to 2017, whereas the temporal holdout test set included the most recently treated patients from 2017 to 2019. Across disorders, the training set contained a greater percentage of improved patients than the test set (63.7% vs 56.7%), and both pretreatment and posttreatment scores were higher in the test set, suggesting a temporal shift toward greater severity.

In nested CV, models with elastic net regularization outperformed those using internal variable selection across all methods except RF, where the omission of redundant predictors may have negatively affected tree splits by reducing ensemble diversity. We therefore proceeded with elastic net-regularized LR, XGBoost, and SVM and a nonregularized RF.

After assessing the calibration of the final models retrained on the entire training set, all models were recalibrated using a 5-fold CV. The recalibration method was selected based on the visual inspection of training set out-of-fold reliability diagrams: Platt scaling (sigmoid) was used when miscalibration showed relatively smooth deviations (benchmark model, LR, and SVM), and isotonic regression was used when deviations were less regular (RF, XGBoost, and both ensembles). The recalibrated models were then evaluated once in the temporally separated holdout test set. Consistent with the lower improvement rate observed in the holdout test set, all models showed some overprediction of the positive class, as reflected by mean predicted probabilities above the observed event rate and negative calibration-in-the-large values. Calibration metrics and curves are provided in Multimedia Appendix 7.

Nested CV and holdout test set performance are shown in Table 2, and receiver operating characteristic curves are presented in Figure 2. The benchmark model yielded the lowest performance (AUC_test 0.695, 95% CI 0.637-0.748). All full phenotypic models showed comparable performance with overlapping CIs, ranging from AUC_test 0.732 (95% CI 0.679-0.784) for LR to 0.749 (95% CI 0.698-0.797) for RF, followed closely by both ensembles. All 21 pairwise DeLong comparisons among the 7 models showed broadly similar performance, except for the benchmark model, which had lower discrimination than RF (P=.04), soft voting (P=.02), and stacking ensemble (P=.03; Multimedia Appendix 8). For threshold-dependent metrics, model-specific decision thresholds were selected on the training set using 10-fold out-of-fold predicted probabilities to maximize the Youden’s J index and were then applied unchanged to the holdout test set. Confusion matrices are shown in Multimedia Appendix 9. Secondary metrics (F₁-score, log loss, and Matthews correlation coefficient) and per-class operating characteristics (sensitivity, specificity, and positive and negative predictive values) are provided in Multimedia Appendix 10.

Genetic models were trained on the full PGS set (615 scores), the preselected reduced PGS set (293 scores), and the psychiatric disorder subset (15 scores). The reduced PGS set performed best across algorithms and was used thereafter. Comparing the difference in AUC (ΔAUC) between PGS-only and PGS + covariates (sex, age, and the first 20 principal components) showed no significant difference; therefore, covariates were excluded from the final models in the interest of parsimony.

Nested CV performance on the reduced multi-PGS set was close to chance for both algorithms. XGBoost performed slightly better (AUC_cv 0.534, 95% CI 0.516-0.552) than LASSO-penalized LR (AUC_cv 0.508, 95% CI 0.480-0.537) and was therefore selected for holdout evaluation, but it did not generalize to the holdout test set (AUC_test 0.512, 95% CI 0.452-0.573), indicating no meaningful predictive signal in this cohort (Multimedia Appendix 14).

The best base learners selected from the 10-fold nested CV were a soft-voting ensemble trained on phenotypic data (AUC_cv 0.703, 95% CI 0.665-0.741) and an XGBoost model trained on multi-PGSs (AUC_cv 0.534, 95% CI 0.516-0.552). The LR meta-learner combining their out-of-sample predicted probabilities yielded an AUC_test of 0.747 (95% CI 0.694-0.797; Multimedia Appendix 15). The DeLong test (P=.97) provided no evidence of improved discrimination by combining genetic and phenotypic predictions over the phenotypic model alone.

Performance metrics of the regressors predicting the continuous posttreatment score on the 0 to 100 outcome scale are presented in Multimedia Appendix 18. All algorithms showed comparable performance; the soft-voting ensemble achieved the lowest RMSE_test of 14.27 (95% CI 12.79-15.81), MAE_test of 10.88 (95% CI 9.92-11.91), and R² of 0.38 (95% CI 0.30-0.46). All algorithms yielded RMSE_test values between 14.27 and 14.52, each below the test set SD (20.1), indicating consistently reduced error relative to a mean-only baseline.

We incorporated multimodal data to predict clinically meaningful improvement following ICBT for depression and anxiety disorders. Briefly, we trained LR, SVM, RF, XGBoost, soft voting, and stacking ensemble models on pretreatment clinical, sociodemographic, and genetic predictors from a large routine-care cohort (N=1790), using MI, elastic net variable selection, nested CV, and a temporal holdout test set. All full phenotypic models achieved moderate predictive performance, and consistent with prior results [28,41], no algorithm was clearly superior. RF and both ensembles surpassed outcome predictions made by therapists at the same clinic [27] and showed higher discrimination than the benchmark model (unadjusted DeLong tests: P=.04, P=.02, and P=.03, respectively). The multi-PGS model showed no predictive power (AUC_test 0.512, 95% CI 0.452-0.573), and the logistic meta-learner ensemble of the best phenotypic and multi-PGS models provided no evidence of a meaningful gain from PGS (P=.97). The holdout performance of all full phenotypic models (AUC _test 0.732‐0.749) compares favorably to other adequately powered and methodologically sound studies that predict CBT outcomes for depression and anxiety using baseline data only (AUCs 0.57‐0.71) [11,64-66]. Our results are thus promising and provide a strong foundation for a future prospective trial to ascertain the model’s real-world utility in the corresponding ICBT context.

Most high-value predictors were those collected during clinical screening, highlighting the utility of easily obtainable data. The screening-based benchmark nevertheless showed lower discrimination than several full phenotypic models, although this comparison cannot fully isolate the contribution of register-based predictors from differences in modeling procedures. Further research is thus needed to determine whether self-reported equivalents of informative register-based predictors (eg, prior psychopathology, medication use, income, and unemployment) can serve as suitable proxies in settings where register linkage is not feasible. PD diagnosis was the strongest phenotypic predictor, reflecting its high response rate (81%) in this transdiagnostic sample. Moreover, the PD PGS was the strongest genetic predictor and the only PGS to survive FDR correction, likely acting as a weak genetic proxy for the clinical phenotype. This finding raised the question of whether model performance was primarily driven by predicting the subgroup for whom ICBT is the most effective treatment. To investigate this, we conducted a post hoc sensitivity analysis, retraining all full phenotypic models on a dataset excluding patients with PD (n=1321; Multimedia Appendix 19). Nested CV discrimination was equivalent for LR and soft-voting ensemble; thus, LR was taken forward to the holdout evaluation on the grounds of parsimony and interpretability, where it retained moderate discrimination (AUC_test 0.693, 95% CI 0.628-0.755). The modest performance drop compared to models trained on the all-disorder dataset is similar to findings from analogous analyses [66] and can be partially attributed to the smaller sample size. This result indicates that model performance was not solely driven by the PD subgroup. To isolate the predictive signal beyond diagnosis more directly, we additionally compared the full phenotypic models with a diagnosis-only model (diagnosis as the sole predictor) in the holdout test set. The full phenotypic models discriminated substantially better than the diagnosis-only model (AUC_test 0.625, 95% CI 0.569-0.679, DeLong P<.001), indicating that they captured individual-level predictive signal beyond between-disorder differences in ICBT outcomes. In practice, this suggests that the models may be most useful as an intake-stage risk stratification tool combining diagnosis with other baseline predictors.

PGSs did not add predictive value in the present cohort and the modeling framework. However, these results should not be interpreted as definitive evidence that genetic data are generally uninformative for psychotherapy outcome prediction. Rather, several characteristics of the current setup likely confined the detectable genetic signal. The multi-PGS framework we adopted, inspired by promising results from Albiñana et al [42], did not prove effective in the present setup, possibly due to several important differences between the studies: the authors trained their models on a much larger cohort (n>100,000), used them to predict psychiatric diagnoses with high established heritability, and included behavioral PGSs that are strongly genetically correlated with those diagnoses. Another important consideration is that our sample contains only mild-to-moderate cases with a rather high level of everyday functioning for whom ICBT is an appropriate treatment modality, whereas more severe patients are referred to specialized outpatient care. This possibly implies less enriched phenotypes and consequently, a weaker genetic signal due to an established correlation between genetic burden and disorder severity [67,68]. The complexity and polygenicity of the studied phenomenon may also play a role in the challenge to capture its genetic underpinnings. Psychotherapeutic treatment outcome is a composite, multifactorial phenotype in which genetic effects are propagated through multiple mediators related to personality traits, learning, executive function, cognition, and other factors, likely translating into diffuse associations of small effect sizes. Finally, treatment outcomes are strongly affected by environmental factors, such as adherence, homework completion, alliance, therapist effects, and idiosyncratic life events, that affect the course of therapy and cannot be accounted for in baseline predictions. Consistent with this, previous work in the same clinical setting has shown that models incorporating predictors collected during ICBT achieve substantially higher predictive accuracy [41].

In clinical practice, the present models could potentially be used for risk stratification at intake by helping to identify patients at elevated risk of nonimprovement who may benefit from additional clinical attention. Plausible downstream actions may include closer monitoring, intensified therapist support, and blended or face-to-face CBT, with the ultimate aspiration to match patients with the most suitable treatment modalities. Such adaptive treatment strategies are particularly relevant in services where ICBT is delivered at scale and clinical resources need to be allocated efficiently. Importantly, while the present models prognostically predict the probability of improvement following ICBT, they do not establish that a patient predicted not to improve would benefit more from a different intervention. Matching patients to alternative treatments is therefore a future possibility that would require prospective evaluation.

A particular advantage of baseline prediction is that it may support the matching of patients to an appropriate level of treatment intensity before treatment initiation. This could help mitigate a drawback of stepped-care approaches informed by within-treatment prediction, in which a patient first undergoes a standard intervention. If this initial treatment format is not beneficial, it may prolong patient suffering, use health care resources inefficiently, and lead to patient dissatisfaction, disengagement, and delayed or reduced future help-seeking.

Temporal validation provides a useful first test of forward-looking performance within the same clinical setting, as the holdout test set consisted of more recently treated patients with higher symptom severity and a higher proportion of nonimprovement than the earlier patients. This severity drift is consistent with changes in referral and admission patterns: as ICBT became more established as an effective treatment modality for common mental health disorders, the clinic accepted a broader and more severe patient population. The models retained moderate performance in this later and more severe sample, supporting their potential relevance for near-future patients in the same clinic. Nevertheless, calibration results indicated some overall overprediction of improvement in the holdout set, consistent with the observed temporal shift in outcome prevalence. Prospective deployment should therefore include calibration monitoring and model updating to reflect current outcome prevalence.

Importantly, statistical prediction is not equivalent to clinical utility. Although the models showed moderate discrimination, their contribution to improved clinical decisions and patient outcomes will require prospective evaluation. The planning of a validation trial in new patients at the same clinic is currently underway, which will test the feasibility of baseline risk prediction under contemporary clinical conditions. Conditional on successful validation, as the next step, a prospective decision-support study will evaluate the impact of early risk stratification and model-informed care on treatment planning and patient outcomes. Finally, external validation in independent ICBT settings in other geographic locations will be required to assess model portability.

For the prospective validation trial, the soft-voting ensemble would be the primary candidate model because it achieved the highest discrimination in nested CV, while showing statistically indistinguishable holdout performance from other full phenotypic models. Before the evaluation in new incoming patients, the model will be retrained on the full available retrospective dataset.

This study has several methodological strengths. First, we used a real-world clinical cohort with high-quality data collection procedures and leveraged national register linkages with excellent coverage. This allowed for the integration of multimodal data (clinical, sociodemographic, and genetic) to build a uniquely comprehensive set of predictors. The PGSs were constructed using large GWAS discovery sets, mitigating the risk that relevant SNPs go undetected in underpowered GWAS. Our models were built on a sufficiently large sample that meets criteria proposed in the literature [69-71] and exceeds the sample sizes of most prior studies [28-31]. We adhered to published guidelines on best practices of predictive model development and reporting [28,33,34]. To derive unbiased performance estimates, we used a nested CV procedure, strictly separating hyperparameter tuning from performance evaluation, and robustly handled missing data using MI nested within the CV folds. Finally, performance was validated using a temporally separated test set, providing a robust real-world estimate of generalizability to future patients.

There are several limitations relevant to the interpretation of this study. First, while efforts were made to ensure model generalizability, the current sample is somewhat clinically and sociodemographically skewed, which may hinder portability to other settings. Specifically, the ICBT clinic from which the sample was drawn is aimed at patients with mild-to-moderate symptom severity, the cohort is on average more educated than the general patient population, may have relatively higher everyday functioning, and be more motivated to undergo treatment because the majority are self-referred. These selection and self-selection biases should be considered when validating the models in new contexts. Moreover, while temporal validation evaluates robustness against data drift, external validation in independent samples will be required to ensure model portability. Second, disorder pooling introduces certain heterogeneity. However, our exploratory comparison showed that combined-data models outperformed single-disorder models, corroborating prior reports [41,72]. This suggests that the added heterogeneity may improve generalizability and mitigate overfitting, in addition to reducing maintenance costs. An additional consideration pertaining to disorder pooling is different symptom reduction thresholds used for operationalization of clinically meaningful improvement (30% for SAD and 40% for depression and PD). Symptom reduction thresholds were defined using disorder-specific criteria suggested by the literature and reflect differences in meaningful symptom reduction across instruments [48-50]. This asymmetry can in principle influence behavior of the transdiagnostic model. However, it is unlikely to account for the main diagnosis-related signal in the present pooled model because SAD still showed the lowest improvement rate and was not a prominent predictor. Instead, the lower relative change threshold for social anxiety may have attenuated between-disorder differences by making disorder subgroups more comparable. Third, the exclusion of observations with missing outcome (~20%) may introduce bias if missingness is not completely random. However, outcome imputation was not feasible as it would cause data leakage during resampling through the necessary combination of predictors with the outcome. Fourth, the reliable quantification of clinically meaningful improvement is challenging, with potential for measurement error. In the absence of objective biomarkers, improvement is a latent construct that cannot be measured directly and is thus operationalized via observable proxies (eg, psychometric scales), which are susceptible to subjectivity and stochastic noise. Construct validity is further compromised by outcome dichotomization. Finally, while substantial for a phenotypic model, the current sample is likely underpowered for a multi-PGS model. Furthermore, an inherent limitation of PGSs is that they are constructed from GWAS of common genetic variants, conferring a risk of missing potentially relevant rare variants. Importantly, the field still lacks a specific PGS for treatment outcomes.

The present baseline predictive models of clinically meaningful improvement after ICBT for depression and anxiety disorders achieved moderate performance. Models incorporating register data generally showed higher discrimination than the benchmark model based on screening data, whereas PGSs added no predictive gain in this cohort and modeling setup. Our models may facilitate early identification of patients at risk of not achieving clinically meaningful improvement following ICBT, thereby informing more targeted clinical decision-making. Next steps include a prospective validation and impact evaluation in new patients at the same facility with subsequent external validation to assess the benefit of model-informed personalized care in other settings.