Comparing lagged impacts of mobility changes and environmental factors on COVID-19 waves in rural and urban India: a Bayesian spatiotemporal modelling study

Previous research in India has identified urbanisation, human mobility and population demographics as key variables associated with higher district level COVID-19 incidence. However, the spatiotemporal dynamics of mobility patterns in rural and urban areas in India, in conjunction with other drivers of COVID-19 transmission, have not been fully investigated. We explored travel networks within India during two pandemic waves using aggregated and anonymized weekly human movement datasets obtained from Google, and quantified changes in mobility before and during the pandemic compared with the mean baseline mobility for the 8-week time period at the beginning of 2020. We fit Bayesian spatiotemporal hierarchical models coupled with distributed lag non-linear models (DLNM) within the integrated nested Laplace approximate (INLA) package in R to examine the lag-response associations of drivers of COVID-19 transmission in urban, suburban, and rural districts in India during two pandemic waves in 2020-2021. Model results demonstrate that recovery of mobility to 99% that of pre-pandemic levels was associated with an increase in relative risk of COVID-19 transmission during the Delta wave of transmission. This increased mobility, coupled with reduced stringency in public intervention policy and the emergence of the Delta variant, were the main contributors to the high COVID-19 transmission peak in India in April 2021. During both pandemic waves in India, reduction in human mobility, higher stringency of interventions, and climate factors (temperature and precipitation) had 2-week lag-response impacts on the Rt of COVID-19 transmission, with variations in drivers of COVID-19 transmission observed across urban, rural and suburban areas. With the increased likelihood of emergent novel infections and disease outbreaks under a changing global climate, providing a framework for understanding the lagged impact of spatiotemporal drivers of infection transmission will be crucial for informing interventions.

second lockdown period (8 weeks) from January 31 to March 27, 2021; 4) Second lockdown (6 weeks) for the Delta wave, from April 18 to May 29, 2021; 5) post-second lockdown period (8 weeks), from November 7 to December 31, 2021, after travel restrictions for COVID-19 had been lifted in India.Table S1.Summary Statistics for data used for wave 1 and Delta wave spatiotemporal models SI Delta Wave Table S2.Wave 2: Adequacy results for models with DLNMs and increasing complexity.Table S3.Wave 2: Adequacy results for models (without DLNMs) using 2-week lag covariates with increasing complexity.Table S4.Model hyperparameters using a range of prior distributions in best fit model 4.1 for Delta Wave Figure S4.Relative intra-district mobility during the Delta wave in India, standardised by pre-pandemic mean baseline levels of mobility for the first eight weeks of 2020 (December 29, 2019-February 22, 2020) for each district.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.Areas shaded in grey are areas for which no data is available.S2), using 2-week lag covariates and leave-one-week-out cross-validation approach.States are ordered by their geographical location.covariates and leave-one-state-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

SI Wave 1
Table S5.Wave 1: Adequacy results for models with DLNMs and increasing complexity.Table S6.Wave 1: Adequacy results for models (without DLNMs) using 2-week lag covariates with increasing complexity.Areas shaded in grey are areas for which no data is available.Areas shaded in grey are areas for which no data is available.12 leave-one-week-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

Figure S48.
Posterior predictive mean Rt during the wave 1 in India, 2020, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag covariates and leave-one-districtout cross-validation approach.Areas shaded in grey are areas for which no data is available. .

Fig S3 .
Fig S3.Relative changes of outbound travel from districts across India during the pandemic compared with average pre-pandemic levels during the 12 weeks from November 10, 2019, to February 22, 2020.(A) Reductions of outbound flows under the first lockdown during the 6-week period from March 22 to May 2, 2020.(B) Changes in outflow during the 8-week period from January 31 to March 27, 2021, before the second lockdown.(C) Reductions of outflows during the 6-week second lockdown from April 18 to May 29, 2021.(D) Changes in outflow during the 8-week period from November 7 to December 31, 2021.Sub-division maps at administrative level I (state) and II (district) were obtained from the GADM version 3.6 (https://gadm.org/).Regions in which outflow data are not available are those represented in green.Areas shaded in grey are areas for which no data is available.

Figure S5 .
Figure S5.Stringency Index of COVID-19 intervention policy implemented during the Delta wave in India.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S6 .
Figure S6.Mean temperature at 2m above the surface during the Delta wave in India.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S7 .
Figure S7.Accumulated weekly precipitation (metres) during the Delta wave in India.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S8 .
Figure S8.Relative humidity during the Delta wave in India.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S9 .
Figure S9.Downward ultraviolet (UV) radiation (KJ/m2 per hour) during the Delta wave in India.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S10 .
Figure S10.Weekly Rt derived from COVID-19 cases reported during the Delta wave in India.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S11 .
Figure S11.Pairwise Pearson correlations between weekly means of variables at district level during the Delta wave in India, 2021.R0: basic reproduction number.Rt: instantaneous reproduction number.ln_R:log(Rt/R0).Cases_rate: new COVID-19 cases reported per 1000 people.Cases_accu_rate: cumulative cases per 1000 people reported since the first week of the wave.mean_intra: intra-district relative mobility.d2m: relative humidity.t2m: mean temperature of air (°C at 2m above the surface of land, sea or inland waters).tp: precipitation (metres).uv: downward ultraviolet radiation.Stringency: index of COVID-19 intervention stringency.Holiday: days of public holidays in a week.pop_sum: total population of each district.pop_density: population number per km2 of each district.

Figure S12 .
Figure S12.Kendall rank correlations between weekly means of variables at district level during the Delta wave in India, 2021.R0: basic reproduction number.Rt: instantaneous reproduction number.ln_R:log(Rt/R0).Cases_rate: new COVID-19 cases reported per 1000 people.Cases_accu_rate: cumulative cases per 1000 people reported since the first week of the wave.mean_intra: intra-district relative mobility.d2m: relative humidity.t2m: mean temperature of air (°C at 2m above the surface of land, sea or inland waters).tp: precipitation (metres).uv: downward ultraviolet radiation.Stringency: index of COVID-19 intervention stringency.Holiday: days of public holidays in a week.pop_sum: total population of each district.pop_density: population number per km2 of each district.

Figure S13 .
Figure S13.Posterior predictive mean Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1) at country level using leave-one-week-out cross-validation approach.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S14 .
Figure S14.Standard deviation (SD) of posterior predictive Rt during the Delta wave in India, 2021,derived from the best fitting model (model 4.1 without DLNMs) at country level using a leave-oneweek-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

