Analyses

Statistical Models All analyses were conducted in R (RStudio; R v4.5.2) using an intention-to-analyse approach on available data, with two-tailed tests and α = .05. Models were estimated using complete-case data per analysis. Unstandardised coefficients (B), standard errors (SE), standardised coefficients (β), test statistics, and 95% confidence intervals (CIs) are reported. Partial η² was derived for ANCOVA-style models where appropriate. H1 (Treatment → depression) Baseline-adjusted ANCOVA-style linear regression models were fitted separately for each follow-up (V6, V7), with change in MADRS from baseline (V2) as the dependent variable. Treatment allocation (psilocybin vs placebo) was entered as the primary predictor, with baseline MADRS included as a covariate: ΔMADRS_wk = β₀ + β₁(arm) + β₂(MADRS_baseline) + ε Model fit was summarised using F-tests and adjusted R². Partial η² was derived to quantify the incremental variance explained by treatment. H2 (OCB → depression change) To test whether acute Oceanic Boundlessness (OCB; V4) predicted depressive symptom change, two complementary approaches were used: Repeated-measures ANCOVA examining whether OCB moderated trajectories of MADRS across V2, V6, and V7 (Time × OCB interaction). Baseline-adjusted ANCOVA-style regressions at each follow-up, with MADRS change as the outcome and OCB, treatment arm, and baseline MADRS entered as predictors: ΔMADRS_wk = β₀ + β₁(OCB) + β₂(arm) + β₃(MADRS_baseline) + ε Partial η² was calculated for OCB where appropriate. H3 (Specific OCB subscales → depression change) Multiple linear regression models were conducted at V6 and V7 with MADRS change as the outcome. The five OCB subscales (Unity, Spiritual Experience, Bliss, Insight, Disembodiment) were entered simultaneously to estimate their unique contributions while accounting for shared variance, alongside treatment arm and baseline MADRS: ΔMADRS_wk = β₀ + Σβᵢ(ASC_subscaleᵢ) + β₆(arm) + β₇(MADRS_baseline) + ε Standardised coefficients and 95% CIs were computed to compare relative effect sizes. This multivariable approach was used to address collinearity and isolate unique experiential predictors. H4 (Treatment → connectedness to nature) Connectedness to Nature Scale (CNS) scores at V2, V6, and V7 were analysed using linear mixed-effects modelling (random intercept for participant). Fixed effects included time, treatment arm, and their interaction: CNS_it = γ₀ + γ₁(time) + γ₂(arm) + γ₃(time × arm) + uᵢ + ε_it Models were estimated using lme4 with Satterthwaite degrees of freedom via lmerTest. The primary inferential test was the time × arm interaction. Estimated marginal means and model-predicted trajectories were examined for interpretation. H5 (Mediation: Treatment → OCB → connectedness to nature) Mediation analyses were conducted using lavaan with nonparametric bootstrapping (5,000 resamples). Separate models were fitted for CNS change at V6 and V7 and for candidate mediators (e.g., OCB total, Bliss, Spiritual Experience). Model structure: a-path: mediator ~ treatment arm b- and c′-paths: CNS change ~ mediator + treatment arm + baseline CNS Indirect effect: a × b Total effect: c′ + (a × b) Baseline CNS was included as a covariate (ANCOVA principle). Indirect, direct, and total effects were reported with bootstrapped 95% CIs and standardised estimates. Effect Size Effect sizes included partial η² (ANCOVA-style models), standardised β coefficients (regressions), and standardised indirect effects (mediation models), all accompanied by 95% CIs. Statistical Power Data were derived from a feasibility phase 2 trial and were not powered to detect a pre-specified effect size. Findings are interpreted with appropriate caution. Inference Criteria Statistical significance was inferred where 95% CIs did not include zero and/or p < .05. Greater emphasis was placed on effect sizes, confidence intervals, and consistency across models rather than single p-values. Assumption Checks Model assumptions (normality, homoscedasticity, influential observations) were assessed visually and statistically. Where necessary, alternative specifications or robust interpretations were considered and documented. Reliability and Robustness Testing Sensitivity analyses examined alternative operationalisations of symptom change (e.g., raw vs percentage change). Exploratory robustness models additionally tested whether OCB effects remained significant after controlling for: Global ASC intensity ASC-DED (challenging/dread experiences) Extended multivariable models were used to determine whether specific subscales (e.g., Spiritual Experience) retained unique predictive effects after accounting for broader experiential intensity. These analyses were treated as hypothesis-generating. Exploratory Analyses Exploratory analyses examined: • Whether OCB differentially predicted short-term (V6) versus longer-term (V7) changes in depression or connectedness. • Associations between changes in connectedness to nature (V2 → V6/V7) and changes in depressive symptoms. All exploratory analyses were clearly labelled as such.