Data were collected via repeated cross-sectional surveys between May 8 and December 1, 2020. A total of 6 waves of data were collected using a web-based panel administered by the research and data collection company Delvinia [14]. The sampling waves occurred from May 8 to 12 (n=1005), May 29 to June 1 (n=1002), June 19 to 23 (n=1005), July 10 to 14 (n=1003), September 18 to 22 (n=1003), and November 27 to December 1 (n=1003). Participants were sampled independently in each wave. The overall response rate was 16.1% (6021/38,987).

Quota sampling based on age, gender, and region was used to obtain a sample that is proportional to the English-speaking population of Canada. Canadians aged ≥18 years were eligible. Respondents provided written informed consent electronically prior to participation. Research ethics approval was obtained from the Centre for Addiction and Mental Health research ethics board.

Survey questions included information on demographics, anxiety, depression, substance use, employment changes, perceived risks, and worries related to COVID-19. A full list of variables considered for analysis is included in Table S1 in Multimedia Appendix 1. Respondents’ anxiety levels were captured using the Generalized Anxiety Disorder–7 questionnaire (GAD-7) [15], a validated inventory measuring the frequency of anxiety symptoms over the past 2 weeks. Raw item scores were used, ranging from 0 (not at all) to 3 (nearly every day). Depressive symptoms were measured using 3 modified questions from the Centre for Epidemiologic Studies Depression Scale (CES-D) [16]. Feelings of depression, loneliness, and hopefulness over the past week were reported on a Likert-style scale between 0 (“Rarely or none of the time [less than 1 day]”) and 3 (“Most or all of the time [5-7 days]”). For further details on the included explanatory variables, see Text S1 in Multimedia Appendix 1.

Given the strong comorbidity of population-level depression and anxiety symptoms [17], as well as their shared neurobiological underpinnings [18], we first examined pairwise correlations among all anxiety questionnaire items (GAD-7) and the 3 available mood/depression (CES-D) questionnaire items, and we applied principal component analysis (PCA). The questions from both scales were similar in scale. Correlations were high between the GAD-7 and CES-D items, allowing for the combination of items from both questionnaires into a single measure of emotional distress using symptoms of both anxiety and depression.

PCA with varimax rotation was used to reduce the number of mood and anxiety variables needed for modeling while retaining as much information as possible from all anxiety and depression variables. PCA is widely used to identify principal axes of variation in psychometric questionnaires [19,20]; in the absence of a validated method of combining items across the GAD-7 and CES-D, PCA was used to retain the largest amount of useful information in a single score incorporating both scales. To properly account for the ordinal nature of the GAD-7 and CES-D variables, polychoric correlations—measuring associations between ordinal variables assumed to be realizations of underlying latent Gaussian distributions [21]—were used. PCA was applied to our outcomes on all observations prior to train and test spitting.

For our data-inclusive approach, all survey questions were considered for inclusion as model predictors, excluding mood and anxiety variables used in our outcome measure. Questions not asked in all 6 survey waves were excluded. Categorical variables were one-hot encoded (1=yes, 0=no for category membership). “Prefer not to answer” responses were treated as missing (see the Predictive Modeling section below).

A total of 6021 respondents provided complete surveys for analysis. The characteristics of the respondents are summarized in Table 1. Demographic distributions of age, sex, and region were representative of the English-speaking Canadian adult population [27,28].

The PCA was initially fit using 10 survey items: 7 from the GAD-7 (anxiety) and 3 from the CES-D (depressive symptoms). The polychoric correlations between the anxiety and depression variables are presented in Figure 1A. All items were moderately to strongly positively correlated (r=0.55-0.89), with the exception of hopefulness (CES-D item; correlation coefficients –0.18 to –0.27). All items except for hopefulness were loaded strongly onto principal component 1 (PC1) (loadings 0.74-0.93), with hopefulness loading weakly onto PC1 (loading –0.28) and strongly onto principal component 2 (PC2) (loading 0.96); this finding supported our choice to model a single latent outcome combining both scales. Given the weak negative correlation between hopefulness and the remaining mood and anxiety variables, hopefulness was dropped from the PCA and instead included as an explanatory variable in downstream modeling. Distributions of the GAD and CES-D scores and PC1 are included in Figure S1 (Multimedia Appendix 1).