Figure S15 .
Figure S15.Posterior predictive mean Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1) at country level using leave-one-state-out cross-validation approach.The weeks in 2021 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S16 .
Figure S16.Standard deviation (SD) of posterior predictive Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1 without DLNMs) at country level using a leave-onestate-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

Figure S17 .
Figure S17.Contribution of spatial random effects to estimates of Rt changes in the base model.Areas shaded in grey are areas for which no data is available.

Figure S18 .
Figure S18.Improvement by using the best fitting model across the country, compared to baseline model.Difference between mean absolute error (MAE) for the baseline model (weekly random effects, spatial random effects and population density) and MAE for the best fitting model (model 4.1 with DLNMs).Districts with positive values (pink) suggest that capturing the nonlinear and delayed impacts of mobility, climate information and intervention stringency, improves the model in these areas.Districts with negative values (blue) suggest that mobility, intervention and climate information did not improve the model fit and other unexplained factors might dominate space-time dynamics in these areas.The MAE of the selected model was smaller than the baseline model for 385 of the 665 (57.9%) districts in India, with the results of model performance provided by geo-political regions in the Table.

Figure S19 .
Figure S19.Observed versus posterior fitted Rt in the capital district of each state using the best fitting model (model 4.1 with DLNMs) at country level.Graphs with a log scale at y-axis show the observed Rt derived from reported case data, and corresponding mean and 95% confidence interval (CI, shaded pink area) of fitted Rt, derived from the best fitting model (model 4.1 with DLNMs) at country level.States are ordered by their geographical location.

Figure S20 .
Figure S20.Observed versus posterior predictive Rt in the capital district of each state, using leaveone-week-out cross-validation approach.Graphs with a log scale at y-axis show the observed Rt derived from reported case data, and corresponding posterior predictive mean and 95% prediction interval (CI, shaded pink area), derived from the best fitting model (model 4.1 with DLNMs) at country level.States are ordered by their geographical location.

Figure S21 .
Figure S21.Contribution of spatial random effects to estimates of Rt changes in the base model.Areas shaded in grey are areas for which no data is available.

Figure S22 .
Figure S22.Improvement of using the best fitting model with 2-week lag covariates (no DLNMs), compared to baseline model with the same lag.Difference between mean absolute error (MAE) for the baseline model and MAE for the best fitting model (Model 4.1).Districts with positive values (pink) suggest that capturing the 2-week lag impacts of mobility, temperature, UV and intervention stringency, improves the model in these areas.Districts with negative values (blue) suggest that mobility, intervention and climate information did not improve the model fit and other unexplained factors might dominate space-time dynamics in these areas.The MAE of the selected model was smaller than the baseline model for 428 of the 665 (64.4%) districts in India, and further improved the best fitting model with DLNMs (Figure S12).Results of model performance are provided by geo-political regions in theTable.Areas shaded in grey are areas for which no data is available.

Figure S23 .
Figure S23.Posterior predictive mean Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag covariates and leave-oneweek-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

Figure S24 .
Figure S24.Standard deviation (SD) of posterior predictive Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag covariates and leave-one-week-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

Fig S25 .
Fig S25.Observed versus posterior predictive Rt in the capital district of each state.Graphs with a log scale at y-axis show the observed Rt derived from reported case data, and corresponding posterior predictive mean and 95% prediction interval (CI, shaded pink area), derived from the best fitting model without DLNMs at country level (model 4.1: base model + mobility + temperature + UV + intervention policy; see SI TableS2), using 2-week lag covariates and leave-one-week-out cross-validation

Figure S26 .
Figure S26.Posterior predictive mean Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag covariates and leave-onestate-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

Figure S27 .
Figure S27.Standard deviation (SD) of posterior predictive Rt during the Delta wave in India, 2021, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag

