We use a combination of Geographical Information Systems (GIS) in ArcGIS Pro 3.1 and statistics in Stata 18 and R. The combination of GIS and statistics enables us to analyze correlations between environmental parameters and COVID-19 death incidence per 100 000 inhabitants for each county, and to present the parameters and resulting correlations spatially on maps.
Data description
We use COVID-19 death incidence data22, the continentality index23, demographic data with percentage of ages 65 and older19, the index of urbanization, the SES index, and a county boundary dataset for all counties of the US24, and integrate these data in ArcGIS Pro 3.1 for preparation and visualization, and for statistical analysis.
Johns Hopkins cumulative total county-level COVID-19 deaths from Jan 22, 2020 through Mar 31, 2022 were obtained from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University22. The data display the 7-day average death numbers. Totals were calculated from the Johns-Hopkins time-series data and converted to death incidence/100 K. For visualization, the resulting table was joined with the attribute table for the administrative borders24 of all counties (Fig. 1).
Total (accumulated) county-wise death incidence for COVID-19 per 100 000 inhabitants, from Jan 22, 2020 to March 31, 2022. Data from JHU CSSE COVID-19 Data, https://github.com/CSSEGISandData/COVID-1922. The State of Nebraska, with a deviating death incidence to neighboring states, is indicated with thicker borders.
In our calculations we used Continentality index data from the CMCC-BioClimInd data set23. This dataset provides data of 35 bioclimatic indicators from historical and future climate simulations. The historical data (1960–1999) employed in our analyses have a spatial resolution of 0.5° and are based on WATCH reanalyses. In this dataset the continentality index is calculated in its simplified form23, this differs from the one developed by Driscoll and Yee Fong (1992) and is calculated as the difference (in °C) between the mean temperature of warmest and of the coldest months of the year. The simplified continentality index used here differs from Driscoll and Yee Fong’s formula by measuring only the difference between the mean temperatures of the warmest and coldest months, without additional weighting factors. This straightforward method relies on monthly temperature data, ensuring robustness and consistency across regions. By directly capturing the seasonal temperature amplitude, it provides a clear and intuitive measure of continental influence. Its simplicity facilitates large-scale spatial comparisons and reduces calculation errors. Therefore, it is well suited for broad regional analyses where data comparability has a strong weight.
The percentage of over 65 year olds in the population, obtained from the U.S. Census Bureau’s American Community Survey from 201919 is a measure that gives the share of the over 65 year olds in the population for each county. We adjusted for the percentage of over 65 year olds as older people were overrepresented in the COVID-19 death statistics9 and also tend to live in coastal areas, like Florida19, where the continentality index is low, to a greater extent. Emerging evidence also suggests that depressive symptoms in middle-aged and elderly individuals significantly elevate the risk of cognitive decline and dementia, which may indirectly exacerbate COVID-19 related mortality risks through compounded health vulnerabilities25.
The degree of urbanization provides the relationship between the population living in urban (and rural) areas and the total population of the municipality26, with other words, population accumulations, or areas of crowding, are identified27 (Supplementary file 1).
The SES index combines the five factors poverty, unemployment, the housing cost burden, no high school diploma and no health insurance28 (Supplementary file 2).
GIS visualization and data preparation
The average continentality value for each county was obtained by spatial join of the county administrative borders file with the continentality raster file23 (Fig. 2).
Average continentality index for each US county (modified from Noce et al.23).
The spatial distribution of the average index of urbanization value for each county was obtained in the same way, from the raster provided by Eurostat29.
The percentage of the population aged over 65y for each country was obtained in tabular form from the US Census Bureau and joined to the administrative border attribute table for visualization (Fig. 3).
The SES for each country was obtained in tabular form from the center of disease control (CDC)28 (Supplementary file 2).
The resulting table of all joins, containing COVID-19 death incidence per 100.000 citizens, the average continentality value, the average index of urbanization value, the percentage of over 65-year-olds of total population and the SES value for each country was downloaded and used for the statistical analyses. The resulting coefficient was joined back to the attribute table in ArcGIS for visualization.
Statistical analysis
We fit linear regression models with regional fixed effects in order to estimate the relationship between continentality index and COVID-19 death incidence (Fig. 4) and visualized the values in Fig. 5 (Table: Supplementary file 3).
Coefficient plot of Continentality index and COVID-19 death incidence with Continentality index 11 as the reference category. Dots represent coefficient estimates, and the bars reflect 95% confidence intervals. Continentality values 39 and 41, that deviate from the trend, only occur in one county each and cannot be regarded as statistically robust. See visualization of coefficient across the US in Fig. 5.
Spatial visualization of Coefficient per US county of Continentality index and COVID-19 death incidence. The coefficient shows how differences in the continentality index are associated with changes in COVID-19 death incidence within the same U.S. Census region, after adjusting for age structure, SES, and crowding. The lower the coefficient, the lower the relative probability of dying in COVID-19.
To estimate the relationship between the continentality index and COVID-19 death incidence, we specify a fixed effects regression model that accounts for the nested data structure, where counties are grouped within broader geographic areas. The dependent variable, \({y}_{ic}\), represents the COVID-19 death incidence measured at the county level c, within area i. The main explanatory variable is the continentality index, also measured at the county level. In addition, we include a vector of control variables \({X}_{ic}\), which capture other observable county-level characteristics that may be associated with the outcome.
The regression model is specified as follows:
$${y}_{ic}={\beta }_{1}{Continentality}_{ic}+{X}_{ic}^{\prime}{\beta }_{2}+{a}_{i}+{\varepsilon }_{ic}$$
In this equation, \({Continentality}_{ic}\) is the variable of interest, \({X}_{ic}\) is a vector of county-level controls, and \({\beta }_{2}\) is a corresponding vector of coefficients. The term \({a}_{i}\) denotes area fixed effects, which absorb all shared unobserved heterogeneity at the area level that may confound the relationship between continentality and COVID-19 death incidence. The error term \({\varepsilon }_{ic}\) captures unobserved variation in the outcome at the county level.
The inclusion of area fixed effects \({a}_{i}\) controls for any persistent characteristics that are common to all counties within a given area—such as regional institutions, historical development patterns, or baseline policy environments—thereby ensuring that the estimated effect of the continental index is identified solely from variation within areas. That is, the fixed effects remove between-area variation from the estimation, allowing the model to compare counties with different values of the continental index but within the same broader area.
The coefficient of interest, \({\beta }_{1}\), thus captures the average within-area association between the continental index and the outcome. It reflects how differences in the continental index across counties within the same area are associated with differences in the outcome, conditional on the included control variables and unobserved area-level characteristics. In other words, \({\beta }_{1}\) estimates how much the outcome changes, on average, when the continental index differs between two counties in the same area, all else being equal.
To account for potential spatial or institutional correlation in the error structure, we cluster standard errors at the area level. This corrects for arbitrary forms of heteroskedasticity and autocorrelation within areas and ensures valid statistical inference in the presence of within-area dependencies.
In addition, we assessed multicollinearity among covariates using the Variance Inflation Factor (VIF) and examined the normality of residuals. Neither diagnostic indicated any concerns.
The inclusion of regional fixed effects allows us to adjust for all factors shared within a region; as such the parameters are estimated exploiting within region variation. This approach is ideal since we adjust for all unobserved factors that are fixed within each region (as defined by the US Census Bureau as South, Midwest, West, and Northeast), such as labor market differences and regional variation in political ideologies. We include a parsimonious set of confounders such as degree of urbanization (crowding), SES index, and the share of the population over 65 years of age in each county.
All analyses used continentality index 11 as the reference category (arbitrary, as the reference category is needed but has no influence on the outcome), so all coefficients in the continentality index correspond to deviations from counties that experience an 11-degree difference between the highest and lowest temperatures in a year. We also conducted robustness checks where we exclude the state of Nebraska, as there appear to be inconsistencies in reporting of COVID-19 death incidence. Nebraska has a noticeably lower death incidence than the surrounding states of Kansas and South Dakota, and in line with California, a state with one of the lowest death incidences, as shown also in other diagrams of the Johns Hopkins data for the states of the US, for example on the webpage 91-divoc.
We did a robustness check and calculated the coefficient by omitting Nebraska. The changes were marginal; for the results presented in this paper, we used the data as they are, including Nebraska.
We chose to normalize our data against the most known factors for the increase in COVID-19 death incidence: the socio-economic status, the age group 65 and older, and crowding. These are all factors within the Social Vulnerability Index (SVI)28. We chose against normalizing for the entire SVI, as the SVI also includes measures that are not known to influence COVID-19 death incidence, such as the population aged 17 and younger. However, we tested normalizing against the entire SVI, and the statistical differences in the spatial pattern were not significant (Supplementary file 4), confirming a lower COVID-19 death incidence in areas with lower continentality and a non-linear rise, with highest death incidence in the areas with highest continentality.
We did not normalize our data against other parameters that likely had local impacts on COVID-19 death incidence, such as distancing rules and vaccinations. These are parameters that varied in space and time across the US states during the years of the pandemic30, and suitable spatial data for the regional-scale perspective of our study does not exist. However, in normalizing our data against the SES, we covered some of these parameters, as socio-economic status, which includes access to information, access to health care, and housing conditions, shows a positive correlation with follow stay-at-home orders31 and getting vaccinated32.