The final PCA included 9 items. In the remainder of this analysis, we refer to PC1 as “emotional distress.” Figure 1B displays a scree plot for this 9-item PCA; the point of inflection occurs at the second PC, indicating that the PCs after the first do not add substantial additional information. The loadings on PC1, emotional distress, are presented in Figure 1C. The item loadings ranged from 0.74 to 0.93. Emotional distress explained 76% of the variance in all 9 variables.

A total of 50 predictor variables were included in our primary XGBoost model. We included 4819 respondents in the training data set; out-of-sample prediction was tested on the remaining 1202 respondents. For the results of the hyperparameter selection, see Text S3 in Multimedia Appendix 1. The final model explained 38.7% of the variance in distress in the out-of-sample prediction (r²=0.387). A scatterplot of the predicted and realized values of emotional distress is presented in Figure S2, Multimedia Appendix 1. In training, the r² value of the model was 0.394, which is close to the r² value for our holdout test set; this indicates that there was no substantial overfitting. For a LASSO model fit to the full test set, r² was 0.354 in the withheld validation set, which is slightly lower than the value for XGBoost in predictive accuracy.

A graphical depiction of model training using gradient-boosted trees [29] is presented in Figure 2A. A visual representation of the first 6 gradient-boosted trees in the fitted model is included in Figure S3 (Multimedia Appendix 1). A partial dependence plot overlaid with COVID-19 positivity rates [30] is presented in Figure S4 (Multimedia Appendix 1).

To determine if the model fit was sensitive to the XGBoost imputation algorithm for variables with missing values, XGBoost was refit with reoptimized hyperparameters and validated using only observations with complete data (n=3689). On the withheld data, this model explained 36.9% of the variance in emotional distress (r²=0.369). Results of a LASSO model fit to the same set of complete observations indicated slightly lower performance compared to the XGBoost fit (LASSO r²=0.346), indicating that distress was reasonably approximated by a linear fit, although at a loss of approximately 2% of the explained variance when nonlinear and interaction relationships were not considered.

The SHAP values are presented in Figure 2B. To aid the direct interpretability of our SHAP value analysis, we modeled each questionnaire item as a function of our latent emotional distress outcome using linear regression. Each 1-unit change in our outcome corresponded to a difference in the question response value of between 0.68 and 0.84 (mean 0.78) across all 9 questions, meaning that a SHAP value of 0.61 (the largest value reported for individuals with high levels of reported financial worry) would translate to an average increase of 0.61*0.78=0.48 across the original question scales. Given that each question was measured as integers ranging from 1 to 4—representing an underlying quantitative scale mapping onto the number of recent days when symptoms were experienced—this maximum SHAP value would represent a predicted marginal change of 16% of the entire spectrum of symptom burden across all questions in the holdout test population: 0.48/(4-1)*100.

The variables with the greatest importance, in descending order, were worry about personal finances due to COVID-19 (mean absolute SHAP value=0.168), worry that oneself or loved ones will become ill with COVID-19 (0.162), age group (0.130), hopefulness (0.115), change in alcohol use due to the pandemic (0.088), and female gender (0.087).

For both financial and illness-related worries, low to moderate levels of worry were associated with decreased or average predicted distress, while severe worries increased predicted distress (0.41 for “very worried” about finances; 0.34 for “very worried” about illness). To a lesser extent, greater hopefulness was associated with decreased predicted distress. The relationship between age group and predicted emotional distress was somewhat linear: membership in the youngest age group (18 to 39 years) increased predicted distress, while membership in the age groups of 40 to 59 years and ≥60 years decreased predicted distress. Female gender increased predicted distress, and both increased and decreased alcohol intake compared to prepandemic intake had higher mean SHAP values compared to no change (mean SHAP values of 0.19 for increased use 0.12 for decreased use; –0.07 no change), indicating a nonlinear effect. Change in cannabis use due to the pandemic was directionally similar (mean SHAP values of 0.25 for increased use, 0.11 for decreased use, and –0.02 for no change), although its overall importance was lower.

Following the top 6 ranked variables, we found a heuristic “elbow point” separating the most important variables from those with lesser importance. For all remaining variables, the mean SHAP values were at or below 0.05. Of note, some important variables demonstrated low mean SHAP values, but their range of values was large. Notably, only a small number of respondents answered “other” (nonbinary) for gender identity (n=47); although these responses had a strong influence on their individual predicted values, mean SHAP values for the feature as a whole remained relatively low (mean absolute SHAP value 0.006, range –0.04 to 0.38).

