Trade, uneven development and people in motion: Used territories and the initial spread of COVID-19 in Mesoamerica and the Caribbean

Mesoamerica and the Caribbean form a region comprised by middle- and low-income countries affected by the COVID-19 pandemic differently. Here, we ask whether the spread of COVID-19, measured using early epidemic growth rates (r), reproduction numbers (Rt), accumulated cases, and deaths, is influenced by how the ‘used territories’ across the regions have been differently shaped by uneven development, human movement and trade differences. Using an econometric approach, we found that trade openness increased cases and deaths, while the number of international cities connected at main airports increased r, cases and deaths. Similarly, increases in concentration of imports, a sign of uneven development, coincided with increases in early epidemic growth and deaths. These results suggest that countries whose used territory was defined by a less uneven development were less likely to show exacerbated COVID-19 patterns of transmission. Health outcomes were worst in more trade-dependent countries, even after controlling for the impact of transmission prevention and mitigation policies, highlighting how structural effects of economic integration in used territories were associated with the initial COVID-19 spread in Mesoamerica and the Caribbean.


Introduction 37
COVID-19 is a pandemic disease caused by SARS-CoV-2, a novel coronavirus of zoonotic origin 38 [1]. In September 2020, half a year into the pandemic, as many COVID-19 cases had been reported in 39 the Western Hemisphere as in the rest of the planet [2] , with the United States reporting the largest 40 number of deaths globally, around 200,000 [3]. Patterns of high transmission were similar to those 41 previously observed in Europe and East Asia after the first global wave of transmission in early 2020 [3, 42 4]. For COVID-19, a previous analysis, restricted to Costa Rica, Panamá and Uruguay, suggested that 43 early epidemic growth may plateau faster in more equitable societies [5]. Thus, transmission patterns 44 may vary across space, associated with a region's uneven development in interrelated social, economic, 45 environmental, and health realms. This uneven development is partially reflected in differences in socio-46 economic indicators among the countries and territories of the region. Indeed, income inequality has 47 been widely recognized as a major driving force of adverse health outcomes [6]. Associated mechanisms 48 of action include the reduced funding of public services [7], such as health care and education. This, in 49 turn, can lead to the erosion of social organizations (which can offer mutual benefit via social 50 cooperation) and of "social capital" [7]. High levels of social cooperation and solidarity can lead to the 51 participation of diverse, often excluded, groups in inclusive policy and decision making [8]; more 52 vulnerable groups are more likely to actively engage in health and social campaigns for health 53 promotion when their views are considered [9] and those views are more likely to be noticed when 54 groups are politically organized [10]. 55 However, beyond income inequality and the erosion of social cooperation, solidarity and 56 political engagement, one of the largest problems that societies face is lack of equitable access to 57 essential resources for the human development of those societies or for their social, cultural and 58 biological well-being [11], despite the abundance of resources that could be leveraged toward such 59 J o u r n a l P r e -p r o o f exponential growth curve; and when 0 < p < 1, the model has sub-exponential growth. The equation 150 presented in (1) has the following solution [21]: 151 ( ) = ( + ) (2) 152 Where A depends on the initial condition C(0) and m > 1 is related to p as follows: 153 The parameter m is expected to decrease as COVID-19 cases decreased in response to transmission 155 reduction policies, while parameters r and A likely reflect structural differences across health systems in 156 Mesoamerica and the Caribbean. We fitted the model presented in (2) by maximum likelihood 157 estimation. We used the mle2 command from the R package bbmle [25]. We fitted the model using a 158 Poisson maximum likelihood function, where the initial parameter search was done using a Genetic 159 Algorithm [26]. Then a local optimization was performed using the . We 160 estimated confidence intervals by profiling the likelihood for each parameter [25]. Estimated 161 parameters for each country are presented in supplementary online Table S3. 162

Rt estimation based on the serial interval 163
We examined the time-varying reproduction number (Rt), which, similar to R0, measures the 164 expected number of new cases generated by a given case, but does not assume that the whole 165 population is susceptible to the disease. To estimate Rt given the above, the only additional information 166 necessary is the serial interval of the pathogen under study, which is the time distribution between 167 symptom onset in an initial case and symptom onset in secondary cases linked by contact tracing to that 168 initial case [22]. We based our serial interval on an estimate from China [28], which was based on the 169 largest sample size available for our period of study and assumes a 10 percent pre-symptomatic 170 9 transmission rate before the originating infected person shows any symptoms [28]. To account for 171 potential uncertainties, we took full advantage of the parameter estimates presented by Du et al. 172 [28].These estimates allowed us to set the 95% credible intervals from their estimates as the minimum 173 and maximum values for the mean and SD of the serial interval, while also setting a standard deviation 174 of 0.2 [5] and employing the estimation method developed by Thompson et al [22]. We reduced 175 potential biases in the estimation of Rt by considering imported cases as potential sources of new local 176 cases, but not as result of local transmission [22]. We estimated Rt at 25, 50, 75 and 100 days after 177 detecting the first case, using the EpiEstim R package [22]. As noted, we chose this method given that its 178 estimation is just based on knowing the serial interval of the pathogen [22]. Our approach thus differs 179 from methods that have more assumptions and require additional parameters for the estimation of Rt 180 [29,30]. Rt estimates are presented in the online supplementary Table S4. 181

