-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path02_model_fitting.R
117 lines (81 loc) · 4.53 KB
/
02_model_fitting.R
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
103
104
105
106
107
108
109
110
111
112
113
114
###############################################################
#### ####
#### Script to fit spatio-temporal additive models ####
#### ####
###############################################################
## Load packages and data
source("00_load_packages_data.R")
##### WARNING: Script may take over a day to run #####
#### Fit baseline model (with only spatio-temporal smooths) ####
model_base <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)),
data = df_model, family = binomial, method = "REML")
# summary(model_base)
# gam.check(model_base)
saveRDS(model_base, file = "output/model_base.rds")
#### Fit full model (with urbanisation, REGIC, months suitable and prior outbreak) ####
model_full <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
urban_prpn + regic_comb + months_suitable.both +
prior_outbreak_fix,
data = df_model, family = binomial, method = "REML")
# summary(model_full)
# gam.check(model_full)
saveRDS(model_full, file = "output/model_full.rds")
#### Add the number of extremely wet months ####
model_wet <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
urban_prpn + regic_comb + months_suitable.both +
prior_outbreak_fix + months_wet,
data = df_model, family = binomial, method = "REML")
# summary(model_wet)
# gam.check(model_wet)
saveRDS(model_wet, file = "output/model_wet.rds")
#### Fit full model using aegypti cut-off (sensitivity analysis) ####
model_aeg <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
urban_prpn + regic_comb + months_suitable.aeg + prior_outbreak_fix,
data = df_model, family = binomial, method = "REML")
# summary(model_aeg)
# gam.check(model_aeg)
saveRDS(model_aeg, file = "output/model_aeg.rds")
#### Fit outbreak100 model (sensitivity analysis) ####
model_full100 <- bam(outbreak_fix100 ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
urban_prpn + regic_comb + months_suitable.both + prior_outbreak_fix100,
data = df_model, family = binomial, method = "REML")
# summary(model_full100)
# gam.check(model_full100)
saveRDS(model_full100, file = "output/model_full100.rds")
#### Fit 75th percentile model (sensitivity analysis) ####
model_full_perc75 <- bam(outbreak_perc75 ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
urban_prpn + regic_comb + months_suitable.both + prior_outbreak_perc75,
data = df_model, family = binomial, method = "REML")
# summary(model_full_perc75)
# gam.check(model_full_perc75)
saveRDS(model_full_perc75, file = "output/model_full_perc75.rds")
#### Add each covariate in turn (sensitivity analysis) ####
## Months suitable
model_temp <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
months_suitable.both,
data = df_model, family = binomial, method = "REML")
saveRDS(model_temp, file = "output/model_temp.rds")
## Prior outbreak
model_prior <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
prior_outbreak_fix,
data = df_model, family = binomial, method = "REML")
saveRDS(model_prior, file = "output/model_prior.rds")
## Urbanisation
model_urb <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
urban_prpn,
data = df_model, family = binomial, method = "REML")
saveRDS(model_urb, file = "output/model_urb.rds")
model_regic <- bam(outbreak_fix ~ s(lon, lat, k = 400) + s(year, k = 10) +
ti(lon, lat, year, d = c(2, 1), k = c(100, 8)) +
regic_comb,
data = df_model, family = binomial, method = "REML")
saveRDS(model_regic, file = "output/model_regic.rds")