Can Using Psychedelics Help Meditation? Looking at Personal, Drug, and Meditation Factors

The primary outcome variable was a single-item question that asked participants: “Overall, do you believe that your psychedelic experience(s) have influenced the quality of your regular meditation practice?” Using a 7-point Likert scale, responses ranged from 1 (strong negative influence) to 7 (strong positive influence). All survey items included are available at the Open Science Framework (https://osf.io/56utj/↗).

Individual factors: Self-reported demographic variables included participant age, gender, and education (in years of schooling). Two life satisfaction questions (“Overall, how satisfied are you with your life nowadays?” and “Overall, to what extent do you feel that the things you do in your life are worthwhile?”) were drawn from questions asked as part of the UK Census [73]. Using an 11-point Likert scale, responses ranged from 0 (not at all) to 10 (completely).

Personality traits were measured using the Big Five Inventory (BFI; [74]), which captures the dimensions of openness to experience, conscientiousness, extraversion, agreeableness, and neuroticism using a 5-point Likert scale ranging from 1 (disagree strongly) to 5 (agree strongly). In this study, Cronbach’s α ranged from 0.74 to 0.87 for these five subscales.

We also captured alcohol, tobacco, and cannabis use. Alcohol use was assessed using an adapted item from screening questions used by the National Institute on Alcohol Abuse and Alcoholism (NIAAA, 2011). Participants were asked about the use of alcohol (“During the last 12 months, how often did you usually have any kind of drink containing alcohol?”) and had ten response options ranging from 0 (not at all) to 9 (everyday). Tobacco (“During the last 12 months, how often, on average, have you smoked or vaporized tobacco?”) and cannabis (“During the last 12 months, how often, on average, have you smoked or vaporized cannabis/weed?”) were similarly assessed but with five response options ranging from 0 (not at all) to 4 (everyday).

Psychedelic-related factors: The frequency of psychedelic use (which excluded microdosing, MDMA/ecstasy, and ketamine) was assessed over the last 12 months and over an individual’s lifetime. Participants were also asked their age when they started consuming psychedelic substances.

The types of psychedelic substances that participants had used at any point in their life were captured by providing a list with the following options: psilocybin/magic mushrooms or truffles, LSD or LSD derivatives, ayahuasca, N,N-DMT, 5-MeO-DMT, mescaline, and others (i.e., an open text box in which participants could specify additional psychedelic substances). Given participants could report exposure to multiple psychedelic substances, binary codes (0 =  No, 1 =  Yes) were generated for exposure to a given psychedelic substance.

To assess the context and set and setting in which participants used psychedelics, single-item questions were asked about intention setting (“Over your lifetime, when you used psychedelics how often did you set clear intentions before a psychedelic experience?”), attention to safety (“Over your lifetime, when you used psychedelics how often did you pay attention to ensuring the setting was safe and supportive?”), and being surrounded by trusted individuals (“Over your lifetime, when you used psychedelics how often did you pay attention to ensuring you were with people that you trusted?”). Responses ranged from 0 (never) to 4 (always). Finally, to understand participants’ experience with psychedelic use, participants were asked about their psychedelic framework (“I have an adequate framework or worldview for understanding or making sense of my psychedelic experience(s)”). Responses ranged from 1 (disagree strongly) to 6 (agree strongly).

Meditation-related factors: Frequency of meditation practice (“During the last 12 months, how regularly have you, on average, practiced formal sitting meditation (for at least 20 minutes per session)?”) was assessed with response options ranging from a minimum of 3 (3 times a week) to 8 (twice a day or more) and was included as a continuous variable. Participants were also asked about their number of years of regular meditation practice (i.e., at least three days a week) and their starting age for practice. To understand participants’ experience with meditation practice, participants were asked about the adequacy of their meditation framework (“I have an adequate framework or worldview for understanding or making sense of my meditation experiences.”). Responses ranged from 1 (disagree strongly) to 6 (agree strongly). Participants were also asked whether they had any retreat experience, the number of retreat days in the last twelve months and in their lifetime as well as the length of their longest retreat. Responses were recorded in number of days. Participants were also asked if they first started meditation practice or if they first used a psychedelic.

Finally, to assess meditation practice, participants were invited to qualitatively describe their meditation practice background. Participant responses were then individually coded into 17 practice types by the fourth (masked) and sixth (masked) authors with expertise in modern mindfulness and Buddhist meditation practices. A bottom-up approach was used for coding (inductive derivations of a large number of categories which were then consolidated) while keeping important theoretical considerations in mind for the classification of meditation practices (e.g., [42,61]). This was done in part due to practical considerations given the varying lengths and specificity of responses from participants. Many participants reported studying under specific contemporary teachers (e.g., Culadasa, Rob Burbea, Sam Harris, Shinzen Young), and in cases where these teachers have created meditation systems that draw on multiple traditions and/or sources in innovative ways that could not be consolidated under a single rubric, the meditation style is coded with the name of the contemporary teacher. Individual names of contemporary teachers were not used when the practices could be consolidated under a general rubric (e.g., Western Non-Dual Meditation). In cases where individuals reported practicing in specific traditional lineages (e.g., Burmese Vipassanā), the practice is coded under the tradition’s name. A significant number of participants (n = 193) reported idiosyncratic practices that involved an assorted range of techniques that they themselves invented or were drawn from numerous teachers and traditions. These were coded as “eclectic/other.” The final practice types are listed in the Meditation-Related Factors subsection of Table 1. Response codes ranged from 0 (participant did not mention the practice or teacher) to 1 (participant listed the practice or teacher). A single participant could have reported multiple practice types. For example, one participant responded to the inquiry with “Mindfulness meditation as taught in Sam Harris’ Waking Up app Open awareness practices as taught by Loch Kelly.” In this case, two codes were applied to this response including the teachings of Sam Harris, and Western Non-Dual Meditation. Binary codes (0 = No, 1 = Yes) were created for exposure to each of the 17 categories.

All data, R code, study items, and consent form used in the study are available at the Open Science Framework (https://osf.io/56utj/↗). The analyses were planned and conducted following data collection and were not preregistered. While pre-registration may be preferred, it is worth noting that the use of cross-validation combines exploratory and confirmatory methods. Many model configurations are simultaneously explored in training data, but the performance of the best model is subsequently confirmed in new (held-out) data.

For descriptive purposes, we first examined zero-order correlations between the outcome and the predictor variables using Pearson’s correlation coefficient. We report false discovery rate (FDR; [75]) corrected p-values for these analyses to control type 1 error rates.

Machine learning analyses were conducted using the ‘tidymodels’ ecosystem [76] in R [77]. We used machine learning for three reasons. First, when standard (i.e., ordinary least squares; OLS) linear regression is used with high-dimensional, correlated predictors, it often yields high variance models that overfit the training data such that they do not generalize well and therefore perform poorly with new data [71,72]. Second, the broader class of linear models (e.g., OLS linear regression, ridge, LASSO, and elastic net regressions) may yield biased models that also perform poorly with new data if the true data-generating process for the outcome is non-linear (e.g., non-linear or interactive effects among predictors and outcome; [78]). Third and more generally, when model performance/fit is evaluated with the same data that were used to train the model, the performance will often overestimate how well the model will perform with new data [72]. Machine learning methods address the first two issues by simultaneously considering several statistical algorithms that adopt different approaches to address the problems of model bias and variance (due to overfitting) by optimizing the bias-variance trade-off [72,79]. Machine learning methods address the third issue by evaluating model performance using new (i.e., held-out) data that were not used to train the models.

In this study, we used three statistical algorithms, OLS linear regression, elastic net linear regression, and random forest. We included OLS linear regression because it is a highly interpretable parametric model that is well-understood and frequently used in the social sciences. We included elastic net linear regression because it is another interpretable parametric model, but it uses regularization techniques that penalize the parameter estimates to yield simpler models that may not overfit training data to the same degree as OLS linear regression [80]. We included random forest because it is a flexible, non-parametric model that may outperform linear models when the data-generating process is non-linear [81]. It also is often more robust to overfitting than OLS linear regression. Furthermore, both elastic net and random forest often outperform OLS linear regression when multi-collinearity among predictors is present [82,83].

We fit all models, tuned hyperparameters, and selected the best model configuration using three repeats of 10-fold cross-validation (i.e., 30 held-out folds; [84]). For elastic net, we tuned two hyperparameters, lambda (λ; the penalty parameter) and alpha (⍺; the mixing parameter that dictates the proportion of L1 vs L2 penalties). For random forest, we also tuned two hyperparameters, mtry (the number of predictors randomly selected to consider for each split) and minimum node size (the minimum observations needed for additional splitting). We used 1500 trees for all random forest configurations. OLS linear regression does not allow tuning any hyperparameters. Sensible values for each of the hyperparameters for elastic net and random forest were provided by relevant tidymodels functions (from ‘tune’ package) and hyperparameter plots were reviewed to confirm that an adequate range of values was considered. We used mean R² across the 30 held-out folds to tune these hyperparameters, select the best model configuration, and evaluate the performance of that best model[85].

Several feature engineering steps were undertaken to improve model performance. For all three statistical algorithms, we winsorized the predictors to two standard deviations above and below the mean to minimize the influence of outliers [86]. For all algorithms, we also imputed missing data on the predictors using a k-nearest neighbor (KNN) approach implemented within ‘tidymodels’ using package guidelines and function defaults [76]. KNN identifies the observations that are most similar to the missing values (i.e., its “nearest neighbors”) and uses the mean value across those cases for imputation [87]. We reduced the dimensionality of the predictors somewhat by substituting the first principal component from principal components analyses (PCA) performed with groups of highly correlated predictors that were identified by preliminary exploratory data analysis using the full sample [88]. These groups of variables included: 1) two items examining life satisfaction, 2) 12-month and lifetime psychedelic use, and 3) 12-month and lifetime retreat practice. Finally, we applied the Yeo-Johnson power transformation to all quantitative predictors to normalize their distributions for use in the linear statistical algorithms (OLS linear regression and elastic net; [84,89]).

