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).