The literature search was concluded in January 2024. We searched five databases for eligible published studies (ISI Web of Science, PubMed, PsychINFO, CINAHL, and Scopus). Furthermore, the Open Access Dissertation and Theses database (oadt.org) was searched for grey literature. We used the following search string: (“smok* OR “nicotine”) AND (“suicide*” OR “suicide attempt”) AND (“psychosis” OR “psychotic” OR “bipolar” OR “schizophren*”). In August 2024, the literature search was updated and the search string was extended to include the search term “tobacco use.”

To be included in the present meta-analysis, studies had to meet seven inclusion criteria. They had to (i) be observational, (ii) provide an odds ratio or sufficient statistical parameters to calculate one for smoking and suicide attempts, (iii) include lifetime or past year reports of suicide attempts, (iv) report the respective sample size, (v) report data of people with psychotic disorders (i.e., schizophrenia, schizoaffective disorder, first-episode psychosis, delusional disorder, bipolar disorder, or psychotic depression), (vi) comprise adult participants (mean sample age = 18 years+), and (vii) be published in English. In all, 3,603 study titles and abstracts were screened after duplicate removal (see Figure 1 for a PRISMA flow chart). Study coding was conducted twice independently by one researcher. Discrepancies were resolved by discussion with another independent and experienced researcher. All data are provided on the OSF. ¹

We coded event vs. non-event data (i.e., numbers of smokers and non-smokers who did or did not attempt suicide) into categories and recorded study-specific variables: Study type (cohort vs. cross-sectional vs. case–control); sample type (inpatients vs. outpatients vs. mixed); diagnosis type (schizophrenia spectrum vs. bipolar spectrum vs. mixed; of note, a further differentiation according to type of diagnosis was not feasible due to low study cell counts in terms of reported diagnoses); sample size; country of publication. In cases of missing information, data were requested from corresponding authors by email. If no responses were received, individual entries were coded as missing.