Data on COVID-19 prevention and mitigation policies 182
We also collected data about policies directed to prevent, mitigate and suppress SARS-CoV-2 183 transmission during the study period, as these might have created heterogeneities in transmission 184 patterns. For this, we used the coronavirus government response tracker developed by the University of 185 Oxford (Table S2) to count the number of policies in place before the detection of the first COVID-19 186 case and 25 days after the first case was detected in each territory. We specifically included the 187 following preparation and prevention policies: the presence of public information campaigns, a defined 188 testing protocol and a contact tracing protocol. To account for mitigation and suppression, we included 189 the following: school closing, workplace closing, cancelation of public events, restrictions on gatherings, 190 public transport closing, stay-at-home requirements, restrictions on internal movement and 191 international travel controls. Figure 2 shows the timeline for the different policies and the detection of 192 the first case across the studied territories. 193

Human development, uneven development, and used territory data 194
Ancillary geographical data relevant to all tested social science theories 195 We included population size estimates to account for heterogeneities that might emerge from 196 different territories having differences in population size and density per area unit. Territory surface 197 areas (in km 2 ) were calculated from Natural Earth 1:50m scale Admin 0 geometries (Table S2) Table S2). The HDI is a "summary measure of average achievement in key 203 dimensions of human development: a long and healthy life, being knowledgeable and have a decent 204 standard of living" (Table S2). The HDI is then corrected by discounting a percent representing inequality 205 in the components (%HD) of human development, thus producing the IHDI (Table S2). We also included 206 the Universal Health Coverage (UHC) service index, a metric developed by the World Health 207 Organization (WHO) to measure access to "health services (including prevention, promotion, treatment, 208 rehabilitation and palliation) of sufficient quality to be effective while also ensuring that the use of these 209 services does not expose the user to financial hardship" (Table S2). UHC values range between 0 and 100 210 (Table S2). These variables quantify the role of human development in the initial spread and growth of 211 COVID-19 across Mesoamerica and the Caribbean. 212 Uneven development data 213 We compiled a series of economic indicators developed by the World Bank including the gross 214 domestic product (GDP) per capita, i.e., the average net economic value added by each resident in a 215 given territory, to measure the role of "uneven development" in the initial spread of COVID-19 across 216 Mesoamerica and the Caribbean. The economic indicators were measured in US$ and corrected by 217 purchasing power parity, a normalized value to account for differences in price levels across countries 218 (Table S2). We also included data on net bilateral aid flows from the United States of America and the 219 United Kingdom (Table S2) into the studied territories from Mesoamerica and the Caribbean. The net 220 bilateral aid flows measure the difference between government grants and loans and repayments of 221 those loans between those two pairs of countries, measured in millions of current US$. We consider 222 these aid flow variables could also provide a measure of uneven development since countries of the 223 periphery are more likely to receive aid from wealthier countries [12,31], and to that end, selected the 224 United States and the UK, the two main donors in the study area (Table S2). 225

Used territory data 226
To quantify the "used territory" impact in the initial spread of COVID-19 across Mesoamerica 227 and the Caribbean, we included data on the number of connected countries, connected international 228 cities and passengers traveling (last annual estimate available) through the main airport for each of the 229 studied territories. These details and data sources are presented in Table S5. Other variables were 230 related to international trade. We included merchandise trade openness, i.e., the ratio of exports plus 231 imports in merchandise goods to that of the territory's GDP (Table S2) as well as overall trade openness, 232 the sum of exports and imports of goods including services divided by a country's GDP from the World 233 Bank (Table S2). We also included additional variables on the composition and diversity of trade from 234 the United Nations Conference on Trade and Development (UNCTAD , Table S2). From the 235 "Merchandise: Product concentration and diversification indices of exports and imports, annual," table, 236 we included the number of products imported and exported at the three-digit SITC, Rev. 3 level. We also 237 included countries' trade concentration indices (normalized Herfindahl-Hirschman indices). For the 238 12 latter, values range between 0 and 1, where high values, i.e., close to 1, indicate that imports and/or 239 exports are concentrated in a few products (Table S2). Finally, we included UNCTAD variables related to 240 foreign direct investments (from the table "Foreign direct investment: Inward and outward flows and 241 stock, annual," Table S2), which are standardized as percentages of GDP. These measure the investment 242 of enterprises from a different country within a target country (in the context of our study, the 243 territories in Mesoamerica and the Caribbean) and can be measured as flows or stocks, which, 244 respectively, measure transactions and accumulated values held at the end of an annual period. For the 245 analysis, we employed the inward metric as it measures the impact of foreign investments within the 246 territories of Mesoamerica and the Caribbean (Table S2). Trade and foreign direct investment data are 247 key to understanding the role of the "used territory" as defined by Santos [20], since these variables 248 relate to the structural underpinnings for the movement of human beings [16]. We hypothesize these 249 variables might play a key role beyond the one directly attributable to human movement on the spread 250 of COVID-19. We note that trade relations are understood within several schools of thought to play 251 integral roles within processes of uneven and under-development and the constitution of the used 252 territory. Here, we group them with the latter for analytical convenience but discuss their intrinsic 253 connections in the conclusion. 254

