-
Notifications
You must be signed in to change notification settings - Fork 0
/
matern-on-graph.Rmd
339 lines (297 loc) · 8.47 KB
/
matern-on-graph.Rmd
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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
---
title: "Matern SPDE model on Metric Graph"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Matern SPDE model on Metric Graph}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
references:
- id: BSW2022a
title: "Gaussian Whittle--Matérn fields on metric graphs"
author:
- family: Bolin
given: David
- family: Simas
given: Alexandre B.
- family: Wallin
given: Jonas
container-title: Bernoulli
type: article
issued:
year: 2023
- id: BSW2022b
title: "Statistical properties of Gaussian Whittle--Matérn fields on metric graphs"
author:
- family: Bolin
given: David
- family: Simas
given: Alexandre B.
- family: Wallin
given: Jonas
container-title: arXiv:2304.10372
type: preprint
issued:
year: 2023
- id: Anderes2020
title: "Isotropic covariance functions on graphs and their edges"
author:
- family: Anderes
given: Ethan
- family: Møller
given: Jesper
- family: Rasmussen
given: Jakob G
container-title: Annals of Statistics
type: article
issue: 48
pages: 2478-2503
issued:
year: 2020
- id: Borovitskiy2021
title: "Matérn Gaussian processes on graphs"
author:
- family: Borovitskiy
given: Viacheslav
- family: Azangulov
given: Iskander
- family: Terenin
given: Alexander
- family: Mostowsky
given: Peter
- family: Deisenroth
given: Marc
- family: Durrande
given: Nicolas
container-title: International Conference on Artificial Intelligence and Statistics
type: article
pages: 2593-2601
issued:
year: 2021
- id: LindgrenRue2015
title: Bayesian Spatial Modelling with R-INLA.
author:
- family: Lindgren
given: Finn
- family: Rue
given: Haavard
container-title: Journal of Statistical Software
type: article
issue: 63
pages: 1-25
issued:
year: 2015
- id: inlabru2019
title: inlabru an R package for Bayesian spatial modelling from ecological survey data
author:
- family: Bachl
given: Fabian E.
- family: Lindgren
given: Finn
- family: Borchers
given: David L.
- family: Illian
given: Janine B.
container-title: Methods in Ecology and Evolution
type: article
issue: 10
pages: 760-766
issued:
year: 2019
- id: sppackage
title: Applied spatial data analysis with R
author:
- family: Bivand
given: Roger S.
- family: Pebesma
given: Edzer
- family: Gomez-Rubio
given: Virgilio
publisher: Springer, NY
type: book
issued:
year: 2013
- id: plotlypackage
title: Interactive Web-Based Data Visualization with R, plotly, and shiny
author:
- family: Sievert
given: Carson
publisher: Chapman and Hall/CRC
type: book
issued:
year: 2020
# Bolin D, Wallin J, Simas A (2023). MetricGraph: Random fields on metric graphs. R package version 1.1.2, https://CRAN.R-project.org/package=MetricGraph.
- id: MetricGraph
title: "MetricGraph: Random fields on metric graphs. R package version 1.1.2"
author:
- family: Bolin
given: David
- family: Simas
given: Alexandre B.
- family: Wallin
given: Jonas
publisher: https://CRAN.R-project.org/package=MetricGraph.
# url:
# - https://CRAN.R-project.org/package=MetricGraph
note:
- R package version 1.1.2
issued:
year: 2023
---
# Introduction
Metric graph is a class of graphs that are embedded in a metric space. This type of graph can model some specific spatial structures, including road networks, river networks, etc.
In this vignette, we will show how to model a Matern SPDE model on a metric graph together with the `R` package `MetricGraph` [@MetricGraph](https://CRAN.R-project.org/package=MetricGraph). It contains functions for working with data and
random fields on compact metric graphs. The main functionality is contained in
the `metric_graph` class, which is used for specifying metric graphs, adding data to them, visualization, and other basic functions that are needed for working with data and random fields on metric graphs.
To know more about the theory of Matern SPDE model on metric graphs, please refer to Bolin et al. [@BSW2022a; @BSW2022b].
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
```{r setup}
library(MetricGraph)
library(ngme2)
set.seed(16)
```
# Simulation on a metric graph
First we will build the graph and construct a uniform mesh on each edges.
```{r fig.align = "center",echo=TRUE}
edge1 <- rbind(c(0,0),c(1,0))
edge2 <- rbind(c(0,0),c(0,1))
edge3 <- rbind(c(0,1),c(-1,1))
theta <- seq(from=pi,to=3*pi/2,length.out = 20)
edge4 <- cbind(sin(theta),1+ cos(theta))
edges = list(edge1, edge2, edge3, edge4)
graph <- metric_graph$new(edges = edges)
graph$plot()
# construct the mesh
graph$build_mesh(h = 0.1)
graph$plot(mesh = TRUE)
# Refine the mesh and print how many mesh nodes
graph$build_mesh(h = 0.005)
length(graph$mesh$h)
```
Next, we simulation non-Gaussian Matern field with NIG noise using `simulate` function.
```{r}
simu_nig <- noise_nig(mu=-10, sigma=30, nu=0.5)
matern_graph <- ngme2::f(
model="matern",
theta_K = log(6), #kappa = 8
mesh=graph,
noise=simu_nig
)
matern_graph
W <- simulate(matern_graph, seed=10)[[1]]
graph$plot_function(X=as.numeric(W))
# 3d plot of the random field
graph$plot_function(X=as.numeric(W), plotly=TRUE)
```
With the latent field generated, next we construct observations.
```{r, fig.align="center",echo=TRUE}
# build observation and A matrix
obs.per.edge <- 200 # how many observations per edge
obs.loc <- NULL
for(i in 1:graph$nE) {
obs.loc <- rbind(obs.loc,
cbind(rep(i,obs.per.edge), runif(obs.per.edge)))
}
n.obs <- obs.per.edge*graph$nE
A <- graph$fem_basis(obs.loc)
# using the model for inference
sigma.e <- 0.1
x1 <- obs.loc[,1]
x2 <- obs.loc[,2]
Y <- 2*x1 - 3*x2 + as.vector(A %*% W + sigma.e * rnorm(n.obs))
df_data <- data.frame(y = Y, edge_number = obs.loc[,1],
distance_on_edge = obs.loc[,2],
x1 = x1, x2 = x2)
graph$clear_observations()
graph$add_observations(data = df_data, normalized = TRUE)
graph$plot(data = "y")
```
Next we estimate the model using ngme2:
```{r model-fit, cache=TRUE, fig.align="center", echo=TRUE}
fit_nig <- ngme(
y ~ 0 + x1 + x2 +
f(model="matern", mesh=graph, name="graph", noise=noise_nig()),
data = graph$get_data(),
control_opt = control_opt(
iterations = 1000,
iters_per_check = 100,
rao_blackwellization = TRUE,
n_parallel_chain = 4,
# verbose = TRUE,
seed = 1
)
)
fit_nig
traceplot(fit_nig, "graph")
traceplot(fit_nig)
```
# Graph model with replicates
The model with replicates is similar to other models in ngme2. We simply provide the `replicate` argument in `ngme` function.
```{r}
n_repl <- 5
df_data <- NULL
for (repl in 1:n_repl) {
obs.per.edge <- 200 # how many observations per edge
obs.loc <- NULL
for(i in 1:graph$nE) {
obs.loc <- rbind(obs.loc,
cbind(rep(i,obs.per.edge), runif(obs.per.edge)))
}
n.obs <- obs.per.edge*graph$nE
A <- graph$fem_basis(obs.loc)
# using the model for inference
sigma.e <- 0.1
x1 <- obs.loc[,1]
x2 <- obs.loc[,2]
Y <- 2*x1 - 3*x2 + as.vector(A %*% W + sigma.e * rnorm(n.obs))
df_data_tmp <- data.frame(y = Y, edge_number = obs.loc[,1],
distance_on_edge = obs.loc[,2],
x1 = x1, x2 = x2, repl = repl)
df_data <- rbind(df_data, df_data_tmp)
}
graph$clear_observations()
graph$add_observations(data = df_data, normalized = TRUE, group="repl")
```
Next step is easy, we simply provide extra `replicate` argument to be same as the group from the graph function.
```{r fig.align = "center",echo=TRUE}
fit_repl <- ngme(
y ~ 0 + x1 + x2 +
f(model="matern", mesh=graph, name="graph", noise=noise_nig()),
data = graph$get_data(),
control_opt = control_opt(
iterations = 2000,
rao_blackwellization = TRUE,
n_parallel_chain = 4,
# verbose = TRUE,
seed = 1
),
replicate = ".group"
)
fit_repl
traceplot(fit_repl, "graph")
traceplot(fit_repl)
```
# Kringing
Kringing is a method for predicting the value of a random field at an unobserved location based on observations of the field at other locations.
We can use `predict` function in `ngme2` to do kringing.
```{r fig.align = "center",echo=TRUE}
# create new observations
locs <- graph$get_mesh_locations()
X = as.data.frame(locs)
names(X) <- c("x1", "x2")
preds <- predict(
fit_repl,
data = X,
# data = list(x1=0, x2=0),
map = list(graph = locs)
)
# plot the prediction
graph$plot_function(preds$mean)
# compare with the true data
graph$plot(data="y")
```
# References {-}