To examine the effects of time, we explored predicted emotional distress values across survey waves. The survey wave variable ranked 14th out of 50 variables with respect to variable importance (mean absolute SHAP value 0.024), indicating that the passing of time was not as influential in our model as other time-independent variables. A partial dependence plot is shown in Figure 2C; when all other variables were held constant, predicted emotional distress was highest in wave 1 (May 8-12; adjusted mean 0.06). After wave 1, the values decreased (wave 2 adjusted mean –0.01). The predicted values decreased most by membership in wave 3 (mean –0.04, June 19-23) and wave 4 (mean –0.03, July 10-14), then increased again in waves 5 and 6 (mean 0.01 for September 18-22 and 0.04 for November 27-December 1).

We next performed exploratory analysis of the importance of two-way variable interactions in our model. Mean SHAP values for two-way variable interactions from the top 15 most important variables are presented in Figure 3; all interactions are shown in Figure S5 (Multimedia Appendix 1). Overall, interactions between features were not substantially important in determining model predictions; the mean absolute SHAP values for the interaction terms ranged from 0 to 0.01, which represents a maximum of approximately 1/17 of the mean contribution of the most important variable (worry about finances). A total of 5 interactions had SHAP values above 0.008 (the visual elbow point); these top interactions included change in alcohol use × worry about finances (SHAP 0.010), worry about getting COVID-19 × worry about finances (SHAP 0.010), hopefulness × high risk for COVID-19 (0.009), hopefulness × female gender (0.009), and worry about getting COVID-19 × cannabis use (0.008).

Figure 4 displays relationships between pairs of variables in the top 5 most important interactions. Plotting both COVID-19–related worries (financial and illness-related) against predicted emotional distress, the distressing effects of illness-related worry were stronger at lower levels of financial worry, where individuals with severe financial worry had the greatest distress regardless of illness-related worries. Additionally, greater hopefulness mitigated differences in distress levels between those who were at high risk for COVID-19 (or had loved ones at high risk) versus those who were not, as well as differences in distress between female and nonfemale respondents. Finally, those who used cannabis had a steeper increase in distress as illness-related worries increased.

To validate our findings, unregularized regression analyses were performed to test the statistical significance of pairwise interaction models with realized distress values. Out of 1225 possible interactions between pairs of explanatory variables in bivariate linear regression on emotional distress, 58 interaction terms had significant associations after Benjamini-Hochberg correction. Of these 58 interactions, 10 involved hopefulness. Of all 5 interactions with SHAP variable importance above the elbow point of 0.008, the top 4 were also statistically significant in the regression analysis (see Table S2 in Multimedia Appendix 1).

In this study, we used machine learning to examine factors associated with emotional distress during the COVID-19 pandemic, informed by self-reported levels of anxiety and depression symptoms, using data from a large national survey. We explored relationships between a wide range of sociodemographic characteristics, substance use patterns, and COVID-19–related perceived risks and worries with distress by examining nonlinear patterns in variable importance and by characterizing the importance of variable interactions.

Our findings that the top predictive factors for emotional distress included worries related to COVID-19 are consistent with a recent study of COVID-19–related anxieties [10], suggesting that the participants experienced substantial health and financial concerns. Female gender contributed substantially to increased distress; this finding is consistent with evidence that even before the pandemic, both general anxiety [31] and depression [32] as well as COVID-19–specific anxieties [10] have been shown to be greater in women than in men. Our top predictive features also included change in alcohol use, which showed a nonlinear effect whereby any change (consuming either more or less alcohol compared to before the pandemic) was associated with an increase in emotional distress. This suggests that individual attempts to mitigate above-average levels of pandemic-related distress by drinking less have not been successful and that those who increased alcohol consumption due to the pandemic have experienced subsequent heightened distress.