Cluster Analysis 255
To characterize differences and similarities among the studied countries, we performed a cluster 256 analysis. In the context of this study, cluster analysis finds potential groups of countries with similar 257 transmission patterns or social science parameters. Briefly, a cluster analysis is a technique that looks at 258 the similarities in the characteristics of objects to separate them in groups, with objects being more 259 similar to members of their own cluster than to members of a different cluster. We specifically did an 260 agglomerative nesting hierarchical cluster analysis [32] of both epidemiological variables and variables 261 J o u r n a l P r e -p r o o f 13 associated with uneven and human development, as well as the used territory. In this type of clustering, 262 objects (here, countries) are systematically grouped in a bottom-up manner where more similar objects 263 are put within the same cluster [32]. 264 The cluster for epidemiological variables included parameters r and A estimated from the early 265 epidemic growth for all the territories considered in this study, Cases, deaths and Rt at day 25, 50, 75 266 and 100 after detecting the first infection in each of the studied territories. The clustering for variables 267 dealing with the used territory, human development, and uneven development included: HDI, %HD, 268 IHDI, UHC, GDP, GDP2, AUS, AUK, CN, CC, TO, TM, IM, EX, CE and CI (for explanations of these 269 parameters, see Table 1). We did not include PA, FI, FO, SI and SO due to these having two or more 270 missing values. Before performing the cluster analysis, we normalized all variables by subtracting the 271 mean and dividing by the standard deviation so distances between territories were not dominated by 272 the variables with the largest numerical scale [32]. In the cluster analysis for disease parameters, we set 273 unreliable Rt estimates (Table S4) to zero as wide confidence limits emerged given that no new cases 274 were reported for periods equal to or longer than the serial interval for disease transmission (Figure 1). 275 Setting these values to zero is an appropriate choice as disease transmission was not observed in the 276 days before the dates (t=50, t=75, t=100) for which Rt was estimated [22]. The zero imputation was 277 unlikely to have created artifacts that affected cluster estimation, as Rt values were likely near zero [33]. 278 To represent the results from the cluster analysis, we made dendrograms. We also present the clusters 279 in a bi-dimensional projection, based on a principal component analysis (PCA) of the variables 280 considered in the cluster analysis [34]. To ease the interpretation of clusters we also performed a 281 segmentation analysis of the clusters, where mean values for each parameter and each cluster were 282 estimated [35,36]. We then made parallel coordinates plots with the standardized values, i.e., removing 283 the mean and dividing by the variance of each variable considered in the cluster analysis [37], of the 284 estimated means for each cluster and parameter [38]. 285

