*2.4. Statistical Analysis*

Missing data were handled by mean imputation, but if exceeding 5%, the case was removed from further statistical analysis. In the preliminary analysis, descriptive statistics were calculated for all measures, including the mean (M), 95% CI, standard deviation (SD), median, skewness, and kurtosis. All scores demonstrated good parametric properties. Therefore, further parametric analyses were conducted. The study hypotheses were tested in several ways. Repeated measures one-way ANOVA was used to examine changes in the mean scores of life satisfaction, perceived stress, and coping styles during the first and second waves of the COVID-19 pandemic. Effect sizes were calculated using the partial eta-squared statistic (η*<sup>p</sup>* 2).

The Pearson's correlation matrix was used to test the bivariate correlations among the variables. Multiple linear regression was conducted to assess the association between life satisfaction during W2 as a dependent variable and perceived stress and coping styles during W1 and W2 of the COVID-29 pandemic as a predictor variable. The mediating role of coping styles in the relationship between perceived stress and life satisfaction was examined in a cross-sectional design, separately for the first and second pandemic waves, using structural equation modeling (SEM). Parallel mediation models (simultaneous analysis of all three coping styles as a mediator) were performed based on maximum likelihood (ML) estimation without missing values and including observed variables. The conditional effect was examined based on a bias-corrected bootstrapping procedure with 1000 samples. A bootstrap confidence interval (95% CI) not including 0 signaled a significant effect.

Next, a two-wave cross-lagged panel design with a time lag of a half-year between W1 and W2 was performed for testing the prospective effect of perceived stress and coping styles on life satisfaction [84]. Two waves of data collection are recommended as the optimal approach to longitudinal design because it curtails the cost of a study in terms of time and money [85]. Structural equation modeling (SEM) was used, based on maximum likelihood (ML) estimation, without missing values, including observed variables, and with the bootstrapping method to test the conditional effect. The parallel mediation model specified included ten observed variables. We decided not to use latent variables because of the problems related to measurement invariance for the total 126 items included in 10 variables (W1 and W2). Both cross-lag and autoregressive paths were used to examine stability and change simultaneously [86,87]. All latent structural models were evaluated using several goodness-of-fit criteria, including ML χ2, *df,* and *p*-value (the ratio χ2/*df* < 5 representing a good fit), standardized root mean square residual (SRMR ≤ 0.08 indicates an acceptable fit), root mean square error of approximation (adequate fit if RMSEA ≤ 0.08), comparative fit index (CFI ≥ 0.90 meaning adequate fit), normed fit index (NFI ≥ 0.95 considering adequate fit), and the Tucker–Lewis index (TLI ≥ 0.90 indicating adequate fit) [88–90]. All statistical analyzes were conducted using SPSS 27 (for descriptive statistics, ANOVA, and correlation analysis) and AMOS 22 (for SEM).

### **3. Results**
