-
Notifications
You must be signed in to change notification settings - Fork 11
/
README.Rmd
157 lines (119 loc) · 5.1 KB
/
README.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
---
output:
md_document:
variant: gfm
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "README-figs/",
cache.path = "README-cache/"
)
options(width = 120) # for printing Stan output
```
# glmmfields <img src="inst/logo.png" align="right" />
[![R build status](https://github.com/seananderson/glmmfields/workflows/R-CMD-check/badge.svg)](https://github.com/seananderson/glmmfields/actions)
[![CRAN status](https://www.r-pkg.org/badges/version/glmmfields)](https://cran.r-project.org/package=glmmfields)
<!-- [![codecov](https://codecov.io/github/seananderson/glmmfields/branch/master/graphs/badge.svg)](https://codecov.io/github/seananderson/glmmfields) -->
The glmmfields R package implements Bayesian spatiotemporal models that allow for extreme
spatial deviations through time. It uses a predictive process approach with
random fields implemented through a multivariate-t distribution instead of a
multivariate normal. The models are fit with [Stan](http://mc-stan.org/).
We published a paper describing the model and package in *Ecology*:
Anderson, S. C., Ward, E. J. 2019. Black swans in space: modelling
spatiotemporal processes with extremes. 100(1):e02403.
<https://doi.org/10.1002/ecy.2403>
You can install the [CRAN version](https://cran.r-project.org/package=glmmfields) of the package with:
```{r, eval=FALSE}
install.packages("glmmfields")
```
If you have a C++ compiler installed, you can install the development version of the package with:
```{r, eval=FALSE}
# install.packages("remotes")
remotes::install_github("seananderson/glmmfields", build_vignettes = TRUE)
```
glmmfields can also fit spatial GLMs with Stan. See the vignette:
```{r, eval=FALSE}
vignette("spatial-glms", package = "glmmfields")
```
## An example spatiotemporal model
```{r libraries, cache=FALSE}
library(glmmfields)
library(ggplot2)
```
Simulate data:
```{r simulate, cache=TRUE}
set.seed(42)
s <- sim_glmmfields(
df = 2.8, n_draws = 12, n_knots = 12, gp_theta = 2.5,
gp_sigma = 0.2, sd_obs = 0.1
)
head(s$dat)
```
```{r plot-sim, cache=TRUE, dependson="simulate"}
print(s$plot)
```
Fit the model:
```{r fit, cache=TRUE, warning=FALSE, message=FALSE, results='hide', dependson="simulate", cache.comments=FALSE}
options(mc.cores = parallel::detectCores()) # for parallel processing
m <- glmmfields(y ~ 0,
data = s$dat, time = "time",
lat = "lat", lon = "lon",
nknots = 12, estimate_df = TRUE, iter = 800, seed = 1
)
```
```{r print, cache=FALSE, dependson="fit"}
print(m)
```
Plot:
```{r plot-predictions, dependson="fit"}
plot(m, type = "prediction") + scale_color_gradient2()
plot(m, type = "spatial-residual")
```
Predictions:
```{r get-predictions, dependson="fit"}
# link scale:
p <- predict(m)
head(p)
# posterior predictive intervals on new observations (include observation error):
p <- predictive_interval(m)
head(p)
```
Use the `tidy` method to extract parameter estimates as a data frame:
```{r tidy, dependson="fit"}
x <- tidy(m, conf.int = TRUE)
head(x)
```
Make predictions on a fine-scale spatial grid:
```{r grid-predictions, dependson="fit"}
pred_grid <- expand.grid(
lat = seq(min(s$dat$lat), max(s$dat$lat), length.out = 25),
lon = seq(min(s$dat$lon), max(s$dat$lon), length.out = 25),
time = unique(s$dat$time)
)
pred_grid$prediction <- predict(m,
newdata = pred_grid, type = "response", iter = 100, estimate_method = "median"
)$estimate
ggplot(pred_grid, aes(lon, lat, fill = prediction)) +
facet_wrap(~time) +
geom_raster() +
scale_fill_gradient2()
```
# References
Anderson, S. C., Ward, E. J. 2019. Black swans in space: modelling
spatiotemporal processes with extremes. 100(1):e02403.
<https://doi.org/10.1002/ecy.2403>
Latimer, A. M., S. Banerjee, H. Sang Jr, E. S. Mosher, and J. A. Silander Jr. 2009. Hierarchical models facilitate spatial analysis of large data sets: a case study on invasive plant species in the northeastern United States. Ecology Letters 12:144–154.
Shelton, A. O., J. T. Thorson, E. J. Ward, and B. E. Feist. 2014. Spatial semiparametric models improve estimates of species abundance and distribution. Canadian Journal of Fisheries and Aquatic Sciences 71:1655–1666.
### NOAA Disclaimer
This repository is a scientific product and is not official communication of the National Oceanic and
Atmospheric Administration, or the United States Department of Commerce. All NOAA GitHub project code is
provided on an ‘as is’ basis and the user assumes responsibility for its use. Any claims against the Department of
Commerce or Department of Commerce bureaus stemming from the use of this GitHub project will be governed
by all applicable Federal law. Any reference to specific commercial products, processes, or services by service
mark, trademark, manufacturer, or otherwise, does not constitute or imply their endorsement, recommendation or
favoring by the Department of Commerce. The Department of Commerce seal and logo, or the seal and logo of a
DOC bureau, shall not be used in any manner to imply endorsement of any commercial product or activity by
DOC or the United States Government.