Econometric Modeling 286
We quantified, using econometric linear models, the impact of variables measuring used 287 territories, human development, and uneven development, and policies to control SARS-CoV-2 288 transmission on early transmission COVID-19 epidemiological parameters. Prior to fitting the linear 289 models, we started with an exploratory analysis of our full dataset. We first examined the Pearson 290 correlation [39] among the epidemiological variables to select which ones to analyze ( Figure 3A). That 291 preliminary analysis showed that cumulative cases and deaths at the chosen dates were highly 292 correlated across the territories studied. Thus, we focused our analysis on cases and deaths at t = 100 293 days. Meanwhile, Rt at t = 25 was positively correlated with death and cases at all times and was the 294 only time at which we were able to reliably estimate this parameter for all the studied countries (Table  295 S4). 296 In contrast, we did not analyze Rt with econometric tools. The Rt estimates at t = 50, t = 75 and t 297 = 100 had between 6 and 8 unreliable estimates each (Table S4) and the correlation with other 298 epidemiological variables suggested that their patterns were either similar to those of cases and deaths 299 or parameters from the early epidemic growth model. We chose r as it was not associated with cases or 300 deaths on any selected dates when examining the exponential growth parameters. Moreover, r was 301 estimated for all territories, something not possible with m, which could not be estimated for Dominica, 302 Dominican Republic, St. Kitts and Nevis and Uruguay. Meanwhile, the covariates (Fig. 3B) had a 303 clustering pattern where human development variables, GDP measurements and foreign direct 304 investment stock (inward and outward) were highly correlated. 305 Similarly, variables related to travel through the main airport, CC, CN, and PA, among the studied 306 territories were highly correlated (Fig. 3B). Given that SI and SO were correlated with human 307 development variables and that two or more observations were missing, we did not consider these 308 variables in our analysis. Similarly, PA, FI and FO were not further considered as those variables also had 309 two or more missing observations. We then selected which variables to consider in the initial "full 310 models" based on the correlation of the covariates with r, Ct and Dt at t = 100. 311 For r, we used linear regression as the variables are assumed to have errors, denoted by , with 312 an equal, independent and normal distribution [40]. For r, the initial model was described by the 313 following general equation: 314 The model presented in (4) accounted for heterogeneities shaped by different population sizes by 317 including weights proportional to population size in the error matrix [40]. Since cases, C, and deaths, D, 318 at t = 100 were count variables, we employed generalized linear models for count variables, using a 319 Poisson distribution. For C and D at t = 100, the initial model was described by the following general 320  (4) and (5), is an intercept; for Dev[j], we alternatively 324 tested HD, IHD, GDP or GDP2, as all these variables were highly correlated ( Figure 3B). In these 325 equations, travel was either CC or CN, import[j] either IM or CI and export [j] either EX, if import[j] was 326 IM, or CE otherwise. In models for C and D, described by equation (5), the link() function is a natural 327 logarithm. Models had population size included as an offset, indicated by offset(log(POP)), to account 328 for heterogeneities in disease counts that emerge from the fact that cases in a given territory are 329 constrained by the population size of that territory [41]. The resulting 16 models described by equations 330 J o u r n a l P r e -p r o o f 16 (4) and (5) were simplified through backward elimination [42]. Backward elimination is an algorithm by 331 which nested models (i.e., models with one parameter less than the "full model" within which they are 332 nested) are compared through the Akaike Information Criterion (AIC) [42]. The model minimizing AIC is 333 selected at each step until a "best model" is found which is not significantly different in its explanatory 334 power from the "full model" but becomes significantly different from the full model if further simplified, 335 through the removal of additional covariates. The resulting best models were then compared using the 336 AIC, a metric that trades off model fit and the number of parameters in a regression [42]. The resulting 337 models were then subject to model diagnostics using the variance inflation factor, VIF, and factors above 338 a threshold value of ten (VIF>10) were removed [43]; we also removed variables that were not 339 significant. 340 For each model, we inspected diagnostic plots to assess homoscedasticity (examining the 341 residuals as a function of fitted values), normality (using a Q-Q plot) and assessing influential points on 342 parameter estimates by plotting residuals as a function of the leverage [40]. 343 To further ensure we did not overfit parameters, we also performed a jackknife [44], a method 344 for parameter inference based on refitting the model leaving one observation out at a time (loo). The 345 parameters fit at each iterative step of the jackknife(̂) are used to estimate average model 346 parameters from all the, i.e., n, loo models (̅ ̅̅̅̅ ). With the jackknife, bias-corrected standard errors for 347 each parameter can be estimated using the deviation of parameter estimates from all the loo models 348 according to the following formula equation [44]: 349 Which can then be used for inference against a normal distribution [45]. 351 We also performed a bootstrap [40], given this is considered one of the best methods for linear 352 model cross-validation [46]. For the linear regression presented in equation (4) [34]. We performed a parametric bootstrap for the count regression models represented by equation (5)  357 [47]. The fitted values were simulated as mean values from a count distribution, Poisson or negative 358 binomial. These simulated values were used to refit the model and build confidence intervals based on 359 the distribution of estimated parameters. We specifically used this method for the count regressions to 360 more accurately represent the original distribution of the count data [47], as models assumed that both 361 coronavirus cases and deaths (equation 5) were proportional to population size. All bootstraps were set 362 to 10,000 replications for each one of the best models. 363 Finally, we also evaluated the goodness of fit of the best models by estimating a predictive 364 coefficient of determination, PR 2 [48] which is an R 2 based on the reduction of the variance in model 365 predictions with data not used for model fitting, measured as the mean squared error (MSE, the mean 366 of the square difference between observed and predicted values), in relation to the variance (var) of the 367 response variable, i.e., r, C and D, according to the following equation [49]: 368 This model performance metric is used to assess prediction, in addition to fit [48]. 370

