Effects of Vaccination and the Spatio-Temporal Diffusion of Covid-19 Incidence in Turkey.

This study assesses the spatio-temporal impact of vaccination efforts on Covid-19 incidence growth in Turkey. Incorporating geographical features of SARS-CoV-2 transmission, we adopt a spatial Susceptible-Infected-Recovered (SIR) model that serves as a guide of our empirical specification. Using provincial weekly panel data, we estimate a dynamic spatial autoregressive (SAR) model to elucidate the short- and the long-run impact of vaccination on Covid-19 incidence growth after controlling for temporal and spatio-temporal diffusion, testing capacity, social distancing behavior and unobserved space-varying confounders. Results show that vaccination growth reduces Covid-19 incidence growth rate directly and indirectly by creating a positive externality over space. The significant association between vaccination and Covid-19 incidence is robust to a host of spatial weight matrix specifications. Conspicuous spatial and temporal diffusion effects of Covid-19 incidence growth were found across all specifications: the former being a severer threat to the containment of the pandemic than the latter.

Introduction

In the absence of pharmaceutical interventions (PI), nonpharmaceutical interventions (NPI) were introduced to flatten the epidemiological curve along with public measures to smooth out the fiscal burden. The warfare against Covid‐19 has reached a turning point upon the use of effective and safe vaccines. As Covid‐19 continues to threaten global public health, only about 63.4% of the world population has received at least one dose of vaccine. As of March 2022, 62% of Turkey's population have been fully vaccinated and an additional 5.7% has been partly vaccinated (Ritchie et al., 2022).

This study assesses the spatio‐temporal impact of vaccination efforts on Covid‐19 incidence growth using panel data for 81 provinces of Turkey over a period of 26 weeks. Incorporating geographical features of SARS‐CoV‐2 transmission, we adopt a spatial SIR model that explicitly takes the spatial distribution of units into account and that serves as a guide of our empirical spatio‐temporal specification. We posit that Covid‐19 incidence growth at location and time depends not only on its own past growth but also on the past growth at other locations. We estimate a dynamic spatial autoregressive (SAR) model to elucidate the short‐run (SR) and the long‐run (LR) impact of vaccination on Covid‐19 incidence growth after controlling for temporal and spatio‐temporal diffusion, testing capacity, social distancing behavior, and space‐varying unobservable confounders. Results indicate that vaccination efforts significantly reduce Covid‐19 incidence growth in Turkey. Specifically, a 1% point increase in vaccination growth rate translates into 0.814%–0.825% points decrease in Covid‐19 incidence growth in the SR and 2.34%–2.40% points decrease in the LR. A breakdown of the LR effect indicates significant direct as well as positive externalities or spatial spillover effects induced by vaccination.

Our study contributes to the emerging literature on the regional dimension of Covid‐19 crisis from a number of perspectives. While the global evolution of the pandemic receives interest among policy makers and scholars, less attention is given to local evolution and its long‐run impact on mitigation efforts. Studies show that local differences in geographical conditions (Gupta, Banerjee, and Das, 2020; Shi et al. 2020) and population density (Carozzi, 2020; Baser, 2021) might play a role in the regional evolution of the pandemic. Additionally, local conditions such as labor market outcomes, pre‐Covid‐19 remote work conditions, health care capacity, and institutional quality are crucial in understanding the differential vulnerability to Covid‐19 in the EU regions (Kapitsinis, 2020; Amdaoud, Arcuri, and Levratto, 2021; Christopoulos, Eleftheriou, and Nijkamp, 2021; Rodríguez‐Pose and Burlina, 2021).

With notable exceptions (Guliyev, 2020; Krisztin, Piribauer, and Wögerer, 2020), there has been a widespread adoption of spatial approaches using cross‐sectional data at the subnational level. These studies use exploratory spatial data analysis (ESDA) techniques that suggest a strong clustering of Covid‐19 cases (Kang et al., 2020; Kim and Castro, 2020; Xiong et al., 2020; Hazbavi et al., 2021; Wang et al., 2021) in conjunction with global spatial models (Sun et al., 2020; Amdaoud, Arcuri, and Levratto, 2021; Andersen et al., 2021) or explore spatial heterogeneity (Mollalo, Vahedi, and Rivera, 2020; Dutta, Basu, and Das, 2021; Mansour et al., 2021). However, cross‐sectional data are incapable of capturing the timing of events, temporal dependence, and most importantly, unobservable heterogeneity. In contrast, we adopt a panel study to address the shortcomings of cross‐sectional analysis.

A related aspect of Covid‐19 modeling is calibration of spatial diffusion. Evidence shows that units surrounded by Covid‐19‐burdened neighbors are more likely to be further affected by the pandemic (Dehghan Shabani and Shahnazi, 2020; Krisztin, Piribauer, and Wögerer, 2020; Florida and Mellander, 2021). In prior studies, this spatial diffusion is captured by a spatial lag model (SAR) among others, without further consideration of the timing of diffusion. We posit that the spatial diffusion, which signals out the importance of regional networks, is not contemporaneous and must account for at least the average incubation period of the virus.

Outside of frequentist approaches, there are alternative methods that encode prior probabilities based on past knowledge of epidemics. These methods employ spatio‐temporal mapping of diseases for establishing an early warning system (Knorr‐Held, 2000; Lawson, 2006) and Bayesian outbreak forecasting methods that suggest the existence of distinct local characteristics that influence the temporal variation of viral infections such as dengue (Jaya and Folmer, 2020, 2021b) and SARS‐CoV‐2 (Jaya and Folmer, 2021a). Bayesian approach with spatio‐temporal random effects enables modeling without covariates. In contrast, our study is concerned with the structural influence of vaccination rather than detecting clusters and hot‐spots of the pandemic. A key element of our setup is the ability to define spatial transmission based on an epidemiological compartmental model.

