files for zenodo

laispfreitas · Jul 8, 2020 · 1188a2a · 1188a2a
1 parent fb3c4f0
commit 1188a2a
Show file tree

Hide file tree

Showing 9 changed files with 14,908 additions and 0 deletions.
diff --git a/ICAR_models_chikungunya_rio.r b/ICAR_models_chikungunya_rio.r
@@ -0,0 +1,220 @@
+##################################################################################
+###                Spatio-temporal distribution of Chikungunya in the          ###
+###                         city of Rio de Janeiro, Brazil                     ###
+##################################################################################
+
+## Laís Picinini Freitas and Alexandra M. Schmidt
+
+## Last update: 8 July 2020
+
+## Preprint: https://www.medrxiv.org/content/10.1101/2020.06.08.20125757v1 
+
+## Objective: 
+# To study the spatio-temporal dynamics of the first chikungunya epidemic in Rio de 
+# Janeiro city, exploring the effects of temperature, green area proportion and 
+# sociodevelopment index.
+
+rm(list=ls())
+
+library(rstan)
+options(mc.cores=parallel::detectCores())
+rstan_options(auto_write=TRUE)
+
+### Data ------------------------------------------------------------------------------------------
+
+# Cases data can be downloaded at the source at the Rio de Janeiro city hall website: 
+# http://www.rio.rj.gov.br/web/sms/exibeConteudo?id=4769664
+
+# Population, sociodevelopment index and land use (green area) data were obtained from the Instituto Pereira Passos at www.data.rio. 
+# Temperature data were obtained from the Brazilian National Institute of Meteorology (http://www.inmet.gov.br/portal/index.php?r=estacoes/estacoesAutomaticas), 
+# the Brazilian Airspace Control Department (https://www.redemet.aer.mil.br/?i=produtos&p=consulta-de-mensagens-opmet), 
+# the Rio de Janeiro State Environmental Institute (http://200.20.53.25/qualiar/home/index), 
+# the Rio de Janeiro Municipal Environmental Secretariat (http://www.data.rio/datasets/dados-hor%C3%A1rios-do-monitoramento-da-qualidade-do-ar-monitorar?orderBy=Data) 
+# and the Alerta Rio System (http://alertario.rio.rj.gov.br/download/dados-meteorologicos/).
+
+
+# For the purpose of this script, we will use the number of cases by neighbourhood and week sampled from one chain of the posterior predicted values from the results of model4:
+cases <- read.csv('simulated_cases.csv')
+y <- matrix(cases$cases, nrow=160, byrow=TRUE)
+
+# covariables: population, sociodevelopment index (sdi), green area proportion
+covs_neigh <- read.csv('covariates_rio.csv') 
+
+# covariables matrix
+x <- matrix(ncol=2, nrow=length(covs_neigh$CodBairro))
+x[,1] <- covs_neigh$sdi
+x[,2] <- (covs_neigh$green_area^(1/3))
+
+# temperature data
+temperature <- read.csv('temperature_rio_2016.csv')
+tmin <- matrix(temperature$min_temperature, nrow=160, byrow=TRUE)
+tmincenter <- (tmin-mean(tmin))/sd(tmin)
+
+N <- nrow(y)
+TAM <- ncol(y)
+
+# calculating the offset (expected cases)
+E <- ((sum(y)/sum(covs_bairro$population))*covs_bairro$population) / TAM
+
+
+### Neighbourhood matrix --------------------------------------------------------------------------
+
+# functions to build the neighbourhood matrix in stan
+library(devtools)
+source_url("https://github.com/stan-dev/example-models/blob/master/knitr/car-iar-poisson/mungeCARdata4stan.R?raw=TRUE")
+
+# loading the neighbourhood matrix (number of neighbours and which are the neighbours)
+numneigh <- scan("NumNeigh.csv")
+neigh <- scan("Neigh.csv")
+
+nbs <- mungeCARdata4stan(neigh, numneigh)
+N <- nbs$N
+node1 <- nbs$node1
+node2 <- nbs$node2
+N_edges <- nbs$N_edges
+
+
+### Models ----------------------------------------------------------------------------------------
+
+
+### Model 0 ####
+
+# Model in stan
+model0 <- stan("model0.stan", 
+               data=list(N=N, T=TAM, m0=0, C0=5,
+                         N_edges=N_edges, node1=node1, node2=node2,
+                         y=y, E=E), 
+               chains=4, iter=10000,
+               control=list(adapt_delta=0.95, max_treedepth=15))
+
+# Summary of the model, WAIC
+print(model0, probs=c(0.05, 0.95), digits_summary=3)
+
+require(loo)
+loo(model0)
+loglikGa <- extract_log_lik(model0)
+waic0 <- waic(loglikGa)
+waic0
+
+check_hmc_diagnostics(model0)
+
+# Traceplots of the chains
+
+library(bayesplot)
+
+traceplot(model0, pars=c("beta0","sigma"), inc_warmup=FALSE)
+
+
+### Model 1 ####
+
+# Model in stan
+model1 <- stan("model1.stan", 
+               data=list(N=N, T=TAM, m0=0, C0=5,
+                         y=y, x=x, E=E, K=2), 
+               chains=4, iter=10000,
+               control=list(adapt_delta=0.95, max_treedepth=15))
+
+# Summary of the model, WAIC
+print(model1, probs=c(0.05, 0.95), digits_summary=3)
+
+require(loo)
+loo(model1)
+loglikGa <- extract_log_lik(model1)
+waic0 <- waic(loglikGa)
+waic0
+
+check_hmc_diagnostics(model1)
+
+# Traceplots of the chains
+
+library(bayesplot)
+
+traceplot(model1, pars=c("beta0","W"), inc_warmup=FALSE)
+
+
+### Model 2 ####
+
+# Model in stan
+model2 <- stan("model2.stan", 
+               data=list(N=N, T=TAM, m0=0, C0=5,
+                         N_edges=N_edges, node1=node1, node2=node2,
+                         y=y, x=x, E=E, K=2), 
+               chains=4, iter=10000,
+               control=list(adapt_delta=0.95, max_treedepth=15))
+
+# Summary of the model, WAIC
+print(model2, probs=c(0.05, 0.95), digits_summary=3)
+
+require(loo)
+loo(model2)
+loglikGa <- extract_log_lik(model2)
+waic0 <- waic(loglikGa)
+waic0
+
+check_hmc_diagnostics(modelo)
+
+# Traceplots of the chains
+
+library(bayesplot)
+
+traceplot(model2, pars=c("beta0","sigma"), inc_warmup=FALSE)
+traceplot(model2, pars=c("W"), inc_warmup=FALSE)
+
+
+### Model 3 ####
+
+# Model in stan
+model3 <- stan("model3.stan", 
+               data=list(N=N, T=TAM, m0=0, C0=5,
+                         N_edges=N_edges, node1=node1, node2=node2,
+                         y=y, tmincenter=tmincenter, E=E), 
+               chains=4, iter=10000,
+               control=list(adapt_delta=0.95, max_treedepth=15))
+
+# Summary of the model, WAIC
+print(model3, probs=c(0.05, 0.95), digits_summary=3)
+
+require(loo)
+loo(model3)
+loglikGa <- extract_log_lik(model3)
+waic0 <- waic(loglikGa)
+waic0
+
+check_hmc_diagnostics(modelo)
+
+# Traceplots of the chains
+
+library(bayesplot)
+
+traceplot(model3, pars=c("beta0","sigma"), inc_warmup=FALSE)
+
+
+### Model 4 ---------------------------------------------------------------------------------------
+
+# Model in stan
+model4 <- stan("model4.stan", 
+                         data=list(N=N, T=TAM, m0=0, C0=5,
+                                   N_edges=N_edges, node1=node1, node2=node2,
+                                   y=y, x=x, tmincenter=tmincenter, E=E, K=2), 
+                         chains=4, iter=10000,
+                         control=list(adapt_delta=0.95, max_treedepth=15))
+
+# Summary of the model, WAIC
+print(model4, probs=c(0.05, 0.95), digits_summary=3)
+
+require(loo)
+loo(model4)
+loglikGa <- extract_log_lik(model4)
+waic0 <- waic(loglikGa)
+waic0
+
+check_hmc_diagnostics(modelo)
+
+# Traceplots of the chains
+
+library(bayesplot)
+
+traceplot(model4,pars=c("beta0","sigma"), inc_warmup=FALSE)
+traceplot(model4, pars=c("W"), inc_warmup=FALSE)
+
+
diff --git a/covariates_rio.csv b/covariates_rio.csv
@@ -0,0 +1,161 @@
+neigh_name,neigh_code,population,sdi,green_area
+SAUDE,1,2749,0.583,0.014690970998359
+GAMBOA,2,13108,0.559,0.021500044743323
+SANTO CRISTO,3,12330,0.569,0.008760813829182
+CAJU,4,20477,0.554,0.017527831214039
+CENTRO,5,41142,0.643,0.002739587981173
+CATUMBI,6,12556,0.58,0.09420806864977
+RIO COMPRIDO,7,43764,0.605,0.356610705524256
+CIDADE NOVA,8,5466,0.594,0
+ESTACIO,9,17189,0.586,0.021986291466026
+SAO CRISTOVAO,10,26510,0.615,0.021403343998946
+MANGUEIRA,11,17835,0.537,0.186612517124832
+BENFICA,12,25081,0.57,0
+PAQUETA,13,3361,0.608,0.30901532613133
+SANTA TERESA,14,40926,0.624,0.524727419144311
+FLAMENGO,15,50043,0.752,0.034391158845145
+GLORIA,16,9661,0.705,0.032255569397758
+LARANJEIRAS,17,45554,0.75,0.336106592054309
+CATETE,18,24057,0.69,0.113333515892728
+COSME VELHO,19,7178,0.693,0.362178516445211
+BOTAFOGO,20,82890,0.733,0.22394421141313
+HUMAITA,21,13285,0.761,0.558477508936037
+URCA,22,7061,0.749,0.446433650008137
+LEME,23,14799,0.723,0.375359834101737
+COPACABANA,24,146392,0.731,0.113376793578974
+IPANEMA,25,42743,0.77,0.001189924149599
+LEBLON,26,46044,0.78,0.164732604151538
+LAGOA,27,21198,0.819,0.079090057204596
+JARDIM BOTANICO,28,18009,0.767,0.374949706218315
+GAVEA,29,16003,0.756,0.527715315673601
+VIDIGAL,30,12797,0.565,0.306282173643349
+SAO CONRADO,31,10980,0.779,0.55808589535582
+PRACA DA BANDEIRA,32,8662,0.681,0
+TIJUCA,33,163805,0.706,0.34762040606862
+ALTO DA BOA VISTA,34,9343,0.54,0.903712959245922
+MARACANA,35,25256,0.722,0
+VILA ISABEL,36,86018,0.667,0.138251762104943
+ANDARAI,37,39365,0.666,0.292364405361085
+GRAJAU,38,38671,0.688,0.646253068256266
+MANGUINHOS,39,36160,0.518,0.031673246484017
+BONSUCESSO,40,18711,0.612,0.007978892667244
+RAMOS,41,40792,0.609,0.001617044744381
+OLARIA,42,57514,0.611,0.077368292711989
+PENHA,43,78678,0.584,0.0698206751473
+PENHA CIRCULAR,44,47816,0.6,0.09156985574522
+BRAS DE PINA,45,59222,0.594,0.002547815774913
+CORDOVIL,46,45202,0.57,0.041288873104201
+PARADA DE LUCAS,47,23923,0.553,0.191353032304033
+VIGARIO GERAL,48,41820,0.531,0.010810640577933
+JARDIM AMERICA,49,25226,0.587,0.017313618926815
+HIGIENOPOLIS,50,15734,0.627,0
+JACARE,51,9276,0.571,0.004616095660138
+MARIA DA GRACA,52,7972,0.626,0
+DEL CASTILHO,53,15610,0.598,0
+INHAUMA,54,45698,0.572,0.039304844746196
+ENGENHO DA RAINHA,55,26659,0.589,0.115705046079132
+TOMAS COELHO,56,22676,0.573,0.265459494630105
+SAO FRANCISCO XAVIER,57,8343,0.613,0.150669643516846
+ROCHA,58,8766,0.627,0.242696081129974
+RIACHUELO,59,12653,0.648,0.209372074462742
+SAMPAIO,60,10895,0.578,0.17343592942983
+ENGENHO NOVO,61,42172,0.617,0.171832364587711
+LINS DE VASCONCELOS,62,37487,0.597,0.37107018261225
+MEIER,63,49828,0.687,0.006089396241485
+TODOS OS SANTOS,64,24646,0.685,0
+CACHAMBI,65,42415,0.654,0.00216392997266
+ENGENHO DE DENTRO,66,45540,0.612,0.204717019868063
+AGUA SANTA,67,8756,0.607,0.763052118139184
+ENCANTADO,68,15021,0.619,0.009099419437473
+PIEDADE,69,43378,0.61,0.219933154139474
+ABOLICAO,70,11356,0.633,0
+PILARES,71,27250,0.593,0.04191417067097
+VILA KOSMOS,72,18274,0.617,0.319105688325958
+VICENTE DE CARVALHO,73,24964,0.579,0.163098707843049
+VILA DA PENHA,74,25465,0.658,0
+VISTA ALEGRE,75,8622,0.636,0
+IRAJA,76,96382,0.615,0.039182164909094
+COLEGIO,77,29245,0.565,0.006268689072128
+CAMPINHO,78,10156,0.61,0.122995583248383
+QUINTINO BOCAIUVA,79,31185,0.603,0.418740226861129
+CAVALCANTI,80,16141,0.566,0.317906209581767
+ENGENHEIRO LEAL,81,6113,0.563,0.173740891867987
+CASCADURA,82,34456,0.589,0.185474334011247
+MADUREIRA,83,50106,0.597,0.1060210456482
+VAZ LOBO,84,15167,0.58,0.076111998918561
+TURIACU,85,17246,0.579,0.123160268399047
+ROCHA MIRANDA,86,44188,0.585,0
+HONORIO GURGEL,87,21989,0.574,0
+OSVALDO CRUZ,88,34040,0.592,0.004204190434114
+BENTO RIBEIRO,89,43707,0.61,0
+MARECHAL HERMES,90,48061,0.577,0.000226266573977
+RIBEIRA,91,3528,0.672,0.069514807468124
+ZUMBI,92,2016,0.676,0.146913542658773
+CACUIA,93,11013,0.602,0.471834074245483
+PITANGUEIRAS,94,11756,0.589,0.073103733632462
+PRAIA DA BANDEIRA,95,5948,0.648,0.139565304984348
+COCOTA,96,4877,0.655,0.2003569580166
+BANCARIOS,97,12512,0.6,0.013050129429196
+FREGUESIA,98,19437,0.61,0.464539627315501
+JARDIM GUANABARA,99,32213,0.72,0.082826998643577
+JARDIM CARIOCA,100,24848,0.61,0.003519216956309
+TAUA,101,29567,0.589,0.020028244134894
+MONERO,102,6476,0.687,0.024439220159813
+PORTUGUESA,103,23856,0.635,0.044727636192725
+GALEAO,104,22971,0.571,0.395794696689267
+CIDADE UNIVERSITARIA,105,1556,0.563,0.343948387826403
+GUADALUPE,106,47144,0.58,0.055952546443109
+ANCHIETA,107,55652,0.566,0.014094250064169
+PARQUE ANCHIETA,108,26212,0.589,0.640092477979658
+RICARDO DE ALBUQUERQUE,109,29310,0.569,0.063215207528372
+COELHO NETO,110,32423,0.573,0.000312711203756
+ACARI,111,27347,0.526,0.135738477402519
+BARROS FILHO,112,14049,0.527,0.040721231877291
+COSTA BARROS,113,28442,0.535,0.080594705036395
+PAVUNA,114,97350,0.563,0.014716897688817
+JACAREPAGUA,115,157326,0.554,0.7272854452225
+ANIL,116,24172,0.632,0.201689418600887
+GARDENIA AZUL,117,17715,0.57,0.02679216546413
+CIDADE DE DEUS,118,36515,0.559,0.02098997656309
+CURICICA,119,31189,0.58,0.130159711482686
+FREGUESIA (JACAREPAGUA),120,70511,0.64,0.387632072645149
+PECHINCHA,121,34709,0.653,0.109360496199858
+TAQUARA,122,102126,0.612,0.23467923370182
+TANQUE,123,37856,0.595,0.415821264729655
+PRACA SECA,124,64147,0.607,0.294255895484776
+VILA VALQUEIRE,125,32279,0.647,0.286117170102319
+JOA,126,818,0.76,0.391813746700496
+ITANHANGA,127,38415,0.527,0.629090109050926
+BARRA DA TIJUCA,128,135924,0.77,0.193005374803183
+CAMORIM,129,1970,0.518,0.779155854011444
+VARGEM PEQUENA,130,27250,0.519,0.646549079048662
+VARGEM GRANDE,131,14039,0.453,0.8411343450598
+RECREIO DOS BANDEIRANTES,132,82240,0.659,0.426417930488442
+GRUMARI,133,167,0.282,0.690459473698816
+DEODORO,134,10842,0.554,0.433126189030173
+VILA MILITAR,135,13184,0.583,0.425999401510555
+CAMPO DOS AFONSOS,136,1365,0.701,0.122309071795856
+JARDIM SULACAP,137,13062,0.641,0.673332674679851
+MAGALHAES BASTOS,138,24430,0.575,0.008565473925822
+REALENGO,139,180123,0.574,0.527090557540516
+PADRE MIGUEL,140,64228,0.582,0.10775660562844
+BANGU,141,243125,0.57,0.521255940239241
+SENADOR CAMARA,142,105515,0.554,0.584448695996251
+SANTISSIMO,143,41458,0.558,0.404179624642746
+CAMPO GRANDE,144,328370,0.572,0.539762120692714
+SENADOR VASCONCELOS,145,30600,0.556,0.479352243625703
+INHOAIBA,146,64649,0.543,0.272995062101142
+COSMOS,147,77007,0.542,0.327173948679956
+PACIENCIA,148,94626,0.536,0.502615175374195
+SANTA CRUZ,149,217333,0.527,0.58913544234935
+SEPETIBA,150,56575,0.517,0.219959290182041
+GUARATIBA,151,110049,0.487,0.697224612725608
+BARRA DE GUARATIBA,152,3577,0.502,0.693142195152656
+PEDRA DE GUARATIBA,153,9488,0.559,0.216846539216702
+ROCINHA,154,69356,0.533,0.169149439479229
+JACAREZINHO,155,37839,0.534,0.002850987739388
+COMPLEXO DO ALEMAO,156,69143,0.532,0.1188597100481
+MARE,157,129770,0.547,0.026309141585875
+PARQUE COLUMBIA,158,9202,0.549,0.034668835576057
+VASCO DA GAMA,159,15482,0.593,0
+GERICINO,160,15167,0.545,0.357699294514938