Results 371
Cluster Analysis 372 The cluster analysis results are presented in Figure 4. At least one cluster, formed by Panamá 373 and the Dominican Republic, was observed when examining the epidemiological parameters ( Figure 4A  374 and 4B). The rest of the clusters had some overlapping member territories as indicated by the color-375 coding. Epidemiological parameters ( Figure 4A) allowed us to define three clear clusters of transmission 376 patterns. Uruguay, as expected, served as an outgroup by not clustering with the rest of the studied 377 territories. In the cluster analysis, Panamá and the Dominican Republic, the two countries with the 378 highest total cases, as shown in Figure 1, formed a cluster. Reflected by the cluster branch's height, the 379 other two clusters showed more similarities. One of the clusters included Nicaragua plus territories that 380 have been relatively unaffected by the coronavirus (Figure 1) Guatemala, Honduras, Bahamas, Grenada and Jamaica. Here it is key to note that these patterns seem 385 to have arisen independently of the COVID-19 related health policies ( Figure 2) and are not an artifact of 386 the epidemiological parameters having different numeric magnitudes, given that all variables were 387 standardized. The goodness of fit of this cluster analysis was high, as indicated by an agglomerative 388 The same clustering patterns were observed when the epidemiological parameter clusters were 390 projected over a PCA-based bi-dimensional projection of the studied epidemiological variables (Figure  391 S2A). PCA loadings (Table S6) suggest that cases and deaths were larger in Panama and the Dominican 392 Republic. In contrast, the other two clusters were mainly separated by differences in Rt and epidemic 393 J o u r n a l P r e -p r o o f 19 growth parameters as we discussed above. This pattern can be more easily observed in a parallel 394 coordinates plot ( Figure 4C) which clearly shows that cluster segmentation was driven by the Dominican 395 Republic and Panamá having a larger burden of cases and deaths through time. In contrast, Uruguay and 396 the third cluster (including most of Mesoamerica as well as Cuba and most of the Greater Antilles 397 territories) differed on Rt and the early epidemic growth parameters, which in these two clusters were 398 larger than in the second cluster (including Belize and Nicaragua and several Caribbean island countries). 399 The mean values for each epidemiological variable and cluster are presented in Table S7. 400 Meanwhile, the cluster of variables associated with human development, uneven development 401 and used territory ( Figure 4B) showed patterns that revealed similarities and differences in compared 402 epidemiological patterns. Four territories did not cluster with any other: Haiti (the most different from 403 all other territories as revealed by the height of its divergence from the others), Cuba, Jamaica, and 404 Puerto Rico. The lack of clustering partially explains the agglomerative coefficient of 0.59, which is lower 405 than observed for the epidemiological parameters. The rest of the territories belonged to four clusters. 406 One of these clusters, related with Puerto Rico, included Panamá and the Dominican Republic, just as 407 observed for the epidemiological parameters. A second group was formed by territories that can be 408 seen as more developed since it included Uruguay (one of the most developed nations in Latin America), 409 Costa Rica, Barbados, Trinidad and Tobago, and The Bahamas. This second cluster was related to a third 410 cluster containing the remaining Caribbean territories that we studied. A fourth cluster included the rest 411 of the Mesoamerican territories. Inspection of the clusters over a PCA-based bi-dimensional projection 412 of the studied variables ( Figure S2B) confirmed this pattern. Based on the PCA loadings (Table S8), the 413 major factor separating the clusters was human development, which was high in the cluster containing 414 Uruguay, but also for Cuba and Puerto Rico, and was smaller for the cluster of Mesoamerican territories, 415 where human development was, nevertheless, higher than for Haiti. This pattern can be observed in the 416 parallel coordinates plot of Figure 4D. The lower connectivity with other places separated the Caribbean 417 J o u r n a l P r e -p r o o f 20 territories from the other clusters. The connectivity also appeared to be a dominant factor in separating 418 Panama and the Dominican Republic into a different cluster ( Figure 4D). The connectivity was measured 419 by the number of connected nations (CN) and connected international cities (CC) at the main airport of 420 each territory. Overall, the cluster analysis is suggestive that clusters based on epidemiological patterns 421 are similar to clusters based on variables related to development and the used territory, even if the 422 match was not perfect. The mean values for the variable and clusters of Figure 4D are presented in 423 Table S9. 424