The bulk of related studies focuses on environmental (Han et al., 2021), socioeconomic (Sannigrahi et al., 2020; Ehlert, 2021), sociodemographic (Raymundo et al., 2021) or geographical (Gupta, Banerjee, and Das, 2020) determinants of Covid‐19. However, a theoretical basis and/or causal reasoning to guide the empirical specification is lacking. Such ad hoc analyses run the risk of including bad controls ‐ a covariate or a collection of covariates whose inclusion increase the bias on the quantity of interest ‐ or excluding good ones, potentially leading to biased estimates and poor inference. In contrast, our empirical strategy is grounded in an epidemiological SIR model that we present in section “A spatial SIR model with vaccination and testing” and that explicitly takes the spatial distribution of units, vaccination, and testing into account. To add transparency, we invoke a graphical model with spatio‐temporal dependence to guide us through the empirical specification in section “Empirical strategy”.

One of the most overlooked issues in the literature is how Covid‐19 burden is defined. A disease containment policy aims to slow down the propagation of the virus (i.e., flattening “the curve”). From a policy‐relevant or epidemiological standpoint, the appropriate measure of Covid‐19 cases is therefore incidence and not prevalence. While the latter is a stock variable measuring the total number of outstanding or cumulative cases, the former is a flow variable measuring new or marginal cases at a given point in time. Our screening of the literature indicates only less than half employ incidence as the outcome variable of interest.

Studies point out the importance of local factors or specific containment strategies such as lockdowns (Bourdin et al., 2020; Borri et al., 2021; Orea and Álvarez, 2021); however, they do not consider how vaccination can play a role in disease containment at the local level, possibly due to being an effective treatment only recently. Our study is one of the first attempts to examine the impact of vaccination at the regional level that hints about how equal access to vaccination and compliance with an effective treatment could be ensured. These aspects are central, especially for isolated, less developed, and less educated regions and even more so for developing countries with sizable regional disparities.

We offer several policy insights. First, a regionally blind policy in a severely Covid‐19‐affected country with sizable local differences is less likely to be effective in controlling the pandemic. While there is inertia in Covid‐19 incidence growth, the spatial SIR model is indicative of a spatial memory in the diffusion of the pandemic, backed up by data. Regions spatially proximate to those with higher incidence growth rates in the previous week are likely to be profoundly affected by the surge. Second, vaccination is the current major policy tool to control the diffusion of the pandemic. Given the spatio‐temporal nature of incidence growth, PI in a region generates spatial spillover to surrounding spatial units. Therefore, policy implementations should be regionally cohesive to contain the pandemic. Incongruent policies across regions amid inadequate vaccination efforts will result in a new wave and will force the government to reintroduce drastic NPIs in Turkey.

A spatial SIR model with vaccination and testing

Inspired by Camacho et al. (1996); Chernozhukov, Kasahara, and Schrimpf (2021); Ghostine et al. (2021); Xu, Wu, and Topcu (2021), we develop a simple Susceptible–Infected–Recovered (SIR) model of Covid‐19, taking into account vaccination, testing, and the distribution of units in space. Variants of SIR model for Covid‐19 have been studied under Susceptible–Asymptomatic–Infected–Recovered (SAIR) with mobility (Aràndiga et al., 2020) and Susceptible–Exposed–Infected–Recovered (SEIR) (Huang, Liu, and Ding, 2020; Danon et al., 2021) models. Here, our aim is to augment existing models by incorporating a spatial dimension within the modeling framework of infections and use the implications of the epidemiological model as a guide for our empirical specification to avoid an ad hoc analysis and enhance the policy relevancy of the empirical findings.

Let be the coordinates, , be the time, , , and be the fraction of susceptible, infected, and recovered individuals respectively at time , location such that and therefore . consists of fraction of individuals who developed immunity via coming through the illness or vaccination. Fig. 1 shows these three compartments and the disease propagation.

Figure 1

Covid‐19 propagation flow for the spatial SIR model. Note:

The laws of motion are given by: where is the spatially weighted average of infections at neighboring units, and indicates the relative connectivity between location and , such that . The infection is transmitted within units at rate (within‐infection rate), and between units at a rate (between‐infection rate). Therefore, the change in the number of infections at location depends not only on its own infections but also on the number of infections of neighboring locations defined by the spatial weighting. This dimension of the model represents the core idea of incorporating the spatiality of the diffusion of Covid‐19. A fraction of susceptible receives vaccination at a rate and move from compartment to . Those who are infected move from compartment to , recover from the disease at a rate and then move from compartment to .

Our outcome of interest is Covid‐19 incidence, . Covid‐19 prevalence, , is governed by the following law of motion: where is the time‐varying rate at which infections are detected. Differentiating equation (4) with respect to time, we obtain the law of motion of :

Dividing equation (5) by equation (4) and using equation (2), we obtain:

We do not observe but only and its spatial average, . Eliminating using equation (4), we obtain:

The left‐hand side (LHS) of equation (7) is the approximate growth rate of Covid‐19 incidence. The first term on the right‐hand side (RHS) is the spatial spillover with denoting the weighted average new Covid‐19 cases at neighboring locations per Covid‐19 positive case at location , the second and the third terms are confounders other than detection and the last term is detection growth rate. Equation (7) is the backbone of the empirical specification shown in Spatio‐temporal models and explains how spatial spillovers affect the evolution of the pandemic.

The basic within () and between () reproduction rates are respectively defined as the ability of the disease to spread within and between provinces in the absence of PI, , . Hence, the total basic reproduction rate is:

Differentiating equation (8) with respect to time and dividing all sides by yield the growth rate of basic reproduction number:

The LHS of equation (9) is the approximate growth of reproduction rate and the second term on the RHS is the recovery growth rate. For , and , if we assume that the growth rate of between infection is greater than that of within infection, that is, , then it follows that . Hence, the first term on the RHS of equation (9) is some weighted average of between and within infection growth.

In the absence of PI, a constant or decreasing basic reproduction rate, , implies:

Equation (10) suggests that without an effective Covid‐19 treatment, the recovery growth must surpass the weighted average infection growth rate in order to prevent basic reproduction rate from growing.

With the introduction of PIs, we define the total effective reproduction rate (), which is the average number of new secondary cases produced by a single Covid‐19 positive case:

where is vaccine efficiency, and is the nonimmunization rate. If , then the pandemic cannot be contained. If , then the disease can be eradicated. Covid‐19 can be eliminated from the population if and only if the total infection rate is less than the ratio of recovered to nonimmunized, that is, . Hence, the less effective the vaccine or the lower the vaccination rates, the more difficult to keep below 1. Note that vaccine efficiency, , is a function of vaccine production technology as well as virus variants. We discuss these issues and how they might affect inference in section “Conclusion”.

Next, we assess how containment efforts evolve in the presence of vaccination. Differentiating equation (11) with respect to time yields the law of motion of :

Dividing all sides by yields the growth rate of effective reproduction number: where the last term on the RHS of equation (13) is the nonimmunization growth. Although equation (13) does not directly answer whether the pandemic can be eradicated or not, it tells us under what conditions the reproduction growth can be slowed down with vaccination. A constant or decreasing effective reproduction rate, implies:

Since effective Covid‐19 treatment helps keep the reproduction number under control, the threshold below which the weighted average infection growth should lie in order to prevent from growing should be lower following the introduction of PI, relative to its absence (see equation(10)). Notice that this holds as long as the proportion of nonimmunized in the population is decreasing, that is .

Our objective is to arrive at an expression that shows the necessary vaccination effort to keep the Covid‐19 spread under control. For this purpose, we further elaborate in equation (14). For simplicity and interpretability, we assume that vaccines are fully efficient (i.e., ) and hence . Then, equation (14) can be written in terms of the vaccination rate necessary to maintain a nongrowing effective reproduction number:

Equation (15) tells us that the change in vaccination rate must be at least as large as the difference between the weighted infection growth among nonimmunized population and the recovery growth among nonvaccinated in order to prevent the effective reproduction number from growing.

Empirical strategy

Data and sample

Our research uses provincial (Nomenclature of Units for Territorial Statistics 3—NUTS‐3) weekly panel data that come from the Ministry of Health of Turkey (MoH, 2021), TURCOVID19 (Uçar et al., 2020), Turkish Statistical Office (Turkstat, 2020) and the Google mobility report (Google, 2021), covering 81 provinces in Turkey for a period of 26 weeks (from 20 February 2021 to 20 August 2021). A detailed description of our dataset is provided in Table 1.

Table 1

Descriptive Statistics

Variable	Description	mean (SD)	Min	Max	Boxplot	Histogram	Source	# of obs.
Covid‐19 burden
Covid‐19 incidence	Weekly # of new positive cases per 100,000	147.06 (147.25)	2.11	1139.94			Uçar, Arslan, and Yapalak (2020) MoH (2021)	2106
Covid‐19 incidence growth rate (%)	‐	4.71 (38.13)	−174.49	190.66			Authors' calculation	2025
Pharmaceutical intervention
Vaccination rate (%, ≥2 doses)	Cumulative weekly # of people with ≥2 doses as % of total province population	18.84 (15.80)	0.63	82.88			Uçar, Arslan, and Yapalak (2020) Turkstat (2020) MoH (2021)	2106
Vaccination growth rate (%)	‐	13.14 (11.88)	1.24	64.19			Authors' calculation	2025
Testing capacity
Testing growth rate	Growth rate of total # of PCR tests	3.52 (8.91)	−13.41	25.36			Uçar, Arslan, and Yapalak (2020) MoH (2021)	26
Social distancing
Workplace mobility		−18.81 (14.28)	−75	17
Retail and recreation mobility		−1.55 (33.80)	−74	143
Grocery and pharmacy mobility	Median % change	46.50 (31.37)	−32	190			Google (2021)	2106
Park mobility	from the baseline	35.33 (65.15)	−67	410
Transit station mobility		3.49 (42.84)	−84	184
Residential mobility		2.87 (6.61)	−11	24

Note: The unit of observation is the province‐week, comprising of a sample period of 26 weeks (20 Feb 2021–20 Aug 2021) and 81 provinces. Changes in mobility for each day are compared with a baseline day before the pandemic outbreak, where the baseline day is the median value over the five‐week period from Jan 3 to Feb 6, 2020 that represents a normal value for that day of the week.

Covid‐19 incidence per 100,000 population is defined as the weekly provincial number of new positive cases. It ranges between 2.11 and 1,140 cases per 100,000 population in the sample period with a mean incidence of 147 cases. The PI captured by vaccination rates is the cumulative weekly number of individuals with at least two doses, as a percentage of the total provincial population. Testing capacity is measured by the growth rate of the total number of PCR tests, which is available only at the national level.