Hopefulness was only weakly to moderately correlated with other anxiety and depression questionnaire items, and it was loaded onto its own principal component in PCA. These findings suggest that hopefulness may not measure the same latent construct captured by the remaining anxiety and depression survey items. This relative lack of cohesion between hopefulness and the other CES-D questions (feelings of loneliness and depression) in latent variable analysis is supported by previous work finding structural inconsistencies among CES-D items [33,34] Given the sudden and temporary nature of imposed pandemic restrictions, we hypothesized that the hopefulness responses represent a more trait-like positive affect rather than more situationally influenced responses to social isolation, such as lonely and depressed feelings. In predicting emotional distress, hopefulness was the fourth most important contributor, with higher hopefulness decreasing overall predicted distress. Hopefulness has been linked to lower emotional distress [35], and it may indicate greater resilience to adversities experienced during the COVID-19 pandemic.

One main aim was to assess changes in predicted emotional distress over time. The survey wave was not among the strongest contributors to emotional distress, with an importance approximately one seventh of that of the top variable. Despite this, variation in predicted values was present across time; the highest predicted emotional distress was observed during survey wave 1 in May 2020, amid nationwide lockdowns and just after the peak 7-day average case count [30] in Canada’s “first wave” of COVID cases. Predicted values declined through waves 2 to 4, when case counts were decreasing and lockdown policies were relaxed. Waves 5 to 6 saw a second increase in predicted emotional distress levels through the fall and winter of 2020, as Canada began its second wave of cases. This pattern mirrored the trajectory of COVID-19 case counts nationwide [30]. Our findings are consistent with results from the Canadian Community Health Survey, which found that mental health worsened due to the pandemic, increasing from September to December in 2020 [36]. These results suggest that Canadians experienced a spike of emotional distress at the start of the pandemic, amid fearful public health messaging and great uncertainty. Following this, as case counts decreased and lockdown measures were lifted, we speculate that increased optimism or a reduction in the perceived threat of COVID-19 may have lowered collective distress. This summer period was followed by the fall and winter months, when seasonal changes and increasing case counts again led to increased distress.

Finally, we conducted an exploratory analysis to examine the importance of variable interactions in predicting emotional distress. Although the importance of these interactions was relatively low (1/17th that of the most important single variable), they played a role in determining the predicted values. In particular, the distressing effects of severe COVID-19–related worries were most pronounced when other worries were not present. Notably, high hopefulness also mitigated the effect of several other factors that increased predicted distress, including female gender and high risk for COVID-19.

This study has several strengths. First, it was conducted on a large national survey sample that was designed to be representative of the Canadian population in age, gender, and region. Second, few other studies have examined changes in mental health outcomes across time throughout the first 10 months of the COVID-19 pandemic. Third, we employed flexible, interpretable machine learning methods to detect nonlinear and interactive relationships of many predictors of anxiety and depression symptoms.

A limitation of this study was its cross-sectional design, which meant that we were not able to track changes in mental health outcomes within the same individuals longitudinally, nor were we able to determine the temporality or direction of associations. However, the repeated cross-sectional study design allowed us to track patterns in mental health outcomes over time. A second limitation was that the survey was not administered prior to the COVID-19 pandemic, so a direct comparison with prepandemic mental health was not possible. A third limitation was that the survey was administered on the web, meaning that Canadians who are not comfortable with technology may have been less likely to participate; furthermore, the complete response rate was low (16.1%), indicating potential selection bias. However, quota sampling techniques were used to represent the Canadian adult population as accurately as possible in the complete analysis data set, and this response rate is similar to that expected for population surveys of this length administered on the web and without financial incentive [37]. Finally, information on mental health histories or prior clinical diagnoses was not available. However, these findings provide insight into emotional distress of the Canadian population at large, independent of clinical diagnoses.

Demographic and COVID-19–related factors were associated with a substantial amount of the variation in emotional distress during the global COVID-19 pandemic. These associations were most strongly driven by COVID-19–related fears, namely worries about personal finance and worries about contracting the illness, and were partly mitigated by high levels of hopefulness.

Rates of negative mental health outcomes such as eating disorders [38], substance use [39], overdose [40], and suicide attempts [40] have risen over the course of the COVID-19 pandemic. Although public health policy has been and continues to be vital to curb the spread of SARS-CoV-2, policy makers should also prioritize the provision of population-level supports to address elevated depression and anxiety symptoms among certain groups. Although we cannot infer causation in our study design, our results indicate that initiatives to mitigate financial worry, alleviate illness-related fear, and promote hopefulness may be effective against symptoms of anxiety and depression in the wake of this and potential future pandemics.

All study data have been made publicly available [41].