Econometric modelling 425
A different way to look at the relationship between disease and social science theory 426 parameters is by performing an econometric analysis. We explored the sixteen potential models 427 explaining the early epidemic growth rate r. We found all candidates for the best model (Table S10) 428 considered concentrations of imports (CI), the number of international cities connected at the main 429 airport (CC) and the number of prevention (PP0) and mitigation (MP0) policies already in place when the 430 first case was detected locally. All VIF diagnostics suggested that variables were unidentifiable, implying 431 there was no reason to worry about problems with parameter estimation due to high correlation in the 432 predictors. The best model was not significantly different from the full model F9, 8=1.101, p>0. 451, 433 meaning fewer parameters explained as much of the variability in r as a model with more parameters. 434 The best model (Table 2) showed that the early epidemic growth increased with all the variables 435 considered (CC, CI, PP0 and MP0). Overall, this model successfully explained the sources of variability, 436 having a coefficient of determination above 60%. Specifically, the model had an ̂2 = 0.61, which was 437 reduced by one-third when using the model to predict observations not considered in the fit (̂2= 438 0.40). All parameters in this model had VIF below 10 (Table S10). The inspection of model diagnostic 439 plots suggests that there were neither outliers, significant deviations from normality, nor influential 440 J o u r n a l P r e -p r o o f 21 points that affected parameter estimates (Fig. S3). As expected, parameter estimates from the jackknife 441 were similar to the estimates obtained by ordinary least squares (Table S11). Confidence intervals from 442 the nonparametric bootstrap also suggest none of the parameters was overfitted (Table S11). 443 Meanwhile, the best model (Table 3) explaining the number of cases, C, at t = 100 days after 444 detecting the first local case dropped a large number of variables after backward elimination (Table S12) (Table S12). The jackknife for this model also found all significant and similar parameters 454 to those found with the negative binomial generalized linear model (Table S13). The parametric 455 bootstrap did not confirm the significance of ASK and TO, a result that might be an artifact of the 456 bootstrap asymptotically approximating the distribution of parameters at the cost of increasing type II 457 error, i.e., failing to reject the null hypothesis when false [40,44,47]. 458 The best model (Table 4) explaining the number of deaths, D, at t = 100 days after the first 459 locally detected case likewise dropped a large number of variables after backward elimination and 460 diagnostic checks (Table S14) previous models, assumptions were met ( Figure S5), so were VIF diagnostics (Table S12). The jackknife 466 for the deaths model confirmed the significance of CC and provided similar estimates for all other 467 parameters (Table S15). The parametric bootstrap, meanwhile, suggested that both CC and USAID were 468 significant while the lower 95% CI of TO was nearly significant, a result also suggested by the jackknife 469 (P<0.08) (Table S15). However, further simplifications to the model resulted in a significant loss of model 470 fit and different overdispersion parameters, a problem suggesting that all parameters presented in Table  471 4 were not overfitted [34]. These results highlight the limitation in negative binomial model selection 472 when overdispersion is not fixed. Overdispersion is a factor that makes reliable model comparison 473 impossible [42,47]. In addition, the satisfactory VIF (Table S14) and model diagnostics ( Figure S5) results 474 suggest that the model was not overfitted. 475

Relevance of social science theory to describing disease spread 477
The main common feature of our results is that variables related to social science theories about 478 used territory, human development, and uneven development were shown to be highly instrumental in 479 explaining the initial spread of COVID-19 in Mesoamerica and the Caribbean. As summarized in a 480 cartographically explicit way in Figure

The feasibility of used territories driving disease spread 496
These patterns can be expected under Santos' conception of the used territory as "made up of 497 objects and actions, and … a synonym for human space, inhabited space," unlike the traditional 498 definition of territory as the land of a state [20]. Santos insightfully argued spatial patterns can go 499 beyond "regional" by considering "horizontalities" as territorial links that emerge from absolute space 500 contiguity. Thus, horizontalities are akin to the rich observations of Tobler on spatiality [51,52], which 501 have been condensed and reinterpreted by others as Tobler's First Law of Geography, the principle 502 stating that 'nearby' phenomena are likely more interrelated [53,54]. In that sense, we can expect the 503 overlapping clusters of territories in Mesoamerica, the Lesser and Greater Antilles, that we described 504 across a spatially contiguous geography as illustrated in Figure 5. More generally, Uruguay illustrates 505 horizontalities, an outgroup, not clustering with Mesoamerica and the Caribbean when examining its 506 disease patterns. However, "verticalities", defined by Santos as forces connecting what is not contiguous 507 in absolute space, also shape used territories [20]. Verticalities can explain clusters like the one observed 508 for the Dominican Republic and Panamá, which, although not contiguous in absolute space, are 509 J o u r n a l P r e -p r o o f 24 connected by having similar social, political and economic relations. Similarly, when looking at indicators 510 related to people in motion, economic structure and trade, the presence of Uruguay in a cluster with 511 countries of Mesoamerica and the Caribbean can be partially explained by human development metrics 512 based on the 'capabilities approach' developed by Amartya Sen [11,55]. This clustering pattern can also 513 be articulated by Amin's theory of uneven development [12]. 514