We planned to use three distinct methods to interpret the best model and identify its important predictors depending on the algorithm that performed best. First, we planned to calculate Shapley Additive Explanations (SHAP; [90]) values to assess feature importance associated with each predictor in the best model. Using a concept drawn from cooperative game theory, SHAP values index the average marginal contribution of a given predictor to the predicted outcome for each observation across all possible combinations of predictors (i.e., local SHAPs [91]). SHAP values are model agnostic (i.e., can be calculated for any statistical algorithm) and have useful properties including additivity (i.e., values for each observation are computed independently and then summed), symmetry (i.e., average values of two features should be equal if they contribute equally to the overall model), and null contribution (i.e., a feature that does not contribute to the model will have a value of 0). Plots of the relationship between predictor values and local SHAP values can aid in the interpretation of the direction of predictor effects. Moreover, the global importance (across all observations) for any predictor can be quantified by averaging the absolute value of the local SHAP values across all observations. Larger global SHAP values suggest greater relative importance of a given predictor to the predictive capacity of the model.

Second, we used a feature ablation procedure to assess the incremental contribution of each predictor to the best model [92]. Specifically, for each predictor we compared two nested models, a full model that contained all predictors and a reduced model that removed (i.e., ablated) the relevant predictor. We used an implementation of Bayesian estimation provided by the ‘tidyposterior’ package [93] to estimate the increase in performance for the full model relative to reduced models using R² from the best 30 held-out folds. Predictors are considered important if the posterior probability that R² is greater for the full vs. reduced model exceeds 0.95. Like SHAP, this feature ablation procedure can be used with any statistical algorithm. We sequentially evaluated predictors (ordered by their global SHAP values) until we reached a predictor that did not have a probability > .95 of increasing R² from the reduced model.

