## Data Analysis - Coelho et al. (2020) - PeerJ
## paper title : "Global expansion of covid19 pandemic is driven by population size and airport connections"
## Data used in this script is available as Supporting information in Xls.
## In the data table, the is one table with data and another with metadata
## save the data in .txt, or .csv and import it to R

#packages
library(fields)
library(ape)
library(car)
library(letsR)
library(QuantPsyc)

#import data
Data <- read.table('Data.txt', head = T, stringsAsFactors = FALSE)


#great-circle distances
geokm <- rdist.earth(as.matrix(Data[, c("Long", "Lat")]),
                     miles = FALSE,
                     R = NULL) 

##model 0  - Assuming only the effect of the environment

mod <- lm(Data[,'b'] ~ 
            Data[, 'MeanTemp'] +
            log(Data[,'MeanPrec']))

summary(mod)
vif(mod)
res <- mod$residuals
res.std <- as.vector((res - mean(res)) / sd(res))
# inspect the spatial autocorrelation in model's residuals
correl <- lets.correl(res.std, geokm, 10, T) # Correlogram


## model 1 - Here assuming only environmental and socioeconomic variables

mod <- lm(Data[,'b'] ~ 
            Data[,'GNI_2018'] + 
            log(Data[,'Population']) + 
            Data[, 'Annual_pop_growth'] + 
            sqrt(Data[, 'Health_investments']) +
            Data[, 'MeanTemp'] +
            log(Data[,'MeanPrec']))

summary(mod)
vif(mod)
res <- mod$residuals
res.std <- as.vector((res - mean(res)) / sd(res))
# inspect the spatial autocorrelation in model's residuals
correl <- lets.correl(res.std, geokm, 10, T) # Correlogram

#standirdized coefficients
lm.beta(mod)



## Model 2 - Adding network centrality
mod <- lm(Data[,'b'] ~ 
            log(Data[, 'Eigenvector_Centrality']) +
            Data[,'GNI_2018'] + 
            log(Data[,'Population']) + 
            Data[, 'Annual_pop_growth'] + 
            sqrt(Data[, 'Health_investments']) +
            Data[, 'MeanTemp'] +
            log(Data[,'MeanPrec']))

summary(mod)
vif(mod)
res <- mod$residuals
res.std <- as.vector((res - mean(res)) / sd(res))
#inspect the spatial autocorrelation in model's residuals
correl <- lets.correl(res.std, geokm, 10, T) # Correlogram
#std coefficients
lm.beta(mod)