Higher human development does not prevent early epidemic growth 515
Following Sen [11,55], we would expect that a higher level of human development would have 516 implied substantial protection against the spread of COVID-19 in Mesoamerica and the Caribbean. Yet 517 we did not find any human development metrics associated with the epidemiological parameters 518 considered a response in our econometric models. We found that r grew with the concentration of 519 imports, a variable positively correlated with human development indices ( Figure 3B). Indeed, the raw 520 correlations indicated a positive association of r with both the human development index (correlation = 521 0.21) and its inequality-adjusted version (correlation = 0.14). This positive association might reflect that 522 countries with higher human development were better able to diagnose SARS-CoV-2 infections. In the 523 studied region this might be the case for countries like Nicaragua, with low human development, that 524 have taken no action to contain, or even diagnose, . This, for instance, contrasts with the 525 laudable epidemic surveillance that Panamá has had in the region, being the country with the best 526 surveillance system in the studied region, the country most successful in diagnosing the virus at the 527 beginning of the pandemic, and the only one that was using molecular information to trace SARS-CoV-2 528 patterns of spread and biological evolution [5,57]. 529 On the other hand, this counterintuitive association, where more human development was 530 associated with larger epidemics, might reflect biases inherent to the human development index 531 construction. Given its components, human development is strongly associated with the gross domestic 532 J o u r n a l P r e -p r o o f 25 product per capita, as shown in Figure 3B. Although able to account for individual wellbeing at a 533 particular scale, the human development index was unable to reflect important societal-scale 534 vulnerabilities (even in an 'inequality-adjusted' form). The data shows that more human development 535 implied less resistance to a pandemic infectious disease such as COVID-19's spread across the studied 536 territories. This raises concern about the formulation of development metrics. It is desirable to reduce 537 the gap in material wellbeing across and within societies and this measure may correlate with some 538 dimensions of human possibility, as originally intended by Sen [11,55] impact on transmission as observed in early epidemic growth rates to having a negative impact on the 543 number of cases 100 days after the detection of the first case. And these association patterns are 544 unlikely to be spurious, given the robustness of our cross-validated results, which come from models 545 explaining at least 60% of the data variability, a value well above what is commonly reported in many 546 areas of the social sciences in which models explaining 35% of the variance may be exceptional [40, 42, 547 59]. 548

Uneven development articulated by used territories 549
Given the limitations of measurements for "human development," and after controlling for the 550 impact of COVID-19 prevention and mitigation policies, we can try to understand the verticalities 551 defining used territories by resorting to the principles of uneven development developed by Samir Amin 552 Similarly, other aspects of unequal trade were associated with increased transmission and 576 burden of COVID-19 in Mesoamerica and the Caribbean. We observed that the concentration of imports 577 increased early epidemic growth and deaths during the study period, likely, at least in part, reflecting 578 emergent structures of trade between core economies and regional blocks within a global political 579 economy of vast differences in power. These dependency mechanisms might explain observed patterns 580 in Caribbean countries like Dominica, whose major imports come from Japan, the major market for its 581 exports [72], but beyond bilateral patterns of trade; they might reflect realities observed in 582 Mesoamerica and the Dominican Republic following the Central American Free Trade Agreement 583 (CAFTA), where the elimination of tax barriers has eased importations from the United States in the 584 countries covered by that treaty [73]; and similarly, with Panamá, which has an independent agreement 585 with similar terms [74]. This promotion of economic specialization [75] has destroyed the ability of 586 peripheral states to satisfy needs for essential goods [76], including those essential to deal with 587 pandemics of infectious diseases through domestic production, as observed for COVID-19 and the 588 shortage of ventilators and other medical supplies in Latin America [77,78]. 589 Meanwhile, exports were associated with a higher number of cases; but when the concentration 590 of exports appeared in tandem with the concentration of imports, it reduced the number of deaths. An 591 explanation for this changing influence may be rooted in the role that exports play in the economy and 592 the shape of the relationship between export concentration and the wealth of nations, as measured by 593 the GDP per capita, which is roughly "U" shaped [79], and the fact that diversification does not always 594 lead to an economy autonomous from market shocks [80]. The non-monotonic relationship between 595 concentration of exports and GDP is, indeed, quite interesting as it might explain the positive effect in 596 the regression analysis if the following are captured: an increased susceptibility in countries whose 597 exports are concentrated by way of primarily extracting natural resources [81] and limitations to 598 diversification in states with small populations or territories [82]. But at the same time, the exports 599 effect on deaths might become negative when it accounts for a national specialization in producing 600 export goods that reflects not the typical core-led development pattern observed but a successful 601 societal mobilization around alternative development strategies that Amin described as 'delinking'. 602 Delinking is not a matter of a country pursuing autarky or rejecting the outside, but the country 603 moderating the terms of its external interactions as part of a reorientation of its social and economic 604 relations to serve the needs of its peoples [83]. A complete investigation of these issues is beyond the 605 goals of this study given numerous unresolved debates around how exports and their trade 606 concentration relate to the growth of wealth [84]. However, when not paired with the concentration of 607 imports, concentrations of exports strongly support the explanation that elements of uneven 608 development are verticalities constructing used territories that accelerated the spread of COVID-19 in 609 Mesoamerica and the Caribbean. 610 The relevance of used territories for the development of spatial theories based on relational 611 geographies 612 As discussed, the used territory offers unique advantages to study phenomena that couple 613 natural and human systems. However, its origin within the discipline of geography enables it to go one 614 step beyond by framing those connections within the spatiality of landscapes. These are landscapes 615 whose attributes, length and time scales, and even appropriate spatial coordinate systems for analysis 616 coevolve with coupled dynamical systems [85]. In that sense, beyond what is possible in human 617 development theory, and as complementary to human and uneven development theories, the used 618 territory is a theoretical tool that eases the understanding of relational geographies [86][87][88]. It 619 contextualizes the placement of people in motion. As measured by the number of international cities 620 connected through the main airport, the motion factor was positively associated with increases in the 621 early exponential growth rate of cases, in cases and in deaths by COVID-19 in Mesoamerica and the 622 Caribbean. Thus, our results illustrate the power of the used territory and of attention to relational 623 geographies generally, as conceptual tools to pose problems within large enough frameworks to find 624 meaningful solutions [89]. 625