Finally, if a linear model (OLS regression or elastic net) was selected as the best-performing algorithm, we planned to review and report its standardized parameter estimates (i.e., βs) to rank order the magnitude of the effects of its predictors on the outcome. This third method was not possible for random forest because that algorithm is non-parametric.

We used mean R² values across the 30 held-out folds to select the best statistical algorithm and associated hyperparameters. The elastic net linear regression (with λ = 0.19 and α = 0.20) had the highest mean R² (R² = 0.266, SE = 0.013) among the three algorithms. The mean R² was 0.252 (SE = 0.011) and 0.251 (SE = 0.012), respectively, for random forest (with mtry = 5 and minimum node size = 30) and OLS linear regression. Therefore, all subsequent analyses were conducted with the elastic net linear regression model. Across the 30 held-out folds, the root mean squared error for this model was 1.04, and the mean absolute error was 0.848 (with the outcome measured on a 7-point scale).

As described earlier, we used three complementary methods to examine the feature importance of the predictors in the elastic net model: 1) SHAP values, 2) feature ablation, and 3) standardized coefficients. Global SHAP values for the predictors are presented in Fig 1. This method suggested that the most important features were psychedelic factors: psychedelic use (lifetime/12 month) and setting intention during psychedelic use. The global SHAP values associated with these two predictors (0.25 and 0.15, respectively) were 2 – 4 times greater than any other predictors in the model. The next three important variables were agreeableness (0.07), exposure to N,N-DMT (0.07), and longest retreat (0.07).

We next used feature ablation to sequentially evaluate the incremental contribution of each of the top predictors (ordered by Global SHAP values) to the model. See Table 2 for complete results. For each predictor, we estimated the posterior distribution for the increase in R² when that predictor was added to the reduced model that contained all other predictors. From this distribution we calculated the mean increase in R², the 95% credible interval for this increase, and the posterior probability that this increase was > 0 (i.e., the probability that mean R²_full> mean R²_reduced for that predictor among the held-out folds). Psychedelic use was associated with mean increase in R² of.048 (95% CI = [.032−.064]), with a posterior probability > .999 that this predictor increased model performance. Setting intentions during psychedelic use was associated with a mean increase in R² of.019 (95% CI = [.013−.024]), with a posterior probability > .999 that this predictor increased model performance. Agreeableness was associated with a mean increase in R² of 0.006 (95% CI = [.002,.011]) with a posterior probability of.985 that this predictor increased model performance. N,N-DMT was associated with a mean increase in R² of 0.005 (95% CI = [.001,.009]) with a posterior probability of.970 that this predictor increased model performance. Longest retreat was associated with a mean increase in R² of.001 (95% CI = [−.002, −.002]) and the probability that this predictor increased model performance was low (posterior probability = .553). As the posterior probability was > .95, no further predictors were evaluated.

To provide a third perspective on feature importance, we display the standardized coefficients for the predictors from the elastic net linear regression in Fig 2. Consistent with the two prior feature importance methods, the four most important predictors, based on the absolute value of their standardized coefficients, were psychedelic use (lifetime/12 month, standardized estimate β = 0.30), setting intentions during psychedelic use (β = 0.18), agreeableness (β = 0.09) and N,N-DMT (β = 0.08). See Fig 2 for the list of variables deemed most important.