Muylaert_et_al_2019/script_muylaert_et_al_2019.R at master · renatamuy/Muylaert_et_al_2019 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# Script Muylaert et al. (2019)

  ####################################
  setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
  ls()
  rm(list = ls())
  gc()
  options(digits=7, scipen=999)
  memory.limit(size= 1.75e13)
  # packages
  require(dplyr)
  require(INLA)
  require(coefINLA)
  # data
  load("data_9316.rda") # data
  load("g_models.rda")  # graph
  # casos
  data_9316 %>% group_by(year) %>% dplyr::summarise(sum(cases)) %>% as.data.frame()
  # casos transformados em zero e um
  data_9316 %>% group_by(year) %>% dplyr::summarise(sum(casesbin)) %>% as.data.frame()

  ##########################################
  # Model with spatial and temporal effects#
  ##########################################
 # Subsetting data
d_0016 <- data_9316[data_9316$variable %in% 2000:2016,] # for binomial
# Creating df from positive counts
d_pos_0016 <- subset(d_0016, cases > 0 ) # for truncated poisson

# Set data from 2015 and 2016 to NA for binomial models
d_0016[d_0016$variable %in% 2015:2016, "casesbin"] <- NA
# Set count data from 2015 and 2016 to NA
d_pos_0016[d_pos_0016$variable %in% 2015:2016, "cases"] <- NA

# binomial
f_mb <-  casesbin ~ scale(s_hosts) +
  scale(rural_workers_approxExtrap) +
  scale(mplu) +
  scale(t) +
  scale(pfor) +
  scale(pasture) +
  scale(sug) +
  scale(peu) +
  scale(mz) +
  f(re_u, model="besag", graph= g, adjust.for.con.comp= FALSE,
    constr= TRUE, scale.model= TRUE) +
  f(year, model = "rw2")

mb_NA1516 <- inla(formula = f_mb,
                        family = "binomial",
                        data = d_0016,
                        verbose = TRUE,
                        control.compute = list(config = TRUE, dic = TRUE, cpo=TRUE))

# Save model
save(mb_NA1516, file = "mb_NA1516.rda")

require(coefINLA)
coefINLA(mb_NA1516, intercept = TRUE, "Greys") +
  labs(title = "Posterior distributions",
       subtitle = "with 95% credible intervals") +
  ggsave(filename = "coefINLA_mb_NA1516.png", dpi = 600,
         width = 18, height = 14, units= "cm")


# Truncated Poisson model
# a very negative and fixed hyperpar approximates pi i to zero
Theta_neg_fixed <-  list(hyper = list(theta = list(initial = -10, fixed = TRUE)))

f_mpos <-  cases ~  scale(s_hosts) +
  scale(rural_workers_approxExtrap) +
  scale(mplu) +
  scale(t) +
  scale(pfor) +
  scale(pasture) +
  scale(sug) +
  scale(peu) +
  scale(mz) +
  f(re_u, model="besag", graph= g, adjust.for.con.comp= FALSE,
    constr= TRUE, scale.model= TRUE) +
  f(year, model = "rw2")

mpos_NA1516 <- inla(formula = f_mpos,
             family = "zeroinflatedpoisson0",
             control.family = list(list(Theta_neg_fixed)),
             data = d_pos_0016,
             verbose = TRUE,
             control.compute = list(config = TRUE, dic = TRUE, cpo = TRUE))

# Improving model
mpos_NA1516_improved <- inla.cpo(mpos_NA1516)

# Save model
save(mpos_NA1516, file = "mb_NA1516.rda")
save(mpos_NA1516_improved, file = "mb_NA1516_improved.rda")

# Coef plot
coefINLA(mpos_NA1516, intercept = TRUE, "Greys") +
  labs(title = "Posterior distributions",
       subtitle = "with 95% credible intervals") +
  ggsave(filename = "coefINLA_mpos_NA1516.png", dpi = 600,
         width = 18, height = 14, units= "cm")