Methods
We analysed data from the birth certificate (CeDAP), which is filled in at birth for each delivery.
Ten Italian Regions and one Autonomous Province (Piedmont, Lombardy, Veneto, Emilia-Romagna, Friuli-Venezia Giulia, and the Province of Trento in Northern Italy; Tuscany and Lazio in Central Italy; Apulia, Campania, and Sicily in Southern Italy) agreed to participate. These Regions cover 84.3% of all the births in Italy.10
We defined March 1st, 2020 - March 31st, 2021 as the pandemic period: this time covered the first two waves of COVID-19 in Italy, corresponding to several restrictions measures. The historical period included the three previous years, from January 2017 to February 2020. For the Campania Region the comparison period started from January 2018 because 2017 data were not available.
The primary outcomes were PTB (live births between 22 and 36 weeks’ GA) and stillbirths, both in singleton and multiple pregnancies. Secondary outcomes were late PTB (32-36 weeks’ GA), very PTB (<32 weeks’ GA), and extremely PTB (<28 weeks’ GA). GA at birth was calculated in completed weeks.
We used univariable and multivariable logistic regression models to examine the association between birth period (pandemic vs historical) and percentage of preterm births estimating odds ratios (ORs) of each outcome. Adjusted analysis included the following variables: maternal country of birth (foreigner vs Italian), maternal age at index birth (continuous), parity (yes/no), maternal education (none or elementary school or primary school diploma; secondary school diploma; University degree), maternal employment (yes/no), pregnancy conceived with assisted reproductive technology (ART, yes/no), sex of the child (female/male). Most of these variables could be in fact considered as mediators or effect modifiers rather than true confounders:11 i.e. mitigation strategies due to COVID-19 could have influenced maternal lifestyle in pregnancy in a different way according to mother’s origin, age, education, employment, and parity. Pregnancy conceived with ART, on the other hand, is an intermediate variable. We therefore chose as the main analysis the unadjusted one. In the adjusted analysis the Lazio Region was not considered because the information on ART was not available for the whole study period.
We further analysed separately singletons and multiples, which could be differently affected by the pandemic restrictions.
Due to privacy regulations, individual data were not shared, and logistic regression analyses were run in each Region. A meta-analysis was performed centrally at the Regional Health Agency of Tuscany.
To estimate the heterogeneity of effects in different Regions, the I² index was calculated.12 When there was no evidence of heterogeneity, the pooled estimate of the effect (OR) was calculated using the inverse variance method (fixed effect model); otherwise, the DerSimonian–Laird weights (random effects model) were used.13 Forest plots were provided to graphically illustrate the effect size estimates for each study as well as the pooled estimate.
We also studied the monthly trend of PTB rates from 2017 to the end of March 2021 using an interrupted time series regression analysis,14 with March 1st 2020 as the date of interruption. In this quasi-experimental technique, one looks for a sharp change in outcomes following public health interventions (the interruption corresponding to the implementation date of COVID-19 mitigation measures). After a visual check of data points, we carried out a log-linear regression analysis of the (log of the) monthly prevalence of PTB over calendar months and introduced a term estimating the level of discontinuity (gap) on March 1st 2020 [at month 38/39], i.e. at the “interruption”. As the prevalence of PTB showed a seasonal trend, we modelled it using Fourier terms (2 pairs of sine and cosine functions).14 In addition, we used the robust Newey-West standard errors for effect estimates to account for residual autocorrelation in the data (“newey” command in Stata).
In all models, we also tested whether the slope had been altered by mitigation measures by running a model including a statistical interaction between slope and period (historical vs pandemic period). As there was no evidence of interaction and the models without interaction had better fit to the data (lower residual MSE), we always used models without interaction.
As a further check of the overall effect of mitigation measures, we carried out a regression of the log of the monthly frequency of PTB over calendar months until February 2020 (just before the pandemic), correcting for seasonal trend and autocorrelation as above. We then computed the expected frequencies for the months following the lockdown (i.e., under the counterfactual scenario of no intervention), and compared them to actual frequencies.
All the analyses were performed using STATA version 15.1 (StataCorp LP, College Station, TX, USA).