Social distancing is captured by daily provincial Google mobility data for six place categories: Grocery and pharmacy, parks, transit stations, retail and recreation, residential, and workplace. Each category shows the change in the length of stay compared with a baseline, measured before the outbreak. This baseline is the median for the corresponding day of the week during the five weeks between January 3 and February 6, 2020. For each week and category, we first calculated the median value of mobility if there is at least one nonmissing value for that week. There were a number of missing values under each category except workplace. Of 2106 observations, 12.2%, 13.99%, 22.97%, 6.7%, and 5.26% were missing, respectively, for grocery and pharmacy, parks, transit stations, retail and recreation, and residential mobility. In order to impute missing values, we first identified the NUTS‐2 (26 regions) of each province. Then, for each place category, the missing value of that week for the province in NUTS‐2 region is imputed using the median value of all provinces in the same NUTS‐2 region and week for which data were nonmissing. There was only one NUTS‐2 region (Ağrı subregion—TRA2) in which data for all provinces were missing but one for the entire sample period of 26 weeks. In that case, we calculated the median value by borrowing from the neighboring province that belonged to the adjacent NUTS‐2 region. Imputing mobility this way is a legitimate approach because NUTS‐2 regions are composed of relatively homogeneous provinces with respect to their economic conditions, labor market outcomes, and production structures.

In the sample period, mean workplace mobility was 19% lower relative to the baseline, indicating an increase in social distancing in workplace. In contrast, grocery and pharmacy, and park mobility were 35% and 36% higher, respectively, relative to the baseline, while retail and recreation and residential mobility were at par with their baseline values. All social distancing measures, except workplace and residential mobility, varied wildly in the sample period, as shown by the range in Table 1.

Fig. 2 displays the ridgeline plot (also known as joy plot) of the provincial Covid‐19 weekly cases using viridis colors that are robust to various forms of colorblindness. The distribution is represented using lowess smoothing on the normalized values with a slight overlap. The plots clearly show the extent and magnitude of the first wave through March–April 2021, followed by the second wave in July–August 2021.

Figure 2

