Association between the Wuhan city travel ban and the arrival time of COVID-19 in other cities

Here we quantify the association between the Wuhan travel shut- down (23 January 2020) and COVID-19 spread.

rm(list = ls())
library(readxl)
library(car)                                                

## Loading required package: carData

library(lmtest)                                                        

## Loading requiredpackage: zoo ##

## Attaching package: 1zoo1

## The following objects are maskedfrom 1package:base1: ##

##    as.Date,as.Date.numeric

library(sandwich)
library(fitdistrplus)                                    

## Loadingrequired package: MASS

## Loading requiredpackage: survival 
## Loading requiredpackage: npsurv 
## Loadingrequired package: lsei

 library(R330)                                             

## Loading requiredpackage: s20x

##

## Attaching package: 1s20x1

## The following object is masked from 1package:car1:

##

##    levene.test

## Loading requiredpackage: leaps ##Loading required package: rgl

## Loading required package: lattice

 infected=read_xlsx("./Travel_ban_model.xlsx",sheet="1")                                                        

Processing the data

infectedW=infected[infected$arrivalday>=2,]
infectedW$Totalflow=infectedW$Train2018+infectedW$Road2018+infectedW$Air2018+1 

Regression analysis

The association between distance, human movement, interventions and timing of COVID-19 spread was assessed by a linear regression.

fit <- lm(arrivalday~lat+lon+log10(Pop2018)+log10(Distance)+ log10(Train2018+1)+ log10(Air2018+1)+ log10(Road2018+1)+log10(Totalflow)+aftershutdown,
data=infectedW)
lm.step<-step(fit)
 

## Start: AIC=113.73

## arrivalday ~ lat + lon + log10(Pop2018) + log10(Distance) + log10(Train2018 + 
##    1) + log10(Air2018 + 1) + log10(Road2018 + 1) + log10(Totalflow) +

## aftershutdown 
##
 

## Step:  AIC=111.91

## arrivalday ~ lat + lon + log10(Pop2018) + log10(Distance) + log10(Air2018 + 
## 1) + log10(Road2018 + 1) + log10(Totalflow) + aftershutdown

##
 

##

## Step:  AIC=110.63

## arrivalday ~ lat + lon + log10(Pop2018) + log10(Distance) + log10(Air2018 + 
## 1) + log10(Totalflow) + aftershutdown

##

##

## Step:  AIC=109.83

## arrivalday ~ lat + lon + log10(Pop2018) + log10(Distance) + log10(Totalflow) + 
##   aftershutdown

##

##

## Step:  AIC=108.57

## arrivalday ~ lat + lon + log10(Pop2018) + log10(Totalflow) + 
##      aftershutdown

##

Selection model

Among five possible regression models examined, the model judged best by the Akaike Information Criterion.

step.model <- stepAIC(fit, direction = "backward", trace = FALSE)

summary(step.model)

##

## Call:

## lm(formula = arrivalday ~ lat + lon + log10(Pop2018) + log10(Totalflow) + 
##      aftershutdown, data = infectedW)

##

## Residuals:

##    Min   1Q  Median    3Q   Max 
## -4.5232 -0.8999 0.0366 0.8347 2.7326 
##

## Coefficients:

## Signif. codes: 0 1***10.001 1**1 0.01 1*1 0.051.1 0.1 1 1 1 
##

## Residual standarderror: 1.217 on 255 degreesof freedom 
## Multiple R-squared:           0.589,         Adjusted R-squared:            0.581 
## F-statistic:             73.1 on 5 and 255 DF,               p-value: <2.2e-16

There appears to be five influential points.
 

 influencePlot(step.model)                                                        

##         StudRes      Hat      CookD 
## 64 -3.4440481 0.02149690 0.041656769

## 72  -3.8878442 0.03598356 0.089102090

## 255-2.3402866 0.08700970 0.085492844

## 257-0.7385126 0.12764904 0.013324957

## 258  0.7114217 0.10448601 0.009861219

Removing the five influential points does not affect the conclusions.

step.model.1 <- lm(arrivalday~lat+lon+log10(Pop2018)+log10(Totalflow)+aftershutdown, data=infectedW[-c(64,72,255,257,258),])

summary(step.model.1)
 

##

## Call:

## lm(formula = arrivalday ~ lat + lon + log10(Pop2018) + log10(Totalflow) + 
##      aftershutdown, data = infectedW[-c(64, 72, 255, 257,258),

##    ])

##

## Residuals:

##    Min   1Q  Median    3Q   Max 
## -4.0551 -0.8616 0.0285 0.7634 2.7182 
##

## Coefficients:

## Signif. codes: 0 1***1 0.0011**1 0.01 1*1 0.05 1.1 0.1 1 1 1 
##

## Residual standarderror: 1.15 on 250 degreesof freedom 
## Multiple R-squared:             0.5949, Adjusted R-squared:  0.5868 
##F-statistic: 73.42 on 5 and 250 DF,                  p-value: <2.2e-16

There are still influential points.

 influencePlot(step.model.1)                                                        

The Cook’s distances do not indicate presence of outliers. Here are a number of observations indicated as having large hat values. However, they are not too far away from they cutoff, and their hat values are quite

similar to each other. So it is hard to justigy removing only a subset of those points, but there are quite a few of them, so we leave them in the model.

par(mfrow = c(2,4)) 
influenceplots(step.model.1)
 

We did not detect issuesregarding heteroscedasticity (studentized Breuch-Pagan test P=0.20)nor normality (Shapiro-Wilk test P=0.06).
 

 lmtest::bptest(step.model.1)                                                        

##

## studentized Breusch-Pagan test 
##

## data:  step.model.1

## BP = 7.2846, df = 5, p-value = 0.2003

 shapiro.test(step.model.1$residuals)                                                        

##

## Shapiro-Wilk normality test 
##

## data: step.model.1$residuals 
## W = 0.98961, p-value= 0.06359

The plot of studentised residuals does not indicate evident non-linearity in the residuals.

par(mfrow = c(2,2)) 
plot(step.model.1)

Diagnostics

Check for spatial correlation in the residuals

Here we performed additional checks to investigate whether there was spatial autocorrelation among the residuals, which could contribute to data dependency. If there were spatial autocorrelation in the residuals, then the pairwise residual differences of cities closer together would tend to be smaller than those of cities far apart. In other words, we would expect a high correlation between the pairwise differences in residuals and pairwise distances of cities. Such a correlation was not evident (r=0.03).

fit.res.dist.mat = as.matrix(dist(residuals(step.model, type = "pearson"))) 
gcs.dist.mat=as.matrix(dist(data.frame(infectedW$lon,infectedW$lat))) 
cor.test(gcs.dist.mat,fit.res.dist.mat)
 

##

##    Pearson1s product-moment correlation 
##

## data:         gcs.dist.mat and fit.res.dist.mat 
## t = 8.5113, df = 68119,p-value < 2.2e-16

## alternative hypothesis: true correlation is not equalto 0 
## 95 percent confidence interval:

##  0.025090330.04009330

## sample estimates:

##       cor

## 0.03259365

Differential accessibility from Wuhan to other cities may contribute to data-dependence. Therefore, we examined the timing of peak inflow from Wuhan City and the arrival time of COVID-19 in each city. The result shows no evident correlation between them (r=-0.07, P=0.25).
 

 cor.test(infectedW$peak_time2,infectedW$arrivalday)                                                        

##

##    Pearson1s product-moment correlation 
##

## data: infectedW$peak_time2 and infectedW$arrivalday 
## t = -1.1663,df = 246, p-value = 0.2446

## alternative hypothesis: true correlation is not equalto 0 
## 95 percent confidence interval:

##  -0.19690241 0.05088306

## sample estimates:

##        cor

## -0.07415409

In addition, we calculated the correlation between the residuals and peak time of Wuhan inflow during 11 to 23 January (15 days before Chinese New Year to the Wuhan shutdown). Again, we found no evident correlation (r=-0.08, P=0.20).
 

 cor.test(step.model$residuals,infectedW$peak_time1)                                                        

##

##    Pearson1s product-moment correlation 
##

## data: step.model$residuals and infectedW$peak_time1 
## t = -1.2948,df = 246, p-value = 0.1966

## alternative hypothesis: true correlation is not equalto 0 
## 95 percent confidence interval:

##  -0.2047447 0.0427291

## sample estimates: 
## cor

##-0.08227597