This study was conducted in Galicia (north-west Spain), a region with 2.7 million inhabitants (25.7% over 65 years of age), practically all of whom (98%) are covered by the Public Health System (PHS). The range of services provided to citizens includes preventive, diagnostic, therapeutic and rehabilitative care, as well as health promotion and maintenance activities. Visits to the doctor are free of charge but out-of-hospital pharmaceutical services are subject to a financial contribution (co-payment). Since patients’ contributions are income-linked, co-payment for medication ranges from 0 to 60%. In Spain, medication in the outpatient setting is mostly dispensed by community pharmacies.

The Galician PHS has a clinical record system which allows access to a unified register of all patients’ clinical information (diagnostic tests, drug prescriptions, International Classification of Primary Care codes, hospitalizations, etc.) for both primary and hospital care.

We used a population-based case–control design [13]. This design is characterized by using data from a representative sample of all cases (in our instance, exhaustive sampling) in a precisely defined and identified population (in our instance, the population attended by the PHS in Galicia), and comparing these to data on persons (controls) randomly extracted from the same population as the cases (population-based case controls), something that could be assumed to give a valid estimate of the prevalence of exposure and covariates in the population of origin [14]. According to Rothman et al., this design can be considered the most desirable option for a case–control study [13].

To study risk of hospitalization, ‘case’ was defined as any patient admitted due to COVID-19, with PCR confirmation, to any public hospital in Galicia since the onset of the pandemic, whose clinical course ended before 1 January 2021. As controls, we selected a random sample of subjects who had no positive PCR (as a result of not having done any PCR or for being PCR negative) during the same period. With the aim of enhancing the efficiency of our analysis of risk of hospitalization [15], controls were randomly selected and matched with cases by age, sex and primary healthcare centre. We selected up to 20 controls for each case.

To assess the risk of progression to severe COVID-19, understood as the risk of requiring admission among COVID-19 positive subjects, we used the same cases as those used to assess the risk of hospitalization (all patients admitted due to COVID-19 with PCR confirmation). As controls, we used all patients with diagnosis of COVID-19, who did not require hospitalization. While the use of these controls unmatched with cases of progression reduces the efficiency of estimates for this model [16], it does not generate confounding bias because (1) the controls were drawn from the same population as the cases [13]; (2) they were selected regardless of exposure [13]; and (3) in the statistical analysis, adjustment was made for available potential confusing variables [15].

To assess susceptibility, cases were defined as all persons diagnosed with COVID-19 confirmed by a positive PCR (hospitalized and not hospitalized) across the study period in Galicia [15]. As controls, we used the same persons as those used to assess risk of hospitalization (subjects who had no positive PCR). As in the case of progression to severe COVID-19, the cases were unmatched with controls, which does not produce any type of bias, only a decrease in efficiency [13, 16].

The study was approved by the Galician Clinical Research Ethics Committee (Comité de Ética de Investigación de Galicia), reference 2020/349, certified by the Spanish Agency of Medicines and Medical Devices (Agencia Española del Medicamento y Productos Sanitarios/AEMPS), and conducted in accordance with the Helsinki Declaration principles and current biomedical research legislation. The study protocol is registered in the EU electronic Register of Post-Authorisation Studies (EUPAS44587) and is available online at https://www.encepp.eu/encepp/viewResource.htm?id=44588.

Data extraction was automated, and the data were anonymized to ensure that the subjects concerned could not be identified.

Automated data extraction was performed by an independent information technology services (IT) company, based on SERGAS’ Complex Data-Analysis System (Sistemas de Información y Análisis Complejos/SIAC). SIAC acts as a data warehouse which stores information for the management of different systems (dispensing of medications, diagnoses and hospitalizations, among others).

The variables of exposure were defined as any of the following anticoagulant medications: direct-acting oral anticoagulants (apixaban, edoxaban, dabigatran and rivaroxaban); vitamin K antagonists (acenocumarol and warfarin); LMWH (enoxaparin, bemiparin, dalteparin, fraxiparine, nadroparin and tinzaparin); and other anticoagulants (fondaparinux and protein C) (see electronic supplementary material [ESM]). Consumption of one type of anticoagulant does not exclude the fact that the patient may consume another type of anticoagulant. We recorded those prescribed and dispensed to cases and controls alike across the study period, in the 6 months leading up to the index date. The index date was defined as the tenth day prior to the PCR+ date, or, for non-PCR+, the same index day as their matched case.

As study covariates, we recorded demographic and anthropometric variables, clinical COVID-19 variables (where applicable), and data on hospitalization, comorbidities (hypertension, diabetes mellitus, chronic obstructive pulmonary disease [COPD], obesity, ischaemic heart disease, cerebrovascular accident, heart failure, atrial fibrillation, chronic renal failure, cancer, asthma, current smoker) and exposure to the other medications prescribed for and dispensed to each of the subjects (antihypertensives, diuretics, non-steroidal anti-inflammatory drugs, hypolipidemic agents, other anticoagulants, antiplatelet agents and glucocorticoids). Polypharmacy (number of different medications prescribed and dispensed for chronic diseases in the last 6 months before the index day) was used as a proxy for the degree of chronicity of the patients [17]. All covariates were recorded for the 6 months prior to the index date.

Generalized linear mixed models were subjected to statistical analysis, with the binomial link function. These models were used because of the structure of the data and because they have many advantages over conditional regression [18,19,20]: (1) they allow the analysis of matched and unmatched models; (2) they permit the introduction of random terms to control for heterogeneity of initial clusters and time periods; (3) strata in which cases coincide in exposures with controls continue to count as events for the calculation and for the estimates.

To construct the models, the following four levels were considered: patient; case–control strata (case and matched control); health centre; and pandemic wave. We used random effects to assess the effect of the pandemic wave, and nested random effects for patients, case–control strata, and health centre. Adjustments were made for potential confounding variables, including sex, age, comorbidities, smoking habit, and each additional pharmacological treatment. Results were expressed in adjusted odds ratio (aOR) with its 95% confidence interval (CI). Adjusted estimates were obtained for the effect of anticoagulant therapy dispensed compared with the absence of any anticoagulant drug therapy. The models were estimated using the glmer function of the lme4 R package [21] (R Software version 4.1.0).