-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathsample_code_pseudocode.txt
More file actions
162 lines (131 loc) · 4.57 KB
/
Copy pathsample_code_pseudocode.txt
File metadata and controls
162 lines (131 loc) · 4.57 KB
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
START
Load required libraries:
sf, terra, dplyr, spdep, INLA, ggplot2
Set random seed
--------------------------------------------------
1. CREATE FAKE DISTRICTS
--------------------------------------------------
Define a Singapore-like bounding box
Convert bounding box into a spatial object
Create a regular grid of square polygons over the box
Assign each grid cell a district ID
Set coordinate reference system
Make geometries valid
--------------------------------------------------
2. CREATE FAKE ENVIRONMENTAL RASTERS
--------------------------------------------------
Create an empty raster template covering the same area
Get x-y coordinates of all raster cells
Simulate rainfall values:
rainfall depends on east-west gradient
rainfall depends on north-south gradient
add random noise
Simulate soil values:
soil depends on sinusoidal spatial pattern
add random noise
Store these simulated values in raster layers:
rain raster
soil raster
--------------------------------------------------
3. SIMULATE CASE LOCATIONS
--------------------------------------------------
Create 2 artificial spatial hotspots
Construct a linear predictor using:
rainfall effect
soil effect
hotspot 1 effect
hotspot 2 effect
random noise
Convert linear predictor into risk values
Convert risk values into sampling probabilities
Choose total number of cases to simulate
Sample raster cells according to the risk probabilities
Jitter sampled points within each raster cell
Create spatial point object for case locations
--------------------------------------------------
4. AGGREGATE CASES TO DISTRICTS
--------------------------------------------------
Spatially join each case point to a district polygon
If any point falls exactly on a boundary and is unmatched:
assign it to the nearest district
Count number of cases in each district
Join district case counts back to district polygons
Replace missing case counts with zero
--------------------------------------------------
5. EXTRACT DISTRICT-LEVEL COVARIATES
--------------------------------------------------
For each district:
calculate mean rainfall from raster
calculate mean soil value from raster
Generate a fake population for each district
Compute expected number of cases E:
total cases × district population / total population
Create numeric region ID for modelling
--------------------------------------------------
6. BUILD SPATIAL ADJACENCY
--------------------------------------------------
Determine which districts share borders with each other
Convert neighbourhood structure into INLA adjacency graph
Read graph into INLA format
--------------------------------------------------
7. PREPARE MODEL DATA
--------------------------------------------------
Drop spatial geometry to create model table
Standardize rainfall
Standardize soil
--------------------------------------------------
8. FIT INLA BYM2 MODEL
--------------------------------------------------
Define model formula:
observed cases ~ rainfall + soil + spatial random effect
Use Poisson likelihood
Use expected counts E as offset/exposure
Fit BYM2 spatial model in INLA
Request:
fitted values
DIC
WAIC
CPO
configuration for posterior summaries
--------------------------------------------------
9. ATTACH MODEL OUTPUTS
--------------------------------------------------
Extract fitted mean counts and credible intervals
For each district:
compute relative risk = fitted mean / expected cases
compute lower and upper RR interval
Extract linear predictor summaries
Approximate posterior probability that RR > 1
Store this probability for each district
--------------------------------------------------
10. CHECK OUTPUTS
--------------------------------------------------
Print fixed effects
Print hyperparameters
Print model fit statistics
Print summary of:
observed cases
relative risk
posterior probability RR > 1
Create district risk ranking table
Sort districts from highest to lowest relative risk
Print ranking
--------------------------------------------------
11. MAKE MAPS
--------------------------------------------------
Map 1:
district boundaries + simulated case points
Map 2:
district-level estimated relative risk
Map 3:
district-level posterior probability that RR > 1
Display all maps
--------------------------------------------------
12. SAVE OPTIONAL OUTPUTS
--------------------------------------------------
Save case coordinates to text file
Save district polygons to shapefile
Save rainfall raster to GeoTIFF
Save soil raster to GeoTIFF
Save district risk results to GeoJSON
END