Conclusion 626
Our cross-national and relational geographic approach to studying the COVID-19 pandemic in 627 Mesoamerica and the Caribbean showed that global political-economic roles of countries are also 628 extremely important to understand how diseases spread within societies. This approach went beyond 629 standard epidemiological analyses focused on public health measures to prevent and mitigate disease 630 transmission, which estimate transmission parameters and alienate human movement from structural 631 factors shaping space as defined by used territories. Some proxies for spatial connections examined in 632 this study were directly and positively associated with SARS-CoV-2 transmission (measured in terms of 633 epidemic growth, cases and deaths), such as the size of the network of cities connected across major 634 airports (parameter CC in Table 1). However, our results also support the epidemiological relevance of 635 social theories that stress a broader spectrum of power-laden global connections. We found that 636 variables key to the formation of space in used territories, e.g., concentration of imports, trade 637 openness and international aid, were positively associated with increases in disease transmission and 638 deaths. These variables, conceived within social science theory, might remain important over much 639 longer timescales than the unfolding of the COVID-19 pandemic and deserve further study in developing 640 mathematical models dealing with spatial and spatio-temporal disease spread. 641    Table 1. The scale at the bottom is for the Pearson correlation coefficient, square size is likewise proportional to the correlation value.

Figure 4
Cluster analysis results. Dendrogram representation of the agglomerative clusters (A) Epidemiological parameters (B) Used territory and human development and uneven development parameters. Parallel coordinates plots for cluster segmentation (C) Epidemiological parameters (D) Used territory and human development and uneven development parameters. In panels A & B branches are colored to indicate clusters. In panel B Haiti, Cuba, Puerto Rico and Jamaica did not belong to any cluster but have the same color to indicate this property. In panels C & D clusters are indicated, respectively, using the same colors used to represent the clusters in panels A & B, and different type lines are used to represent Cuba, Haiti, Jamaica and Puerto Rico in panel D, which are indicated in the inset legend. The xaxis of panels C & D indicate the different parameters we used in the cluster analysis.

Figure 5
Map of Mesoamerica and the Caribbean displaying relations between disease and used territory, human development and uneven development parameters. Weavings are used to depict the two patterns of clustering. In the map we also included Uruguay, since data for that country was used as an outgroup for the underlying analyses.

Supplementary Online Figure Legends
Figure S1 COVID-19 timeline in Mesoamerica the Caribbean and Uruguay, covering March 2020. In the plot the date for the first case and the implementation of the different prevention, mitigation and suppression policies is indicated by symbols along time (x axis). Figure 2 in the main manuscript shows the timeline for the first 100 days of transmission.

Figure S2
Principal component Analysis representation of the agglomerative cluster analysis (A) Epidemiological Parameters (B) Used territory and human and uneven development parameters. In both panels territories are presented by their ISO 3166-1 alpha-2 two letter code presented in Figure 1 from the main text. In each plot clusters are represented with the same colors as in Figure 4 from the main text. Axes in each plot indicate the 1 st (Dim1) and 2 nd (Dim2) principal component, and the percent of variability explained by each principal component is presented within parenthesis. In both panels colors used to indicate clusters are the same as those presented in Figure 4.

Figure S3
Diagnostics for the model presented in Table 2. Each plot has self-explanatory titles and labels.

Figure S4
Diagnostics for the model presented in Table 3. Each plot has self-explanatory titles and labels.

Figure S5
Diagnostics for the model presented in Table 4. Each plot has self-explanatory titles and labels.   Table 3 Parameter estimates for the best negative binomial generalized linear model explaining total COVID-19 Cases 100 days after the detection of the first case across nations from Mesoamerica and Caribbean. In the model estimates for population size in 2020 for each country were used as offset to account for differences in the size of the population at risk. J o u r n a l P r e -p r o o f Table 4 Parameter estimates for the best negative binomial generalized linear model explaining total COVID-19 Deaths 100 days after the detection of the first case across nations from Mesoamerica and Caribbean. In the model estimates for population size in 2020 for each country were used as offset to account for differences in the size of the population at risk.