We used the Newcastle-Ottawa Quality Assessment Scale (Wells et al., 2021) to assess study quality regarding risk of bias assessments according to an adapted approach for cross-sectional studies (Dürlinger and Pietschnig, 2022). Primary study quality was assessed by an experienced researcher [J.P.] and rated according to four criteria (ratings are provided at https://osf.io/brzs4). Rating sum scores were subsequently correlated with primary study effect sizes to assess potential influences of study quality on effect size estimation. No statistically significant influences of primary study quality ratings on effect estimation were identified in the present study in a precision-weighted meta-regression (Q = 1.78; p = 0.182).

We used random-effects models to calculate summary effects. For our calculations, odds ratios were transformed into log odds ratios and, after analysis, back-transformed for ease of interpretation. All analyses were conducted with the open source software R 4.10 (R Core Team, 2023) [packages: metafor (Viechtbauer, 2010), data.table (Dowle et al., 2019), esc (Lüdecke et al., 2019), puniform (van Aert et al., 2020), readxl (Wickham et al., 2019), tidyverse (Wickham and Wickham, 2017), and weightr (Coburn et al., 2019)]. Between-studies heterogeneity was interpreted according to well-established thresholds of the I²-index [i.e., 25, 50, and 75% representing lower thresholds of small, moderate, and large heterogeneity, respectively; (Higgins et al., 2003)]. Following well-established benchmarks, Odds Ratios of 1.44, 2.48, and 4.27 (corresponding to logOR = 0.36, 0.91, and 1.45, respectively) were interpreted as the lower thresholds of small, medium, and large effects (Cohen, 1988).

We used mixed-effects subgroup analyses to assess potential influences of categorical moderators, specifically, diagnosis type (bipolar vs. schizophrenia vs. mixed), sample type (inpatients vs. outpatients vs. mixed), and study type (case–control vs. cohort vs. cross-sectional). We used precision-weighted mixed-effects meta-regressions to assess the potential effects of the continuous variables sex (percentage of women within samples) and publication year.

To detect potential influences of confounding dissemination bias, we used several different bias detection approaches (i.e., visual funnel plot inspection, formal funnel plot asymmetry-based assessments, selection model assessments, and p-value-based methods), following current recommendations (Siegel et al., 2022). First, we visually inspected funnel plots for potential asymmetry. Second, Sterne and Egger’s regression test (Egger et al., 1997) and Trim-and-Fill (Duval and Tweedie, 2000) were applied. The Sterne and Egger regression regresses effect standard deviates (i.e., standardized effects) on study precision. A significant intercept can be interpreted as evidence of publication bias. In Trim-and-Fill, funnel plot asymmetry is assessed and excessive effects in terms of strength according to the estimated summary effect are trimmed. Subsequently, the summary effect is re-estimated according to the trimmed data set. This iterative procedure is repeated until no further asymmetry is detected. Then, trimmed studies are reinserted, and presumably missing effects are symmetrically imputed on the opposite side of the funnel. Notably, the resulting adjusted effect size should not be seen as a corrected estimate but should be rather interpreted in the sense of a sensitivity analysis.

Third, we examined the effects of a moderate one-tailed selection based on Vevea and Woods (2005) selection model approach. In this method, effect sizes are weighted according to the publication likelihood of specific assumed p-values (we followed a standard weighting scenario as suggested by the authors). Their distribution is then compared to an unweighted model. The resulting effect estimates should be similar in the absence of publication bias.

Finally, we used p-uniform (Van Assen et al., 2015) and p-curve (Simonsohn et al., 2014) to assess potential bias. Both methods are based on the idea that study p-values are uniformly distributed in the presence of a null effect. With increasing non-null effect strength, the p-value distribution becomes increasingly right-skewed. An identical phenomenon should be observable when focusing exclusively on p-values that fall below the typically assumed significance threshold of 0.05. In p-curve, formal tests allow an assessment (i) of the evidential value of a given set of study effects as well as (ii) effect-distorting effects of p-hacking (i.e., extensive use of questionable research practices) based on the observed p-value distributions.

In p-uniform, conditional distributions of significant p-values based on the corresponding population effect size are compared with a uniform distribution to assess evidence for p-hacking. Moreover, this approach allows a summary effect and confidence interval estimation based only on the observed significant p-values and their associated degrees of freedom. Following current recommendations, we interpreted p-values <0.10 or differences between estimated effects exceeding 20% as indicative of publication bias in all detection methods (Siegel et al., 2022).

We included 22 studies comprising k = 27 independent effect sizes (N = 11,452 patients; mean sample age = 39.5 years). Most participants were schizophrenic spectrum (n = 4,797) and bipolar spectrum patients (n = 4,020), with the remaining samples comprising both diagnosis types (n = 2,635). The majority of samples comprised outpatients (n = 6,269), with the remaining samples consisting of inpatients (n = 2,566) and mixed patient groups (n = 2,590). Primary study characteristics are detailed in Table 1.

Effect sizes were differentiated according to study type (Q = 12.141, p < 0.001), with cohort studies yielding a moderate effect (OR = 2.76, p < 0.001) that was significantly larger than effects from cross-sectional (OR = 1.58, p < 0.001) and case–control studies (OR = 1.28, p = 0.106).

Groups with different diagnoses showed predominantly small non-trivial effects that did not significantly differ between groups (Q = 0.864, p = 0.353). However, patients with bipolar disorder showed the numerically largest and non-trivial effect (OR = 2.02, p < 0.001), followed by patients with psychosis (OR = 1.64, p < 0.001), and the mixed group (OR = 1.44, p < 0.001).

Similarly, no significant differences between sample types were observed (Q = 1.614, p = 0.204), yielding mostly non-trivial but small effects. The outpatient group showed numerically somewhat larger effects (OR = 1.92, p < 0.001) compared to the inpatient (OR = 1.67, p = 0.008) and the mixed groups (OR = 1.42, p = 0.002). Numerical results of subgroup summary effects are detailed in Table 2.

A precision-weighted meta-regression of publication year on effect sizes indicated significant decreases in effect strength over time (b = −0.097, Q = 17.61, p < 0.001; Figure 3). Cumulative meta-analyses supported this finding, indicating almost exclusively continuously decreasing summary effect size estimates when individual study effects were added in turn according to study publication years (Figure 4).

In another meta-regression, the percentage of women within samples showed no significant influences on effect sizes (b = 0.004, Q = 1.05, p = 0.305; see Figure 5).

Visual inspection of the funnel plot showed substantial asymmetry, suggesting inflation of the observed summary effect (Figure 6). Formal tests in terms of Sterne and Egger’s regression method yielded significant results, thus supporting this interpretation (Z = 3.81, p < 0.001). Results from Trim-and-Fill were consistent with findings from Sterne and Egger’s regression, indicating nine missing effect sizes on the left side of the funnel plot and yielding an adjusted effect of OR = 1.51 (LogOR = 0.42).

In contrast, the selection model approach did not show evidence for publication bias, indicating an adjusted summary effect that remained within 20% of the value of the observed effect (i.e., OR = 1.69 vs. 1.70, respectively).

There was no significant evidence for p-hacking according to p-curve analysis in binomial and continuous tests. Moreover, both binomial and continuous tests indicated evidence for an appropriate evidential value of our summary effect calculations (see Figure 7). Similarly, p-uniform did not yield evidence for p-hacking either (p = 0.989). p-uniform-based summary effect estimations showed a non-trivial effect, yielding an OR = 1.63 (95% CI [1.52; 1.75], p < 0.001).

Some limitations of this meta-analysis need to be acknowledged. First, due to the correlational design of the present study, it was not possible to empirically establish a causal mechanism in regard to the smoking status and suicide attempt link of individuals with psychosis and bipolar disorder. However, the observed non-trivial link provides a useful means for practical care recommendations. Second, some statistical noise may have been introduced by different design characteristics between primary studies that could not be entirely accounted for in our present meta-analysis (e.g., in terms of different smoking operationalizations or psychosis and bipolar disorder assessments). However, this is common in meta-analytical investigations and was accounted for by modeling between-studies heterogeneity in terms of random-effects calculations. Finally, self-reports of lifetime suicide attempts are susceptible to misclassification (Castelein et al., 2015) and prone to survival bias (only reports by survivors). Nonetheless, suicide attempts have been frequently been observed to predict suicide better than suicide ideation (Millner et al., 2015; May and Klonsky, 2016; Ribeiro et al., 2016). There is no reason to suspect that any of these points may have introduced any systematic influences on summary effect estimations, thus providing evidence for the salience of our observed summary effect.