Figure S29 .
Figure S29.Relative intra-district mobility during wave 1 in India, standardised by pre-pandemic mean baseline levels of mobility for the first eight weeks of 2020(December 29, 2019 -February 22, 2020)    for each district.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S30 .
Figure S30.Stringency Index of COVID-19 intervention policy implemented during wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S31 .
Figure S31.Mean temperature at 2m above the surface during wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S32 .
Figure S32.Accumulated weekly precipitation (metres) during wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S33 .
Figure S33.Relative humidity during wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S34 .
Figure S34.Downward ultraviolet (UV) radiation (KJ/m2 per hour) during wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S35 .
Figure S35.Weekly Rt derived from COVID-19 cases reported during the wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S36 .
Figure S36.Pairwise Pearson correlations between weekly means of variables at district level during the wave 1 in India, 2020.R0: basic reproduction number.Rt: instantaneous reproduction number.ln_R: log(Rt/R0).Cases_rate: new COVID-19 cases reported per 1000 people.Cases_accu_rate: cumulative cases per 1000 people reported since the first week of the wave.mean_intra: intra-district relative mobility.d2m: relative humidity.t2m: mean temperature of air (°C at 2m above the surface of land, sea or inland waters).tp: precipitation (metres).uv: downward ultraviolet radiation.Stringency: index of COVID-19 intervention stringency.Holiday: days of public holidays in a week.pop_sum: total population of each district.pop_density: population number per km2 of each district.

Figure S37 .
Figure S37.Kendall rank correlations between weekly means of variables at district level during the wave 1 in India, 2020.R0: basic reproduction number.Rt: instantaneous reproduction number.ln_R: log(Rt/R0).Cases_rate: new COVID-19 cases reported per 1000 people.Cases_accu_rate: cumulative cases per 1000 people reported since the first week of the wave.mean_intra: intra-district relative mobility.d2m: relative humidity.t2m: mean temperature of air (°C at 2m above the surface of land, sea or inland waters).tp: precipitation (metres).uv: downward ultraviolet radiation.Stringency: index of COVID-19 intervention stringency.Holiday: days of public holidays in a week.pop_sum: total population of each district.pop_density: population number per km2 of each district.

Figure S38 .
Figure S38.Posterior predictive mean Rt during wave 1 in India, 2020, derived from the best fitting model (model 4.1) at country level using leave-one-week-out cross-validation approach.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S40 ..Fig S41 .
Figure S40.Posterior predictive mean Rt during wave 1 in India, 2020, derived from the best fitting model (model 4.1) at country level using leave-one-district-out cross-validation approach.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.

Figure S42 .
Figure S42.Contribution of spatial random effects to estimates of Rt changes in the base model.Areas shaded in grey are areas for which no data is available.

Figure S43 .
Figure S43.Improvement by using the best fitting model across the country, compared to baseline model.Difference between mean absolute error (MAE) for the baseline model (weekly random effects, spatial random effects and population density) and MAE for the best fitting model (model 4.1 with DLNMs).Districts with positive values (pink) suggest that capturing the nonlinear and delayed impacts of mobility, climate information and intervention stringency, improves the model in these areas.Districts with negative values (blue) suggest that mobility, intervention and climate information did not improve the model fit and other unexplained factors might dominate space-time dynamics in these areas.The MAE of the selected model was smaller than the baseline model for 430 of the 661 (65.17%) districts in India, with the results of model performance provided by geo-political regions in the Table.

Figure S44 .
Figure S44.Observed versus posterior fitted Rt in the capital district of each state using the best fitting model (model 4.1 with DLNMs) at country level.Graphs with a log scale at y-axis show the observed Rt derived from reported case data, and corresponding mean and 95% confidence interval (CI, shaded pink area) of fitted Rt, derived from the best fitting model (model 4.1 with DLNMs) at country level.States are ordered by their geographical location.

Figure S45 .
Figure S45.Observed versus posterior predictive Rt in the capital district of each state, using leaveone-week-out cross-validation approach.Graphs with a log scale at y-axis show the observed Rt derived from reported case data, and corresponding posterior predictive mean and 95% prediction interval (CI, shaded pink area), derived from the best fitting model (model 4.1 with DLNMs) at country level.States are ordered by their geographical location.

Figure S46 .
Figure S46.Posterior predictive mean Rt during the wave 1 in India, 2020, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag covariates and leave-one-weekout cross-validation approach.Areas shaded in grey are areas for which no data is available.

Fig S49 .
Fig S49.Standard deviation (SD) of posterior predictive Rt during wave 1 in India, 2020, derived from the best fitting model (model 4.1 without DLNMs) at country level using 2-week lag covariates and leave-one-district-out cross-validation approach.Areas shaded in grey are areas for which no data is available.

which was not certified by peer review)
International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(

Table S7 . Model hyperparameters using a range of prior distributions in best fit model 4.1 for Wave 1 Figure S28.
COVID-19 cases reported by district each week during wave 1 in India.The weeks in 2020 investigated are numbered in maps.Areas shaded in grey are areas for which no data is available.