Provincial Covid‐19 weekly cases in Turkey. Note: Provinces are given from bottom to top, left to right, in an ascending order of cumulative Covid‐19 cases for a period of 26 weeks (from 20 Feb. to 20 Aug., 2021). [Colour figure can be viewed at

Exploratory spatial data analysis (ESDA)

Choropleth maps are commonly‐used visual tools to inspect the geographical change in variables of interest defined by administrative boundaries where the distribution is captured by shading colors. However, they show abrupt changes at the boundaries of spatial units when in fact such changes occur as a continuum. In the context of disease mapping, methods that capture continuous change in space may be more appropriate (Jaya and Folmer, 2020). For this reason, we use kernel density estimates of the intensity function of two‐dimensional spatial point patterns and isopleths in the mapping of Covid‐19 incidence and vaccination growth. The intensity function, at each grid point , , defines the expected number of objects per unit area at location in the study region and is given by where is the kernel function, is the Euclidean distance between data point and grid point , is the kernel bandwidth at grid point , is the area of the subregion of R over which the kernel function is evaluated at grid point , is a constant of proportionality, and is the number of objects located at data point . We use a quartic kernel function and the minimum number of data points to specify the method to be used for setting the value of the kernel bandwidth.

The geographical distribution of the intensity function of Covid‐19 incidence and vaccination growth rates are displayed respectively in Figs. 3 and 4, eight weeks apart. The administrative provincial boundaries and province names are superimposed for informativeness. In the first three maps in Fig. 3, corresponding to weeks 2, 10, and 18, there is limited change in the intensity of local growth of incidence, which shifts from the north to the east. This period corresponds to an increase in the speed at which people are vaccinated in the northern regions. An increase in vaccination growth rate in the north leads to a significant fall in Covid‐19 incidence growth. These observations point to the effectiveness and the spatial dimension of vaccination.

Figure 3

Intensity of Covid‐19 incidence growth: (a) Week 2 (27 Feb–5 Mar, 2021); (b) Week 10 (24 Apr–30 Apr, 2021); (c) Week 18 (19 Jun–25 Jun, 2021); (d) Week 26 (14 Aug–20 Aug, 2021). [Colour figure can be viewed at wileyonlinelibrary.com].

Figure 4

Intensity of vaccination growth: (a) Week 2 (27 Feb–5 Mar, 2021); (b) Week 10 (24 Apr–30 Apr, 2021); (c) Week 18 (19 Jun–25 Jun, 2021); (d) Week 26 (14 Aug–20 Aug, 2021). [Colour figure can be viewed at wileyonlinelibrary.com].

Fig. 5 plots the evolution of the strength of global spatial dependence in Covid‐19 incidence growth and vaccination growth over the sample period using the Moran's I statistic along with a 99% confidence interval. In both subfigures, we employ a distance‐based and a socioeconomic spatial weight matrix. Distance‐based spatial weight matrix is specified as an inverse distance whose elements contain the inverse of the distance between the centroid of provinces and . For the fact that a row‐normalization renders the inverse distance weight matrix asymmetric and that the impact of unit on unit is not the same as that of unit on unit , its economic interpretation in terms of distance decay is no longer valid (Elhorst, 2001; Anselin, 2013). Elhorst (2001) and Kelejian and Prucha (2010) suggest a normalization procedure where each element of the prenormalized inverse distance matrix is divided by its largest eigenvalue.

Figure 5

Global spatial autocorrelation: (a) Covid‐19 incidence growth; (b) vaccination growth. [Colour figure can be viewed at wileyonlinelibrary.com].

We specify a row‐normalized socioeconomic spatial weight matrix with elements where is the average years of schooling between the ages of 25 and 64 in 2020 and is the absolute value of this difference between provinces and . The underlying idea behind the use of a socioeconomic spatial weight matrix is to capture the cognitive nonphysical connectivity between provinces.

A plethora of other spatial weighting schemes could be used by empirical researchers in the assessment of spatial associations. Our choice is primarily driven by human mobility and by our expectation regarding the nature of the spatial associations under Covid‐19 burden and PI. First, although provinces surrounded by neighbors with high infection rates are more likely to be affected, it is unrealistic to assume that disease propagation does not exert a direct impact beyond neighboring units determined by administrative boundaries, regardless of how such adjacency might be defined. Therefore, threshold distance or contiguity‐based spatial weighting schemes might have less representative power of spatial associations in our context. In contrast, inverse distance and socioeconomic weight matrices provide information about the spatial linkages, defined by nonadministrative boundaries and structural factors that affect mobility.

Fig. 6 displays the intensity plots and the connectivity graph to gain further insight. Fig. 6a,b shows the extent of dissimilarity of the spatial weight matrices with all elements partitioned into 15 bins of equal length in gray scale with darker tones indicating stronger spatial association. Fig. 6c further demonstrates how each province is connected to the remaining provinces. Notice that both spatial weight matrices exhibit the same connectivity pattern, yet they are conspicuously different from one another. For one, the diagonal elements of an inverse distance are zero by definition, yet the diagonal elements of the socioeconomic weight matrix are positive.

Figure 6

Spatial weight matrix (a,b) intensity plots (inverse distance and socioeconomic, respectively) and (c) connectivity graph. [Colour figure can be viewed at wileyonlinelibrary.com].

The Moran's I calculated for each time period does not well inform us about the overall strength of spatial autocorrelation in incidence and vaccination growth. We therefore compute the overall global Moran's I along with the associated ‐value in Fig. 5, which can be obtained by averaging across time periods, . The test statistic for the spatial autocorrelation in panel data follows a normal distribution, where , is the number of cross sections, is the number of time periods, and is the element of the spatial weight matrix connecting provinces and . (Beenstock and Felsenstein, 2019).

Albeit of small magnitude, the overall global Moran's I in Fig. 5 indicates that the null hypothesis of spatial randomness can be rejected in favor of spatial clustering for both variables and spatial weighting schemes, pointing out that the assumption of independence of observations is unjustified. The socioeconomic weight matrix, having the lowest dependence and highest stability, reflects the importance of geographical connectivity, which might overemphasize the impact of nonphysical cognitive ties across provinces under certain circumstances.

Fig. 7 shows the local indicators for spatial association (LISA) map of the dominant regime or quadrant for each province using an inverse distance spatial weighting scheme. In Fig. 7a,b, two types of dominant regimes are determined for the Covid‐19 incidence growth and vaccination growth. In the left panel in both subfigures, for any regime to be dominant relative to every other regime (), the number of weeks the province falls into a quadrant of the Moran scatterplot () must be greater than the sum of the number of weeks during which the province falls into all other quadrants: (labeled as Dominant regime I). In the right panel, for any regime to be dominant, the number of weeks the province falls into a quadrant of the Moran scatterplot () must be greater than the number of weeks during which the province falls into each individual quadrant other than : (labeled as Dominant regime II). The former approach being more restrictive considers the frequency of being detected in a given spatial regime by comparing the frequency of one spatial regime with the summation of the frequency of the remaining regimes. The latter approach directly compares the frequency of a spatial regime with others. Hence, under the first measure (dominant regime I), one is more likely to observe an abundance of no dominant spatial regime.

Figure 7

LISA maps of dominant regimes: (a) Covid‐19 incidence growth; (b) vaccination growth. Note: LISA maps show the dominant regime (i.e.,

In terms of Covid‐19 incidence growth, southeastern provinces form a cluster of high incidence growth for both dominant regimes. On the other hand, a limited number of northwestern regions are located in a low incidence growth cluster for dominant regime I. Western, southern, and northern regions form a low incidence cluster for the dominant regime II. Central regions act as transition regions with higher incidence growth for this second dominant regime. LISA maps for vaccination growth highlight that southeastern regions form the cluster of high vaccination, regardless of the dominant regime. On the other hand, western regions are mostly composed of provinces in lower vaccination growth cluster.

Spatio‐temporal models

We estimate the discrete‐time analogue of equation (7). Using matrix notation, we have a spatio‐temporal model with provinces over weeks: where is a vector indicating the growth rate of Covid‐19 incidence, are the province‐fixed effects, is an vector of ones, is the temporal diffusion rate, is a spatial weight matrix whose elements signify the relative connectedness between province and , is the spatio‐temporally weighted average growth rate of Covid‐19 incidence with being the spatio‐temporal diffusion rate, is pharmaceutical intervention (PI) captured by vaccination growth, is the covariate on social distancing measures captured by Google mobility trends, is the contemporaneous testing growth rate, and are the spatially uncorrelated stochastic errors.

Unlike prior attempts to model the spatial diffusion of the pandemic (Guliyev, 2020; Krisztin, Piribauer, and Wögerer, 2020), equation (16) excludes contemporaneous spatial lag of the dependent variable and the associated endogenous interaction effect (i.e., ) but instead incorporates a spatio‐temporal lag of the dependent variable along with an endogenous interaction effect (i.e., ) on the grounds that the incubation period of the virus ranges between 2 and 10–14 days with an average incubation period of about 4–5 days (WHO, 2020; CDC, 2021). Failure to account for this lag might result in a misspecification.

Research shows that individuals become fully immunized two weeks after the second dose, and the optimum gap between the doses is about eight weeks. The gap between the first and the second dose has shortened in Turkey down to 4–6 weeks following the emergence of Delta variant. We therefore conjecture that the lag between incidence and vaccination growth should be at least four weeks in order to detect this effect. Similarly, Google mobility trends are lagged by two weeks, which is equivalent to the upper end of the incubation period.

For the fact that the spatial units are not randomly sampled and are exhaustive of the population, spatial fixed‐effects model is more appropriate than random effects (Elhorst, 2012). Province fixed effects capture time‐invariant but space‐varying confounders such as population density, education, income, health literacy, age composition, hospital/ICU/health workforce capacity that are unlikely to change over our sample period.

Graphical models are great tools for tractability and transparency. They convey prior knowledge and help gain insight into the problem. The graphical model leading to equation (16) is given in Fig. 8. For simplicity, the diagram depicts spatio‐temporal dependence with only three spatial units along with signed edges. Covid‐19 incidence growth of each spatial unit depends on its past (purple edges representing the temporal diffusion rate), as well as other spatial unit's temporal lag (red edges representing the spatio‐temporal diffusion rate). The diffusion in both terms is incidence‐aggravating. Notice that the interaction between spatial units does not directly occur through vaccination or social distancing. In our proposed model, the channel through which vaccination growth in a spatial unit at time affect Covid‐19 incidence growth of other spatial units at time is through vaccination affecting its own Covid‐19 incidence growth, which spills over to other spatial units in the following period.

Figure 8

Path diagram among three spatial units with spatio‐temporal dependence. Note:

Within each spatial unit, the spatial fixed effects, shown by the gray nodes, are correlated with observed covariates and therefore act as confounders. The exposure variable, lagged vaccination growth, shown by the teal nodes, directly affects the contemporaneous Covid‐19 incidence growth. Graphically, we are interested in the (direct) effect of vaccination growth (thick teal edges). This pharmaceutical effect is expected to be negative. Similarly, increased social distancing, shown by the yellow nodes, is expected to reduce incidence growth. Anecdotal evidence in Turkey suggests a loosening of social distancing behavior or mask wearing following full vaccination despite having the same virus load as nonvaccinated. This is shown by the negative sign of the edge . Therefore, social distancing acts as a mediator between vaccination and Covid‐19 incidence for each spatial unit. The indirect effect of vaccination on incidence, through social distancing, is the product of the sign of the edges that constitute the path , which is positive. On the other hand, the direct effect of vaccination is negative. This suggests that the total effect of vaccination, direct and mediated through social distancing, is theoretically ambiguous. The direction of the total effect therefore depends on the relative size of the pharmaceutical and the indirect effect. The pharmaceutical effect of vaccination should be very large relative to moral hazard induced by vaccination so that the total effect is negative. What this reasoning shows is that lack of social distancing induced by gaining immunity undermines vaccination efforts. Finally, based upon the spatial SIR model, we include time‐varying but space‐invariant testing growth rate as the major predictor of Covid‐19 burden, shown by the light blue node.

The interpretation of point estimates in spatial models may lead to erroneous conclusions. Therefore, drawing inferences based on the partial derivative of the effect of a regressor on the outcome variable is a more valid approach (LeSage and Pace, 2009). Following Elhorst (2014), equation (16) can be rewritten to assess direct and indirect effects as: where is the identity matrix and captures the province‐fixed effects and the errors.

Our ultimate interests are of several folds. First, we would like to identify the direct effects of vaccination growth on Covid‐19 incidence growth. In the SR, defined as the period in which both the temporal and the spatio‐temporal diffusion is absent, this direct effect is the diagonal elements of the following matrix of partial derivatives, which is simply the parameter estimate of the impact of vaccination growth ( ) in equation (16): where is an identity matrix of order . Similarly, in the LR, defined as the period in which both the temporal and the spatio‐temporal diffusion are accounted for, this direct effect of vaccination growth on Covid‐19 incidence growth is the diagonal elements of the following matrix of partial derivatives:

Second, we would like to assess empirically whether vaccination exerts a positive externality over space, that is, whether an increase in vaccination growth rate in a particular province reduces Covid‐19 incidence growth over time in all other provinces. This indirect effect can be assessed in the SR and the LR. Since equation (16) does not contain both the contemporaneous spatial lag of Covid‐19 incidence and the contemporaneous spatial lag of vaccination growth, SR indirect effects do not occur (Elhorst, 2014). However, the LR indirect effect of vaccination growth on Covid‐19 incidence growth is the average of either the row or the column sums of the off‐diagonal elements of the RHS of equation (19). In the LR, the sum of the direct and indirect effects constitutes the LR total effect of vaccination growth on Covid‐19 burden.

In the following section, we estimate equation (16) using distance‐based and socioeconomic spatial weight matrices; provide evidence that the SAR with time and space‐time lagged dependent variable is the appropriate specification and compute SR and LR effects of vaccination.

Results

The estimation results of equation (16) are reported using a distance‐based and a socioeconomic spatial weight matrix, respectively in columns (1) and (2) of Table 2. As argued in Leenders (2002) and Elhorst (2010), we can only have a subjective view, based on theory at hand, as to which spatial weighting scheme is the most appropriate because the spatial data generating process and spatial associations are unknown. Instead, as often conducted by empirical researchers, we leverage the impact of a host of spatial weighting schemes on our results by conducting a robustness check with several well‐known specifications. These robustness checks, explored in section A of the Appendix, show that the choice of the spatial weighting does not significantly alter our inferences. Additionally, an extension replacing vaccination growth () by the natural log of vaccination rates () is assessed in section B of the Appendix. Given a fairly large time dimension, estimating equation (16) requires that is stationary over time. We show, through a battery of panel unit root tests, that specification in equation (16) is a stationary process. Due to space limits, these tests are relegated to section C of the Appendix.

Table 2

Effects of vaccination growth on Covid‐19 incidence growth, one‐way fixed‐effects SAR with time and space‐time lagged dependent variable

Outcome variable	Growth rate of Covid‐19 incidence
	(1)	(2)
Spatio‐temporal dependence (one‐week lag)
Temporal lag of incidence growth (ϕ)	0.152*** (0.056)	0.153** (0.061)
Spatio‐temporal lag of incidence growth (ψ)	0.506*** (0.053)	0.504*** (0.065)
Spatial lag of incidence growth (ρ)	Constrained to zero
Pharmaceutical intervention (four‐week lag)
Vaccination growth rate	−0.814*** (0.090)	−0.825*** (0.088)
Social distancing behavior (two‐week lag)
Workplace mobility	0.360*** (0.056)	0.329*** (0.059)
Retail and recreation mobility	0.137* (0.073)	0.154** (0.078)
Grocery and pharmacy mobility	−0.162*** (0.051)	−0.153*** (0.052)
Park mobility	0.073** (0.033)	0.050 (0.038)
Transit station mobility	−0.175*** (0.052)	−0.163*** (0.058)
Testing capacity
Testing growth rate, national	1.848*** (0.094)	1.847*** (0.095)
Log‐pseudolikelihood	−193.40	−196.34
R‐squared
Overall	0.4688	0.4666
Between	0.0961	0.0662
Within	0.4805	0.4786
Akaike information criterion	397.77	403.56
Chi‐square test (ϕ=ψ=0)	381.83***	523.77***
LR test
Time versus time and space‐time lag	131.1***	125.3***
Space‐time versus time and space‐time lag	4.06**	3.98**
Residual Moran's I	0.0143***	0.0093**
Spatial weight matrix	Inverse distance	Socioeconomic

Note: Effective number of observations is 1620. The unit of observation is the province‐week. Both specifications include spatial fixed effects. All mobility measures are defined as percentage change from the baseline. See section “Exploratory spatial data analysis” for the details of spatial weight matrices. Standard errors in parentheses are clustered at the NUTS‐2 level (26 regions). LR tests respectively report the chi‐square for the null hypothesis that a SAR with time and space‐time lagged dependent variable can be constrained to a SAR with time‐lagged or space‐time lagged dependent variable.

*, **, and ***denote statistical significance at the 10%, 5%, and 1% levels, respectively.

For both spatial weighting schemes in Table 2, the coefficient on the contemporaneous spatial lag of incidence growth () is constrained to zero on the grounds that the transmission of the virus takes an average incubation period of at most about a week. However, both the temporal and the spatio‐temporal diffusion rate (i.e., and ) estimates show that the spatial spillover of Covid‐19 incidence growth and the intra‐province inertia in growth increase the speed at which Covid‐19 incidence propagates in the current period. In Table 2, the temporal diffusion rate of Covid‐19 incidence growth is 0.152 and 0.153 while the spatio‐temporal diffusion rate is 0.506 and 0.504, respectively, for the inverse distance and socioeconomic spatial weighting schemes.

A model comparison can be performed using the maximum log‐likelihood criterion (Stakhovych and Bijmolt, 2009; Elhorst, 2010). The log‐pseudolikelihood for both weighting schemes in Table 2 are significantly larger than those reported for alternative weight matrices (See Table A.1 of the Appendix), indicating that the spatial associations are accurately captured by the inverse distance and socioeconomic spatial weighting, not only theoretically but also statistically.

We test the null hypothesis that and are jointly statistically indistinguishable from zero. Rejection of the null indicates that the SAR model with time‐ and space‐time lagged dependent variable is an appropriate specification. Otherwise, the model boils down to a possibly nonspatial and/or a dynamic model. For both models in Table 2, the Chi‐square test indicates that the temporal and the spatio‐temporal diffusion rates are jointly statistically significant. Further, the likelihood ratio (LR) test suggests that a SAR with time and space‐time lagged dependent variable cannot be constrained to a SAR with either time or space‐time lagged dependent variable at conventional test levels.

In line with our expectations, four‐week‐lagged vaccination growth has a negative impact on incidence growth across both specifications. With respect to social distancing behavior captured by two‐week lags, higher workplace mobility (along with retail and recreation and parks) translates into faster growth of incidence, despite being the area most easily monitored. On the contrary, grocery and pharmacy and transit station mobility are negatively associated with Covid‐19 incidence growth. The negative effect of transit stations on Covid‐19 owes to mask requirements on transportation and the personal code (HES code) introduced by the Ministry of Health to reduce the presence of positive‐tested individuals in congested public places. The negative effect of grocery and pharmacy is attributed to significant shifts from on‐site to online grocery shopping. Finally, testing capacity, as the major predictor, increases Covid‐19 incidence growth.

In order to assess the impact of PI on Covid‐19 burden, we hone in on the SR and LR effects, displayed in Fig. 9 along with a 95% confidence interval. Since neither specification contains both the contemporaneous spatial lag of Covid‐19 incidence and the contemporaneous spatial lag of vaccination growth, SR indirect effects are zero by construction. However, SR direct effects do occur. The effect of % point increase in vaccination growth at a particular province on Covid‐19 incidence growth in the same province is and in the SR. In the LR, this direct effect is and , respectively for the inverse distance and socioeconomic spatial weighting schemes.

Figure 9

Marginal effects of vaccination growth. Note: All effects are derived from those reported in Table 

The LR indirect effect is of special importance for it shows whether vaccination exerts positive externalities over space. For both spatial weight matrices employed in Table 2, vaccination growth curbs Covid‐19 incidence growth in the LR through creating positive externality with an LR effect of and . The LR total effects are displayed in the far right in each plot in Fig. 9. It captures the impact of vaccination growth at a particular province on Covid‐19 incidence growth in the same province as well as its impact via spatial spillovers. The LR total effect is and , depending on the spatial weight matrix.

Fig. 10 shows the extent of spatial autocorrelation in the residuals, respectively for columns (1) and (2) of Table 2 along with a 99% confidence interval. The overall Moran's I statistic indicates that the null hypothesis of spatial randomness can be rejected in favor of spatial clustering for the inverse distance but cannot be rejected for the socioeconomic weighting scheme at a 1% significance level. Nevertheless, the residual spatial autocorrelation is extremely weak in both cases (0.0143 for inverse distance and 0.0093 for the socioeconomic weight matrix).

Figure 10

Residual spatial autocorrelation. [Colour figure can be viewed at wileyonlinelibrary.com].

The weak but persisting spatial autocorrelation in the residuals might be an indication of spatial heterogeneity and call for methods such as geographically weighted regression (GWR). Although GWR is able to exploit spatial heterogeneity and may reduce residual spatial autocorrelation over global models, it is likely that the comprehensiveness and complexity of the model considered in this article and the longitudinal nature of our data leave very little room for improvement through GWR‐type approaches.

Conclusion

This study uses provincial panel data for 81 provinces in Turkey over the period of February–August 2021 to leverage the spatio‐temporal impact of vaccination on Covid‐19 incidence. Our empirical strategy is grounded in a spatial SIR model with vaccination and testing. The epidemiological model shows that Covid‐19 incidence growth depends on the weighted average number of new Covid‐19 cases at spatially proximate locations per every positive case at a given location, testing capacity and confounders. Using our prior knowledge, a path diagram is developed to add transparency and tractability to the analysis. Together with the spatial SIR model, it serves as a basis for the estimation of an SAR model with time and space‐time lagged Covid‐19 incidence growth.

There is significant inertia and spatio‐temporal diffusion in Covid‐incidence growth in Turkey, even after accounting for vaccination, testing capacity, social distancing behavior, and unobservable confounding. These two distinct sources of Covid‐19 diffusion act as a double‐edged sword. On one edge, a steady diffusion in space and time prolongs the pandemic and exacerbates its detrimental public health effects. On the other edge, it conduces positive externalities of vaccination over space due to spatial spillover of Covid‐19. The combined size of temporal and spatio‐temporal diffusion falls short of that of pharmaceutical effects. This justifies rapid vaccination to outrun the burden.

The major takeaway of this study is that the speed at which individuals are vaccinated reduces the speed at which Covid‐19 propagates, directly through its pharmaceutical effects via reduced risk of contraction and duration of infectious period, and indirectly through beneficial spatial spillovers. These effects persist not only for spatially proximate locations but also for distant ones. Conversely, the slower the vaccination efforts, the greater the threat to surrounding regions. Therefore, policy actions should be regionally cohesive and harmonious and a tremendous effort should be spent on boosting the speed of vaccination across all regions of Turkey despite occasional fluctuations in vaccine supply.

Several factors create barriers to a full‐fledged containment policy. Covid‐19 crisis has been a global concern since the late 2019. While maintaining public health is the primary objective, vaccination efforts are threatened by unequal access, lack of compliance with effective treatment, and the emergence of virus variants. For instance, Gamma variant is 1.7–2.4 times more contagious than wild type (Faria et al., 2021) and Delta variant is the most contagious of all its predecessors (Mahase, 2021). Table 3 displays the efficacy of authorized vaccines in Turkey by variant. Although a direct comparison of these figures is not meaningful due to differences in methods, samples, and lack of transportability of the effects to other geographies, clinical evidence from several countries shows that inactive vaccines are less effective than mRNA. The only phase 3 trial conducted in Turkey reports that the efficacy of Pfizer‐BioNTech (95%) is higher than that of CoronaVac (83.5%) against the wild type (original SARS‐CoV‐2) (Tanrıöver et al., 2021). Similarly, the efficacy of a given vaccine is lower relative to the wild type for emerging variants such as Delta. This suggests that the vaccination rate necessary to keep the effective reproduction number close to zero and end the pandemic in the pre‐Delta is no longer enough in the post‐Delta period.

Table 3

Variants and Vaccine Efficiency

Verified variants in Turkey by WHO	Efficacy of authorized vaccines (%)
	Pfizer‐BioNTech	CoronaVac
Alpha (B.1.1.7)	93.7 *, 90 †	71−91 ‡
Beta (B.1.351)	72−75 ‡, 85 †	65.9 †
Delta (B.1.617.2)	88 ‡	N/A
Gamma (P.1)	85 †	65.9 †

Bernal et al. (2021).

Cevik et al. (2021).

Mahase (2021).

The effectiveness of vaccines or the minimum effective scale thereof is lower due to production technology, anti‐vaccination groups, and virus variants. This will keep the effective reproduction rate above the level needed to contain the pandemic, at least in the short run, during which neither of these threats can be addressed promptly. Therefore, the virus keeps spreading despite vaccination efforts. Consistent with this insight, our short‐run estimates suggest attenuated direct effects of vaccination growth relative to those of long run. As a result, drastic NPIs may have to be reintroduced in the short run to compensate for transitory setbacks.

This study is not devoid of limitations. The first is the imputation of a limited number of missing values of social distancing measures. Relying on the homogeneous nature of provinces within NUTS‐2 regions in terms of economic conditions, labor market outcomes, and production structures, missing data on provincial Google mobility measures are imputed by the spatial median of the remaining provinces in the same region and period. Although other strategies might yield different imputed datasets, our imputation method is simple, contextually admissible, and cost‐justified relative to complex ones. Second, Covid‐19 modeling efforts are undercut due to underreporting of cases and testing data (Uçar et al., 2020). Entrenched data problems exist on provincial Covid‐19 deaths, which are not publicly available. Although attempts are being made to extract this information via excess mortality measures for a small fraction of provinces in Turkey (Yaman, 2021), structural issues remain as to how excess mortality should be measured and attributed to Covid‐19 pandemic.

Appendix S1 Supplementary Material

Click here for additional data file.